Q: VarCorr argument in makeLmer for more compicated random effects

[b1azk0]

I'm trying to get models of the following form to work:

makeLmer(y ~ x1 + x2 + (1|id) + (0 + x1 | id), fixef=b, VarCorr=V2, sigma=s, data=X) makeLmer(y ~ x1 + x2 + (1|id) + (0 + x1 | id) + (0 + x2 | id), fixef=b, VarCorr=V2, sigma=s, data=X) makeLmer(y ~ x1 + x2 + (x1 + x2 |id), fixef=b, VarCorr=V2, sigma=s, data=X) makeLmer(y ~ x1 + x2 + *(x1 x2 |id)**, fixef=b, VarCorr=V2, sigma=s, data=X)

but can't figure out how the VarCorr matrix should look like.

The data I'm playing with is coded as:

id <- 50 obs <- 14 x1<-mvrnorm(id, mu=1, Sigma=0.1, empirical = TRUE)

x2<-mvrnorm(id*obs, mu=1, Sigma=0.1, empirical = TRUE)

X1 <- data.frame(id=1:id, obs=1, x1=x1) X <- extend(X1, along="obs", n=obs) X$x2 <- x2 b = c(0.1, 0.8, 0.3) s=1

V1=1 # works for random intercept V2 = matrix(c(0.5, 0.01, 0.01, 0.5), 2, 2) # works for 1 random effect V3= matrix(c(1, 1))

[pitakakariki]

Short answer: you want respectively a list with 2 elements, a list with three elements, a 3x3 matrix, and a 4x4 matrix.

Longer answer: I'm working on a way to make this simpler to work out. In the mean time I'd recommend the following:

Make some placeholder data for y, fit a model and extract the random effect variances:

X$y <- rnorm(700)
fm <- lmer(y ~ x1 +x2 + (1|id) + (0+x1|id), data=X)
vc <- VarCorr(fm)

Have a look at the structure:

str(vc)
# List of 2
#  $ id  : num [1, 1] 0
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr "(Intercept)"
#   .. ..$ : chr "(Intercept)"
#   ..- attr(*, "stddev")= Named num 0
#   .. ..- attr(*, "names")= chr "(Intercept)"
#   ..- attr(*, "correlation")= num [1, 1] 1
#   .. ..- attr(*, "dimnames")=List of 2
#   .. .. ..$ : chr "(Intercept)"
#   .. .. ..$ : chr "(Intercept)"
#  $ id.1: num [1, 1] 0.00791
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr "x1"
#   .. ..$ : chr "x1"
#   ..- attr(*, "stddev")= Named num 0.0889
#   .. ..- attr(*, "names")= chr "x1"
#   ..- attr(*, "correlation")= num [1, 1] 1
#   .. ..- attr(*, "dimnames")=List of 2
#   .. .. ..$ : chr "x1"
#   .. .. ..$ : chr "x1"
#  - attr(*, "sc")= num 1.03
#  - attr(*, "useSc")= logi TRUE
#  - attr(*, "class")= chr "VarCorr.merMod"

Don't worry about the attributes, just note that you've got a list with two numbers.

V <- list(4, 9)

Now:

makeLmer(y ~ x1 + x2 + (1|id) + (0 + x1 | id), fixef=b, VarCorr=V, sigma=s, data=X)
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x1 + x2 + (1 | id) + (0 + x1 | id)
#    Data: X
# REML criterion at convergence: 2268.812
# Random effects:
#  Groups   Name        Std.Dev.
#  id       (Intercept) 2       
#  id.1     x1          3       
#  Residual             1       
# Number of obs: 700, groups:  id, 50
# Fixed Effects:
# (Intercept)           x1           x2  
#         0.1          0.8          0.3 

Don't forget to double check that the variances went to the correct effects - I haven't implemented proper checks for this yet in makeLmer.

[b1azk0]

Great! Thank you very much

[paulwarren]

Kia ora Pete As part of doing a power analysis for our own prospective study I'm trying to implement a model based on published data (Table 2 of this paper: https://link.springer.com/article/10.3758/s13428-016-0779-0). The table is pasted below. The model includes two random effects with slopes: (congruent | participant) and (congruent | stimulus). The table of model results gives the SDs for intercepts and slopes, which I can convert to Variances for the VarCorr argument in makeLmer, but it also gives correlations and not covariances. I'm not sure how to proceed, as VarCorr requires covariances.

Apologies if this is a dumb statistical question - I am a linguist and not a statistician.

Ngā mihi Paul

[pitakakariki]

Kia ora Paul.

Just multiply all three numbers together to get the covariance.

i.e. sd1 x sd2 x correlation

[ChloeYing2022]

Short answer: you want respectively a list with 2 elements, a list with three elements, a 3x3 matrix, and a 4x4 matrix.

Longer answer: I'm working on a way to make this simpler to work out. In the mean time I'd recommend the following:

Make some placeholder data for y, fit a model and extract the random effect variances:

X$y <- rnorm(700)
fm <- lmer(y ~ x1 +x2 + (1|id) + (0+x1|id), data=X)
vc <- VarCorr(fm)

Have a look at the structure:

str(vc)
# List of 2
#  $ id  : num [1, 1] 0
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr "(Intercept)"
#   .. ..$ : chr "(Intercept)"
#   ..- attr(*, "stddev")= Named num 0
#   .. ..- attr(*, "names")= chr "(Intercept)"
#   ..- attr(*, "correlation")= num [1, 1] 1
#   .. ..- attr(*, "dimnames")=List of 2
#   .. .. ..$ : chr "(Intercept)"
#   .. .. ..$ : chr "(Intercept)"
#  $ id.1: num [1, 1] 0.00791
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : chr "x1"
#   .. ..$ : chr "x1"
#   ..- attr(*, "stddev")= Named num 0.0889
#   .. ..- attr(*, "names")= chr "x1"
#   ..- attr(*, "correlation")= num [1, 1] 1
#   .. ..- attr(*, "dimnames")=List of 2
#   .. .. ..$ : chr "x1"
#   .. .. ..$ : chr "x1"
#  - attr(*, "sc")= num 1.03
#  - attr(*, "useSc")= logi TRUE
#  - attr(*, "class")= chr "VarCorr.merMod"

Don't worry about the attributes, just note that you've got a list with two numbers.

V <- list(4, 9)

Now:

makeLmer(y ~ x1 + x2 + (1|id) + (0 + x1 | id), fixef=b, VarCorr=V, sigma=s, data=X)
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x1 + x2 + (1 | id) + (0 + x1 | id)
#    Data: X
# REML criterion at convergence: 2268.812
# Random effects:
#  Groups   Name        Std.Dev.
#  id       (Intercept) 2       
#  id.1     x1          3       
#  Residual             1       
# Number of obs: 700, groups:  id, 50
# Fixed Effects:
# (Intercept)           x1           x2  
#         0.1          0.8          0.3 

Don't forget to double check that the variances went to the correct effects - I haven't implemented proper checks for this yet in makeLmer.

Hi may I ask how to set up "V <- list(4, 9)" based on the "str(V)" information? I am also wondering whether I can fit a model based on 3 participants and use the random effect structure from this model to generate new model for power simulation? Does the sample size of my pilot data matter in the power analysis? Thanks!

[pitakakariki]

As long as the parameter values you get from your pilot study are reasonable it shouldn't stop you from doing a power analysis. With only three participants you might want to check your estimates against similar studies though.

I suspect that if you're trying to fit complicated random effect structures then three participants won't give you much information.

[ChloeYing2022]

As long as the parameter values you get from your pilot study are reasonable it shouldn't stop you from doing a power analysis. With only three participants you might want to check your estimates against similar studies though.

I suspect that if you're trying to fit complicated random effect structures then three participants won't give you much information.

Hi, If I want to add random slopes for the fixed effects, how should I create the variance-covariance matrix for the model like this: model1 <- makeLmer(y ~ (x1 * x2 + (x1+ x2|subject) + (x1+x2|item)? both X1 and X2 have three levels using sliding difference contrast so it will need a 5x5 matrix? But not sure about how to set it up. Any assistance is much appreciated.

[pitakakariki]

It looks like you'll need a list of two 5x5 matrices, one for subject and one for item. This does seem like a slightly unusual model to me though, with the 9 fixed effects for x1 * x2 but then only five random slopes. I would recommend talking to a statistician to help with this.

Also please note that this issue tracker is here for bug reports, and unfortunately I do not have a lot of time to assist people with their analyses.