2
$\begingroup$

TL:DR version - I am trying to simulate data from a gamma distribution and then fit a Generalized Linear Mixed Model (GLMM) to recover the parameters. The parameter recovery for the fixed effects is agreeable but the random intercept parameter doesn't seem correct.

Specific version

I am simulating some data from a scenario where there are,

  • 40 subjects, who are participating in a longitudinal study with data collected at,
  • 3 timepoints (T0, T1, T2) with,
  • 30 replicates at each timepoint.

In the code below, this is done using a function with a number of different arguments to specify the subjects, repeats, fixed effects values as well as the random intercept SD.

# Loading packages ####

library(tidyverse)
library(lme4)

# Setting up arguments for function 'GenData' ####

subjects <- 40 # number of subjects
repeats  <- 30 # number of replicates per time (3*repeats per subject)
beta_T0   <-  0.2  # Intercept = exp(beta_0)
beta_T1  <-  -0.1  # Effect at Time 1 relative to base-line effect
beta_T2  <-  0.15  # Effect at Time 2 relative to base-line effect 
sigmaI <- 0.2 # between-subject SD

# 'GenData' function ####

GenData <- function(subjects, repeats, beta_0, beta_T1, beta_T2, sigmaI) {

  # Setting seed

  set.seed(248)

  # Creating data frame

  df_data <- expand.grid(I = 1:subjects, r = 1:repeats, T = 0:2) %>% 
    arrange(I, T, r)

  # Random intercept for each individual

  df_I <- tibble(I = unique(df_data$I))
  df_I$deltaI <- rnorm(nrow(df_I), 0, sigmaI)

  # Beta intercepts for Time factor

  df_beta <- tibble(T = c(0,1,2),
                    beta = c(beta_T0, beta_T0 + beta_T1, beta_T0 + beta_T2))

  # Combining it together and taking exponential of the mean

  df_data <- df_data %>% 
    left_join(df_I, by = "I") %>% 
    left_join(df_beta, by = "T") %>%
    mutate(
      ln_mu  = beta + deltaI,
      mu     = exp(ln_mu)
      ) 

  # Calculating the shape parameter using the formula, mean = shape x scale and generating individual values from a gamma distribution

  df_data$shape <- df_data$mu / 0.15 
  df_data$y <- 0.0
  for(i in 1:nrow(df_data)) {
    df_data$y[i] <- rgamma(1,shape = df_data$shape[i], scale = 0.15)
    }

  # Setting as factors

  df_data$I   <- factor(df_data$I)
  df_data$T   <- factor(df_data$T)

  return (df_data)
}

# Generating data ####

df_data <- GenData(subjects, repeats, beta_T0, beta_T1, beta_T2, sigmaI)

# Specifying model and getting model summary ####

fit <- glmer(y ~ 1 + T + (1|I), 
             family = Gamma(link = "log"),
             data = df_data)

summary(fit)

From the output below, one can see that the intercept (0.20283) is close to what it was set (0.2), the beta_T1 (-0.09364) is also close to what it was set (-0.1), and beta_T2 (0.14496) is also close to what it was set (0.15). Thus, all fixed effect parameters are recovered. However, the random intercept SD (0.1002) is almost half of what it was set to (0.2).

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Gamma  ( log )
Formula: y ~ 1 + T + (1 | I)
   Data: df_data

     AIC      BIC   logLik deviance df.resid 
  4099.8   4130.8  -2044.9   4089.8     3595 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3413 -0.7125 -0.0964  0.6089  4.4542 

Random effects:
 Groups   Name        Variance Std.Dev.
 I        (Intercept) 0.01005  0.1002  
 Residual             0.12574  0.3546  
Number of obs: 3600, groups:  I, 40

Fixed effects:
            Estimate Std. Error t value Pr(>|z|)    
(Intercept)  0.20283    0.03594   5.644 1.66e-08 ***
T1          -0.09364    0.01439  -6.509 7.56e-11 ***
T2           0.14496    0.01438  10.078  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
   (Intr) T1    
T1 -0.200       
T2 -0.200  0.500

UPDATE (31/01/2020)

I've done this simulation 1000 times (i.e., set.seed(1) to set.seed(1000)) and below are the histograms and means for the effects.

  • T0 fixed effect (mean from 1000 simulations is 0.1978955; originally set to 0.2)

enter image description here

  • T1 fixed effect (mean from 1000 simulations is -0.100678; originally set to -0.1)

enter image description here

  • T2 fixed effect (mean from 1000 simulations is 0.1497004; originally set to 0.15)

enter image description here

  • Random intercept effect (mean from 1000 simulations is 0.1068291; SD originally set to 0.2)

    UPDATE on 19/02/2020 - Mean value and histogram was that of the variance. This has now been changed to SD for both. Note - this doesn't change the issue on hand.

enter image description here

Am I missing something obvious here?

$\endgroup$
6
$\begingroup$

This is an interesting observation that I am able to reproduce in lme4 (version 1.1-21). I have also implemented the Gamma mixed model in my GLMMadaptive package (version 0.6-9 available currently on GitHub), which seems to recover the correct value for the standard deviation of the random intercepts. The following code illustrates the issue and compares the parameter estimates from the two packages:

simulate_gamma_mixed <- function (n, seed = 1) {
    if (!exists(".Random.seed", envir = .GlobalEnv)) 
        runif(1)
    RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    on.exit(assign(".Random.seed", RNGstate, envir = .GlobalEnv))
    K <- 30 # number of measurements per subject
    # we constuct a data frame with the design: 
    DF <- data.frame(id = rep(seq_len(n), each = K),
                     time = gl(3, 10, length = n * K, labels = paste0("T", 0:2)))
    # design matrix for the fixed effects
    X <- model.matrix(~ time, data = DF)
    betas <- c(0.2, -0.1, 0.15) # fixed effects coefficients
    sigma_b <- 0.2 # standard deviation of random intercepts
    scale <- 0.15 # scale of the Gamma distribution
    # we simulate random effects
    b <- rnorm(n, sd = sigma_b)
    # linear predictor
    eta <- c(X %*% betas + b[DF$id])
    # we simulate Gamma longitudinal data with log link
    DF$y <- rgamma(n * K, shape = exp(eta) / scale, scale = scale)
    DF
}

##########################################################################################
##########################################################################################

library("GLMMadaptive") # you need version >= 0.6-9
library("lme4")
M <- 100 # number of simulations
betas_glmer <- betas_mixmod <- matrix(as.numeric(NA), M, 3)
sigma_b_glmer <- sigma_b_mixmod <- numeric(M)
for (m in seq_len(M)) {
    DF_m <- simulate_gamma_mixed(n = 50, seed = 2020 + m)
    fm_m <- glmer(y ~ time + (1 | id), data = DF_m, 
                  family = Gamma(link = "log"), nAGQ = 15)
    betas_glmer[m, ] <- fixef(fm_m)
    sigma_b_glmer[m] <- sqrt(VarCorr(fm_m)$id[1, 1])
    #####
    gm_m <- mixed_model(y ~ time, random = ~ 1 | id, data = DF_m, 
                        family = Gamma(link = "log"), nAGQ = 15)
    betas_mixmod[m, ] <- fixef(gm_m)
    sigma_b_mixmod[m] <- sqrt(gm_m$D)
}

colMeans(betas_glmer)
mean(sigma_b_glmer)
####
colMeans(betas_mixmod)
mean(sigma_b_mixmod)
$\endgroup$
  • $\begingroup$ @BenBolker could perhaps comment on this. $\endgroup$ – Dimitris Rizopoulos Feb 2 at 20:47
  • 2
    $\begingroup$ sigh, will take a look when I get a chance, have posted this as an issue. Wonder if there might be multiple optima? Should try with glmmTMB too ... github.com/lme4/lme4/issues/557 $\endgroup$ – Ben Bolker Feb 4 at 22:44
  • $\begingroup$ @BenBolker Ran the above code (673 simulations) with glmmTMB too and it seems to be able to recover the parameters (note, this is the mean of all simulations). T0: 0.2012063 (set to 0.2 originally); T1: -0.1012887 (set to -0.1 originally); T2: 0.1499098 (set to 0.15 originally); Random intercept: 0.1960435 (set to 0.2 originally). $\endgroup$ – nahorp Feb 28 at 0:10
  • $\begingroup$ Ran the above code (209 simulations) with brms too and it seems to be able to recover the parameters (note, this is the mean of all simulations). T0: 0.2025091 (set to 0.2 originally); T1: -0.1024149 (set to -0.1 originally); T2: 0.1506764 (set to 0.15 originally); Random intercept: 0.2051822 (set to 0.2 originally) $\endgroup$ – nahorp Mar 1 at 6:22

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.