Parameter recovery (Gamma distribution and model with a random intercept)

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.

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_datay[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),
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) • T1 fixed effect (mean from 1000 simulations is -0.100678; originally set to -0.1) • T2 fixed effect (mean from 1000 simulations is 0.1497004; originally set to 0.15) • 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. Am I missing something obvious here?

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)
• @BenBolker could perhaps comment on this. Feb 2 '20 at 20:47
• 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 Feb 4 '20 at 22:44
• @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). Feb 28 '20 at 0:10
• 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) Mar 1 '20 at 6:22