# R/t3way.R In WRS2: A Collection of Robust Statistical Methods

```t3way <- function(formula, data, tr = 0.2){

if (missing(data)) {
mf <- model.frame(formula)
} else {
mf <- model.frame(formula, data)
}
mcl <- match.call()

J <- nlevels(mf[,2])
K <- nlevels(mf[,3])
L <- nlevels(mf[,4])
p <- J*K*L
grp <- c(1:p)
lev.col <- 2:4
var.col <- 1
alpha <- 0.05

data=as.matrix(mf)

temp=selby2(data,lev.col,var.col)
lev1=length(unique(temp\$grpn[,1]))
lev2=length(unique(temp\$grpn[,2]))
lev3=length(unique(temp\$grpn[,3]))
gv=apply(temp\$grpn,2,rank)
J=lev1
K=lev2
L=lev3
data=temp\$x

if(is.matrix(data))data=listm(data)
data <- lapply(data, as.numeric)
tmeans<-0
h<-0
v<-0

if(length(grp) != p) stop("Incomplete design! It needs to be full factorial!")

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]]]))
v[i]<-(length(data[[grp[i]]])-1)*winvar(data[[grp[i]]],tr)/(h[i]*(h[i]-1))
}
v<-diag(v,p,p)   # Put squared standard errors in a diag matrix.
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)
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
alval<-c(1:999)/1000
#  Do test for factor A
cmat<-kron(cj,kron(ik,il))  # Contrast matrix for factor A
Qa<-johan(cmat,tmeans,v,h,alpha)
A.p.value=t3pval(cmat,tmeans,v,h)
# Do test for factor B
cmat<-kron(ij,kron(ck,il))  # Contrast matrix for factor B
Qb<-johan(cmat,tmeans,v,h,alpha)
B.p.value=t3pval(cmat,tmeans,v,h)
# Do test for factor C
cmat<-kron(ij,kron(ik,cl))  # Contrast matrix for factor C
#Qc<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qc<-johan(cmat,tmeans,v,h,alval[i])
if(Qc\$teststat>Qc\$crit)break
}
C.p.value=irem/1000
# Do test for factor A by B interaction
cmat<-kron(cj,kron(ck,il))  # Contrast matrix for factor A by B
for(i in 1:999){
irem<-i
Qab<-johan(cmat,tmeans,v,h,alval[i])
if(Qab\$teststat>Qab\$crit)break
}
AB.p.value=irem/1000
# Do test for factor A by C interaction
cmat<-kron(cj,kron(ik,cl))  # Contrast matrix for factor A by C
for(i in 1:999){
irem<-i
Qac<-johan(cmat,tmeans,v,h,alval[i])
if(Qac\$teststat>Qac\$crit)break
}
AC.p.value=irem/1000
#Qac<-johan(cmat,tmeans,v,h,alpha)
# Do test for factor B by C interaction
cmat<-kron(ij,kron(ck,cl))  # Contrast matrix for factor B by C
#Qbc<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qbc<-johan(cmat,tmeans,v,h,alval[i])
if(Qbc\$teststat>Qbc\$crit)break
}
BC.p.value=irem/1000
# 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<-johan(cmat,tmeans,v,h,alpha)
for(i in 1:999){
irem<-i
Qabc<-johan(cmat,tmeans,v,h,alval[i])
if(Qabc\$teststat>Qabc\$crit)break
}
ABC.p.value=irem/1000
result <- list(Qa=Qa\$teststat,A.p.value=A.p.value,
Qb=Qb\$teststat,B.p.value=B.p.value,
Qc=Qc\$teststat,C.p.value=C.p.value,
Qab=Qab\$teststat,AB.p.value=AB.p.value,
Qac=Qac\$teststat,AC.p.value=AC.p.value,
Qbc=Qbc\$teststat,BC.p.value=BC.p.value,
Qabc=Qabc\$teststat,ABC.p.value=ABC.p.value,
call = mcl, varnames = colnames(mf))
class(result) <- c("t3way")
result
}
```

## Try the WRS2 package in your browser

Any scripts or data that you put into this service are public.

WRS2 documentation built on May 2, 2019, 4:46 p.m.