1
\$\begingroup\$

I recently read an article in which the authors present Algorithm 1 to compute a simultaneous confidence interval. I was interested in implementing Algorithm 1 in R here. The algorithm is freely available at the link.

Please review this to verify that I've implemented the algorithm correctly. It uses a toy example with the built-in R datasets "trees." I tried to comment my code and name variables in a way similar to the published algorithm.

library(MASS)
alpha<-.05
B=100000
#Use toy example from R's trees dataset
myfit<-lm(Height~Girth+Volume, data=trees)
myfit
param<-coef(myfit)
coefficients<-length(param)
mycov<-vcov(myfit)
#Simulate 10K bootstrap samples using the estimated mean vector and covariance matrix
#from the regression model
#Step 1.
theata.tilde.b<-mvrnorm(B, mu=param, Sigma=mycov)
head(theata.tilde.b)

#Step 2
rbj<-apply(theata.tilde.b, 2, rank)

#Step 3.
sample.b.rank.upper<-apply(rbj, 1, max)
sample.b.rank.lower<-apply(rbj, 1, min)

boot.df<-data.frame(cbind(theata.tilde.b, rbj, sample.b.rank.upper, sample.b.rank.lower))
names(boot.df)<-c("B0", "B1", "B2", "Rank0", "Rank1", "Rank2", "SampleBRankUpper", "SampleBRankLower")

upper.index<-B*(1-alpha/2)
lower.index<-B*(alpha/2)
#Sort datasets by max and then min rank of the sample-b rank
upper<-(boot.df[order(boot.df$SampleBRankUpper),])
lower<-(boot.df[order(boot.df$SampleBRankLower),])

#Step 5
upper[upper.index,1:coefficients]
lower[lower.index,1:coefficients]
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.