2
$\begingroup$

I struggle to simulate linear mixed model data for the following within-subject within-item factorial design (with crossed participant and item effects) where the categorical independent variable has three levels:

\begin{aligned} &Y_{c,s,i}\ =\ \beta_0\ +\ S_{0,s}\ +\ I_{0,i}\ +\ C_{0,s,i}\ +\ \left(\beta_1\ +\ S_{1,s}\ +\ I_{1,i}\right)X_{1,c}\ +\left(\beta_2\ +\ S_{2,s}\ +\ I_{2,i}\right)X_{2,c}\ +\ \epsilon_{c,s,i},\\ &\begin{pmatrix} S_{0,s}\\ S_{1,s}\\ S_{2,s} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & \rho_{s_{01}}\tau_{00}\tau_{11} & \rho_{s_{02}}\tau_{00}\tau_{22}\\ \rho_{s_{01}}\tau_{00}\tau_{11} & \tau_{11}^2 & \rho_{s_{12}}\tau_{11}\tau_{22}\\ \rho_{s_{02}}\tau_{00}\tau_{22} & \rho_{s_{12}}\tau_{11}\tau_{22} & \tau_{22}^2 \end{pmatrix} \end{bmatrix},\\ &\begin{pmatrix} I_{0,i}\\ I_{1,i}\\ I_{2,i} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \omega_{00}^2 & \rho_{i_{01}}\omega_{00}\omega_{11} & \rho_{i_{02}}\omega_{00}\omega_{22}\\ \rho_{i_{01}}\omega_{00}\omega_{11} & \omega_{11}^2 & \rho_{i_{12}}\omega_{11}\omega_{22}\\ \rho_{i_{02}}\omega_{00}\omega_{22} & \rho_{i_{12}}\omega_{11}\omega_{22} & \omega_{22}^2 \end{pmatrix} \end{bmatrix},\\ &C_{0,s,i}\ \sim\ N \left(0,\ \phi^2\right),\\ &\epsilon_{c,s,i}\ \sim\ N \left(0,\ \sigma^2\right) \end{aligned}

$Y_{c,s,i}$ refers to the dependent variable, the subscripts stand for condition (c), subject (s) and item (i), the by-subject, by-item and by-subject-item-combination intercepts are $S_{0,s}$, $I_{0,i}$ and $C_{0,s,i}$. The corresponding by-subject and by-item slopes are $S_{1,s}$, $S_{2,s}$ and $I_{1,i}$, $I_{2,i}$, respectively. The fixed effects are $\beta_0$, $\beta_1$ and $\beta_2$ and the residual error is $\epsilon_{c,s,i}$.

I hope the mapping between this model formulation and simulata_data() is clear:

  • $\beta_1$ -> effect1
  • $\tau_{00}$ -> subj_interc_sd
  • $\rho_{s_{12}}$ -> subj_cor_s1_s2
  • $\omega_{11}$ -> item_slope1_sd
  • $\phi$ -> subj_item_interc_sd
  • ...

The default parameters I used correspond to the estimates I got from an analysis of previous data.
As I am using the mvrnorm() function with empirical = TRUE I expected that the estimated parameters from lmer and my specified parameters for the simulation would match closely.
However, they differ more than I expected:

  • effect1: should be 12, estimated as 13.771 (Fixed effects table: contrast1 in summary below)
  • item_cor_s2_i: should be 0.3, estimated as 0.11 (Random effects table: Corr(item, contrast2) in summary below))
  • item_cor_s1_s2: should be -0.7, estimated as -0.88 (Random effects table: Corr(contrast1, contrast2) in summary below)

Since this is my first time programming a simulation for mixed models, I'm not sure if the (for me) unexpected results are caused by a programming error or if the default parameters are difficult to estimate for lmer.
I would be grateful if somebody more experienced in mixed model simulations could comment on this.

Running this code ...

library("MASS")
library("lme4")    
set.seed(123)

simulate_data <- function(n_items = 60, n_subj = 27, n_subj_item_rep = 30,
                          mu = 308, 
                          contrast1 = c(2/3, -1/3, -1/3), contrast2 = c(-1/3, 2/3, -1/3), 
                          effect1 = 12, effect2 = -4.5, 
                          subj_interc_sd = 43, 
                          subj_slope1_sd = 17, 
                          subj_slope2_sd = 112, 
                          subj_cor_s1_i = 0.4,
                          subj_cor_s2_i = -0.65,
                          subj_cor_s1_s2 = -0.55,
                          item_interc_sd = 18, 
                          item_slope1_sd = 6, 
                          item_slope2_sd = 6, 
                          item_cor_s1_i = 0.4,
                          item_cor_s2_i = 0.3,
                          item_cor_s1_s2 = -0.7,
                          subj_item_interc_sd = 20,
                          residual_sd = 180) {

  # generate design
  dat <- expand.grid(item = seq_len(n_items), subject = seq_len(n_subj))
  dat$item_group <- dat$item %% 3 + 1    
  dat$subject_group <- dat$subject %% 3 + 1
  dat$condition <- (dat$item_group + dat$subject_group) %% 3 + 1

  # generate appropriate contrasts
  # R requires the generalized inverse of the contrast matrix
  contrast_mat <- ginv(rbind(contrast1, contrast2))
  dat$contrast1 <- contrast_mat[, 1][dat$condition]
  dat$contrast2 <- contrast_mat[, 2][dat$condition]

  # fixed grand mean
  dat$mean_response <- mu

  # condition effects
  dat$contrast_effect1 <- effect1 * dat$contrast1
  dat$contrast_effect2 <- effect2 * dat$contrast2

  # random effects subjects: intercept, slope1, slope2
  subject_effect <- mvrnorm(n_subj, 
                            mu = c(0, 0, 0), 
                            Sigma = matrix(c(subj_interc_sd^2, 
                                             subj_cor_s1_i * subj_slope1_sd * subj_interc_sd,
                                             subj_cor_s2_i * subj_slope2_sd * subj_interc_sd,
                                             subj_cor_s1_i * subj_slope1_sd * subj_interc_sd,
                                             subj_slope1_sd^2,
                                             subj_cor_s1_s2 * subj_slope1_sd * subj_slope2_sd,
                                             subj_cor_s2_i * subj_slope2_sd * subj_interc_sd,
                                             subj_cor_s1_s2 * subj_slope1_sd * subj_slope2_sd,
                                             subj_slope2_sd^2), nrow = 3), 
                            empirical = TRUE)  # mu & sigma as empirical (not population!) parameters for testing

  dat$subject_intercept <- subject_effect[dat$subject, 1] 
  dat$subject_slope1 <- subject_effect[dat$subject, 2] * dat$contrast1
  dat$subject_slope2 <- subject_effect[dat$subject, 3] * dat$contrast2

  # random effects items: intercept, slope1, slope2
  item_effect <- mvrnorm(n_items, 
                         mu = c(0, 0, 0), 
                         Sigma = matrix(c(item_interc_sd^2, 
                                          item_cor_s1_i * item_slope1_sd * item_interc_sd,
                                          item_cor_s2_i * item_slope2_sd * item_interc_sd,
                                          item_cor_s1_i * item_slope1_sd * item_interc_sd,
                                          item_slope1_sd^2,
                                          item_cor_s1_s2 * item_slope1_sd * item_slope2_sd,
                                          item_cor_s2_i * item_slope2_sd * item_interc_sd,
                                          item_cor_s1_s2 * item_slope1_sd * item_slope2_sd,
                                          item_slope2_sd^2), nrow = 3), 
                         empirical = TRUE)

  dat$item_intercept <- item_effect[dat$item, 1] 
  dat$item_slope1 <- item_effect[dat$item, 2] * dat$contrast1
  dat$item_slope2 <- item_effect[dat$item, 3] * dat$contrast2

  # random effects subject-item combinations
  dat$subj_item_intercept <- 0 
  if (n_subj_item_rep > 1) {  # if more than 1 observations of the same subject-item combination
    dat$subj_item_intercept <- mvrnorm(n_subj * n_items, mu = 0, Sigma = subj_item_interc_sd^2, empirical = TRUE)[, 1]
  }

  # generate subject-item replications according to number of observations of the same subject-item combination
  dat <- do.call("rbind", replicate(n_subj_item_rep, dat, simplify = FALSE))

  # residual variation
  dat$residual_sd <- mvrnorm(nrow(dat), mu = 0, Sigma = residual_sd^2, empirical = TRUE)[, 1] 

  # calculate response
  dat$response <- dat$mean_response + dat$contrast_effect1 + dat$contrast_effect2 + 
    dat$item_intercept + dat$item_slope1 + dat$item_slope2 + 
    dat$subject_intercept + dat$subject_slope1 + dat$subject_slope2 + 
    dat$subj_item_intercept +
    dat$residual_sd

  # convert categorical variables to factors
  dat[, 1:5] <- lapply(dat[, 1:5], factor)

  # set specified contrasts
  contrasts(dat$condition) <- contrast_mat

  dat
}

d <- simulate_data()

# if more than 1 observations of the same subject-item combination ..
summary(lmer(response ~ contrast1 + contrast2 + (contrast1 + contrast2|item) +
               (contrast1 + contrast2|subject) + (1|subject:item), d))
# .. otherwise
# summary(lmer(response ~ contrast1 + contrast2 + (contrast1 + contrast2|item) +
#                (contrast1 + contrast2|subject), d))

... I get the following output

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ contrast1 + contrast2 + (contrast1 + contrast2 | item) +  
    (contrast1 + contrast2 | subject) + (1 | subject:item)
   Data: d

REML criterion at convergence: 643678

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1695 -0.6709 -0.0006  0.6728  4.2132 

Random effects:
 Groups       Name        Variance Std.Dev. Corr       
 subject:item (Intercept)   423.44  20.578             
 item         (Intercept)   337.69  18.376             
              contrast1      36.02   6.002   0.38      
              contrast2      66.85   8.176   0.11 -0.88
 subject      (Intercept)  1819.74  42.658             
              contrast1     270.89  16.459   0.45      
              contrast2   12677.59 112.595  -0.65 -0.60
 Residual                 32431.73 180.088             
Number of obs: 48600, groups:  subject:item, 1620; item, 60; subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept)  308.000      8.600   35.82
contrast1     13.771      3.534    3.90
contrast2     -4.708     21.737   -0.22

Correlation of Fixed Effects:
          (Intr) cntrs1
contrast1  0.410       
contrast2 -0.619 -0.557
$\endgroup$
  • $\begingroup$ What is your question? $\endgroup$ – AdamO Mar 20 '18 at 15:18
  • 1
    $\begingroup$ Thank you for providing runnable code. Unfortunately, I don't see the output you refer to in the R output. Which estimates are you unhappy about? That said, have you tried running the simulation with different seeds and looking how much the estimates vary? $\endgroup$ – Stephan Kolassa Mar 21 '18 at 11:23
  • $\begingroup$ @Stephan Kolassa: I hope my edit makes this clear. The estimates do indeed vary. However, I expected that they would match closer with my specified parameters for the simulation. $\endgroup$ – statmerkur Mar 22 '18 at 16:57
2
$\begingroup$

Here is my attempt to simulate the data. I've tried to simplify it.

library(MASS)
library(mvtnorm)
library(dplyr)
simulate_data <- function(n_items = 60, n_subj = 27, n_subj_item_rep = 30,
                          mu = 308, 
                          contrast1 = c(2/3, -1/3, -1/3), contrast2 = c(-1/3, 2/3, -1/3), 
                          effect1 = 12, effect2 = -4.5, 
                          subj_interc_sd = 43, 
                          subj_slope1_sd = 17, 
                          subj_slope2_sd = 112, 
                          subj_cor_s1_i = 0.4,
                          subj_cor_s2_i = -0.65,
                          subj_cor_s1_s2 = -0.55,
                          item_interc_sd = 18, 
                          item_slope1_sd = 6, 
                          item_slope2_sd = 6, 
                          item_cor_s1_i = 0.4,
                          item_cor_s2_i = 0.3,
                          item_cor_s1_s2 = -0.7,
                          subj_item_interc_sd = 20,
                          residual_sd = 180) {
  contrast_mat <- ginv(rbind(contrast1, contrast2))
  m_s <- diag(c(subj_interc_sd, subj_slope1_sd, subj_slope2_sd))
  m_s[lower.tri(m_s)] <- c(subj_cor_s1_i, subj_cor_s2_i, subj_cor_s1_s2)
  m_s[upper.tri(m_s)] <- c(subj_cor_s1_i, subj_cor_s2_i, subj_cor_s1_s2)
  subject_effects <- rmvnorm(n_subj, sigma = sdcor2cov(m_s))
  m_i <- diag(c(item_interc_sd, item_slope1_sd, item_slope2_sd))
  m_i[lower.tri(m_i)] <- c(item_cor_s1_i, item_cor_s2_i, item_cor_s1_s2)
  m_i[upper.tri(m_i)] <- c(item_cor_s1_i, item_cor_s2_i, item_cor_s1_s2)
  item_effects <- rmvnorm(n_items, sigma = sdcor2cov(m_i))

  expand.grid(
    item = seq_len(n_items), 
    subject = seq_len(n_subj),
    replicate = seq_len(n_subj_item_rep)
  ) %>%
    mutate(
      item_group = item %% 3 + 1,
      subject_group = subject %% 3 + 1,
      condition = (item_group + subject_group) %% 3 + 1,
      contrast1 = contrast_mat[, 1][condition],
      contrast2 = contrast_mat[, 2][condition],
      fixed = mu + effect1 * contrast1 + effect2 * contrast2,
      random_subject = subject_effects[subject, 1] + 
        subject_effects[subject, 2] * contrast1 + 
        subject_effects[subject, 3] * contrast2,
      random_item = item_effects[item, 1] + 
        item_effects[item, 2] * contrast1 + 
        item_effects[item, 3] * contrast2,
      item_subject = interaction(item, subject),
      random_interaction = rnorm(n_items * n_subj, sd = subj_item_interc_sd)[item_subject],
      noise = rnorm(n(), sd = residual_sd),
      response = fixed + random_subject + random_item + random_interaction + noise
    ) %>%
    mutate_at(c("item_group", "subject_group", "item", "subject", "condition"), factor) -> dat
  contrasts(dat$condition) <- contrast_mat
  return(dat)
}

Next we can simulate the data several times and model these datasets. In the example below I only extract the fixed effect parameters.

set.seed(123)
replicate(
  10, {
  simulate_data() %>%
    lmer(
      formula = response ~ contrast1 + contrast2 + 
        (contrast1 + contrast2|item) +
        (contrast1 + contrast2|subject) + 
        (1|subject:item)
    ) %>%
    fixef()
}) %>%
  t() -> sim
summary(sim)

The output of the summary (below) indicates quite strong differences in estimates among the simulations. This isn't surprising given that several standard deviations are very high ($\sigma$ = 180, $\tau_{22} = 112$). As a result the signal-to-noise ratio is quite low.

  (Intercept)      contrast1       contrast2      
 Min.   :306.9   Min.   :12.32   Min.   :-25.122  
 1st Qu.:309.4   1st Qu.:13.59   1st Qu.:-19.757  
 Median :310.3   Median :16.05   Median : -7.593  
 Mean   :311.3   Mean   :15.70   Mean   : -6.242  
 3rd Qu.:310.6   3rd Qu.:17.61   3rd Qu.:  2.968  
 Max.   :320.1   Max.   :18.36   Max.   : 18.226  
$\endgroup$
  • $\begingroup$ Thanks for simplifying the code. AFAICS, there is no difference in how you calculate the response variable compared to my approach. So, do you thinks my coding is, apart from being overly complicated, correct? I explicitly used MASS::mvrnorm() with empirical = TRUE to specify empirical not population parameters for the multivariate normal distributions. As a result my estimates vary not as strongly as yours but still quite strongly. I think this is a result of the low signal-to-noise ratio you mentioned. Can you add these two issues to your answer? $\endgroup$ – statmerkur Mar 31 '18 at 11:59

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.