I want to fit a GLMM model using the glmm() function in R. Its documentation is given here. I'm struggling to understand how to fit it given the way it's described in the documentation. Specifically, I'm unsure how to set up the formula for the random effects (which is to be given in the form of a list).

The model I'm trying to fit is of the form

$$logit(\pi_{it}) = \beta^T x_{it} + u_i^T z_{it}$$

where $\beta = [\beta_1, \beta_2, \beta_3]^T$, $x_{it}=[1,x_{i1},x_{i2}]^T$, $u_i \sim N(0,\Sigma)$ where $\Sigma$ is a diagonal matrix with elements $\sigma_i^2,\sigma_{t1}^2,\sigma_{t2}^2,\sigma_{t3}^2,\sigma_{t4}^2,\sigma_{t5}^2$.

The response variable is binary and $\pi_{it}$ is the probability of success for individual $i$ at time $t$. My corresponding dataframe is of the form

       y x1         x2 t1 t2 t3 t4 t5 i
1  FALSE   0 -0.56166863  1  0  0  0  0 1
2  FALSE   0 -0.83614911  0  1  0  0  0 1
3  FALSE   1 -0.40625873  0  0  1  0  0 1
4  FALSE   1  1.19090616  0  0  0  1  0 1
5  FALSE   1  0.76887964  0  0  0  0  1 1
6  FALSE   1  0.58393470  1  0  0  0  0 2
7  FALSE   0 -0.98137470  0  1  0  0  0 2
8   TRUE   1 -2.19152501  0  0  1  0  0 2
9   TRUE   1 -0.96884116  0  0  0  1  0 2
10  TRUE   1 -0.02810732  0  0  0  0  1 2

Following the documentation for glmm() given above, I tried the following (among others)

  model <- glmm(y ~ x1 + x2, random=list(y~0+i,~0+t1,~0+t2,~0+t3,~0+t4,~0+t5),varcomps.names=c("u","t1","t2","t3","t4","t5"),data=dat,family.glm=bernoulli.glmm)

But the result is clearly incorrect. For example this gives summary

Link is: "logit (log odds)"

Fixed Effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.73635    0.06090 -12.092  < 2e-16 ***
xt1          0.23423    0.08412   2.785  0.00536 ** 
xt2         -0.34128    0.04329  -7.883  3.2e-15 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Variance Components for Random Effects (P-values are one-tailed):
    Estimate Std. Error z value Pr(>|z|)/2
u  1.235e-07  1.751e-07   0.706      0.240
t1 1.181e-02  1.673e-02   0.706      0.240
t2 3.733e-03  5.298e-03   0.705      0.240
t3 2.117e-04  3.012e-04   0.703      0.241
t4 8.407e-03  1.191e-02   0.706      0.240
t5 2.628e-02  3.726e-02   0.705      0.240

I know this is incorrect because the mixed effects were generated from a $N(0,I_5)$ distribution, and so the variance estimates should be close to $1$.

What is the correct way to write the mixed effect formula in this instance?


After messing around I realized I should be using interactions, I think like this

  model <- glmm(y ~ xt1 + xt2, random=list(y~0+i,~0+i:t1,~0+i:t2,~0+i:t3,~0+i:t4,~0+i:t5),varcomps.names=c("subject","t1","t2","t3","t4","t5"),

Can someone more knowledgeable with glmm() please confirm this is correct, or offer some feedback? The results are close enough it could be correct, but also further than I would expect.

  • $\begingroup$ Reading through old questions asked here for glmm() and glmer() I still feel lost. I want to model a random subject intercept, and also random intercepts for subject $i$ at each time $t=1,2,\dots,5$. Every model I construct vastly under estimates the true variance. Any help at all is appreciated. $\endgroup$ – Xiaomi Sep 28 '18 at 15:35
  • $\begingroup$ Do you have a complete working example including the data or how to generate/simulate it? $\endgroup$ – Sextus Empiricus Sep 30 '18 at 8:27

The error from the quote below has two reasons

I know this is incorrect because the mixed effects were generated from a $N(0,I_5)$ distribution, and so the variance estimates should be close to $1$.

1 Data generation

I was going trough your code http://pastebin.com/RjWUaMVE and something seems wrong with this part...

Sig <- diag(1,T+1,T+1)  #Variance covariance matrix


for(row in 1:n){
  for(col in 1:T){
    Ui <- mvrnorm(1,mu=rep(0,T+1),Sigma=Sig)  #Generate random intercepts U_i
    z1 <- beta[1] + beta[2]*Xit1[row,col] + beta[3]*Xit2[row,col] + sum(Z[col,]*Ui)  #The RHS of model equation
    pi <- exp(z1)/(1+exp(z1))  #Predicted probability given covariates
    Y[row,col] <- runif(1) < pi  #Simulate binary response with predicted probability

The term $u^T_{i}z_{it}$ is not implemented correctly in the code.

  1. It is a term that is different for every $i$, but you have it also inside the loop for(col in 1:T) which means that for every $t$ you are recomputing the $u_i$ and you will not get the same $u_i$ for every $z_{it}$ with the same $i$ but different $t$.

  2. You have only a single case for each combination of $i$ and $t$. This makes it very hard to estimate the different $\sigma$ (they will only relate to the heteroskedasticity but if the sigma are not different then there is no heteroskedasticity). When you have multiple measurements for the same $i$ and $t$, then the sigma can be related to the correlation between samples within a single class of $i$ and $t$ (although this will be very difficult for Bernoulli distributed data).

The code below changes this:


n <- 100 # Sample size
T <- 5 # Number of time points
rit <- 10 # number of repetitions within each combination i and t

beta <- c(-1,1,-0.2)  #True beta vector
Sig <- diag(c(1,0.5,0.5,0.5,1,1),T+1,T+1)  #Variance covariance matrix
Z <- matrix(c(1,1,0,0,0,0,
              1,0,0,0,0,1),ncol=6,byrow=TRUE)  #Z Design matrix

Xit1 <- cbind(rbinom(n,1,0.5),rbinom(n,1,0.5),rbinom(n,1,0.5),rbinom(n,1,0.5),rbinom(n,1,0.5))  #Generate Bernoulli covariates
Xit2 <- cbind(rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1))  #Generate normal covariates
Y <- matrix(rep(NA,n*T),ncol=T)  #Matrix to store simulated measurements. Each individual has 5 measurements in 1 row.

for(row in 1:(n/rit)){
  Ui <- mvrnorm(1,mu=rep(0,T+1),Sigma=Sig)  #Generate random intercepts U_i
  for(col in 1:T){
    for (i in 0:rit) {
    z1 <- beta[1] + beta[2]*Xit1[row*rit-i,col] + beta[3]*Xit2[row*rit-i,col] + sum(Z[col,]*Ui)  #The RHS of model equation
    pi <- exp(z1)/(1+exp(z1))  #Predicted probability given covariates
    Y[row*rit-i,col] <- runif(1) < pi  #Simulate binary response with predicted probability

y <- as.vector(t(Y))  #Transforms Y matrix into vector
xit1 <- as.vector(t(Xit1)) #Transforms Xit1 matrix into vector
xit2 <- as.vector(t(Xit2)) #Transforms Xit2 matrix into vector
i <- rep(1:(n/rit),each=rit*T)  #Subject label
t1 <- rep(c(1,0,0,0,0),n)  #Binary variable, = 1 if measurement was at time 1
t2 <- rep(c(0,1,0,0,0),n)  #Binary variable, = 1 if measurement was at time 2
t3 <- rep(c(0,0,1,0,0),n)  #Binary variable, = 1 if measurement was at time 2
t4 <- rep(c(0,0,0,1,0),n)  #Binary variable, = 1 if measurement was at time 2
t5 <- rep(c(0,0,0,0,1),n)  #Binary variable, = 1 if measurement was at time 2
t <- rep(c(1:T),n)   #Stores time in non-binary form

dat <- data.frame(y=y,xt1=xit1,xt2=xit2,t1=t1,t2=t2,t3=t3,t4=t4,t5=t5,i=factor(i),t=factor(t))

notice especially (1) the changes in the position of the call Ui <- mvrnorm... which makes that the term Ui doesn not continually change for the same i (2) and the addition of the variable rit which specifies how often the combination of the same $i$ and $t$ are repeated.

2 The call to the glmm function.

The glmm function seems to be able to only specify random intercept term and the parameter that you use is the group for which you want this to change. For instance when you call


then you have a random intercept per $i$ and per $t$. But you are not looking for a random intercept per time category, you are looking for a random slope / main effect of $t$ per individual $i$.

You function that you mention later


does it indeed correctly. But, the interpretation as interaction might not be so right.

GLMM alternatives.

There are alternatives that are (1) faster (2) more intuitive (3) and more general (ie you can have random slope coefficient besides a random intercept term and you can have correlation/dependency between different slope/intercept terms)

MASS::glmmPQL(y ~ xt1 + xt2, random= list(i = pdDiag(~1+t1+t2+t3+t4+t5)),
             data=dat,family = binomial())

lme4::glmer( y ~ xt1 + xt2 + (1+t1+t2+t3+t4+t5 || i) , 
             data=dat,family = binomial())

You could try the GLMMadaptive package that will fit the model using the adaptive Gaussian quadrature, e.g.,

fm <- mixed_model(y ~ x1 + x2, random = ~ timeF || id, data = dat, family = binomial())

where timeF is the factor variable for the different time points. The use of the || symbol in the random argument specifies that the covariance matrix of the random effects is diagonal.

  • 1
    $\begingroup$ Thanks for your answer, your code does work and I'll be using it if I can't get glmm() to work :) Unfortunately I have been instructed to use glmm() specifically $\endgroup$ – Xiaomi Sep 28 '18 at 5:25

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.