Dear List,
I'm trying to create a logistic linear mixed effect model with 3 fixed
factors, 2 of which are partially crossed (I think this is the right
term?). Factor A has 2 levels, factor B has 2 levels that only occur
within one level of A, and factor C has 2 levels. There are also 2 crossed
random factors, subject and item. I've provided an example data set with
the same structure at the end of this email.
My question is, how do I code the partially crossed fixed effects (A and B)
in order to get all the fixed effects of interest?
Here's what I'd like to get from the model:
1. main effect of A
2. effect of B (which occurs within one level of A)
3. main effect of C
4. A:C interaction
5. B:C interaction (which occurs within one level of A)
So far I have tried the following:
1. Code factor B as two levels, which are only valid within level 1 of A.
Within level 2 of A, factor B = NA. This gives me the following
combinations:
A B
1 1
1 2
2 NA
I get an error using this formula (random intercepts only, for
convenience):
glmer(dv ~ 1 + A + B + C + A:C + B:C +
(1|subject) +
(1|item),
data = mydata,
family = binomial)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
2. Use dummy coding for factor B, where A1B1 is coded as 1 and everything
else is 2, and then look at the effect of B within A1. Recoding the B
factor in this way gives me:
A BnoNA
1 1
1 2
2 2
glmer(dv ~ 1 + A + BnoNA + C + A:C + BnoNA:C +
(1|subject) +
(1|item),
data = mydata,
family = binomial)
This model converges without errors/warnings, but I don't know how to get
the effect of B within A1. Also, it seems like I don't want to model the
main effects and interactions of 'BnoNA' (i.e. the effect of A1B1), but I
don't know how else to write the formula...
3. Remove factor A from the model, recode B with 3 levels (i.e. 1 = A1B1, 2
= A1B2, 3 = A2), then use contrasts of B levels to get the main effect of A
(i.e. (B = 1 and B = 2) vs B = 3).
A B B3levels
1 1 1
1 2 2
2 NA 3
glmer(dv ~ 1 + B3levels + C + B3levels:C +
(1|subject) +
(1|item),
data = mydata,
family = binomial)
This model converges without errors/warnings, but I don't know how to add
the contrasts for (B1 and B2) vs B3, or the interaction between this
contrast coding (i.e. main effect of A) and C.
4. Run 2 separate models. The first model uses the full data set to get the
effects of A, C and the interaction. The second model uses a subset of the
data (only level 1 of A) to get the effect of B within A1, and the B:C
interaction. My concern with this solution is that the random effects
structure would differ between the two models, because the random effects
for subjects and items would be based on data from both levels of A in the
first model and only A1 in the second model.
Apologies for the novice question, messy data simulation code and possible
confusion of terms! I searched the list archives and I googled extensively
but couldn't find an answer.
Thanks in advance.
Becky
### Here is code to produce some example data of the same structure.
A <- c(rep(rep(c(1,1,1,1,0,0), each = 13), 13),
rep(rep(c(0,1,1,1,1,0), each = 13), 13),
rep(rep(c(0,0,1,1,1,1), each = 13), 13),
rep(rep(c(1,0,0,1,1,1), each = 13), 13),
rep(rep(c(1,1,0,0,1,1), each = 13), 13),
rep(rep(c(1,1,1,0,0,1), each = 13), 13))
B <- c(rep(rep(c(1,1,0,0,NA,NA), each = 13), 13),
rep(rep(c(NA,1,1,0,0,NA), each = 13), 13),
rep(rep(c(NA,NA,1,1,0,0), each = 13), 13),
rep(rep(c(0,NA,NA,1,1,0), each = 13), 13),
rep(rep(c(0,0,NA,NA,1,1), each = 13), 13),
rep(rep(c(1,0,0,NA,NA,1), each = 13), 13))
C <- c(rep(rep(c(1,0,1,0,1,0), each = 13), 13),
rep(rep(c(0,1,0,1,0,1), each = 13), 13),
rep(rep(c(1,0,1,0,1,0), each = 13), 13),
rep(rep(c(0,1,0,1,0,1), each = 13), 13),
rep(rep(c(1,0,1,0,1,0), each = 13), 13),
rep(rep(c(0,1,0,1,0,1), each = 13), 13))
mydata <- data.frame(A = factor(A),
B = factor(B),
C = factor(C),
subject = factor(rep(1:78, each = 78)),
item = factor(rep(1:78, 78)),
dv = factor(sample(0:1, 78*78, replace = TRUE)))
# recode B factor - change NAs into 0s
mydata$BnoNA <- mydata$B
mydata$BnoNA[which(is.na(mydata$BnoNA))] <- 0
# recode B factor - change NAs into a 3rd level
mydata$B3levels <- as.numeric(mydata$B)
mydata$B3levels[which(is.na(mydata$B3levels))] <- 3
mydata$B3levels <- as.factor(mydata$B3levels)