R/bwwtrim.R

bwwtrim <-
function(J,K,L,data,tr=.2,grp=c(1:p),alpha=.05,p=J*K*L){
#  Perform a between by within by within (three-way) anova
#  on trimmed means where
#
#  J independent groups, KL dependent groups
#
#  The variable data is assumed to contain the raw
#  data stored in list mode. data[[1]] contains the data
#  for the first level of all three factors: level 1,1,1.
#  data][2]] is assumed to contain the data for level 1 of the
#  first two factors and level 2 of the third factor: level 1,1,2
#  data[[L]] is the data for level 1,1,L
#  data[[L+1]] is the data for level 1,2,1. data[[2L]] is level 1,2,L.
#  data[[KL+1]] is level 2,1,1, etc.
#
#  The default amount of trimming is tr=.2
#
#  It is assumed that data has length JKL, the total number of
#  groups being tested.
#
if(is.data.frame(data))data=as.matrix(data)
if(is.list(data))data=bwwna(J,K,L,data) # remove missing values
if(is.matrix(data))data=bwwmatna(J,K,L,data) #remove missing values
#                                     and convert to list mode
if(!is.list(data))stop("The 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")
}
tmeans<-0
h<-0
v<-0
for (i in 1:p){
tmeans[i]<-mean(data[[grp[i]]],tr)
h[i]<-length(data[[grp[i]]])-2*floor(tr*length(data[[grp[i]]]))
#    h is the effective sample size
}
v=bwwcovm(J,K,L,data,tr=tr)
ij<-matrix(c(rep(1,J)),1,J)
ik<-matrix(c(rep(1,K)),1,K)
il<-matrix(c(rep(1,L)),1,L)
jm1<-J-1
cj<-diag(1,jm1,J)
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
lm1<-L-1
cl<-diag(1,lm1,L)
for (i in 1:lm1)cl[i,i+1]<-0-1
#  Do test for factor A
cmat<-kron(cj,kron(ik,il))  # Contrast matrix for factor A
Qa=bwwtrim.sub(cmat, tmeans, v, h,p)
Qa.siglevel <- 1 - pf(Qa, J - 1, 999)
# Do test for factor B
cmat<-kron(ij,kron(ck,il))  # Contrast matrix for factor B
Qb=bwwtrim.sub(cmat, tmeans, v, h,p)
 Qb.siglevel <- 1 - pf(Qb, K - 1, 999)
# Do test for factor C
cmat<-kron(ij,kron(ik,cl))  # Contrast matrix for factor C
Qc<-bwwtrim.sub(cmat, tmeans, v, h,p)
Qc.siglevel <- 1 - pf(Qc, L - 1, 999)
# Do test for factor A by B interaction
cmat<-kron(cj,kron(ck,il))  # Contrast matrix for factor A by B
Qab<-bwwtrim.sub(cmat, tmeans, v, h,p)
Qab.siglevel <- 1 - pf(Qab, (J - 1) * (K - 1), 999)
# Do test for factor A by C interaction
cmat<-kron(cj,kron(ik,cl))  # Contrast matrix for factor A by C
Qac<-bwwtrim.sub(cmat, tmeans, v, h,p)
Qac.siglevel <- 1 - pf(Qac, (J - 1) * (L - 1), 999)
# Do test for factor B by C interaction
cmat<-kron(ij,kron(ck,cl))  # Contrast matrix for factor B by C
Qbc<-bwwtrim.sub(cmat, tmeans, v, h,p)
Qbc.siglevel <- 1 - pf(Qbc, (K - 1) * (L - 1), 999)
# Do test for factor A by B by C interaction
cmat<-kron(cj,kron(ck,cl))  # Contrast matrix for factor A by B by C
Qabc<-bwwtrim.sub(cmat, tmeans, v, h,p)
Qabc.siglevel <-1-pf(Qabc,(J-1)*(K-1)*(L-1), 999)
list(Qa=Qa,Qa.p.value=Qa.siglevel,Qb=Qb,Qb.p.value=Qb.siglevel,
Qc=Qc,Qc.p.value=Qc.siglevel,Qab=Qab,Qab.p.value=Qab.siglevel,
Qac=Qac,Qac.p.value=Qac.siglevel,Qbc=Qbc,Qbc.p.value=Qbc.siglevel,
Qabc=Qabc,Qabc.p.value=Qabc.siglevel)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.