2
$\begingroup$

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 
$\endgroup$
3
$\begingroup$

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 !

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Does this answer your question ? If so, please consider marking it as accepted. Thanks ! $\endgroup$ – Robert Long Jul 19 at 4:46

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.