R/trimww.R

trimww <-
function(J,K,x,grp=c(1:p),p=J*K,tr=.2){
#
#  Perform a J by K anova using trimmed means with
#  repeated measures on both factors.
#
#  tr=.2 is default trimming
#
#  The R variable data is assumed to contain the raw
#  data stored 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.
#
#  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.list(x))x<-elimna(matl(x))
if(is.matrix(x))x<-elimna(x)
data<-x
if(is.matrix(data))data<-listm(data)
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 stored in x is")
print(length(data))
print("Warning: These two values are not equal")
}
if(p!=length(grp))stop("Apparently a subset of the groups was specified that does not match the total number of groups indicated by the values for J and K.")
tmeans<-0
h<-length(data[[grp[1]]])
v<-matrix(0,p,p)
for (i in 1:p)tmeans[i]<-mean(data[[grp[i]]],tr=tr,na.rm=TRUE)
v<-covmtrim(data,tr=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)
Qa<-trimww.sub(cmat,tmeans,v,h,J,K)
#Qa.siglevel<-1-pf(Qa$teststat,J-1,999)
Qa.siglevel<-1-pf(Qa,J-1,999)
# Do test for factor B
cmat<-kron(ij,ck)  # Contrast matrix for factor B
#Qb<-johansp(cmat,tmeans,v,h,J,K)
Qb<-trimww.sub(cmat,tmeans,v,h,J,K)
Qb.siglevel<-1-pf(Qb,K-1,999)
# 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)
Qab<-trimww.sub(cmat,tmeans,v,h,J,K)
Qab.siglevel<-1-pf(Qab,(J-1)*(K-1),999)
list(Qa=Qa,Qa.siglevel=Qa.siglevel,
Qb=Qb,Qb.siglevel=Qb.siglevel,
Qab=Qab,Qab.siglevel=Qab.siglevel)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.