R/tsplitbt.R

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)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.