# How to understand the calculations behind a robust mixed factorial function in R, and the implications of running it? Example provided

In "Discovering Statistics Using IBM SPSS Statistics" Andy Field recommends using Robust versions of existing statistics in situations where not all assumptions are met. I'm running a 2x3(bxw) mixed factorial anova where the between groups are n=17 and n=21 and where the group data aren't normal.

In R, I've run bwtrimbt() from the Rallfun-v40.txt file downloadable from Rand Wilcox's website. It's provided some interesting results, but I feel uncomfortable using them without understanding how the calculation differs from a usual mixed factorial ANOVA, and in what senses it is more robust. It uses bootstrap sampling and trimmed means, but you could perform both of these procedures in the context of a usual ANOVA, so there must be more to its robustness.

On investigation, bwtrimbt() performs the bootstrapping, but beyond that relies upon the tsplit() function. The mean trimming is done in the usual way using the standard R mean function with a tr argument).

tsplit(): • works out the effective sample size (the sample size after the trimming) • creates a covariance matrix for the sample trimmed means • creates contrast matrices for the main effects and the interaction

However the main effect and interaction are calculated using a further function: johansp(), and I'm finding it hard to understand the calculations, how they compare to those in a standard mixed factorial ANOVA, and what makes them 'robust.'

Here's some example data to input into it:

source(file.choose()) #then open Rallfun-v40.txt (from Wilcox's website, see above) from your working directory

# Example data:
b1w1 <- c(-11.1, 5.6, 11.1, 16.7, -5.6, 0, -22.2, -5.6, 0, 0, 16.7, 5.6, 55.6, 11.1, 0,  0, 11.1)
b1w2 <- c(0, 5.6, 11.1, 16.7, 5.6, 0, 16.7, 5.6, 11.1, -5.6, 22.2, 11.1, 44.4, 11.1, -22.2, 11.1, 5.6)
b1w3 <- c(-5.6, 5.6, -5.6, 27.8, 5.6, 5.6, 5.6, 22.2, 16.7, -11.1, 5.6, 16.7, -11.1,  44.4, 22.2, 5.6, 22.2)
b2w1 <- c(16.7, 50, 33.33, 22.2, 16.7, 44.4, -16.7, 44.4, 11.1, 44.4, -16.7, 16.7, 16.7, 11.1, 50, 0, 27.8, 16.7, 16.7, 22.2, 16.7)
b2w2 <- c(11.1, 50, 33.3, 0, 11.1, 22.2, 5.6, 0, -5.6, 11.1, -5.6, 0, 16.7, 0, 16.7, 5.6, 11.1, -5.6, 16.7, 44.4, 44.4)
b2w3 <- c(11.1, 11.1, 16.7, 27.8, -16.7, 11.1, 11.1, 5.6, 11.1, -22.2, 5.6, -9.9, 11.1,  11.1,16.7, 27.8, 5.6, 11.1, 27.8, 44.4, 72.2)
complete_data = list(b1w1, b1w2, b1w3, b2w1, b2w2, b2w3)

# Some variables to feed through:
J = 2 # number of between subjects factor levels
K = 3 # number of within subjects factor levels
tr = 0.2 # 20% trimm
p=J*K
grp=c(1:p)
data = complete_data
tmeans<-0
for (i in 1:p)tmeans[i]<-mean(data[[grp[i]]],tr,na.rm=TRUE) # vector of trimmed means

# A covariance matrix for the sample trimmed means
h<-0
v<-matrix(0,p,p)
klow<-1-K  # for indexing, see below
kup<-0 # for indexing, see below
for (j in 1:J){
h[j]<-length(data[[grp[j]]])-2*floor(tr*length(data[[grp[j]]])) # this works out the sample size after the trimming ("the effective sample size")
klow<-klow+K
kup<-kup+K
sel<-c(klow:kup)
v[sel,sel]<-covmtrim(data[grp[klow:kup]],tr)
}
#         [,1]      [,2]       [,3]      [,4]      [,5]      [,6]
#[1,] 6.789829 3.0038021  1.2657005  0.000000  0.000000  0.000000
#[2,] 3.003802 5.6854332  0.6213583  0.000000  0.000000  0.000000
#[3,] 1.265701 0.6213583 17.3762032  0.000000  0.000000  0.000000
#[4,] 0.000000 0.0000000  0.0000000 21.750385  6.485926 -2.788295
#[5,] 0.000000 0.0000000  0.0000000  6.485926 10.386490  3.281377
#[6,] 0.000000 0.0000000  0.0000000 -2.788295  3.281377  9.342241

# Make a contrast matrix for the within subjects main effect
ij<-matrix(c(rep(1,J)),1,J)
km1<-K-1
ck<-diag(1,km1,K)
for (i in 1:km1)ck[i,i+1]<-0-1
cmat<-kron(ij,ck)
#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]    1   -1    0    1   -1    0
#[2,]    0    1   -1    0    1   -1

# We can then feed these into johansp() to get the main effect as follows

johansp(cmat,tmeans,v,h,J,K)
# i.e., the contrast matrix, the trimmed means, the covariance matrix, the effective sample size, the number of between subjects levels, and the number of within subjects levels



johansp() takes in the input as follows:

cmat <- cmat # the contrast matrix
vmean <- tmeans # the trimmed means
vsqse <- v # the covariance matrix
h <- h # a vector containing the effective sample size for each of the between subjects groups
J <- J # the number of between subjects levels
K <- K # the number of within subjects levels



Then it performs the test. Can anybody explain or decipher what the calculations are achieving here? And how they compare to a standard mixed ANOVA? And what makes them more 'robust' (other than the fact we've inputted trimmed means and used bootstrapping)?

Here is how the test runs:


# the trimmed means are put into a single column matrix:
yvec<-matrix(vmean,length(vmean),1)

#          [,1]
#[1,]  3.536364
#[2,]  8.600000
#[3,]  9.618182
#[4,] 21.379231
#[5,]  9.838462
#[6,] 12.400000

# the contrast matrix is multiplied by the covariance matrix:
m1 = cmat%*%vsqse
#         [,1]      [,2]        [,3]      [,4]      [,5]      [,6]
#[1,] 3.786027 -2.681631   0.6443422 15.264459 -3.900563 -6.069671
#[2,] 1.738102  5.064075 -16.7548449  9.274221  7.105113 -6.060864

#Then this is multiplied by a transposition of the the contrast matrix
test <- m1%*%t(cmat)
#          [,1]      [,2]
#[1,] 25.632680 -1.156865
#[2,] -1.156865 34.984897

invc<-solve(test) # this gives us the inverse of the test matrix
#            [,1]        [,2]
#[1,] 0.039071009 0.001291983
#[2,] 0.001291983 0.028626486

m1 <- t(yvec)%*%t(cmat)
#         [,1]     [,2]
#[1,] 6.477133 -3.57972

m2 <- m1%*%invc
#          [,1]        [,2]
#[1,] 0.2484432 -0.09410647

m3 <- m2%*%cmat
#          [,1]       [,2]       [,3]      [,4]       [,5]       [,6]
#[1,] 0.2484432 -0.3425496 0.09410647 0.2484432 -0.3425496 0.09410647

test <- m3%*%yvec
#         [,1]
#[1,] 1.946074



So, to start with, I have no idea what all of that matrix multiplication has actually achieved? What does 1.946074 tell me?

The next bit contains further matrix multiplication, and again I don't understand what it is achieving. It generates a test statistic that can be checked against an F-distribution.


temp<-0
klim<-1-K  # for indexing, see below
kup<-0  # for indexing, see below
for (j in 1:J){
klim<-klim+K
kup<-kup+K
Q<-matrix(0,p,p)
for (k in klim:kup)Q[k,k]<-1
mtem<-vsqse%*%t(cmat)%*%invc%*%cmat%*%Q
temp[j]<-(sum(diag(mtem%*%mtem))+(sum(diag(mtem)))^2)/(h[j]-1)
}

# looking at the matrix multiplication. Imagine that we are on the first iteration of the loop (so j = 1). klim will be 1 (the lower bound for indexing on this iteration), while kup will be 3 (the upper bound for indexing).

Q<-matrix(0,p,p) #  create Q sub j
for (k in klim:kup)Q[k,k]<-1
#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]    1    0    0    0    0    0
#[2,]    0    1    0    0    0    0
#[3,]    0    0    1    0    0    0
#[4,]    0    0    0    0    0    0
#[5,]    0    0    0    0    0    0
#[6,]    0    0    0    0    0    0

m1 <- -vsqse%*%t(cmat)
#            [,1]      [,2]
#[1,]  -3.7860267 -1.738102
#[2,]   2.6816310 -5.064075
#[3,]  -0.6443422 16.754845
#[4,] -15.2644588 -9.274221
#[5,]   3.9005632 -7.105113
#[6,]   6.0696712  6.060864

m2 <- m1%*%invc
#            [,1]        [,2]
#[1,] -0.150169481 -0.05464722
#[2,]  0.098231330 -0.14150205
#[3,] -0.003528127  0.47879985
#[4,] -0.608379936 -0.28520978
#[5,]  0.143219253 -0.19835495
#[6,]  0.244978710  0.18134315

m3 <- m2%*%cmat
#             [,1]        [,2]        [,3]         [,4]        [,5]        [,6]
#[1,] -0.150169481  0.09552226  0.05464722 -0.150169481  0.09552226  0.05464722
#[2,]  0.098231330 -0.23973338  0.14150205  0.098231330 -0.23973338  0.14150205
#[3,] -0.003528127  0.48232798 -0.47879985 -0.003528127  0.48232798 -0.47879985
#[4,] -0.608379936  0.32317016  0.28520978 -0.608379936  0.32317016  0.28520978
#[5,]  0.143219253 -0.34157421  0.19835495  0.143219253 -0.34157421  0.19835495
#[6,]  0.244978710 -0.06363556 -0.18134315  0.244978710 -0.06363556 -0.18134315

mtem<- m3%*%Q
#             [,1]        [,2]        [,3] [,4] [,5] [,6]
#[1,] -0.150169481  0.09552226  0.05464722    0    0    0
#[2,]  0.098231330 -0.23973338  0.14150205    0    0    0
#[3,] -0.003528127  0.48232798 -0.47879985    0    0    0
#[4,] -0.608379936  0.32317016  0.28520978    0    0    0
#[5,]  0.143219253 -0.34157421  0.19835495    0    0    0
#[6,]  0.244978710 -0.06363556 -0.18134315    0    0    0

m1 = mtem%*%mtem
#           [,1]         [,2]        [,3] [,4] [,5] [,6]
#[1,]  0.03174135 -0.010886516 -0.02085483    0    0    0
#[2,] -0.03879991  0.135105766 -0.09630585    0    0    0
#[3,]  0.04959880 -0.346905696  0.29730689    0    0    0
#[4,]  0.12209928  0.001976158 -0.12407544    0    0    0
#[5,] -0.05576027  0.191239509 -0.13547924    0    0    0
#[6,] -0.04239953 -0.048810385  0.09120992    0    0    0

matrix_diagonal = diag(m1)
# [0.03174135, 0.13510577, 0.29730689, 0.00000000, 0.00000000, 0.00000000]

num1 <- sum(matrix_diagonal)  # 0.464154
num2 <- (sum(diag(mtem)))^2 # 0.7546444
den <- h[j]-1 #10 i.e., the effective sample size for the first between subjects group, minus 1
temp[j] <- (num1 + num2)/den  # 0.1218798

# after running the loop to completion, temp = [0.1218798 , 0.2006582]

# it then calculates and returns the following details:

A<-.5*sum(temp) # 0.161269 #for use in calculating cval and the df2
df1<-nrow(cmat) # the number of rows in the contrast matrix, i.e., 2
df2<-nrow(cmat)*(nrow(cmat)+2)/(3*A) # so, 16.53552 ...I don't understand this calculation
cval<-nrow(cmat)+2*A-6*A/(nrow(cmat)+2) #2.080635  ...I don't understand this calculation
test<-test/cval  #the test statistic
#          [,1]
#[1,] 0.9353273
sig<- 1-pf(test,df1,df2) # calculates the significance by checking the test statistic against an F-distribution  # 0.4122603



when called from the tsplit() function, "test" and "sig" (from above) are reported as "Qb" and "Qb.p.value"

I'm trying to work out how the above calculation compares to the standard calculation for a mixed factorial ANOVA. When I started investigating this I was hoping to see analogues to parts of the usual formulae (the various sums of squares etc), but I'm finding the actual calculation behind the robust function very difficult to interpret. Any help appreciated! How is this working? What makes it robust? In what ways does it differ from a standard mixed factorial ANOVA?

Field seems to vastly prefer using robust methods to transforming the data, or performing more straightforward trimming and bootstrapping. I'd like to get a handle on what it is about them that is preferable, and what the implications of using them are.