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

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

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

mf1 <- match.call()
m <- match(c("formula", "data", "id"), names(mf1), 0L)
mf1 <- mf1[c(1L, m)]
mf1\$drop.unused.levels <- TRUE
mf1[[1L]] <- quote(stats::model.frame)
mf1 <- eval(mf1, parent.frame())

random1 <- mf1[, "(id)"]
depvar <- colnames(mf)[1]

## check which one is the within subjects factor
if (all(length(table(random1)) == table(mf[,3]))) {
ranvar <- colnames(mf)[3]
fixvar <- colnames(mf)[2]
} else {
ranvar <- colnames(mf)[2]
fixvar <- colnames(mf)[3]
}

K <- length(table(mf[, ranvar]))  ## number of repeated measurements
J <- length(table(mf[, fixvar]))  ## number of levels
p <- J*K
grp <- 1:p

fixsplit <- split(mf[,depvar], mf[,fixvar])
indsplit <- split(mf[,ranvar], mf[,fixvar])
dattemp <- mapply(split, fixsplit, indsplit, SIMPLIFY = FALSE)
data <- do.call(c, dattemp)

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)

result <- list(Qa=Qa\$teststat, A.p.value=as.numeric(Qa\$siglevel), A.df = Qa\$df, Qb=Qb\$teststat, B.p.value=as.numeric(Qb\$siglevel), B.df = Qb\$df,
Qab=Qab\$teststat, AB.p.value=as.numeric(Qab\$siglevel), AB.df = Qab\$df,
call = cl, varnames = c(depvar, fixvar, ranvar))
class(result) <- c("bwtrim")
result
}

tsplit <- bwtrim
```

## 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.