BMLR

Misc

• Beyond Multiple Linear Regression, Chapter 8: Introduction to Multilevel Models
• Variable Types
  • Level One: variables measured at the most frequently occurring observational unit
    • i.e. Variables that (for the most part) have different values for each row
    • i.e. Vary for each repeated measure of a subject and vary between subjects
  • Level Two: variables measured on larger observational units
    • i.e. Constant for each repeated measure of a subject but vary between each subject
• Data Description in the Examples
  • na (outcome) - Negative Affect (i.e. Performance Anxiety)
  • large - binary
    • large = 1 means the performance type is a Large Ensemble
    • large = 0 means the performance type is either a Solo or Small Ensemble
  • orch - binary
    • orch = 1 means the instrument is Orchestral
    • orch = 0 means the instrument is either a Keyboard or it’s Vocal
• Two-Stage Model
  • The 2-stage model shows a clearer depiction of how a multi-level model is fit in principle
    • 2-stage model uses Ordinary Least Squares to fit a system of regression equations in stages
    • The multi-level model fits a composite of that system of equations using Maximum Likelihood Estimation
  • Process
    • Level 1 models are fitted for each subject/unit
    • Using the estimated level 1 intercepts and slopes as outcome variables, Level 2 intercept and slope models respectively are fit
    • The errors from the Level 2 models are the random effects
  • Issues
    • Weights every subject the same regardless of the number of repeated observations
    • Responds to missing individual slopes (i.e. subjects never exposed to treatment) by simply dropping those subjects
    • Does not share strength effectively across subjects
  • Specification

    • Level 1

      \[ Y_{ij} = a_i + b_i \text{large}_{ij} + \epsilon_{ij} \]

      • \(b_i\) is the fixed effect for subject/unit \(i\) and observation \(j\) (i.e. each subject has repeated observations)

    • Level 2

      \[ \begin{align} a_i &= \alpha_0 + \alpha_1 \text{orch}_i + u_i \\ b_i &= \beta_0 + \beta_1 \text{orch}_i + v_i \end{align} \]

      • \(u_i\) and \(v_i\) are the random effects for subject/unit \(i\)
      • \(a_i\) is the true (i.e. not estimated since no hat) mean of performance anxiety (outcome) when subject plays solos or small ensembles (large = 0)
      • \(b_i\) is the true mean difference in performance anxiety (outcome) for subject_i between large ensembles and other performance types (large contrast)

Chapter 8: Multilevel Models

Multi-Level

• Process

  • Level 1 and Level 2 equations have been combined through substitution then reduced into a composite model
  • Parameters estimated through Maximum Likelihood Estimation

• Misc

  • Error Distribution
    • Correlation between parameters is accounted for using a multivariate normal distribution estimated through a variance-covariance matrix of the error terms (aka random effects) of each Level (see link for details)
  • Adding Level 1 variables to a model formula should change within-person variability (\(\hat\sigma^2\) and \(\hat\sigma\))
    • Even without adding a Level 1 variable, small changes could occur due to numerical estimation procedures used in likelihood-based parameter estimates
  • Optimization methods
    • Usually very little difference between ML and REML parameter estimates
    • Restricted Maximum Likelihood (REML) is preferable when the number of parameters is large or the primary interest is obtaining estimates of model parameters, either fixed effects or variance components associated with random effects
    • Maximum Likelihood (MLE) should be used if nested fixed effects models are being compared using a likelihood ratio test, although REML is fine for nested models of random effects that have the same fixed effects.
  • P-Values
    • Can’t be calculated because the exact distribution of the test statistics under the null hypothesis (no fixed effect) is unknown, primarily because the exact degrees of freedom is not known
    • Rule of Thumb: t-values (ratios of parameter estimates to estimated standard errors) with absolute value above 2 indicate significant evidence that a particular model parameter is different than 0
    • Packages that do report p-values of fixed effects typically using conservative assumptions, large-sample results, or approximate degrees of freedom for a t-distribution
      • Using {lmerTest} will produce a summary of {lme4} with pvals for coefficients
    • Parametric Bootstrap can be used to approximate the distribution of the likelihood test statistic and produce more accurate p-values by simulating data under the null hypothesis
  • Model Comparison
    • Pseudo-R²
      • See below in Simple to Complex Model Building Workflow >>
        • Random Slopes and Intercepts Model (with 1 covariate)
        • Random Slopes and Intercepts Model (with 2 covariates)
    • AIC, BIC
      • See Random Intercepts (with 2 covariates)
    • Likelihood Ratio Test (LR-Test)
      • See Random Slopes and Intercepts Model (with 3 covariates)

• Composite Specification

  \[ Y_{ij} = [\alpha_0 + \alpha_1 \text{orch}_i + \beta_0 \text{large}_{ij} + \beta_1 \text{orch}_i \text{large}_{ij}] + [u_i + v_i \text{large}_{ij} + \epsilon_{ij}] \]

  • Composite model of level 1 and level 2 equations
  • Random Slopes and Intercepts Model (with 2 covariates)
  • Fixed Effects: \(\alpha_0\), \(\alpha_1\), \(\beta_0\) and \(\beta_1\)
  • \(u_i + v_i \cdot \text{large}_{ij} + \epsilon_{ij}\) is the interesting part.
    • This part has all the error terms and any variables in the Level 2 equations
    • The first part is a typical regression with interaction terms.

• Interpretation

  • See Two-stage model for definitions for \(a_i\) and \(b_i\)

    • Keyboardists and Vocalists (orch_i = 0)

      \[ \begin{align} a_i &= \alpha_0 + u_i \\ b_i &= \beta_0 + v_i \end{align} \]

    • Orchestral Instrumentalists (orch_i = 1)

      \[ \begin{align} a_i = (\alpha_0 + \alpha_1) + u_i \\ b_i = (\beta_0 + \beta_1) + v_i \end{align} \]

  • Mean Performance Anxiety (outcome) when:
    • Keyboardists or Vocalists (orch = 0) play solos or small ensembles (large = 0): \(\alpha_0\)
    • Keyboardists or vocalists (orch = 0) play large ensembles (large = 1): \(\alpha_0 + \beta_0\)
    • Orchestral instrumentalists (orch = 1) play solos or small ensembles (large = 0): \(\alpha_0 + \alpha_1\alpha_0 + \alpha_1\)
    • Orchestral instrumentalists (orch = 1) play large ensembles (large = 1): \(\alpha_0 + \alpha_1 + \beta_0 + \beta_1\)

• Model Summary (using lmer4::lmer())

  #>     Linear mixed model fit by REML ['lmerMod'] 
  #> A)  Formula: na ~ orch + large + orch:large + (large | id) 
  #>         Data: music 
  #> B)  REML criterion at convergence: 2987 
  
  #> B2)      AIC      BIC  logLik deviance df.resid 
  #>         3007    3041    -1496    2991      489 
  
  #>     Random effects: 
  #>       Groups  Name        Variance Std.Dev. Corr 
  #> C)  id      (Intercept)    5.655   2.378         
  #> D)            large        0.452   0.672   -0.63 
  #> E)  Residual              21.807   4.670         
  #> F)  Number of obs: 497, groups:  id, 37 
  
  #>     Fixed effects: 
  #>                 Estimate Std. Error t value 
  #> G)  (Intercept)  15.930      0.641  24.83 
  #> H)  orch          1.693      0.945   1.79 
  #> I)  large        -0.911      0.845  -1.08 
  #> J)  orch:large   -1.424      1.099  -1.30

  • Definitions
    • A: How our multilevel model is written in R, based on the composite model formulation.
    • B: Measures of model performance. Since this model was fit using REML, this line only contains the REML criterion.
    • B2: If the model is fit with ML instead of REML, the measures of performance will contain AIC, BIC, deviance, and the log-likelihood.
    • C: Estimated variance components (\(\hat \sigma^2_u\) and \(\hat \sigma_u\)) associated with the intercept equation in Level Two. (between-unit variability)
    • D: Estimated variance components (\(\hat \sigma^2_v\) and \(\hat \sigma_v\)) associated with the large ensemble (large = 1) effect equation in Level Two. (Also between-unit variability but just for the slope)
      • Also, in the “Corr” column, the estimated correlation (\(\hat \rho_{uv}\)) between the two Level Two error terms.
    • E: Estimated variance components (\(\hat \sigma^2\) and \(\hat \sigma\) associated with the Level One equation. (within-unit variability)
    • F: Total number of performances where data was collected (Level One observations = 497) and total number of subjects (Level Two observations = 37).
    • G: Estimated fixed effect (\(\hat \alpha_0\)) for the intercept term, along with its standard error and t-value (which is the ratio of the estimated coefficient to its standard error).
    • H: Estimated fixed effect (\(\hat \alpha_1\)) for the orchestral instrument (orch = 1) effect, along with its standard error and t-value.
    • I: Estimated fixed effect (\(\hat \beta_0\)) for the large ensemble (large = 1) effect, along with its standard error and t-value.
    • J: Estimated fixed effect (\(\hat \beta_1\)) for the interaction between orchestral instruments (orch = 1) and large ensembles (large = 1), along with its standard error and t-value.
  • Interpretations
    • Fixed effects:
      • \(\hat \alpha_0 = 15.9\) — The estimated mean performance anxiety (outcome) for solos and small ensembles (large = 0) for keyboard players and vocalists (orch = 0) is 15.9.
      • \(\hat \alpha_1 = 1.7\) — Orchestral instrumentalists (orch = 1) have an estimated mean performance anxiety (outcome) for solos and small ensembles (large = 0) which is 1.7 points higher than keyboard players and vocalists (orch = 0).
      • \(\hat \beta_0 = −0.9\) — Keyboard players and vocalists (orch = 0) have an estimated mean decrease in performance anxiety (outcome) of 0.9 points when playing in large ensembles (large = 1) instead of solos or small ensembles (large = 0).
      • \(\hat \beta_1 = −1.4\) — Orchestral instrumentalists (orch = 1) have an estimated mean decrease in performance anxiety (outcome) of 2.3 points when playing in large ensembles (large = 1) instead of solos or small ensembles (large = 0), 1.4 points greater than the mean decrease among keyboard players and vocalists (orch = 0).
        • He’s calculating simple slope/marginal effect (See Regression, Interactions) and choosing to use that as the interpretation for the interaction term
          • From Ch. 1 (link): “We interpret the coefficient for the interaction term by comparing slopes under fast and non-fast conditions; this produces a much more understandable interpretation for a reader than attempting to interpret the -0.011 (interaction coef) directly”
            • Fast and non-fast are levels of a main effect and a term in the interaction.
        • Calculation: The interaction is the difference in slopes for one interaction variable (e.g. large) at different values of the other interaction variable (e.g. orch)
          • Regression eq (ignoring the random effects part):
            • \(\hat Y_i = 15.93 + 1.693\text{orch}_i - 0.911\text{large} − 1.424\text{orch}_i \times \text{large}_i\)
          • Examining the large slope so only dealing with terms that contain “large”
            • large = 1 is reference category in the interaction that the coefficient is associated with so it remains constant
          • When orch = 0 and large = 1, the slope for large is -0.911 = \(\hat \beta_0\)
            • \(- 0.911\text{large} − 1.424\text{orch}_i \times \text{large}_i = -0.911 \times 1 - 1.424 \times \boldsymbol{0} \times 1 = -0.911 = \hat \beta_0\)
          • When orch = 1 and large = 1, the slope for large is -2.335 (uses this one for his interpretation)
            • \(- 0.911\text{large} − 1.424\text{orch}_i \times \text{large}_i = -0.911 \times 1 - 1.424 \times \boldsymbol{1} \times 1 = -0.911 - 1.424 = -2.335\) (i.e. “decrease… of 2.3 points”)
          • The difference in these slopes is the interaction coefficient: \(-2.335 - (-0.911) = -1.424\)
    • Variance components
      • \(\hat \sigma_u = 2.4\) — The estimated standard deviation of performance anxiety (outcome) for solos and small ensembles (large = 0) is 2.4 points, after controlling for instrument (orch) played.
      • \(\hat \sigma_v = 0.7\) — The estimated standard deviation of differences in performance anxiety (outcome) between large ensembles (large = 1) and other performance types (large = 0) is 0.7 points, after controlling for instrument (orch) played.
      • \(\hat \rho_{uv} = −0.64\) — The estimated correlation between performance anxiety scores (outcome) for solos and small ensembles (large = 0) and increases in performance anxiety (outcome) for large ensembles (large = 1) is -0.64, after controlling for instrument (orch) played.
        • Those subjects with higher performance anxiety scores for solos and small ensembles tend to have greater decreases in performance anxiety for large ensemble performances.
      • \(\hat \sigma = 4.7\) — The estimated standard deviation in residuals for the individual regression models is 4.7 points.

Model Building Workflow: Simple to Complex

• Workflow should include:
  • Fixed effects that allow one to address primary research questions
  • Fixed effects that control for important covariates at all levels
  • Investigation of potential interactions
  • Variables that are centered where interpretations can be enhanced
  • Important variance components
  • Removal of unnecessary terms
  • A final model that tells a “persuasive story parsimoniously”
• Model Comparison (Also see multi-level >> misc >> p-values & model comparison)
  • LR-Test for nested models
  • AIC, BIC for unnested models
  • Consider models that involve removing terms that have t-stats < |2|
• Models
  • Unconditional Means (aka Random Intercepts)
  • Random Slopes and Intercepts Model (with 1 covariate)
  • Random Slopes and Intercepts Model (with 2 covariates)
  • Random Intercepts (with 2 covariates)
  • Random Slopes and Intercepts Model (with 3 covariates)
  • Final Model

Unconditional Means (aka Random Intercepts)

• Specification

  • Level 1: \(Y_{ij} = a_i + \epsilon_{ij} \quad \text{where}\: \epsilon \sim \mathcal{N}(0, \sigma^2)\)
    • \(a_i\) is the true mean response of all observations for subject_i
  • Level 2: \(a_i = \alpha_0 + u_i \quad \text{where} \: u_i \sim \mathcal{N}(0, \sigma^2_u)\)
    • Each subject’s intercept is assumed to be a random value from a normal distribution centered at \(\alpha_0\) with variance \(\sigma^2_u\).
  • Composite: \(Y_{ij} = \alpha_0 + u_i + \epsilon_{ij}\)

• Model

  #Model A (Unconditional means model)
  model.a <- lmer(na ~ 1 + (1 | id), REML = T, data = music)
  ##  Groups  Name        Variance Std.Dev.
  ##  id      (Intercept)  4.95    2.22   
  ##  Residual            22.46    4.74
  ##  Number of Level Two groups =  37
  ##            Estimate Std. Error t value
  ## (Intercept)    16.24    0.4279  37.94

• Interpretation

  • \(\hat \alpha_0 = 16.2\) — The estimated mean performance anxiety score across all performances and all subjects.

  • \(\hat \sigma^2 = 22.5\) — The estimated variance in within-person deviations.

  • \(\hat \sigma^2_u = 5.0\) (rounded from 4.95) — The estimated variance in between-person deviations.

  • ICC

    \[ \hat \rho = \frac{\text{between-person variability}}{\text{total variability}} = \frac{\hat \sigma^2_u}{\hat \sigma^2_u + \hat \sigma ^2} = \frac{5.0}{5.0 + 22.5} = 0.182 \]

    • 18.2% of the total variability in performance anxiety scores are attributable to differences among subjects.
      • *For plain random intercepts models only*, you can also say this same number says that the average correlation for any pair of responses from the same individual is a moderately low 0.182.
    • As \(\rho\) approaches 0: responses from an individual are essentially independent and accounting for the multilevel structure of the data becomes less crucial.
    • As \(\rho\) approaches 1: repeated observations from the same individual essentially provide no additional information and accounting for the multilevel structure becomes very important.
    • Effective Sample Size (ESS): The number of independent pieces of information we have for modeling
      • With \(\rho\) near 0: ESS approaches the total number of observations.
      • With \(\rho\) near 1: ESS approaches the number of subjects in the study.

Random Slopes and Intercepts Model (with 1 covariate)

• 1 - Level 1

• Specification

  • Level 1: \(Y_{ij} = a_i + b_i \text{large}_{ij} + \epsilon_{ij}\)

  • Level 2

    \[ a_i = \alpha_0 + u_i \\ b_i = \beta_0 +v_i \]

  • Composite

    \[ \begin{aligned} &Y_{ij} = [\alpha_0 + \beta_0 \text{large}_{ij}] + [u_i + v_i \text{large}_{ij} + \epsilon_{ij}] \\ &\begin{aligned} \text{where}\quad \epsilon &\sim \mathcal{N}(0, \sigma^2) \: \text{and}\\ \left[ \begin{array}{cc} u_i \\ v_i \end{array} \right] &\sim \mathcal{N} \left(\left[\begin{array}{cc} 0\\0 \end{array}\right], \left[\begin{array}{cc} \sigma^2_u\\\rho\sigma_u \sigma_v \quad \sigma^2_v \end{array}\right]\right) \end{aligned} \end{aligned} \]

• Model

  model.b <- lmer(na ~ large + (large | id), data = music)
  ##  Groups  Name        Variance Std.Dev. Corr 
  ##  id      (Intercept)  6.333  2.517         
  ##          large        0.743  0.862    -0.76
  ##  Residual            21.771  4.666
  ##  Number of Level Two groups =  37
  ##            Estimate Std. Error t value
  ## (Intercept)  16.730    0.4908  34.09
  ## large        -1.676    0.5425  -3.09

• Interpretation

  • \(\hat \alpha_0 = 16.7\) — The mean performance anxiety level (outcome) before solos and small ensemble performances (large = 0).
  • \(\hat \beta_0 = −1.7\) — The mean decrease in performance anxiety (outcome) before large ensemble performances (large = 1).
    • Subjects had a performance anxiety level of 16.7 before solos and small ensembles, and their anxiety levels were 1.7 points lower, on average, before large ensembles, producing an average performance anxiety level before large ensembles of 15.0
    • Statistically significant since |t-value| > 2
  • \(\hat \sigma^2 = 21.8\) — The variance in within-person deviations.
  • \(\hat \sigma^2_u = 6.3\) — The variance in between-person deviations in performance anxiety scores (outcome) before solos and small ensembles (large = 0).
  • \(\hat \sigma^2_v = 0.7\) — The variance in between-person deviations in increases (or decreases) (slope) in performance anxiety scores (outcome) before large ensembles (large = 1).
  • \(\hat \rho_{uv} = −0.76\) (Corr column)
    • Slopes and intercepts are negatively correlated
    • A strong negative relationship between a subject’s performance anxiety (outcome) before solos and small ensembles (large = 0) (Intercept) and their (typical) decrease in performance anxiety (Slope) before large ensembles (large = 1) .

• Pseudo-R² (Not always a reliable performance measure)

  • Unconditional Means vs Random Intercepts and Slopes

    \[ \text{Psuedo R}^2_{L1} = \frac{\hat \sigma^2 (\text{Model A}) - \hat\sigma^2(\text{Model B})}{\hat \sigma^2(\text{Model A})} = \frac{22.5 - 21.8}{22.5} = 0.031 \]

    • Where Model A: Unconditional Means; Model B: Random intercepts and slopes
    • FYI
      • \(\text{percent decrease} = \frac{\text{original value} - \text{new value}}{\text{orginal value}}\)
        • pseudo-R² uses this one
        • If negative, it’s a percent increase.
      • \(\text{percent increase} = \frac{\text{new value} - \text{original value}}{\text{original value}}\)
        • If negative, it’s a percent decrease.

  • Positive value means Model B is an improvement over Model A

    • Means the variance in Model B’s error terms is smaller than Model A’s, and therefore better fits the data. (I think)

  • The estimated within-person variance \(\hat \sigma^2\) decreased by 3.1% (from 22.5 to 21.8) from the unconditional means model

  • Implies that only 3.1% of within-person variability in performance anxiety scores can be explained by performance type

  • Values of \(\hat \sigma^2_u\) and \(\hat \sigma^2_v\) from Model B cannot be compared to between-person variability from Model A using pseudo-R², since the inclusion of performance type has changed the interpretation of these values

  • Issues

    • “Because of the complexity of estimating fixed effects and variance components at various levels of a multilevel model, it is not unusual to encounter situations in which covariates in a Level Two equation for the intercept (for example) remain constant (while other aspects of the model change), yet the associated pseudo R-squared values differ or are negative.” (?)

Random Slopes and Intercepts Model (with 2 covariates)

• 1 - Level 1; 1 - Level 2

• Specification, Model, Interpretation (see Multi-Level section above)

  • Says the large effect (slope) varies across subjects/units

• Pseudo-R²

  \[ \begin{align} \text{Psuedo R}^2_{L2_u} &= \frac{\hat \sigma^2_u (\text{Model B}) - \hat\sigma^2_u(\text{Model C})}{\hat \sigma^2_u(\text{Model B})} = \frac{6.33 - 5.66}{6.33} = 0.106 \\ \text{Psuedo R}^2_{L2_v} &= \frac{\hat \sigma^2_v (\text{Model A}) - \hat\sigma^2_v(\text{Model B})}{\hat \sigma^2_v(\text{Model A})} = \frac{0.74 - 0.45}{0.74} = 0.392 \end{align} \]

  • Model B is Random Slopes and Intercepts Model (with 1 covariate)

  • Model C is Random Slopes and Intercepts Model (with 2 covariates)

  • The addition of Level 2 variable, orch, in Model C:

    • (Top) Decreased the between-person variability in mean (intercept) performance anxiety (outcome) before solos and small ensembles (large = 0) by 10.6%
    • (Bottom) Decreased the between-person variability in the effect (slope) of large ensembles (large = 1) on performance anxiety (outcome) by 39.2%

Random Intercepts (with 2 covariates)

• 1 - Level 1; 1 - Level 2

• Instead of assuming that the large ensemble (large) effects, after accounting for instrument played (orch), vary by individual, we are assuming that large ensemble effect is fixed across subjects.

• Specification

  • Level 1: \(Y_{ij} = a_i + b_i \text{large}_{ij} + \epsilon_{ij}\)

  • Level 2

    \[ \begin{align} &a_i = \alpha_0 + \alpha_1 \text{orch}_i + u_i\\ &b_i = \beta_0 + \beta_1 \text{orch}_i \\ &\text{where}\: \epsilon \sim \mathcal{N}(0, \sigma^2)\: \text{and} \: u_i \sim \mathcal{N}(0, \sigma^2_u) \end{align} \]

    • Only difference is the random effect/error term, \(v_i\), isn’t specified

  • Composite

    • Not given, but probably very similar to the previous model just without \(v_i\)

• Model

  model.c2 <- lmer(na ~ orch + large + orch:large + (1|id), data = music)
  ##  Groups  Name        Variance Std.Dev.
  ##  id      (Intercept)  5.13    2.27   
  ##  Residual            21.88    4.68
  ##  Number of Level Two groups =  37
  ##            Estimate Std. Error t value
  ## (Intercept)  15.9026    0.6187  25.703
  ## orch          1.7100    0.9131   1.873
  ## large        -0.8918    0.8415  -1.060
  ## orch:large   -1.4650    1.0880  -1.347

• Interpretation

  • Same as previous model except there’s no between-units variation estimate for the slope (see previous model and multi-level model section)
  • Estimates are similar to the previous model’s. Largest difference seems to be the between-unit variation for the Intercept (5.66 vs 5.13)

• Comparison

          df  AIC
  model.c  8 3003
  model.c2 6 2999
          df  BIC
  model.c  8 3037
  model.c2 6 3025

  • Model C is the previous model and Model C2 is this model
  • Both criterion pick this model to predict best
  • The more complex model (e.g. varying intercepts and varying slopes) doesn’t always perform best

Random Slopes and Intercepts Model (with 3 covariates)

• 1 - Level 1; 2 - Level 2

• MPQnem has been centered (mean = 31.63)

  • Otherwise some interpretations would involve MPQnem = 0 when the minimum score is 11 which would make the interpretation nonsensical

• Specification

  • Level 1: \(Y_{ij} = a_i + b_i \text{large}_{ij} + \epsilon_{ij}\)

  • Level 2

    \[ \begin{align} a_i = \alpha_0 + \alpha_1 \text{orch}_i + \alpha_2 \text{MPQnem}_i + u_i \\ b_i = \beta_0 + \beta_1 \text{orch}_i + \beta_2 \text{MPQnem}_i + v_i \end{align} \]

  • Composite

    \[ Y_{ij} = [\alpha_0 + \alpha_1 \text{orch}_i + \alpha_2 \text{mpqnem}_i + \beta_0 \text{large}_{ij} + \beta_1 \text{orch}_i \text{large}_{ij} + \beta_2 \text{mpqnem}_i \text{large}_{ij}] + [u_i + v_i \text{large}_{ij} + \epsilon_{ij}] \]

• Model

  model.d <- lmer(na ~ orch + mpqnem + large +
                      orch:large + mpqnem:large +
                      (large | id),
                      data = music)
  ##  Groups  Name        Variance Std.Dev. Corr 
  ##  id      (Intercept)  3.286  1.813         
  ##          large        0.557  0.746    -0.38
  ##  Residual            21.811  4.670
  ##  Number of Level Two groups =  37
  ##               Estimate Std. Error t value
  ## (Intercept)   16.25679    0.54756 29.6893
  ## orch           1.00069    0.81713  1.2246
  ## cmpqnem        0.14823    0.03808  3.8925
  ## large         -1.23484    0.84320 -1.4645
  ## orch:large    -0.94927    1.10620 -0.8581
  ## cmpqnem:large -0.03018    0.05246 -0.5753

• Interpretation

  • Compared to Random Slopes and Intercepts Model (with 2 covariates):
    • Directions of the effects of instrument (orch) and performance type (large) are consistent
    • Effect sizes and levels of significance are reduced because of the relative importance of the negative emotionality (mpqnem) term
  • \(\alpha_0 = 16.26\) — The estimated mean performance anxiety (outcome) for solos and small ensembles (large = 0) is 16.26 for keyboard players and vocalists (orch=0) with an average level of negative emotionality at baseline (mpqnem = 31.63)
  • \(\hat \alpha_1 = 1.00\) — Orchestral instrument (orch = 1) players have an estimated mean performance anxiety (outcome) level before solos and small ensembles (large = 0) which is 1.00 point higher than keyboardists and vocalists (orch = 0), controlling for the effects of baseline negative emotionality.
  • \(\hat \sigma^2 = 0.15\) — A one point increase in baseline negative emotionality (mpqnem = 0) is associated with an estimated 0.15 mean increase in performance anxiety (outcome) levels before solos and small ensembles (large = 0), after controlling for instrument (orch).
  • \(\hat \beta_0 = −1.23\) — Keyboard players and vocalists (orch = 0) with an average level of baseline negative emotionality levels (mpqnem = 31.63) have an estimated mean decrease in performance anxiety level (outcome) of 1.23 points before large ensemble performances (large = 1) compared to other performance types (large = 0).
    • See multi-level>> interpretation for explainer on the interaction slope interpretations
  • \(\hat \beta_1 = −0.95\) — After accounting for baseline negative emotionality (mpqnem = 0), orchestral instrument players (orch = 1) have an estimated mean performance anxiety level (outcome) before solos and small ensembles (large = 0) which is 1.00 point higher than keyboardists and vocalists (orch = 0), while the mean performance anxiety (outcome) of orchestral players (orch = 1) is only .05 points higher before large ensembles (large = 1) (a difference of .95 points).
    • See multi-level >> interpretation for explainer on the interaction slope interpretations
  • \(\hat \beta_2 = −0.03\) — After accounting for instrument, a one-point increase in baseline negative emotionality is associated with an estimated 0.15 mean increase in performance anxiety levels before solos and small ensembles, but only an estimated 0.12 increase before large ensembles (a difference of .03 points).

• LR-Test for comparing nested models

  drop_in_dev <- anova(model.d, model.c, test = "Chisq")
  #>         npar  AIC  BIC logLik  dev Chisq Df      pval
  #> model.c    8 3007 3041  -1496 2991    NA NA        NA
  #> model.d   10 2996 3039  -1488 2976 14.73  2 0.0006319

  • MLE must be used instead of REML to estimate the parameters of each model
  • model.d is the current and more complex model, Random Slopes and Intercepts Model (with 3 covariates)
  • model.c is the less complex model, Random Slopes and Intercepts Model (with 2 covariates)
  • Process
    • The likelihood is larger (and the log-likelihood is less negative) under the larger model (Model D);
    • \(\text{Chisq} = 14.734 = -2 \cdot (-1488 - (-1496))\)
    • Using dof = 2, signifying the number of additional terms in Model D, we obtain a p-value of .0006.
  • A p-value < 0.05 indicates the difference in log-likelihoods is significantly different from 0, and therefore the more complex model, model.d, fits the data better.
  • Note: dropping the mpqnem:large term which has a t-stat < |2| (-0.5753) produces a better model than model.d according the LR-test

Final Model

• Description

  • One of other valid final models
  • Performance type categorical is no longer collapsed into the binary “large” ensemble/not large and is has now been collapsed into the binary, solo, not solo
  • Audience categorical has been transformed to dummies: students, juried, public with instructor as the reference category.
  • Varying/random slopes for previous, students, juried, public, and solo
  • mpqnem is in the solo level 2 equation, so the combined model has an interaction between the two in the fixed effects.

• Specification

  • Level 1: \(Y_{ij} = a_i + b_i \text{previous}_{ij} + c_i \text{students}_{ij} + d_i \text{juried}_{ij} + e_i \text{public}_{ij} + f_i \text{solo}_{ij} + \epsilon_{ij}\)

  • Level 2

    \[ \begin{align} a_i &= \alpha_0 + \alpha_1 \text{mpqpem}_i + \alpha_2 \text{mpqab}_i + \alpha_3 \text{orch}_i + \alpha_4 \text{mpqnem}_i + u_i \\ b_i &= \beta_0 + v_i \\ c_i &= \gamma_0 + w_i \\ d_i &= \delta_0 + x_i \\ e_i &= \epsilon_0 + y_i \\ f_i &= \zeta_0 + \zeta_1 \text{mpqnem}_i + z_i \end{align} \]

    • Variance-Covariance matrix

      \[ \left[\begin{array}{cc} u_i\\v_i\\w_i\\x_i\\y_i\\z_i \end{array} \right] \sim \mathcal{N} \left( \left[\begin{array}{cc} 0\\0\\0\\0\\0\\0\end{array} \right], \begin{bmatrix} \sigma_u^2 \\ \sigma_{uv} & \sigma_v^2 \\ \sigma_{uw} & \sigma_{vw} & \sigma^2_w \\ \sigma_{ux} & \sigma_{vx} & \sigma_{wx} & \sigma_x^2 \\ \sigma_{uy} & \sigma_{vy} & \sigma_{wy} & \sigma_{xy} & \sigma_{y}^2 \\ \sigma_{uz} & \sigma_{vz} & \sigma_{wz} & \sigma_{xz} & \sigma_{yz} & \sigma_z^2 \end{bmatrix} \right) \]

    • 6 variance terms and 15 correlation terms at Level Two, along with 1 variance term at Level One.

      • The number of correlation terms is equal to the number of unique pairs among Level Two random effects

• Model

  # Model F (One - of many - reasonable final models)
  model.f <- lmer(na ~ previous + students + juried + 
      public + solo + mpqpem + mpqab + orch + mpqnem + 
      mpqnem:solo + (previous + students + juried + 
      public + solo | id), REML = T, data = music)
  
  ##  Groups  Name        Variance Std.Dev. Corr       
  ##  id      (Intercept) 14.4802  3.805               
  ##          previous     0.0707  0.266    -0.65     
  ##          students     8.2151  2.866    -0.63  0.00
  ##          juried      18.3177  4.280    -0.64 -0.12
  ##          public      12.8094  3.579    -0.83  0.33
  ##          solo         0.7665  0.876    -0.67  0.47
  ##  Residual            15.2844  3.910     
  ##             
  ##  0.84           
  ##  0.66  0.58     
  ##  0.49  0.21  0.90
  ## 
  ##  Number of Level Two groups =  37
  ##            Estimate Std. Error t value
  ## (Intercept)  8.36883    1.91369  4.3731
  ## previous    -0.14303    0.06247 -2.2895
  ## students     3.61115    0.76796  4.7022
  ## juried       4.07332    1.03130  3.9497
  ## public       3.06453    0.89274  3.4327
  ## solo         0.51647    1.39635  0.3699
  ## mpqpem      -0.08312    0.02408 -3.4524
  ## mpqab        0.20377    0.04740  4.2986
  ## orch         1.53138    0.58384  2.6230
  ## mpqnem       0.11465    0.03591  3.1930
  ## solo:mpqnem  0.08296    0.04158  1.9951

• Interpretation

  • In general
    • Performance anxiety (outcome) is higher when a musician is performing in front of students, a jury, or the general public rather than their instructor
      • students, juried, and public are indicator variables created from the audience categorical variable (so that “Instructor” is the reference level in this model)
    • Performance anxiety (outcome) is lower for each additional diary the musician previously filled out
      • previous is the number of previous diary entries filled out by that individual
    • Musicians with lower levels of positive emotionality (mpqpem) and higher levels of absorption (mpqab) tend to experience greater performance anxiety (outcome)
    • Those who play orchestral instruments experience (orch = 1) more performance anxiety than those who play keyboards or sing (orch = 0)
  • Key terms
    • \(\hat \alpha_4 = 0.11\) (mpqnem, fixed) — A one-point increase in baseline level of negative emotionality (mpqnem) is associated with an estimated 0.11 mean increase in performance anxiety (outcome) for musicians performing in an ensemble group (solo=0), after controlling for previous diary entries (previous), audience (dummy vars), positive emotionality (mpqpem), absorption (mpqab), and instrument (orch).
    • \(\hat \zeta_1 = 0.08\) (solo:mpqnem, fixed) — When musicians play solos (solo = 1), a one-point increase in baseline level of negative emotionality (mpqnem) is associated with an estimated 0.19 mean increase in performance anxiety (outcome), 0.08 points (73%) higher than musicians playing in ensemble groups (solo = 0), controlling for the effects of previous diary entries (previous), audience (dummy vars), positive emotionality (mpqpem), absorption (mpqab), and instrument (orch).
      • See multi-level >> interpretation for explainer on the interaction slope interpretations
  • Addressing the researchers’ primary hypothesis:
    • After controlling for all these factors, we have significant evidence that musicians with higher levels of negative emotionality (mpqpem) experience higher levels of performance anxiety (outcome), and that this association is even more pronounced (interaction) when musicians are performing solos (solo = 1) rather than as part of an ensemble group (solo = 0).

Chapter 9: Longitudinal Data

• Beyond Multiple Linear Regression, Chapter 9: Two-Level Longitudinal Data
• Research Questions
  • Which factors most influence a school’s performance in Minnesota Comprehensive Assessment (MCA) testing?
  • How do the average math MCA-II scores for 6th graders enrolled in charter schools differ from scores for students who attend non-charter public schools? Do these differences persist after accounting for differences in student populations?
  • Are there differences in yearly improvement between charter and non-charter public schools?
• Data Answers
  • Within school—changes over time
  • Between schools—effects of school-specific covariates (charter or non-charter, urban or rural, percent free and reduced lunch, percent special education, and percent non-white) on 2008 math scores and rate of change between 2008 and 2010.
• Misc
  • Missing data is a common phenomenon in longitudinal studies