R/tsplit.R

tsplit <-
function(J,K,data,tr=.2,grp=c(1:p),p=J*K){
#  Perform a J by K anova on trimmed means with
#  repeated measures on the second factor. That is, a split-plot design
#  is assumed, with the first factor consisting of independent groups.
#
#  The  variable data is assumed to contain the raw
#  data stored in list mode or a matrix.
#  If in list mode, data[[1]] contains the data
#  for the first level of both factors: level 1,1.
#  data[[2]] is assumed to contain the data for level 1 of the
#  first factor and level 2 of the second: level 1,2
#  data[[K]] is the data for level 1,K
#  data[[K+1]] is the data for level 2,1, data[2K] is level 2,K, etc.
#
#  The default amount of trimming is tr=.2
#
#  It is assumed that data has length JK, the total number of
#  groups being tested, but a subset of the data can be analyzed
#  using grp
#
if(is.data.frame(data))data=as.matrix(data)
x<-data
       if(is.matrix(x)) {
                y <- list()
                for(j in 1:ncol(x))
                        y[[j]] <- x[, j]
                data <- y
        }
if(!is.list(data))stop("Data are not stored in list mode or a matrix")
if(p!=length(data)){
print("The total number of groups, based on the specified levels, is")
print(p)
print("The number of groups in data is")
print(length(data))
print("Warning: These two values are not equal")
}
if(p!=length(grp)){
print("Apparently a subset of the groups was specified")
print("that does not match the total number of groups ")
stop("indicated by the values for J and K.")
}
tmeans<-0
h<-0
v<-matrix(0,p,p)
klow<-1-K
kup<-0
for (i in 1:p)tmeans[i]<-mean(data[[grp[i]]],tr,na.rm=TRUE)
for (j in 1:J){
h[j]<-length(data[[grp[j]]])-2*floor(tr*length(data[[grp[j]]]))
#    h is the effective sample size for the jth level of factor A
#   Use covmtrim to determine blocks of squared standard errors and
#   covariances.
klow<-klow+K
kup<-kup+K
sel<-c(klow:kup)
v[sel,sel]<-covmtrim(data[grp[klow:kup]],tr)
}
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
#  Do test for factor A
cmat<-kron(cj,ik)  # Contrast matrix for factor A
Qa<-johansp(cmat,tmeans,v,h,J,K)
# Do test for factor B
cmat<-kron(ij,ck)  # Contrast matrix for factor B
Qb<-johansp(cmat,tmeans,v,h,J,K)
# Do test for factor A by B interaction
cmat<-kron(cj,ck)  # Contrast matrix for factor A by B
Qab<-johansp(cmat,tmeans,v,h,J,K)
list(Qa=Qa$teststat,Qa.p.value=Qa$p.value,
Qb=Qb$teststat,Qb.p.value=Qb$p.value,
Qab=Qab$teststat,Qab.p.value=Qab$p.value)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.