Module 5: Random Slope Models

Part 1: Run and Interpret Random Intercept Model

1. Treating thk as the DV, and school as the level-2 clustering variable, estimate a conditional random intercept model with prethk, cc, and gprethk as predictors. Interpret the results. Obtain fit statistics (AIC and BIC) for this baseline model.

This model shows ‘cc’ and ‘gprethk’ are statistically significant. The tobacco prevention curriculum seems to have an impact. As Dr. Broda stated, “This suggests that school context may be quite important in preventing tobacco use. Schools with prevention curricula, as well as schools where students on average tend to know more about tobacco prevention, both seem to have a positive impact on tobacco prevention knowledge for individual students.”

The fit statistics are as follows:

AIC = 4712.6 BIC =4744.9

Part 2: Adding Random Slopes

2. Now, add a random slope for prethk at level 2. Stata people: use the default independent error covariance structure. Interpret the fixed estimate for prethk and the random variance. Obtain fit statistics and decide whether or not a random slope for prethk improves the fit of the model.

The random slope for prethk does not improve the fit of the model. AIC increased, which is not a good sign.

The fit statistics are as follows: AIC = 4715.5 BIC = 4758.6

3. Stata people: rerun the same model this time using the unstructured covariance option. How is this model different from the model in #3 above? R folks: this is the model you ran above :).

Part 3: Testing a Cross-Level Interaction

4. Finally, test a cross-level interaction between prethk and gprethk (don’t forget to create new variables representing the centered versions of prethk and gprethk). Interpret the results of the interaction effect. Use margins (Stata) or interplot (R) to obtain predicted slopes for gprethk for students with prethk (centered) scores from -2 to 2. Visualize and explain your results.This week, we are going to use data from Gavin and Hofmann (2002), a study on organizational climate and attitudes published in Leadership Quarterly. Here, we have individuals soldiers nested within companies. This is the same dataset that Garson uses in Chapter 6, so you can recreate his analysis.

The interaction between prethk and gprethk appears to be signifcant. The model improved when this interaction was included. The slope estimate for gprethk is going to have the most impact on thk for someone who has the highest prethk score.

AIC = 4702.9 BIC = 4740.5

Load Some Packages to Help with the Analysis and Data Management:

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(psych))
library(haven)
library(tibble)

Load and Explore the Data

schoolsmoke <- haven::read_dta("schoolsmoke.dta")
glimpse(schoolsmoke)

Rows: 1,600
Columns: 5
$ school  <dbl> 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193...
$ thk     <dbl> 1, 2, 3, 1, 2, 2, 3, 3, 2, 2, 4, 2, 2, 3, 4, 4, 1, 3, 4, 2, 1, 3, 1, 3, 2, 4, ...
$ prethk  <dbl> 5, 3, 4, 0, 1, 3, 1, 2, 3, 2, 2, 2, 3, 3, 2, 3, 3, 3, 2, 1, 3, 0, 1, 0, 1, 3, ...
$ cc      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ gprethk <dbl> 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846...

schoolsmoke.1 <- schoolsmoke %>%
  mutate(., cc.factor = as_factor(cc))

glimpse(schoolsmoke.1)  

Rows: 1,600
Columns: 6
$ school    <dbl> 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 1...
$ thk       <dbl> 1, 2, 3, 1, 2, 2, 3, 3, 2, 2, 4, 2, 2, 3, 4, 4, 1, 3, 4, 2, 1, 3, 1, 3, 2, 4...
$ prethk    <dbl> 5, 3, 4, 0, 1, 3, 1, 2, 3, 2, 2, 2, 3, 3, 2, 3, 3, 3, 2, 1, 3, 0, 1, 0, 1, 3...
$ cc        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ gprethk   <dbl> 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.1538...
$ cc.factor <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

Null Model

Treating “thk” as the DV, and “school” as the level-2 clustering variable (we will ignore classroom clustering for now- stay tuned!), estimate a null model. Compute, report, and interpret the ICC, as well as the likelihood ratio test at the bottom of the output. Is it worth conducting MLM on these data?

model.null <- lmer(thk ~  (1|school), REML = FALSE, data = schoolsmoke.1)
summary(model.null)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ (1 | school)
   Data: schoolsmoke.1

     AIC      BIC   logLik deviance df.resid 
  4833.0   4849.1  -2413.5   4827.0     1597 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9039 -0.7734  0.0693  0.9240  1.7430 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.09438  0.3072  
 Residual             1.16219  1.0780  
Number of obs: 1600, groups:  school, 28

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  2.60750    0.06538 25.83292   39.88   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lmerTest::rand(model.null)

ANOVA-like table for random-effects: Single term deletions

Model:
thk ~ (1 | school)
             npar  logLik    AIC   LRT Df Pr(>Chisq)    
<none>          3 -2413.5 4833.0                        
(1 | school)    2 -2445.6 4895.2 64.16  1  1.147e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

model.1 <- lmer(thk ~ prethk + (1|school), REML = FALSE, data = schoolsmoke.1)
summary(model.1)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ prethk + (1 | school)
   Data: schoolsmoke.1

     AIC      BIC   logLik deviance df.resid 
  4733.5   4755.0  -2362.7   4725.5     1596 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.12273 -0.82727  0.03082  0.83756  2.15165 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.07519  0.2742  
 Residual             1.09321  1.0456  
Number of obs: 1600, groups:  school, 28

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 2.152e+00  7.414e-02 6.142e+01   29.02   <2e-16 ***
prethk      2.175e-01  2.122e-02 1.599e+03   10.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
prethk -0.598

Question 1: Adding Level-2 Covariates

model.2 <- lmer(thk ~ prethk + cc.factor + gprethk + (1|school), REML = FALSE, data = schoolsmoke.1)
summary(model.2)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ prethk + cc.factor + gprethk + (1 | school)
   Data: schoolsmoke.1

     AIC      BIC   logLik deviance df.resid 
  4712.6   4744.9  -2350.3   4700.6     1594 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.23332 -0.79567  0.02374  0.79254  2.16704 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.01734  0.1317  
 Residual             1.09334  1.0456  
Number of obs: 1600, groups:  school, 28

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 1.076e+00  2.548e-01 3.024e+01   4.224 0.000203 ***
prethk      2.126e-01  2.138e-02 1.571e+03   9.944  < 2e-16 ***
cc.factor1  4.308e-01  7.554e-02 2.330e+01   5.702 7.94e-06 ***
gprethk     4.183e-01  1.193e-01 3.387e+01   3.507 0.001301 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) prethk cc.fc1
prethk      0.000              
cc.factor1 -0.288  0.000       
gprethk    -0.963 -0.179  0.148

Using the modelsummary and broom.mixed Packages to Organize Your Results:

library(modelsummary)
library(broom.mixed)

Step Three: Add a Random Slope for prethk

model.3 <- lmer(thk ~ prethk + cc.factor + gprethk + (prethk|school), REML = FALSE, data = schoolsmoke.1)
summary(model.3)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ prethk + cc.factor + gprethk + (prethk | school)
   Data: schoolsmoke.1

     AIC      BIC   logLik deviance df.resid 
  4715.5   4758.6  -2349.8   4699.5     1592 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.32829 -0.81670  0.02938  0.81294  2.14320 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 school   (Intercept) 0.049006 0.22137       
          prethk      0.002063 0.04542  -0.97
 Residual             1.090740 1.04438       
Number of obs: 1600, groups:  school, 28

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.90319    0.24823 25.53020   3.639 0.001216 ** 
prethk       0.21136    0.02323 23.87728   9.097 3.16e-09 ***
cc.factor1   0.42343    0.07337 19.72400   5.771 1.27e-05 ***
gprethk      0.50405    0.11449 24.46798   4.403 0.000183 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) prethk cc.fc1
prethk     -0.065              
cc.factor1 -0.295 -0.004       
gprethk    -0.951 -0.174  0.160

Using lmerTest to Evaluate Random Slope:

lmerTest::rand(model.3)

ANOVA-like table for random-effects: Single term deletions

Model:
thk ~ prethk + cc.factor + gprethk + (prethk | school)
                            npar  logLik    AIC    LRT Df Pr(>Chisq)
<none>                         8 -2349.8 4715.5                     
prethk in (prethk | school)    6 -2350.3 4712.6 1.0762  2     0.5839

Getting Fancy: Modeling a Cross-Level Interaction

This is more advanced stuff, but I want to show you how you might model a cross-level interaction. Last week, we looked at the interaction between LEAD and TSIG, which is an interaction between two level-1 (soldier) covariates. In MLM, it is really easy to look at cross-level interactions, as well. So, let’s say we want to test whether there is an interaction between TSIG (level-1) and GTSIG (level-2). Put another way, is there an interaction between how I feel about the importance of my job, and how my peers feel about the importance of their job? Probably- it is worth testing.

Step One: Created Centered Versions of your Two Predictors

Centering a variable means to subtract the mean from every observation. So, now that variables has a mean of 0 (0 now means average). When we look at interaction effects, it is usually a good idea to center any continuous predictors before analyzing. This facilitates easier interpretation of the results, and it ensures that the interaction is being evaluated near the sample mean.

Centering works like this: you use the mutate function from dplyr to create a new variable that is equal to the original variable minus the mean (which we get using the mean function). Just to check our work, we can call the describe function from psych and make sure that the means of those new variables are 0.

schoolsmoke.2 <- schoolsmoke.1 %>%
  mutate(.,
         prethk_cent = prethk - mean(prethk),
         gprethk_cent = gprethk - mean(gprethk))

psych::describe(schoolsmoke.2, fast = TRUE)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

Step Two: Run the MLM and Interpret the Results

model.4 <- lmer(thk ~ prethk_cent + cc.factor + gprethk_cent + (1|school), REML = FALSE, data = schoolsmoke.2)
summary(model.4)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ prethk_cent + cc.factor + gprethk_cent + (1 | school)
   Data: schoolsmoke.2

     AIC      BIC   logLik deviance df.resid 
  4712.6   4744.9  -2350.3   4700.6     1594 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.23332 -0.79567  0.02374  0.79254  2.16704 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.01734  0.1317  
 Residual             1.09334  1.0456  
Number of obs: 1600, groups:  school, 28

Fixed effects:
              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)  2.382e+00  5.235e-02 2.179e+01  45.499  < 2e-16 ***
prethk_cent  2.126e-01  2.138e-02 1.571e+03   9.944  < 2e-16 ***
cc.factor1   4.308e-01  7.554e-02 2.330e+01   5.702 7.94e-06 ***
gprethk_cent 4.183e-01  1.193e-01 3.387e+01   3.507   0.0013 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) prthk_ cc.fc1
prethk_cent  0.000              
cc.factor1  -0.701  0.000       
gprethk_cnt -0.123 -0.179  0.148

Interaction Effects Model

model.5 <- lmer(thk ~ prethk_cent + cc.factor + gprethk_cent + prethk_cent:gprethk_cent + (1|school), REML = FALSE, data = schoolsmoke.2)
summary(model.5)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: thk ~ prethk_cent + cc.factor + gprethk_cent + prethk_cent:gprethk_cent +  
    (1 | school)
   Data: schoolsmoke.2

     AIC      BIC   logLik deviance df.resid 
  4702.9   4740.5  -2344.4   4688.9     1593 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.24644 -0.81395  0.04367  0.80525  2.07026 

Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.0165   0.1284  
 Residual             1.0857   1.0420  
Number of obs: 1600, groups:  school, 28

Fixed effects:
                          Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)              2.359e+00  5.205e-02 2.243e+01  45.322  < 2e-16 ***
prethk_cent              2.080e-01  2.135e-02 1.572e+03   9.743  < 2e-16 ***
cc.factor1               4.300e-01  7.453e-02 2.337e+01   5.770 6.64e-06 ***
gprethk_cent             3.651e-01  1.189e-01 3.469e+01   3.071 0.004131 ** 
prethk_cent:gprethk_cent 2.290e-01  6.663e-02 1.502e+03   3.437 0.000605 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) prthk_ cc.fc1 gprth_
prethk_cent  0.008                     
cc.factor1  -0.695  0.000              
gprethk_cnt -0.104 -0.171  0.148       
prthk_cnt:_ -0.128 -0.063 -0.003 -0.130

Should We Include That Interaction? Comparing model.4 with model.5:

anova(model.5, model.4)

Data: schoolsmoke.2
Models:
model.4: thk ~ prethk_cent + cc.factor + gprethk_cent + (1 | school)
model.5: thk ~ prethk_cent + cc.factor + gprethk_cent + prethk_cent:gprethk_cent + 
model.5:     (1 | school)
        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
model.4    6 4712.6 4744.9 -2350.3   4700.6                         
model.5    7 4702.9 4740.5 -2344.4   4688.9 11.764  1   0.000604 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step Three: Use interplot to Obtain More Easily Interpretable Results

interplot::interplot(model.5, var1 = "gprethk_cent", var2 = "prethk_cent")

Using the modelsummary and broom.mixed Packages to Organize Your Results:

library(modelsummary)
library(broom.mixed)

models <- list(model.4, model.5)
modelsummary(models, output = "markdown")

	Model 1	Model 2
(Intercept)	2.382	2.359
	(0.052)	(0.052)
cc.factor1	0.431	0.430
	(0.076)	(0.075)
gprethk_cent	0.418	0.365
	(0.119)	(0.119)
prethk_cent	0.213	0.208
	(0.021)	(0.021)
sd__(Intercept)	0.132	0.128
sd__Observation	1.046	1.042
prethk_cent × gprethk_cent		0.229
		(0.067)
AIC	4712.6	4702.9
BIC	4744.9	4740.5
Log.Lik.	-2350.308	-2344.426

HTML Version That You Can Open in Word:

modelsummary(models, output = 'msum.html', title = 'MLM Estimates')

[WARNING] This document format requires a nonempty <title> element.
  Please specify either 'title' or 'pagetitle' in the metadata,
  e.g. by using --metadata pagetitle="..." on the command line.
  Falling back to 'msum'