tsplitbt <-
function(J,K,x,tr=.2,alpha=.05,JK=J*K,grp=c(1:JK),nboot=599,
SEED=TRUE,monitor=FALSE){
#
# A bootstrap-t for performing a split-plot design
# with trimmed means.
# By default, 20% trimming is used with B=599 bootstrap samples.
#
#
# The R variable x is assumed to contain the raw
# data stored in list mode or in a matrix.
# If in list mode, x[[1]] contains the data
# for the first level of both factors: level 1,1.
# x[[2]] is assumed to contain the data for level 1 of the
# first factor and level 2 of the second: level 1,2
# x[[K]] is the data for level 1,K
# x[[K+1]] is the data for level 2,1, x[[2K]] is level 2,K, etc.
#
# If the data are in a matrix, column 1 is assumed to
# correspond to x[[1]], column 2 to x[[2]], etc.
#
# When in list mode x is assumed to have length JK, the total number of
# groups being tested, but a subset of the data can be analyzed
# using grp
#
if(SEED)set.seed(2)
if(is.data.frame(x))x=as.matrix(x)
if(is.matrix(x)) {
y <- list()
ik=0
il=c(1:K)-K
for(j in 1:J){
il=il+K
zz=x[,il]
zz=elimna(zz)
for(k in 1:K){
ik=ik+1
y[[ik]]=zz[,k]
}}
x <- y
}
JK<-J*K
data<-list()
xcen<-list()
for(j in 1:length(x)){
data[[j]]<-x[[grp[j]]] # Now have the groups in proper order.
xcen[[j]]<-data[[j]]-mean(data[[j]],tr) # Centered data
}
x<-data
#
set.seed(2) # set seed of random number generator so that
# results can be duplicated.
# Next determine the n_j values
nvec<-NA
jp<-1-K
for(j in 1:J){
jp<-jp+K
nvec[j]<-length(x[[j]])
}
if(min(nvec)<10){
print("Warning: with small sample sizes, a bootstrap-t method can")
print("result in estimated variances equal to zero, resulting in")
print("this function terminating and giving no results and peculiar error messages.")
}
blist<-list()
print("Taking bootstrap samples. Please wait.")
testmat<-matrix(NA,ncol=3,nrow=nboot)
for(iboot in 1:nboot){
iv<-0
for(j in 1:J){
temp<-sample(nvec[j],replace = T)
for(k in 1:K){
iv<-iv+1
tempx<-xcen[[iv]]
blist[[iv]]<-tempx[temp]
}}
if(monitor)print(paste("Bootstrap iteration" ,iboot, "is complete"))
btest<-tsplit(J,K,blist,tr)
testmat[iboot,1]<-btest$Qa
testmat[iboot,2]<-btest$Qb
testmat[iboot,3]<-btest$Qab
}
lcrit<-round((1-alpha)*nboot)
temp3<-sort(testmat[,1])
crit.Qa<-temp3[lcrit]
temp3<-sort(testmat[,2])
crit.Qb<-temp3[lcrit]
temp3<-sort(testmat[,3])
crit.Qab<-temp3[lcrit]
temp4<-tsplit(J,K,x,tr=tr)
list(Qa=temp4$Qa,Qb=temp4$Qb,Qab=temp4$Qab,crit.Qa=crit.Qa,crit.Qb=crit.Qb,crit.Qab=crit.Qab)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.