# R - Data Simulation with Multiple Random Slopes

I am trying to run the following model:

        I(week^2):mutation_status +
(week + I(week^2) | subject_id) ,
data = sim_dat)


This is the output I get from this model. The correlation between week and I(week)^2 is rather high (-0.95) and I was curious to know how I can change my data simulation code to lower that value. I am also having issues with convergence and would love to get thoughts on how to avoid issues regarding convergence. Thank you!

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]

Formula: fetal_weight ~ week + mutation_status + week:mutation_status + I(week^2) +
I(week^2):mutation_status + (week + I(week^2) | subject_id)
Data: sim_dat

REML criterion at convergence: 114

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.9343 -0.3980 -0.0075  0.4168  1.9928

Random effects:
Groups     Name        Variance  Std.Dev. Corr
subject_id (Intercept)   1.51524  1.2309
week         61.01832  7.8114   0.16
I(week^2)   496.20236 22.2756  -0.19 -0.95
Residual                 0.05892  0.2427
Number of obs: 100, groups:  subject_id, 20

Fixed effects:
Estimate Std. Error       df t value       Pr(>|t|)
(Intercept)                  3.1331     0.4226  17.9961   7.413 0.000000713699 ***
week                        -1.2465     3.5210  18.0053  -0.354          0.727
mutation_statusY             0.5061     0.5977  17.9961   0.847          0.408
I(week^2)                   34.5706    10.8146  18.0084   3.197          0.005 **
week:mutation_statusY       -1.0697     4.9795  18.0053  -0.215          0.832
mutation_statusY:I(week^2) 202.5521    15.2942  18.0084  13.244 0.000000000101 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) week   mttn_Y I(w^2) wk:m_Y
week        -0.160
muttn_sttsY -0.707  0.113
I(week^2)    0.145 -0.966 -0.103
wk:mttn_stY  0.113 -0.707 -0.160  0.683
mtt_Y:I(^2) -0.103  0.683  0.145 -0.707 -0.966


Code for data simulation:

set.seed(123)  # set the seed
J <- 20  # number of individuals (clusters)
cs <- 5  # number of time points (cluster size)
gam <- c(0, 0.75, 1.34)  # fixed effects
G <- matrix(c(1.75, 0, 0,
0, 0.0125, 0,
0,0,0.0625),
nrow = 3)  # random effect variances (G-matrix)
sigma2 <- 0.08

mutation_status <- rep(0:1, each = 5)
mutation_status <- rep(mutation_status, times = 10)

X <- cbind(1, seq_len(cs),
(seq_len(cs))^2)  # for each individual
X <- X[rep(seq_len(cs), J), ]
X <- X[rep(seq_len(cs), J), ]
X[,2] <- X[,2]*0.05
X[,3] <- X[,3]*0.05

# repeat each row cs times
pid <- seq_len(J)  # individual id
pid <- rep(pid, each = cs)

# Generate person-level (lv-2) random effects
uj <- lmf::rmnorm(J, mean = rep(0, 2), varcov = G)

# Generate repeated-measure-level (lv-1) error term
eij <- rnorm(J * cs, sd = sqrt(sigma2))

# Compute beta_j's
betaj <- matrix(gam, nrow = J, ncol = 3, byrow = TRUE) + uj

# Compute outcome:
y <- rowSums(X * betaj[pid, ]) + eij + mutation_status*X[ , 3]*10 + 3

# Form a data frame
sim_dat <- tibble(y, time = X[ , 2], pid, mutation_status)
sim_dat <- sim_dat[,c(3,2,4,1)]
colnames(sim_dat) <- c("subject_id", "week", "mutation_status", "fetal_weight")
sim_dat$$week_2 <-sim_dat$$week*5 + 15

sim_dat$$mutation_status <- as.factor(ifelse(sim_dat$$mutation_status == 0, "N", "Y"))
sim_dat$$subject_id <- as.factor(sim_dat$$subject_id )
#Model Building


First note that your code won't run without the lmf package.

The next thing I see is:

> cor(sim_dat$$week, sim_dat$$week^2)
[1] 0.9811049


So its not surprising that the model has difficulty converging. If you centre the variable it converges without warnings:

> sim_dat$$week0 <- sim_dat$$week - mean(sim_dat\$week)

> m1 <- lmer(fetal_weight ~ week0 + mutation_status + week0:mutation_status + I(week0^2) +
I(week0^2):mutation_status + (week0 + I(week0^2) | subject_id), data = sim_dat)
>  summary(m1)

Random effects:
Groups     Name        Variance  Std.Dev. Corr
subject_id (Intercept)   1.57276  1.2541
week0         2.88079  1.6973  0.11
I(week0^2)  252.60782 15.8936  0.64 0.83
Residual                 0.06426  0.2535
Number of obs: 100, groups:  subject_id, 20

Fixed effects:
Estimate Std. Error t value
(Intercept)                   4.0115     0.4005  10.016
week0                         7.6290     0.7383  10.333
mutation_statusY              4.2372     0.5664   7.481
I(week0^2)                   35.3376     9.9347   3.557
week0:mutation_statusY       61.6479     1.0441  59.042
mutation_statusY:I(week0^2) 190.3985    14.0498  13.552



Now whether these are sensible results or not, I'm not so sure. Personally I always think it's very ambituous to fit random slopes for quadratic terms, but that might be just me !

• Does this answer your question ? If so, please consider marking it as accepted. Thanks ! Jul 19, 2020 at 4:46