# load them
pacman::p_load(mlmRev, tidyverse, lme4, merTools)
# data from mlmRev package
data(Hsb82)
#
str(Hsb82)
## 'data.frame':    7185 obs. of  8 variables:
##  $ school : Ord.factor w/ 160 levels "8367"<"8854"<..: 59 59 59 59 59 59 59 59 59 59 ...
##  $ minrty : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sx     : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 2 1 ...
##  $ ses    : num  -1.528 -0.588 -0.528 -0.668 -0.158 ...
##  $ mAch   : num  5.88 19.71 20.35 8.78 17.9 ...
##  $ meanses: num  -0.434 -0.434 -0.434 -0.434 -0.434 ...
##  $ sector : Factor w/ 2 levels "Public","Catholic": 1 1 1 1 1 1 1 1 1 1 ...
##  $ cses   : num  -1.0936 -0.1536 -0.0936 -0.2336 0.2764 ...
# make a copy of data to the dta object 
dta <- Hsb82  

names(dta) <- c("School", "Minority", "Gender", "SES", "Math",
                "Ave_SES", "Sector", "C_SES")

theme_set(theme_bw())
# Math vs SES by sector
ggplot(dta, aes(C_SES, Math, color = Sector)) +
 geom_smooth(method = "lm") +
 geom_jitter(alpha = 0.2) +
 labs(x = "Centered SES", y = "Math achievement score") +
 theme(legend.position = c(.1, .8))
## `geom_smooth()` using formula 'y ~ x'

to facilitate comparison with sas output

#options(contrasts = c(factor = “contr.SAS”, ordered = “contr.poly”))

# full
summary(m4a <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
        Sector:C_SES + (C_SES | School), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +  
##     (C_SES | School)
##    Data: dta
## 
## REML criterion at convergence: 46503.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.15926 -0.72319  0.01704  0.75444  2.95822 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  School   (Intercept)  2.380   1.5426       
##           C_SES        0.101   0.3179   0.39
##  Residual             36.721   6.0598       
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           12.1279     0.1993  60.856
## Ave_SES                5.3329     0.3692  14.446
## SectorCatholic         1.2266     0.3063   4.005
## C_SES                  2.9450     0.1556  18.928
## Ave_SES:C_SES          1.0393     0.2989   3.477
## SectorCatholic:C_SES  -1.6427     0.2398  -6.851
## 
## Correlation of Fixed Effects:
##             (Intr) Av_SES SctrCt C_SES  A_SES:
## Ave_SES      0.256                            
## SectorCthlc -0.699 -0.356                     
## C_SES        0.075  0.019 -0.053              
## A_SES:C_SES  0.019  0.074 -0.026  0.293       
## SctrC:C_SES -0.052 -0.027  0.077 -0.696 -0.351
# reduced
summary(m4b <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
        Sector:C_SES + (1 | School), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +  
##     (1 | School)
##    Data: dta
## 
## REML criterion at convergence: 46504.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1701 -0.7249  0.0148  0.7542  2.9655 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  School   (Intercept)  2.375   1.541   
##  Residual             36.766   6.064   
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           12.1282     0.1992  60.885
## Ave_SES                5.3367     0.3690  14.463
## SectorCatholic         1.2245     0.3061   4.000
## C_SES                  2.9421     0.1512  19.457
## Ave_SES:C_SES          1.0444     0.2910   3.589
## SectorCatholic:C_SES  -1.6422     0.2331  -7.045
## 
## Correlation of Fixed Effects:
##             (Intr) Av_SES SctrCt C_SES  A_SES:
## Ave_SES      0.256                            
## SectorCthlc -0.699 -0.356                     
## C_SES        0.000  0.000  0.000              
## A_SES:C_SES  0.000  0.000  0.000  0.295       
## SctrC:C_SES  0.000  0.000  0.000 -0.696 -0.351
# will refit be ML estimation
anova(m4a, m4b)
## refitting model(s) with ML (instead of REML)
## Data: dta
## Models:
## m4b: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES + (1 | School)
## m4a: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES + (C_SES | School)
##     npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## m4b    8 46513 46568 -23249    46497                     
## m4a   10 46516 46585 -23248    46496 1.0016  2     0.6061
# residualual plot
plot(m4b, xlab = "Fitted values", ylab = "Pearson residuals", 
     pch = 20, cex = .5, type = c("p", "g"))

# quick summary of model parameter estimates
fastdisp(m4b)
## lmer(formula = Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + 
##     Sector:C_SES + (1 | School), data = dta)
##                      coef.est coef.se
## (Intercept)          12.13     0.20  
## Ave_SES               5.34     0.37  
## SectorCatholic        1.22     0.31  
## C_SES                 2.94     0.15  
## Ave_SES:C_SES         1.04     0.29  
## SectorCatholic:C_SES -1.64     0.23  
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  School   (Intercept) 1.54    
##  Residual             6.06    
## ---
## number of obs: 7185, groups: School, 160
## AIC = 46520.8
# estimate fixed-effects parameters by simulation
m4b_fe <- FEsim(m4b, 1000)
# plot for fixed effects
plotFEsim(m4b_fe) + 
 labs(title = "Coefficient Plot", x = "Median Effect Estimate")

# estimate random-effects parameters by simulation
m4b_re <- REsim(m4b)
# normality plot for random effects
plotREsim(m4b_re)

The end