Using Generalized Mixed Models In R

Question

I'm trying to run a logistic regression but for a within-subject design so I ended up using a genearlized mixed model (glmer and lme functions in R). There is one predictor (with 3 conditions) and a binary outcome coded as 0 or 1.

I was wondering if anyone had advice on which model is better to use and why.

I'm also having trouble interpreting the output, and the summary seems to be leaving out one of the three conditions in my predictor. When I run the code below:

 D1_data <- read.table("data", header = TRUE, sep = ",") 
 print(D1_data)

 library(nlme)
 model <- lme(output ~ condition, data = D1_data, random = ~1|participant)

 anova(model)
 summary(model)


 log.model <- glmer(output ~ condition + (1|participant), data=D1_data, family=binomial(link=logit))

 summary(log.model)
 anova(log.model)

The summary of only includes the first 2 conditions (in both instances) for example:

 Fixed effects:
                Estimate Std. Error z value Pr(>|z|)
 (Intercept)       40.03     427.04   0.094    0.925
 condition1   -39.63     427.04  -0.093    0.926
 condition2   -37.07     427.04  -0.087    0.931

Even though there should be a 3rd condition. Could someone please help me to find out what the issue is and provide some advice on how to get the 3 conditions to be shown! Thanks in advanced I'm very new to R so please excuse me if I'm asking a simple question.

Answer 1

Interesting thread! I'll chime in for the fun of it.

lme()

When you formulate your lme() model, you are implicitly assuming that your output variable is a continuous random variable. For example, output could be a variable such as body weight, which in principle could take any value within a plausible range of values.

The lme() model you stated would be interpreted with this assumption of a continuous output variable in mind:

library(nlme)
model <- lme(output ~ condition, 
             random = ~1|participant, 
             data = D1_data)

The model will tell you how the expected value of the output variable depends on condition for a typical participant (i.e., a participant with a random intercept equal to 0).

Because R uses dummy coding for the predictor condition, it will encode the effect of condition by including two dummy variables in the model. If your condition predictor has 3 levels, denoted as 0, 1 and 2, R will treat the first of these levels as the reference condition and define the dummy variables that will help you compare the remaining non-reference conditions to the reference one. Specifically, the two dummy variables will be defined as:

condition1 = 1 if condition is equal to 1 and 0 else

and

condition2 = 1 if condition is equal to 1 and 0 else

Note that these dummy variables are NOT physically added to your dataset; they are however used by R when it fits the model as if they were included.

When you interpret the summary of your lme() model, you will see R report the expected value of the output variable for the reference condition (i.e., condition 0) for the typical subject: look for it in the fixed effects portion of your summary, at the intersection of the Intercept row and the Estimate column. R will also report the difference in the expected values of the output variable between condition 1 and condition 0: look for it at the intersection of the condition1 row and the Estimate column. The difference in the expected values of the output variable between condition 2 condition 0 will be found at the intersection of the condition2 row and the Estimate column.

If you parametrize the model like this (i.e., if you force R to fit a model without a fixed intercept):

library(nlme)
model <- lme(output ~ -1 + condition, 
             random = ~1|participant, 
             data = D1_data)

this will prompt R to report the expected values of the output variable for the typical subject for all 3 conditions (rather than the expected value for condition 0 and the differences in expected values between conditions 2 and 3, respectively, and condition 0). In other words, expect to see something like this in your model output:

Fixed effects:
                Estimate            Std. Error z value Pr(>|z|)
 condition0       40.03     
 condition1       40.03 + (-39.63)     
 condition2       40.33 + (-37.07)      

glmer()

If your output variable is binary, you need to formulate your glmer() model using a binomial family and a logit link. For example, output could be a variable such as whether or not the study subject is overweight (where 1 = overweight, 0 = not overweight). A binary variable can only take 2 values: 1 or 0.

This type of glmer() model is called a mixed effects binary logistic regression model and can be fitted like this in R:

library(lme4)

model <- glmer(output ~ condition + (1|participant),                 
             data = D1_data)

The model will investigate how the log-odds that the output variable takes the value 1 rather than 0 for a typical participant (i.e., a participant with a random intercept equal to 0) depends on condition.

Answer 2

Regarding the missing condition: This seems to be a result of dummy/effect coding. That is what happens in the background of the regression functions as these are not able to deal with ordinal data directly. It splits your ordinal predictor (your conditions) into three minus one binary predictors (condition 1: 0, 1; condition 2: 0,1) - one for each condition. The estimate that you get for the two presented conditions is to be interpreted relatively to the third ("missing") condition, which is described here as the intercept.

Regarding model selection: You want to conduct a logistic regression. Therefore, you should use the modeling approach accordingly. In your case, this would be your glmer model. lme is a purely linear model which does not suit your assumptions.