R/t2waypbg.R

t2waypbg <-
function(J,K,x,alpha=.05,nboot=NA,grp=NA,est=onestep,...){
#
#   Two-way ANOVA for independent groups based on
#   robust measures of location
#   and a percentile bootstrap method.

#   The data are assumed to be stored in x in list mode or in a matrix.
       #  If grp is unspecified, it is assumed 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 factor: level 1,2
        #  x[[j+1]] is the data for level 2,1, etc.
        #  If the data are in wrong order, grp can be used to rearrange the
        #  groups. For example, for a two by two design, grp<-c(2,4,3,1)
        #  indicates that the second group corresponds to level 1,1;
        #  group 4 corresponds to level 1,2; group 3 is level 2,1;
        #  and group 1 is level 2,2.
#
#   Missing values are automatically removed.
#
if(is.data.frame(x))x=as.matrix(x)
JK<-J*K
if(is.matrix(x))x<-listm(x)
if(!is.na(grp)){
yy<-x
for(j in 1:length(grp))
x[[j]]<-yy[[grp[j]]]
}
if(!is.list(x))stop("Data must be stored in list mode or a matrix.")
for(j in 1:JK){
xx<-x[[j]]
x[[j]]<-xx[!is.na(xx)] # Remove any missing values.
}
#
# Create the three contrast matrices
#
       ij <- matrix(c(rep(1, J)), 1, J)
        ik <- matrix(c(rep(1, K)), 1, K)
       jm1 <- J - 1
        cj <- diag(1, jm1, J)
        for(i in 1:jm1)
                cj[i, i + 1] <- 0 - 1
        km1 <- K - 1
        ck <- diag(1, km1, K)
        for(i in 1:km1)
                ck[i, i + 1] <- 0 - 1
conA<-t(kron(cj,ik))
conB<-t(kron(ij,ck))
conAB<-t(kron(cj,ck))
ncon<-max(nrow(conA),nrow(conB),nrow(conAB))
if(JK!=length(x)){
print("Warning: The number of groups does not match")
print("the number of contrast coefficients.")
}
if(is.na(nboot)){
nboot<-5000
if(ncon<=4)nboot<-2000
}
m1<-matrix(0,nrow=JK,ncol=nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
for(j in 1:JK){
paste("Working on group ",j)
data<-matrix(sample(x[[j]],size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
m1[j,]<-apply(data,1,est,...)
}
bootA<-matrix(0,ncol(conA),nboot)
bootB<-matrix(0,ncol(conB),nboot)
bootAB<-matrix(0,ncol(conAB),nboot)
testA<-NA
testB<-NA
testAB<-NA
testvecA<-NA
testvecB<-NA
testvecAB<-NA
for (d in 1:ncol(conA)){
bootA[d,]<-apply(m1,2,trimpartt,conA[,d])
# A vector of length nboot containing psi hat values
# corresponding to the dth linear contrast
testA[d]<-sum((bootA[d,]>0))/nboot
testA[d]<-min(testA[d],1-testA[d])
}
for (d in 1:ncol(conB)){
bootB[d,]<-apply(m1,2,trimpartt,conB[,d])
# A vector of length nboot containing psi hat values
# corresponding to the dth linear contrast
testB[d]<-sum((bootB[d,]>0))/nboot
testB[d]<-min(testB[d],1-testB[d])
}
for (d in 1:ncol(conAB)){
bootAB[d,]<-apply(m1,2,trimpartt,conAB[,d])
# A vector of length nboot containing psi hat values
# corresponding to the dth linear contrast
testAB[d]<-sum((bootAB[d,]>0))/nboot
testAB[d]<-min(testAB[d],1-testAB[d])
}
#
#  Determine critical value
#
Jm<-J-1
Km<-K-1
JKm<-(J-1)*(K-1)
dvecA <- alpha/c(1:Jm)
dvecB <- alpha/c(1:Km)
dvecAB <- alpha/c(1:JKm)
testA<-(0 - 1) * sort(-2 * testA)
testB<-(0 - 1) * sort(-2 * testB)
testAB<-(0 - 1) * sort(-2 * testAB)
sig <- sum((testA < dvecA[1:Jm]))
if(sig > 0)
print("Significant result obtained for Factor A: Reject")
if(sig == 0)
print("No significant result Factor A: Fail to reject")
sig <- sum((testB < dvecB[1:Km]))
if(sig > 0)
print("Significant result obtained for Factor B: Reject")
if(sig == 0)
print("No significant result Factor B: Fail to reject")
sig <- sum((testAB < dvecAB[1:JKm]))
if(sig > 0)
print("Significant Interaction: Reject")
if(sig == 0)
print("No significant Interaction: Fail to reject")
list(testA=testA,crit.vecA=dvecA,testB=testB,crit.vecB=dvecB,testAB=testAB,crit.vecAB=dvecAB)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.