2
$\begingroup$

This seems like a relatively basic question, but I can't find good pointers after two days of searching.

I am trying to fit a generalized linear mixed model to data obtained in an experiment. The data are repeated measures on subjects over two conditions (one main effect), and I am interested in measuring that main effect.

The data is the following:

subject <- factor(rep(c(1,2,3,4,5,6,7), each=2))
n <- rep(85, 14)
Ncorrect <- c(67, 70, 65, 80, 70, 76, 67, 62, 65, 71, 63, 75, 58, 71)
Out <- rep(c(0,1), 7) #condition encoded as Integer
OutF <- factor(Out) #condition encoded as logical
expData <- data.frame(subject, n, Ncorrect, Out, OutF)

I want to fit a mixed model to take into consideration the variation of the main effect (slope) between subjects. The problem I have is that I get different models fitted depending on whether I encode experimental conditions as factor vs encoding them as 0,1 integers.

fitInt <- glmer(cbind(Ncorrect, n-Ncorrect) ~ Out + (Out+0|subject),
              data=expData, family=binomial) #Out is an integer

fitFac <- glmer(cbind(Ncorrect, n-Ncorrect) ~ OutF + (OutF+0|subject),
             data=expData, family=binomial) #OutF is factor

First the similarity: they both return almost identical fixed effects (Out and Intercept), which is what I'm interested in.

Then, fitFac shows df=9, while fitInt shows df=11, and they differ in AIC and BIC. In this example fitInt seems more in line with the 3 parameters I want: (Intercept), Out (main effect), and the standard deviation of Out among subjects.

The outputs for the random factors for each fit are the following:

summary(fitInt)
Random effects:
 Groups  Name Variance Std.Dev.
 subject Out  0.1466   0.3829  

summary(fitLog)
Random effects:
 Groups  Name      Variance  Std.Dev.  Corr 
 subject OutLFALSE 9.417e-08 0.0003069      
         OutLTRUE  1.466e-01 0.3828849 -0.96

It seems to me like fitLog is fitting a different random effect for each level of Out, and a correlation between those two? That would explain the Random effects summary and the differences in degrees of freedom.

As asked for in a comment, the Z matrices for both fits are the following:

getME(fitFac,"Z")
14 x 14 sparse Matrix of class "dgCMatrix"
   [[ suppressing 14 column names ‘1’, ‘1’, ‘2’ ... ]]

1  1 . . . . . . . . . . . . .
2  . 1 . . . . . . . . . . . .
3  . . 1 . . . . . . . . . . .
4  . . . 1 . . . . . . . . . .
5  . . . . 1 . . . . . . . . .
6  . . . . . 1 . . . . . . . .
7  . . . . . . 1 . . . . . . .
8  . . . . . . . 1 . . . . . .
9  . . . . . . . . 1 . . . . .
10 . . . . . . . . . 1 . . . .
11 . . . . . . . . . . 1 . . .
12 . . . . . . . . . . . 1 . .
13 . . . . . . . . . . . . 1 .
14 . . . . . . . . . . . . . 1


getME(fitInt,"Z")
14 x 7 sparse Matrix of class "dgCMatrix"
   1 2 3 4 5 6 7
1  . . . . . . .
2  1 . . . . . .
3  . . . . . . .
4  . 1 . . . . .
5  . . . . . . .
6  . . 1 . . . .
7  . . . . . . .
8  . . . 1 . . .
9  . . . . . . .
10 . . . . 1 . .
11 . . . . . . .
12 . . . . . 1 .
13 . . . . . . .
14 . . . . . . 1

Since I want just the random effect of subject on the main effect (Out), matrix 2 seems more in line with that.

The problem seems to be that I want one effect for Out, but when it's coded as a factor (even if binary), an effect for each level seems to be fitted, leading to problems with all the variance falling into one of the levels. So having Out as an integer seems to give me the model I want, but it also seems to be messier according to responses/comments and some previous experience of mine? This looks like a basic issue, maybe I'm missing something obvious I have somehow managed to do without so far.

The main effects I'm interested in don't change between the fits, but I have to apply this analysis to other data too, and I want to understand what's happening before publishing the results.

Thanks for any help.

$\endgroup$
  • $\begingroup$ +1. One of your comments is not quite accurate: a "condition encoded as Integer" would be created with the assignment Out <- rep(factor(c(0,1)), 7). Your question amounts to asking what the difference is. $\endgroup$ – whuber Dec 3 '18 at 20:02
  • 1
    $\begingroup$ Could you use getME(fitInt2,”Z”) and 'getME(fitLog2,”Z”)` to get the Z matrix, so that we can compare and find where the difference come from. $\endgroup$ – user158565 Dec 3 '18 at 23:53
  • $\begingroup$ Edited to include the getME output $\endgroup$ – dherrera Dec 4 '18 at 3:40
3
$\begingroup$

What you observe is due to how the formula interface in R works when you give it a factor, and you exclude the intercept. That is, the formulas ~ Out + 0 and ~ OutF + 0 are not equivalent. To see this, compare the output of the model.matrix() function that is used internally to construct the corresponding design matrices, i.e., compare

 model.matrix(~ OutF + 0, data = DF)

with

model.matrix(~ Out + 0, data = DF)

To make it equivalent you need to replace the second by

model.matrix(~ I(1 - Out) + Out + 0, data = DF)

Hence, compare

glmer(cbind(Ncorrect, n - Ncorrect) ~ Out + (I(1 - Out) + Out + 0 | subject),
      data = expData, family = binomial())

with

glmer(cbind(Ncorrect, n - Ncorrect) ~ OutF + (OutF + 0 | subject),
      data = expData, family = binomial())
$\endgroup$
  • $\begingroup$ Thanks, that makes it clearer. That seems to "make sense" when defining a model like in your examples, because the model with only 1 parameter is unreasonable. But still, in my particular problem, it seems like that way of formula interfacing would only allow me to have 2 random effects, and not just 1 random effect for Out, if it is coded as a factor. In this case, the model with 1 random effect seems reasonable, so it is odd that I wouldn't be able to define it. Is there a workaround to define the desired model (Intercept + Out + random(Out)) with Out as factor? $\endgroup$ – dherrera Dec 4 '18 at 15:09
  • $\begingroup$ It is not clear, to me at least to me, how you want the random-effects structure to look like and why. In any case, remember that random effects are actually modeling correlations. Moreover, in the particular case of mixed effects logistic regression that you have, you need to be aware of the fact that the estimated fixed-effects coefficient you will obtain have an interpretation conditional on the random effects. $\endgroup$ – Dimitris Rizopoulos Dec 4 '18 at 16:10
  • $\begingroup$ Thanks for the bit on interpretation of the parameters, it's something I'm looking into right now. Then, I think the random-effect structure I'm referring to is quite common: fixed intercept, random slope structure. Just that in this case the "slope" is the effect of Out=1. But in all those cases I end up with 2 random effects. As to why that structure: I was divided between random intercept or one intercept per subject (fixed subject effect), because I'm only interested in the "slope", and thought maybe random intercepts could interfere with slope estimation. It's something I'm looking into. $\endgroup$ – dherrera Dec 4 '18 at 18:33
  • $\begingroup$ Indeed this is a feature of the formula interface. You will need to create the dummy variables yourself and include them in the formula. $\endgroup$ – Dimitris Rizopoulos Dec 5 '18 at 8:15
2
$\begingroup$

You are probably used to the classification of statistical variables treated as predictor variables into numerical, discrete or categorical.

An example of numerical predictor variable would be body weight (kg). In R, this would have to be converted to numeric using something like

variable <- as.numeric(variable),

where variable is the predictor variable of interest, located in the R working space.

Count variables are examples of discrete predictor variables and they should be converted to integer variables in R using something like

variable <- as.integer(variable).

Binary variables are also examples of discrete predictor variables.

Categorical predictor variables whose categories are not ordered are handled in R using the concept of unordered factor. The conversion code in R would look something like

variable <- as.factor(variable, ordered = FALSE).

Categorical predictor variables whose categories are ordered are handled in R using the concept of ordered factor. The conversion code in R would look something like

variable <- as.factor(variable, ordered = TRUE).

If your predictor variable is a categorical variable with two categories (be them ordered or unordered), you should follow the rules outlined above and declare it as a factor. Since R is a bit finicky about how it handles ordered categorical predictors (see Is it wrong for `R` to use Polynomial contrasts for ordered categorical variables? and Polynomial contrasts for regression), your best bet is to actually declare them as unordered categorical factors.

In your example, Out is a categorical predictor variable with unordered levels (aka a nominal variable), so just convert it to a factor:

Out <- rep(c(0,1),7)

Out <- factor(Out, levels = c(0,1))

If you wanted to re-label the levels from 0/1 to FALSE/TRUE for your predictor variable, then you could use this code:

Out <- rep(c(0,1),7)

Out <- factor(Out, levels = c(0,1), labels = c("FALSE","TRUE"))

The only time you would use a logical variable - rather than a factor variable - is if that variable were a flag keeping track of some important properties of your data. For example, a flag which keeps track of which record of your data contains at least one missing value. You would NOT use a logical variable as a predictor in your model - only to subset your data according to the value of the flag variable.

R is smart enough to allow you to work with a binary predictor variable (taking the integer values 0 or 1) instead of the factor variable I recommended. But you shouldn't force its hand by throwing logical predictor variables in your model. Why not be consistent with how you code predictor variables in all of your models following the suggestions I provided? This will keep you out of trouble and ensure a consistent workflow.

Addendum

What happens when you try the R code below?

fitFac.one.one <- glmer(cbind(Ncorrect, n-Ncorrect) ~ 
         1 + OutF + (1 + OutF|subject),
         data=expData, family=binomial) #OutF is factor

summary(fitFac.one.one)

The default parametrization involves adding a 1 in front of the factor, rather than a 0, and corresponds to fitFac.one.one above. In this parametrization, the logit of the probability p of a correct response for the typical subject is expressed as

logit(p) = beta0 + beta1*OutF.

In other words, logit(p) = beta0 when OutF = 0 and logit(p) = beta0 + beta1 when OutF = 1, where beta1 is the fixed effect coefficient of OutF in the fitFac.one.one model and beta0 is the fixed effect intercept in that model. Thus, beta1 represents the difference between two logits. The model then allows subjects which are not typical to have logits described by coefficients which vary randomly about beta0 and beta1, respectively. So you will see two standard deviations for the random effects which encapsulate the magnitude of this variation.

In general, adding a 0 in front of a factor in the fixed effects part of your model changes how the effect of the first factor present in the model formula is parametrized (even though there may be other factors in the model). Try this:

fitFac.zero.zero <- glmer(cbind(Ncorrect, n-Ncorrect) ~ 
         0 + OutF + (0 + OutF|subject),
         data=expData, family=binomial) #OutF is factor

summary(fitFac.zero.zero)

In this changed parametrization, the logit of the probability p of a correct response for the typical subject is expressed as

logit(p) = delta0*OutF0 + delta1*OutF1 

where OutF0 is a dummy variable set to 1 when OutF = 0 and set to 0 when OutF = 1. Furthermore, OutF1 is a dummy variable set to 1 when OutF = 1 and set to 0 when OutF = 0. In other words, for a typical subject, logit(p) = delta0 when OutF0 = 1 and OutF1 = 0 (i.e., when OutF = 0), while logit(p) = delta1 when OutF0 = 0 and OutF0 = 1 (i.e., when OutF = 1). Here, delta0 and delta1 are the fixed effect coefficients of OutF0 and OutF1 in the fitFac.zero.zero model. The model then allows subjects which are not typical to have logits described by coefficients which vary randomly about delta0 and delta1, respectively. So you will see two standard deviations for the random effects which encapsulate the magnitude of this variation. But this time, delta1 is no longer a difference in logits between the two levels of OutF!

$\endgroup$
  • 1
    $\begingroup$ Thanks for the answer, I'm still learning good practices, so it's great to know I shouldn't use logicals, I'll change my workflow for that. But my problem with glmer still remains. I edited the question to compare integer vs factor now. The integer encoding gives me the model I want, but the factor doesn't. I want 1 main effect and 1 random effect, but the factor encoding gives me 2 random effects, one for each level. I can't see how I would get the model I want encoding Out as a factor, I'm probably missing something obvious. $\endgroup$ – dherrera Dec 4 '18 at 3:47
  • $\begingroup$ @dherrera: Please see my addendum. This is not necessarily an obvious issue - if anything is unclear, just post a comment here. Maybe others on this forum will chime in to hekp you. $\endgroup$ – Isabella Ghement Dec 4 '18 at 5:07
  • $\begingroup$ I meant the Addendum to my original answer... $\endgroup$ – Isabella Ghement Dec 4 '18 at 5:15
  • 1
    $\begingroup$ Thanks, the addendum did help me understand what's going on, and it's nicely complemented by the other answer. But there is still the issue of a seemingly reasonable model (fixed Intercept + slope, random slope only) that I can't see how could be defined with Out being a factor (better explained in my comment to Dimitris). $\endgroup$ – dherrera Dec 4 '18 at 15:16
  • 1
    $\begingroup$ Thank you very much, it does look like it's the exact same issue. I think this is all the information that my question asked for, and I certainly understand what's happening much better now. $\endgroup$ – dherrera Dec 4 '18 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.