Introduction

This exercise gives a basic introduction to a multilevel regression and shows how you can replicate the popularity data multilevel models from the book Multilevel analysis: Techniques and applications, Chapter 2. In this exercise, the software packages lme4 and lmerTest for R (Windows) are used. Results should be very similar to results obtained with other software packages. However, due to convergence and rounding issues, you might notice minor differences.  

Preparation

This exercise expects:

  • Basic knowledge of multilevel analyses (the first two chapters of the book are sufficient).
  • Basic knowledge of coding in R.
  • Installation of R packages lme4, haven, tidyverse, RColorBrewer and lmerTest. This exercise was made using R version 4.3.0

Step 1: Setting up packages

If you have not yet installed all below-mentioned packages, you can install them by the command install.packages("NAMEOFPACKAGE").

library(lme4) # for the analysis
library(haven) # to load the SPSS .sav file
library(tidyverse) # needed for data manipulation.
library(RColorBrewer) # needed for some extra colours in one of the graphs
library(lmerTest) # to get p-value estimations that are not part of the standard lme4 package
library(viridis) # for making nice color plots

Step 2: Data input

The popularity dataset popular2.sav contains characteristics of pupils in different classes. The main goal of this exercise is to find models and test hypotheses about the relation between these characteristics and the popularity of pupils (according to their classmates). This dataset can be downloaded here.

The file can be loaded using the following syntax (make sure the working directory is the one in which this data set is located).

popular2 <- read_sav(file ="popular2.sav")

There are some variables in the dataset that we do not use, so we can select the variables we will use and have a look at the first few observations.

popular2 <- select(popular2, pupil, class, extrav, sex, texp, popular)
head(popular2) 
## # A tibble: 6 × 6
##   pupil class extrav sex        texp popular
##   <dbl> <dbl>  <dbl> <dbl+lbl> <dbl>   <dbl>
## 1     1     1      5 1 [girl]     24     6.3
## 2     2     1      7 0 [boy]      24     4.9
## 3     3     1      4 1 [girl]     24     5.3
## 4     4     1      3 1 [girl]     24     4.7
## 5     5     1      5 1 [girl]     24     6  
## 6     6     1      4 0 [boy]      24     4.7

Step 3: Plotting the data

Before we start the analysis, we can plot the relationship between extraversion and popularity, without taking into consideration the multilevel structure of the data.

ggplot(data  = popular2,
       aes(x = extrav,
           y = popular))+
  geom_point(size = 1.2,
             alpha = .8)+
  theme_minimal()+
  labs(title = "Popularity vs. Extraversion")

Now we can add a regression line to this plot.

ggplot(data  = popular2,
       aes(x = extrav,
           y = popular))+
  geom_point(size     = 1.2,
             alpha    = .8)+
  geom_smooth(method = lm,
              se     = FALSE, 
              col    = "black",
              linewidth   = .5, 
              alpha  = .8)+ # to add regression line
  theme_minimal()+
  labs(title    = "Popularity vs. Extraversion",
       subtitle = "add regression line")
## `geom_smooth()` using formula = 'y ~ x'

There is an increasing relation between popularity and extraversion.

So far we have ignored the nested multilevel structure of the data. We can show this multilevel structure by colour coding the different classes.

ggplot(data    = popular2,
       aes(x   = extrav,
           y   = popular,
           col = class))+ #to add the colours for different classes
  geom_point(size     = 1.2,
             alpha    = .8)+
  theme_minimal()+
  theme(legend.position = "none")+
  scale_color_gradientn(colours = rainbow(100))+
  labs(title    = "Popularity vs. Extraversion",
       subtitle = "add colours for different classes")

Now we can draw different regression lines for the 100 different classes in the dataset

ggplot(data      = popular2,
       aes(x     = extrav,
           y     = popular,
           col   = class,
           group = class))+ #to add the colours for different classes
  geom_point(size     = 1.2,
             alpha    = .8)+ 
  theme_minimal()+
  theme(legend.position = "none")+
  scale_color_gradientn(colours = rainbow(100))+
  geom_smooth(method = lm,
              se     = FALSE,
              linewidth   = .5, 
              alpha  = .8)+ # to add regression line
  labs(title    = "Popularity vs. Extraversion",
       subtitle = "add colours for different classes and regression lines")
## `geom_smooth()` using formula = 'y ~ x'

We clearly see that the relationship between extraversion and popularity is not the same in all classes, but on average there is a clear positive relationship. In this exercise, we will show the estimation of these different slopes (and how to explain these differences). Again, for more information please refer to the book Multilevel analysis: Techniques and applications.

We can also colour code the most extreme regression lines.

# To colour code the extremes, we need to write a small function that calculates the regression lines and adds a column indicating which clusters have the most extreme.
f1 <- function(data, x, y, grouping, n.highest = 3, n.lowest = 3){
  groupinglevel <- data[,grouping]
  res           <- data.frame(coef = rep(NA, length(unique(groupinglevel))), group = unique(groupinglevel))
  names(res)    <- c("coef", grouping)
  for(i in 1:length(unique(groupinglevel))){
    data2    <- as.data.frame(data[data[,grouping] == i,])
    res[i,1] <- as.numeric(lm(data2[, y] ~ data2[, x])$coefficients[2])
  }
  top    <- res %>% top_n(n.highest, coef)
  bottom <- res %>% top_n(-n.lowest, coef)
  res    <- res %>% mutate(high_and_low = ifelse(coef %in% top$coef, "top",  ifelse(coef %in% bottom$coef, "bottom", "none")))
  data3  <- left_join(data, res)
  return(data3)
}

Now we can use this function on the popularity data.

f1(data = as.data.frame(popular2), 
   x    = "extrav",
   y    = "popular",
   grouping = "class",
   n.highest = 3, 
   n.lowest = 3) %>%
  ggplot()+
  geom_point(aes(x     = extrav,
                 y     = popular, 
                 fill  = class, 
                 group = class),
             size     =  1, 
             alpha    = .5, 
             shape    = 21, 
             col      = "white")+
  geom_smooth(aes(x     = extrav,
                  y     = popular,
                  col   = high_and_low,
                  group = class,
                  size  = as.factor(high_and_low),
                  alpha = as.factor(high_and_low)),
              method = lm,
              se     = FALSE)+
  theme_minimal()+
  theme(legend.position = "none")+
  scale_fill_gradientn(colours = rainbow(100))+
  scale_color_manual(values=c("top"      = "blue",
                              "bottom"   = "red",
                              "none"     = "grey40"))+
  scale_size_manual(values=c("top"       = 1.2,
                              "bottom"   = 1.2,
                              "none"     = .5))+
  scale_alpha_manual(values=c("top"      = 1,
                             "bottom"    = 1,
                             "none"      =.3))+
  labs(title="Linear Relationship Between Popularity and Extraversion for 100 Classes",
       subtitle="The 6 classes with the most extreme relationship have been highlighted red and blue")

In the three classes that are drawn in blue the relationship is positive, in the three classes drawn in red the relationship is negative.

Step 4: Analysing the data

Model 1: Intercept only model

The first model that we replicate is the intercept only model. The model at the pupil level is \[Popularity_{ij}=\beta_{0j}+e_{ij}.\] The intercept randomly varies across classes but we do not yet include class-level predictor variables to explain (part of) this variability. The model at the class level is \[\beta_{0j}=\gamma_{00}+u_{0j}.\]
Combined we get: \[Popularity_{ij}= \gamma_{00}+u_{0j}+e_{ij}.\] This model includes a mean \(\gamma_{00}\), a random effect \(u_{0j}\) at the class level related to the intercept and a random effect \(e_{ij}\) at the pupil level.

The R syntax is as follows:

model1 <- lmer(formula = popular ~ 1 + (1|class),
                           data    = popular2, REML = FALSE)

If we look at the different inputs for the lmer function we see the following:

  1. “popular”, which indicates the dependent variable we want to predict.
  2. a “~”, that is comparable to the ‘=’ of the regression equation.
  3. a “1” in the formula which indicates the intercept. Since this is an intercept only model, we do not have any independent variables here.
  4. a random intercept, specified in the random part (within brackets). Again the value 1 is to indicate the intercept and the variables right of the vertical “|” bar are used to indicate grouping variables, in this case the class ID. So the dependent variable ‘popular’ is predicted by an intercept and a random error term for the intercept.
  5. the dataset we want to use is given after the data = command
  6. the model is estimated using full maximum likelihood: REML = FALSE

For more information on the lmer() function see the lme4 manual

The following output displays the results

summary(model1) 
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: popular ~ 1 + (1 | class)
##    Data: popular2
## 
##      AIC      BIC   logLik deviance df.resid 
##   6333.5   6350.3  -3163.7   6327.5     1997 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5662 -0.6983  0.0021  0.6758  3.3173 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.6945   0.8333  
##  Residual             1.2218   1.1053  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  5.07786    0.08696 99.90774    58.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If we look at the output we see under fixed effects the estimate of the intercept is reported, which is 5.078. As there are no predictor variables in this model, the intercept is the mean popularity score. By default, the lmer function does only give test statistics and estimates, but no p-values. However, because we use the lmerTest package we do get p-values for the fixed effects. We observe the intercept is significantly different from zero (p<2e-16).

We see under the random effects that the intercept variance on the class level 0.6945 and residual variance on the pupil level is 1.2218.

Model 2: First level predictors

Now we can add first (student) level predictors. The first level predictors are sex and extraversion. For now we assume their effects are fixed across the classes. The model at the pupil level is \[Popularity_{ij}=\beta_{0j}+\beta_1*sex_{ij}+ \beta_2*extraversion_{ij}+e_{ij}.\] We allow the random intercept to vary across classes but do not yet include class level predictor variables to explain part of the random intercept variability: \[\beta_{0j}=\gamma_{00}+u_{0j}.\] Furthermore, the slopes for sex and extraversion are fixed: \[\beta_1=\gamma_{10},\] \[\beta_2=\gamma_{20}.\] Combined we get: \[Popularity_{ij}= \gamma_{00}+\gamma_{10}*sex_{ij}+\gamma_{20}*extraversion_{ij}+u_{0j}+e_{ij}.\] There are three regression coefficients in this model: an intercept \(\gamma_{00}\), an effect of sex \(\gamma_{10}\) and an effect of extraversion \(\gamma_{20}\). Furthermore, there is a random effect \(u_{0j}\) at the class level and a random effect \(e_{ij}\) at the pupil level.

Before fitting the model we can plot the difference in effect for the two genders. We see there likely is an average difference between genders, but no difference in slopes (regression coefficients).

ggplot(data = popular2, 
       aes(x   = extrav,
           y   = popular, 
           col = as.factor(sex)))+
  geom_point(size     = 1, 
             alpha    = .7)+
  geom_smooth(method   = lm,
              se       = T, 
              size     = 1.5, 
              linetype = 1, 
              alpha    = .7)+
  theme_minimal()+
  labs(title    = "Linear Relationship Between Popularity and Extraversion for the 2 Genders", 
       subtitle = "The linear relationship between the two is similar for both genders, with a clear intercept difference")+
  scale_color_manual(name   =" Gender",
                     labels = c("Boys", "Girls"),
                     values = c("lightblue", "pink"))

The first level predictors sex and extraversion are now added to the fixed part of the model so that 1 in the part outside the brackets to the right of the ~ sign becomes 1+sex+extrav. Th R syntax of this model is

model2 <- lmer(formula = popular ~ 1 + sex + extrav + (1|class), 
               data    = popular2, REML = FALSE)
summary(model2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: popular ~ 1 + sex + extrav + (1 | class)
##    Data: popular2
## 
##      AIC      BIC   logLik deviance df.resid 
##     4944     4972    -2467     4934     1995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2113 -0.6578 -0.0048  0.6739  2.9771 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.6204   0.7876  
##  Residual             0.5915   0.7691  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 2.141e+00  1.170e-01 3.977e+02   18.31   <2e-16 ***
## sex         1.253e+00  3.741e-02 1.929e+03   33.50   <2e-16 ***
## extrav      4.415e-01  1.615e-02 1.959e+03   27.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) sex   
## sex    -0.100       
## extrav -0.706 -0.085

The intercept is now 2.14. This is the mean popularity score for a subject with sex=0 and extraversion=0. The regression coefficient for sex is 1.25, and the regression coefficient for extraversion 0.44 (note on notation: \(4.415e-01\) means \(4.415*10^{-1}\), which is equal to \(0.4415\)). In the last column of the fixed effects table of the output we see the p-values, which indicate all regression coefficients are significantly different from 0. Note that the results of this model are not given in the book.

Model 3: First and second level predictors

We now add the second level predictor variable teacher experience while keeping both sex and extraversion in the model (as both had significant effects). The pupil level model is \[Popularity_{ij}=\beta_{0j}+\beta_1*sex_{ij}+ \beta_2*extraversion_{ij}+e_{ij}.\] We now use teacher experience to explain part of the variability of the random intercept. The model at the class level is \[\beta_{0j}=\gamma_{00}+\gamma_{01}*experience_j+u_{0j}.\]
The slopes for sex and extraversion are fixed: \[\beta_1=\gamma_{10},\] \[\beta_2=\gamma_{20}.\] Combined we get: \[Popularity_{ij}= \gamma_{00}+\gamma_{10}*sex_{ij}+\gamma_{20}*extraversion_{ij}+\gamma_{01}*experience_j+u_{0j}+e_{ij}.\] This model includes four regression coefficients: an intercept \(\gamma_{00}\), an effect of sex \(\gamma_{10}\), an effect of extraversion \(\gamma_{20}\) and an effect of teacher experience \(\gamma_{01}\). Furthermore, there is a random effect \(u_{0j}\) at the class level and a random effect \(e_{ij}\) at the pupil level.

The R syntax for this model is

model3 <- lmer(popular ~ 1 + sex + extrav + texp + (1 | class), data=popular2, REML = FALSE)
summary(model3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: popular ~ 1 + sex + extrav + texp + (1 | class)
##    Data: popular2
## 
##      AIC      BIC   logLik deviance df.resid 
##   4874.3   4907.9  -2431.1   4862.3     1994 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1794 -0.6492 -0.0067  0.6708  3.0103 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.2888   0.5374  
##  Residual             0.5914   0.7690  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 8.093e-01  1.688e-01 2.338e+02   4.794 2.91e-06 ***
## sex         1.254e+00  3.726e-02 1.951e+03  33.653  < 2e-16 ***
## extrav      4.545e-01  1.615e-02 1.958e+03  28.134  < 2e-16 ***
## texp        8.841e-02  8.676e-03 1.038e+02  10.190  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) sex    extrav
## sex    -0.040              
## extrav -0.592 -0.090       
## texp   -0.801 -0.037  0.141

The results for the fixed effects show that both the level 1 and level 2 predictor variables have significant effects.

Model 4a: First and second level predictors with a random slope for sex

Now we also include the random slope for gender. In other words, the effect of gender is allowed to vary across the classes.

The model at the pupil level is \[Popularity_{ij}=\beta_{0j}+\beta_{1j}*sex_{ij}+ \beta_2*extraversion_{ij}+e_{ij}.\] We use teacher experience to explain part of the random intercept variability: \[\beta_{0j}=\gamma_{00}+\gamma_{01}*experience_j+u_{0j}.\]
The slope for sex is random but we do not yet include predictor variables to explain slope variability: \[\beta_{1j}=\gamma_{10}+u_{1j},\] and the model for the slope of extraversion is \[\beta_2=\gamma_{20}.\] Combined we get: \[Popularity_{ij}= \gamma_{00}+\gamma_{10}*sex_{ij}+\gamma_{20}*extraversion_{ij}+\gamma_{01}*experience_j+u_{0j}+u_{1j}*sex_{ij}+e_{ij}.\] This model includes the same four regression coefficients as the previous model. We now have three random effects: \(u_{0j}\) for the random intercept, \(u_{1j}\) for the random slope of sex and \(e_{ij}\) at the pupil level.

To fit the model we add the variables for which we want to add a random slope for sex to the random part of the input. This means that (1|class)becomes (1+sex|class).

model4a <- lmer(formula = popular ~ 1 + sex + extrav + texp + (1 + sex | class),
               data    = popular2, REML = FALSE)
## boundary (singular) fit: see help('isSingular')
summary(model4a)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: popular ~ 1 + sex + extrav + texp + (1 + sex | class)
##    Data: popular2
## 
##      AIC      BIC   logLik deviance df.resid 
##   4877.0   4921.8  -2430.5   4861.0     1992 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1682 -0.6626 -0.0115  0.6675  3.0485 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  class    (Intercept) 0.317532 0.56350       
##           sex         0.002204 0.04695  -1.00
##  Residual             0.590680 0.76856       
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 8.357e-01  1.696e-01 2.362e+02   4.927 1.57e-06 ***
## sex         1.252e+00  3.760e-02 1.027e+03  33.295  < 2e-16 ***
## extrav      4.548e-01  1.614e-02 1.956e+03  28.177  < 2e-16 ***
## texp        8.675e-02  8.633e-03 1.043e+02  10.049  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) sex    extrav
## sex    -0.082              
## extrav -0.593 -0.088       
## texp   -0.797 -0.037  0.146
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

First of all we observe a warning was generated. Apparently, the model with a random slope for sex is too complicated for the data that we have.

Model 4b: First and second level predictors with random slope for extraversion

We now fit a model with a random slope for extraversion. The model at the pupil level is \[Popularity_{ij}=\beta_{0j}+\beta_1*sex_{ij}+ \beta_{2j}*extraversion_{ij}+e_{ij}.\] We use teacher experience to explain part of the random intercept variability: \[\beta_{0j}=\gamma_{00}+\gamma_{01}*experience_j+u_{0j}.\]
The slope for sex is fixed: \[\beta_1=\gamma_{10}.\] The slope of extraversion is random but we do not yet include any class-level predictors to explain slope variability: \[\beta_{2j}=\gamma_{20}+u_{2j}.\] Combined we get: \[Popularity_{ij}= \gamma_{00}+\gamma_{10}*sex_{ij}+\gamma_{20}*extraversion_{ij}+\gamma_{01}*experience_j+u_{0j}+u_{2j}*extraversion_{ij}+e_{ij}.\] This model includes the same four regression coefficients as the previous model. We now have three random effects: \(u_{0j}\) for the random intercept, \(u_{2j}\) for the random slope of extraversion and \(e_{ij}\) at the pupil level. The random part is now (1 + extrav |class).

model4b <- lmer(formula = popular ~ 1 + sex + extrav + texp + (1 + extrav |class), 
               data    = popular2, REML = FALSE)
summary(model4b)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: popular ~ 1 + sex + extrav + texp + (1 + extrav | class)
##    Data: popular2
## 
##      AIC      BIC   logLik deviance df.resid 
##   4828.8   4873.6  -2406.4   4812.8     1992 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1816 -0.6467 -0.0215  0.6640  2.9686 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  class    (Intercept) 1.28051  1.1316        
##           extrav      0.03393  0.1842   -0.89
##  Residual             0.55180  0.7428        
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 7.377e-01  1.952e-01 1.851e+02    3.78 0.000212 ***
## sex         1.252e+00  3.655e-02 1.916e+03   34.26  < 2e-16 ***
## extrav      4.526e-01  2.448e-02 9.850e+01   18.49  < 2e-16 ***
## texp        9.087e-02  8.598e-03 1.038e+02   10.57  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) sex    extrav
## sex    -0.031              
## extrav -0.719 -0.058       
## texp   -0.687 -0.040  0.087

We see that no warning is given for this model.

The fixed effects show that all predictors have significant effects.

Model 5: First and second level predictors with a random slope for extraversion and a cross-level interaction

The previous model showed that the effect of extraversion varies across the classes. As a final step, we can add a cross-level interaction between teacher experience and extraversion. In other words, we want to investigate whether the differences in the relation between extraversion and popularity across the classes can be explained by the experience of the teachers. In this next step to reproduce Model M2 from Table 2.3, we add the cross-level interaction between extraversion and teacher experience. This means we have to add teacher experience as a predictor for the coefficient of extraversion.

The model at the pupil level is: \[Popularity_{ij}=\beta_{0j}+\beta_1*sex_{ij}+ \beta_{2j}*extraversion_{ij}+e_{ij},\] in which we use teacher experience to predict the random intercept \[\beta_{0j}=\gamma_{00}+\gamma_{01}*experience_j+u_{0j},\]
have a fixed slope for sex \[\beta_1=\gamma_{10},\] and use teacher experience to predict the random slope for extraversion. \[\beta_{2j}=\gamma_{20}+\gamma_{21}*experience_j+u_{2j}.\] Combined we get: \[Popularity_{ij}= \gamma_{00}+\gamma_{10}*sex_{ij}+\gamma_{20}*extraversion_{ij}+\gamma_{01}*experience_j+\gamma_{21}*extraversion_{ij}*experience_j+u_{2j}*extraversion_{ij}+u_{0j}+e_{ij}.\] This model includes five regression coefficients: an intercept \(\gamma_{00}\), an effect of sex \(\gamma_{10}\), an effect of extraversion \(\gamma_{20}\), an effect of teacher experience \(\gamma_{01}\), and a cross-level interaction effect \(\gamma_{21}\) for extraversion with teacher experience. We have three random effects: \(u_{0j}\) for the random intercept, \(u_{2j}\) for the random slope of extraversion and \(e_{ij}\) at the pupil level.

The cross-level interaction term between extraversion and teacher experience can be created by the ‘:’ sign or by multiplying the terms. The R syntax for this model is

model5<-lmer(formula = popular ~ 1 + sex + extrav + texp+ extrav:texp + (1 + extrav | class), 
             data    = popular2)
summary(model5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: popular ~ 1 + sex + extrav + texp + extrav:texp + (1 + extrav |  
##     class)
##    Data: popular2
## 
## REML criterion at convergence: 4780.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.12872 -0.63857 -0.01129  0.67916  3.05006 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  class    (Intercept) 0.478639 0.69184       
##           extrav      0.005409 0.07355  -0.64
##  Residual             0.552769 0.74348       
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -1.210e+00  2.719e-01  1.093e+02  -4.449 2.09e-05 ***
## sex          1.241e+00  3.623e-02  1.941e+03  34.243  < 2e-16 ***
## extrav       8.036e-01  4.012e-02  7.207e+01  20.031  < 2e-16 ***
## texp         2.262e-01  1.681e-02  9.851e+01  13.458  < 2e-16 ***
## extrav:texp -2.473e-02  2.555e-03  7.199e+01  -9.679 1.15e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sex    extrav texp  
## sex          0.002                     
## extrav      -0.867 -0.065              
## texp        -0.916 -0.047  0.801       
## extrav:texp  0.773  0.033 -0.901 -0.859

The interaction term is denoted by extrav:texp under Fixed effects and is estimated at -0.0247. All regression coefficients in this model are significant.

To interpret the fixed part of the model we plug in the estimated values of the regresssion coefficients: \[Popularity_{ij}= -1.21+1.241*sex_{ij}+0.8036*extraversion_{ij}+0.2262*experience_j-0.0247*extraversion_{ij}*experience_j.\] This equation can be rewritten such that the effect of extraversion becomes a function of teacher experience: \[Popularity_{ij}= -1.21+1.241*sex_{ij}+0.2262*experience_j+(0.8036-0.0247*experience_j)*extraversion_{ij}.\]

The regression coefficients can now be interpreted as follows:

Intercept = -1.21. This is the average popularity score for boys with an extraversion equal to zero and whose teacher has zero years experience. It should be noted that the latter two predictor variables do not have the value zero in the data set.

Effect of sex=1.241. On the average females have a populatiry score 1.241 higher than the boys (for same values of extraversion and teacher experience).

Main effect of teacher experience=0.2262. This effect is interpreted for extraversion=0: for each extra year of teacher experience the popularity score increases by 0.2262 (keeping sex constant).

Effect of extraversion as a function of teacher experience: 0.8036-0.0247*experience. This effect is equal to 0.8036 if teacher experience=0 and decreases by 0.0247 for each extra year of teacher experience. In other words: the effect of extraversion becomes weaker when teacher experiences increases.

Graphical representation of final model

In a plot we can also clearly see that years of teacher experience has an influence on both the intercept and the regression coefficient of extraversion on popularity.

ggplot(data = popular2,
       aes(x = extrav, 
           y = popular, 
           col = as.factor(texp)))+
  viridis::scale_color_viridis(discrete = TRUE)+
  geom_point(size     = .7,
             alpha    = .8)+
  geom_smooth(method = lm,
              se     = FALSE,
              size   = 2,
              alpha  = .8)+
  theme_minimal()+
  labs(title    = "Linear Relationship for Different Years of Teacher Experience as Observed", 
       subtitle = "The linear relationship between the two is not the same for all classes", 
       col      = "Years of\nTeacher\nExperience")
## `geom_smooth()` using formula = 'y ~ x'

The intercept is higher if teacher experience is higher. The effect of extraversion is lower if teacher experience is higher.