pairdepb <-
function(x,tr=.2,alpha=.05,grp=0,nboot=599){
#
# Using the percentile t bootstrap method,
# compute a .95 confidence interval for all pairwise differences between
# the trimmed means of dependent groups.
# By default, 20% trimming is used with B=599 bootstrap samples.
#
# x can be an n by J matrix or it can have list mode
#
if(is.data.frame(x)) x <- as.matrix(x)
if(!is.list(x) && !is.matrix(x))stop("Data must be stored in a matrix or in list mode.")
if(is.list(x)){
if(sum(grp)==0)grp<-c(1:length(x))
# put the data in an n by J matrix
mat<-matrix(0,length(x[[1]]),length(grp))
for (j in 1:length(grp))mat[,j]<-x[[grp[j]]]
}
if(is.matrix(x)){
if(sum(grp)==0)grp<-c(1:ncol(x))
mat<-x[,grp]
}
if(sum(is.na(mat)>=1))stop("Missing values are not allowed.")
J<-ncol(mat)
connum<-(J^2-J)/2
bvec<-matrix(0,connum,nboot)
set.seed(2) # set seed of random number generator so that
# results can be duplicated.
print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(nrow(mat),size=nrow(mat)*nboot,replace=TRUE),nrow=nboot)
xcen<-matrix(0,nrow(mat),ncol(mat))
for (j in 1:J)xcen[,j]<-mat[,j]-mean(mat[,j],tr) #Center data
it<-0
for (j in 1:J){
for (k in 1:J){
if(j<k){
it<-it+1
bvec[it,]<-apply(data,1,tsub,xcen[,j],xcen[,k],tr)
# bvec is a connum by nboot matrix containing the bootstrap test statistics.
}}}
bvec<-abs(bvec) #Doing two-sided confidence intervals
icrit<-round((1-alpha)*nboot)
critvec<-apply(bvec,2,max)
critvec<-sort(critvec)
crit<-critvec[icrit]
psihat<-matrix(0,connum,5)
dimnames(psihat)<-list(NULL,c("Group","Group","psihat","ci.lower","ci.upper"))
test<-matrix(NA,connum,4)
dimnames(test)<-list(NULL,c("Group","Group","test","se"))
it<-0
for (j in 1:J){
for (k in 1:J){
if(j<k){
it<-it+1
estse<-yuend(mat[,j],mat[,k])$se
dif<-mean(mat[,j],tr)-mean(mat[,k],tr)
psihat[it,1]<-grp[j]
psihat[it,2]<-grp[k]
psihat[it,3]<-dif
psihat[it,4]<-dif-crit*estse
psihat[it,5]<-dif+crit*estse
test[it,1]<-grp[j]
test[it,2]<-grp[k]
test[it,3]<-yuend(mat[,j],mat[,k])$teststat
test[it,4]<-estse
}}}
list(test=test,psihat=psihat,crit=crit)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.