Multilevel Modeling, Module 5: Random Slope Models

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.

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

suppressPackageStartupMessages(library(tidyverse))

package ‘tibble’ was built under R version 3.6.2package ‘dplyr’ was built under R version 3.6.2

suppressPackageStartupMessages(library(lme4))

This week, we are going to continue 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 7, so you can recreate his analysis.

Load and Explore the Data

glimpse(lq2002)

Rows: 2,042
Columns: 28
$ ROWID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22…
$ COMPID   <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3,…
$ SUB      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22…
$ LEAD01   <dbl> 2, 4, 4, 1, 1, 2, 3, 3, 3, 4, 2, 2, 4, 4, 1, 2, 4, 4, 4, 4, 4, 3, 2, 4, 5, 3,…
$ LEAD02   <dbl> 2, 4, 4, 1, 1, 3, 2, 3, 4, 3, 4, 3, 4, 4, 1, 1, 5, 4, 5, 4, 2, 4, 3, 4, 5, 4,…
$ LEAD03   <dbl> 2, 2, 2, 1, 1, 2, 2, 2, 3, 3, 1, 1, 4, 4, 1, 2, 3, 4, 1, 4, 2, 4, 2, 4, 4, 1,…
$ LEAD04   <dbl> 3, 4, 4, 1, 1, 4, 1, 4, 2, 4, 3, 1, 4, 4, 1, 4, 5, 4, 1, 5, 1, 4, 2, 3, 4, 3,…
$ LEAD05   <dbl> 2, 4, 2, 1, 3, 2, 2, 2, 3, 3, 1, 2, 4, 4, 1, 3, 3, 4, 1, 2, 2, 3, 3, 3, 3, 1,…
$ LEAD06   <dbl> 2, 4, 4, 1, 2, 4, 2, 2, 3, 4, 3, 4, 4, 4, 1, 3, 5, 4, 4, 2, 2, 4, 4, 4, 5, 4,…
$ LEAD07   <dbl> 5, 4, 4, 4, 1, 3, 2, 2, 3, 3, 3, 1, 4, 3, 1, 5, 3, 4, 1, 3, 2, 3, 1, 3, 3, 3,…
$ LEAD08   <dbl> 5, 4, 4, 2, 1, 4, 2, 4, 1, 4, 4, 3, 4, 4, 1, 2, 5, 4, 1, 3, 2, 4, 4, 3, 5, 3,…
$ LEAD09   <dbl> 5, 4, 4, 3, 1, 2, 3, 2, 3, 3, 2, 2, 4, 4, 1, 4, 1, 4, 1, 4, 2, 2, 1, 4, 3, 3,…
$ LEAD10   <dbl> 5, 4, 4, 3, 1, 4, 4, 2, 2, 4, 3, 2, 4, 4, 1, 2, 5, 4, 4, 4, 1, 4, 3, 4, 4, 3,…
$ LEAD11   <dbl> 4, 3, 2, 2, 1, 3, 2, 4, 3, 3, 2, 1, 4, 4, 1, 2, 4, 4, 1, 4, 2, 3, 3, 3, 4, 3,…
$ TSIG01   <dbl> 1, 3, 5, 1, 4, 4, 2, 2, 4, 4, 3, 3, 4, 4, 1, 1, 4, 4, 1, 5, 1, 4, 2, 4, 4, 3,…
$ TSIG02   <dbl> 5, 4, 5, 3, 4, 4, 4, 2, 3, 4, 4, 3, 4, 4, 1, 5, 4, 4, 5, 5, 4, 4, 3, 4, 4, 3,…
$ TSIG03   <dbl> 5, 3, 5, 2, 4, 4, 4, 4, 4, 4, 3, 3, 4, 4, 1, 4, 5, 4, 5, 5, 4, 4, 3, 4, 5, 3,…
$ HOSTIL01 <dbl> 0, 3, 1, 2, 2, 0, 2, 0, 0, 1, 1, 2, 0, 0, 4, 3, 1, 0, 4, 3, 4, 1, 3, 0, 2, 0,…
$ HOSTIL02 <dbl> 0, 2, 0, 1, 0, 0, 0, 0, 4, 0, 0, 2, 0, 0, 4, 0, 0, 0, 0, 1, 4, 1, 3, 0, 0, 0,…
$ HOSTIL03 <dbl> 0, 0, 0, 3, 0, 0, 1, 0, 4, 0, 0, 0, 0, 0, 4, 0, 0, 0, 4, 2, 2, 0, 3, 0, 0, 0,…
$ HOSTIL04 <dbl> 0, 0, 0, 3, 0, 0, 1, 0, 4, 0, 0, 1, 0, 0, 4, 0, 0, 0, 4, 3, 2, 0, 1, 0, 0, 0,…
$ HOSTIL05 <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 4, 0, 0, 2, 0, 0, 2, 0, 0, 0, 4, 2, 3, 0, 1, 0, 0, 0,…
$ LEAD     <dbl> 3.363636, 3.727273, 3.454545, 1.818182, 1.272727, 3.000000, 2.272727, 2.72727…
$ TSIG     <dbl> 3.666667, 3.333333, 5.000000, 2.000000, 4.000000, 4.000000, 3.333333, 2.66666…
$ HOSTILE  <dbl> 0.0, 1.2, 0.4, 1.8, 0.4, 0.0, 0.8, 0.0, 3.2, 0.2, 0.2, 1.4, 0.0, 0.0, 3.6, 0.…
$ GLEAD    <dbl> 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.88257…
$ GTSIG    <dbl> 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.54166…
$ GHOSTILE <dbl> 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, …

Data Management Tip of the Week: Calculating Group (Cluster) Means

Here is something that comes up a lot in MLM. Often, you want to take an individual-level measure and calculate the cluster, or group, mean. This is “aggregating up” an individual measure to see if it is a significant predictor at level 2. We did this last week with the task significance measure (TSIG). My own sense of task significance is important, but it is probably also important to consider how my colleagues feel about their work as well.

lq2002.1 <- lq2002 %>%
  group_by(COMPID) %>%
  mutate(.,
            glead = mean(LEAD),
            gtsig = mean(TSIG),
            ghostile = mean(GHOSTILE)) %>%
  ungroup()

glimpse(lq2002.1)  

Rows: 2,042
Columns: 31
$ ROWID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22…
$ COMPID   <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3,…
$ SUB      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22…
$ LEAD01   <dbl> 2, 4, 4, 1, 1, 2, 3, 3, 3, 4, 2, 2, 4, 4, 1, 2, 4, 4, 4, 4, 4, 3, 2, 4, 5, 3,…
$ LEAD02   <dbl> 2, 4, 4, 1, 1, 3, 2, 3, 4, 3, 4, 3, 4, 4, 1, 1, 5, 4, 5, 4, 2, 4, 3, 4, 5, 4,…
$ LEAD03   <dbl> 2, 2, 2, 1, 1, 2, 2, 2, 3, 3, 1, 1, 4, 4, 1, 2, 3, 4, 1, 4, 2, 4, 2, 4, 4, 1,…
$ LEAD04   <dbl> 3, 4, 4, 1, 1, 4, 1, 4, 2, 4, 3, 1, 4, 4, 1, 4, 5, 4, 1, 5, 1, 4, 2, 3, 4, 3,…
$ LEAD05   <dbl> 2, 4, 2, 1, 3, 2, 2, 2, 3, 3, 1, 2, 4, 4, 1, 3, 3, 4, 1, 2, 2, 3, 3, 3, 3, 1,…
$ LEAD06   <dbl> 2, 4, 4, 1, 2, 4, 2, 2, 3, 4, 3, 4, 4, 4, 1, 3, 5, 4, 4, 2, 2, 4, 4, 4, 5, 4,…
$ LEAD07   <dbl> 5, 4, 4, 4, 1, 3, 2, 2, 3, 3, 3, 1, 4, 3, 1, 5, 3, 4, 1, 3, 2, 3, 1, 3, 3, 3,…
$ LEAD08   <dbl> 5, 4, 4, 2, 1, 4, 2, 4, 1, 4, 4, 3, 4, 4, 1, 2, 5, 4, 1, 3, 2, 4, 4, 3, 5, 3,…
$ LEAD09   <dbl> 5, 4, 4, 3, 1, 2, 3, 2, 3, 3, 2, 2, 4, 4, 1, 4, 1, 4, 1, 4, 2, 2, 1, 4, 3, 3,…
$ LEAD10   <dbl> 5, 4, 4, 3, 1, 4, 4, 2, 2, 4, 3, 2, 4, 4, 1, 2, 5, 4, 4, 4, 1, 4, 3, 4, 4, 3,…
$ LEAD11   <dbl> 4, 3, 2, 2, 1, 3, 2, 4, 3, 3, 2, 1, 4, 4, 1, 2, 4, 4, 1, 4, 2, 3, 3, 3, 4, 3,…
$ TSIG01   <dbl> 1, 3, 5, 1, 4, 4, 2, 2, 4, 4, 3, 3, 4, 4, 1, 1, 4, 4, 1, 5, 1, 4, 2, 4, 4, 3,…
$ TSIG02   <dbl> 5, 4, 5, 3, 4, 4, 4, 2, 3, 4, 4, 3, 4, 4, 1, 5, 4, 4, 5, 5, 4, 4, 3, 4, 4, 3,…
$ TSIG03   <dbl> 5, 3, 5, 2, 4, 4, 4, 4, 4, 4, 3, 3, 4, 4, 1, 4, 5, 4, 5, 5, 4, 4, 3, 4, 5, 3,…
$ HOSTIL01 <dbl> 0, 3, 1, 2, 2, 0, 2, 0, 0, 1, 1, 2, 0, 0, 4, 3, 1, 0, 4, 3, 4, 1, 3, 0, 2, 0,…
$ HOSTIL02 <dbl> 0, 2, 0, 1, 0, 0, 0, 0, 4, 0, 0, 2, 0, 0, 4, 0, 0, 0, 0, 1, 4, 1, 3, 0, 0, 0,…
$ HOSTIL03 <dbl> 0, 0, 0, 3, 0, 0, 1, 0, 4, 0, 0, 0, 0, 0, 4, 0, 0, 0, 4, 2, 2, 0, 3, 0, 0, 0,…
$ HOSTIL04 <dbl> 0, 0, 0, 3, 0, 0, 1, 0, 4, 0, 0, 1, 0, 0, 4, 0, 0, 0, 4, 3, 2, 0, 1, 0, 0, 0,…
$ HOSTIL05 <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 4, 0, 0, 2, 0, 0, 2, 0, 0, 0, 4, 2, 3, 0, 1, 0, 0, 0,…
$ LEAD     <dbl> 3.363636, 3.727273, 3.454545, 1.818182, 1.272727, 3.000000, 2.272727, 2.72727…
$ TSIG     <dbl> 3.666667, 3.333333, 5.000000, 2.000000, 4.000000, 4.000000, 3.333333, 2.66666…
$ HOSTILE  <dbl> 0.0, 1.2, 0.4, 1.8, 0.4, 0.0, 0.8, 0.0, 3.2, 0.2, 0.2, 1.4, 0.0, 0.0, 3.6, 0.…
$ GLEAD    <dbl> 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.88257…
$ GTSIG    <dbl> 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.54166…
$ GHOSTILE <dbl> 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, …
$ glead    <dbl> 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.882576, 2.88257…
$ gtsig    <dbl> 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.541667, 3.54166…
$ ghostile <dbl> 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, 1.0416667, …

The above command has two parts. First, we use the group_by command from our friend, the dplyr package, to group all responses in their respective clusters/companies determined by COMPID. This is a first step to declare that all the operations that follow are to be performed at the cluster (company) level. Second, we use the mutate command, which we have used before, to create new variables that are equal to the mean of the original variables. Since we just grouped by COMPID, this means we get the mean of the variables by company.

So, we calculate new variables called glead, gtsig, and ghostile that are equal to the mean of LEAD, TSIG, and HOSTILE by company. This is actually identical to what is already calculated for us in the dataset- the variables GLEAD, GTSIG, and GHOSTILE. I just wanted to show you how to do it so you can try it on your own!

Running Random Slopes Models

Step One: Run Your Preferred Model Without Random Slopes

Let’s assume that we ran a bunch of models with this army data (which we did, last week!), and arrived at a model that we thought best represented our data and matched our research questions. That is the model below, with a random intercept, two level-1 predictors (LEAD and TSIG), and one level-2 predictor (GTSIG). Let’s rerun that model to remind ourselves and also get measures of model fit (AIC and BIC) that we can use as we go along.

model.1 <- lmer(HOSTILE ~ TSIG + LEAD + GTSIG + (1|COMPID), REML = FALSE, data = lq2002.1)
summary(model.1)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: HOSTILE ~ TSIG + LEAD + GTSIG + (1 | COMPID)
   Data: lq2002.1

     AIC      BIC   logLik deviance df.resid 
  5533.7   5567.4  -2760.8   5521.7     2036 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3166 -0.6960 -0.2211  0.5170  3.4998 

Random effects:
 Groups   Name        Variance Std.Dev.
 COMPID   (Intercept) 0.009905 0.09952 
 Residual             0.867050 0.93116 
Number of obs: 2042, groups:  COMPID, 49

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.40718    0.25785  13.214
TSIG        -0.17234    0.02479  -6.953
LEAD        -0.35177    0.02996 -11.741
GTSIG       -0.27729    0.08258  -3.358

Correlation of Fixed Effects:
      (Intr) TSIG   LEAD  
TSIG   0.105              
LEAD  -0.205 -0.511       
GTSIG -0.942 -0.227  0.010

Step Two: Add a Random Slope for TSIG

Now, the next logical step after testing predictors is to see whether or not any predictors should be treated as random at level-2. Here, we are going to test a random slope for TSIG. This is effectively allowing the slope of TSIG (which was overall negative) to vary according to which company a solider is in. In some companies, that slope could now be higher, in some companies, it could be lower.

In terms of coding, to add a random slope, you just list the predictor name again after the colon. Notice now that TSIG appears twice in the code- once to the left of the ~, and once on the right side of the pipe operator (|).

model.2 <- lmer(HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG|COMPID), REML = FALSE, data = lq2002.1)

boundary (singular) fit: see ?isSingular

summary(model.2)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG | COMPID)
   Data: lq2002.1

     AIC      BIC   logLik deviance df.resid 
  5517.3   5562.3  -2750.6   5501.3     2034 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6985 -0.6803 -0.2299  0.4954  3.4722 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 COMPID   (Intercept) 0.21072  0.4590        
          TSIG        0.01234  0.1111   -1.00
 Residual             0.85150  0.9228        
Number of obs: 2042, groups:  COMPID, 49

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.16962    0.25633  12.365
TSIG        -0.15484    0.03037  -5.099
LEAD        -0.34450    0.02964 -11.624
GTSIG       -0.23218    0.07685  -3.021

Correlation of Fixed Effects:
      (Intr) TSIG   LEAD  
TSIG  -0.112              
LEAD  -0.227 -0.406       
GTSIG -0.897 -0.170  0.031
convergence code: 0
boundary (singular) fit: see ?isSingular

Using lmerTest to Evaluate Random Slope:

lmerTest::rand(model.2)

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

Model:
HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG | COMPID)
                        npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>                     8 -2750.6 5517.3                         
TSIG in (TSIG | COMPID)    6 -2760.8 5533.7 20.388  2   3.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step Three: Add a Random Slope for LEAD

Same thing here, except now we are testing a random slope for LEAD:

model.3 <- lmer(HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG + LEAD|COMPID), REML = FALSE, data = lq2002.1)

boundary (singular) fit: see ?isSingular

summary(model.3)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG + LEAD | COMPID)
   Data: lq2002.1

     AIC      BIC   logLik deviance df.resid 
  5521.0   5582.9  -2749.5   5499.0     2031 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5998 -0.6759 -0.2227  0.4955  3.5053 

Random effects:
 Groups   Name        Variance Std.Dev. Corr       
 COMPID   (Intercept) 0.24520  0.4952              
          TSIG        0.01646  0.1283   -0.71      
          LEAD        0.01176  0.1085   -0.31 -0.45
 Residual             0.84526  0.9194              
Number of obs: 2042, groups:  COMPID, 49

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.19270    0.26113  12.226
TSIG        -0.15408    0.03211  -4.799
LEAD        -0.34390    0.03470  -9.910
GTSIG       -0.24041    0.07788  -3.087

Correlation of Fixed Effects:
      (Intr) TSIG   LEAD  
TSIG  -0.085              
LEAD  -0.242 -0.481       
GTSIG -0.893 -0.156  0.024
convergence code: 0
boundary (singular) fit: see ?isSingular

Using lmerTest to Evaluate Random Slope:

lmerTest::rand(model.3)

boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular

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

Model:
HOSTILE ~ TSIG + LEAD + GTSIG + (TSIG + LEAD | COMPID)
                               npar  logLik    AIC     LRT Df Pr(>Chisq)   
<none>                           11 -2749.5 5521.0                         
TSIG in (TSIG + LEAD | COMPID)    8 -2756.6 5529.3 14.2255  3   0.002614 **
LEAD in (TSIG + LEAD | COMPID)    8 -2750.6 5517.3  2.2441  3   0.523305   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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.

Step Two: Run the MLM and Interpret the Results

model.4 <- lmer(HOSTILE ~ tsig_cent + LEAD + gtsig_cent + (1|COMPID), REML = FALSE, data = lq2002.2)
summary(model.4)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: HOSTILE ~ tsig_cent + LEAD + gtsig_cent + (1 | COMPID)
   Data: lq2002.2

     AIC      BIC   logLik deviance df.resid 
  5533.7   5567.4  -2760.8   5521.7     2036 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3166 -0.6960 -0.2211  0.5170  3.4998 

Random effects:
 Groups   Name        Variance Std.Dev.
 COMPID   (Intercept) 0.009905 0.09952 
 Residual             0.867050 0.93116 
Number of obs: 2042, groups:  COMPID, 49

Fixed effects:
            Estimate Std. Error t value
(Intercept)  1.99845    0.09384  21.296
tsig_cent   -0.17234    0.02479  -6.953
LEAD        -0.35177    0.02996 -11.741
gtsig_cent  -0.27729    0.08258  -3.358

Correlation of Fixed Effects:
           (Intr) tsg_cn LEAD  
tsig_cent   0.491              
LEAD       -0.960 -0.511       
gtsig_cent -0.020 -0.227  0.010

Interaction Effects Model

model.5 <- lmer(HOSTILE ~ tsig_cent + LEAD + gtsig_cent + tsig_cent:gtsig_cent + (1|COMPID), REML = FALSE, data = lq2002.2)
summary(model.5)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: HOSTILE ~ tsig_cent + LEAD + gtsig_cent + tsig_cent:gtsig_cent +      (1 | COMPID)
   Data: lq2002.2

     AIC      BIC   logLik deviance df.resid 
  5530.2   5569.6  -2758.1   5516.2     2035 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4241 -0.6923 -0.2381  0.5068  3.4601 

Random effects:
 Groups   Name        Variance Std.Dev.
 COMPID   (Intercept) 0.009531 0.09763 
 Residual             0.864954 0.93003 
Number of obs: 2042, groups:  COMPID, 49

Fixed effects:
                     Estimate Std. Error t value
(Intercept)           1.97511    0.09418  20.971
tsig_cent            -0.16826    0.02481  -6.781
LEAD                 -0.34961    0.02993 -11.682
gtsig_cent           -0.29177    0.08223  -3.548
tsig_cent:gtsig_cent  0.15799    0.06757   2.338

Correlation of Fixed Effects:
            (Intr) tsg_cn LEAD   gtsg_c
tsig_cent    0.479                     
LEAD        -0.958 -0.507              
gtsig_cent  -0.012 -0.232  0.008       
tsg_cnt:gt_ -0.107  0.070  0.032 -0.071

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

anova(model.5, model.4)

Data: lq2002.2
Models:
model.4: HOSTILE ~ tsig_cent + LEAD + gtsig_cent + (1 | COMPID)
model.5: HOSTILE ~ tsig_cent + LEAD + gtsig_cent + tsig_cent:gtsig_cent + 
model.5:     (1 | COMPID)
        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
model.4  6 5533.7 5567.4 -2760.8   5521.7                           
model.5  7 5530.2 5569.6 -2758.1   5516.2 5.4564      1     0.0195 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Step Three: Use interplot to Obtain More Easily Interpretable Results

Interpreting significant interactions is hard, there’s no way around it. But, the margins command can be super useful for this. We will use margins to compare hypothetical soldiers with different levels of TSIG who are in companies that have different levels of GTSIG. And then we will plot the results to create a nice visual.

interplot::interplot(model.5, var1 = "gtsig_cent", var2 = "tsig_cent")

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

modelsummary(models, output = "markdown")

`data_frame()` is deprecated as of tibble 1.1.0.
Please use `tibble()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

	Model 1	Model 2
(Intercept)	1.998	1.975
	(0.094)	(0.094)
gtsig_cent	-0.277	-0.292
	(0.083)	(0.082)
LEAD	-0.352	-0.350
	(0.030)	(0.030)
sd__(Intercept)	0.100	0.098
sd__Observation	0.931	0.930
tsig_cent	-0.172	-0.168
	(0.025)	(0.025)
tsig_cent × gtsig_cent		0.158
		(0.068)
AIC	5533.7	5530.2
BIC	5567.4	5569.6
Log.Lik.	-2760.833	-2758.105

HTML Version That You Can Open in Word:

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