> print("R Descriptive Statistics from School-Level-2 Data")
[1] "R Descriptive Statistics from School-Level-2 Data"
> univSummary(data = Example2Means, vars = c("NperschoolID", "CM_math", 
+ "size", "private"))

 Variable       N     Mean      SD   Variance     Min      Max
 NperschoolID 160   44.906  11.855    140.538  14.000   67.000
 CM_math      160   12.621   3.118      9.720   4.240   19.719
 size         160 1097.825 629.506 396278.347 100.000 2713.000
 private      160    0.438   0.498      0.248   0.000    1.000

> print("R Descriptive Statistics from Student-Level-1 Data")
[1] "R Descriptive Statistics from Student-Level-1 Data"
> c("NperschoolID", "CM_math", "size", "private")
[1] "NperschoolID" "CM_math"      "size"         "private"     
> univSummary(data = Example2, vars = c("math", "size", "private"))

 Variable    N     Mean      SD   Variance     Min      Max
 math     7185   12.748   6.878     47.310  -2.832   24.993
 size     7185 1056.862 604.172 365024.409 100.000 2713.000
 private  7185    0.493   0.500      0.250   0.000    1.000

> Example2$size100 = (Example2$size - 1000)/100
> print("Model 0: Single-Level Empty Means, NO Random Intercept")
[1] "Model 0: Single-Level Empty Means, NO Random Intercept"
> Model0 = gls(data = Example2, model = math ~ 1)
> glsSummary(Model0, CI = TRUE, digits = 3)

Model Fit for Predicting math:
       AIC       BIC     logLik -2*logLik
 48106.917 48120.676 -24051.459 48102.917

Residual Variance (sigma^2):
47.310 

Fixed Effects (with 95% CI):
 Term          Est    SE   DF t.value p.value  Lower  Upper
 Intercept  12.748 0.081 7184 157.099   <.001 12.589 12.907

> print("Model 1: Two-Level Empty Means, WITH Random Intercept")
[1] "Model 1: Two-Level Empty Means, WITH Random Intercept"
> Model1 = lmer(data = Example2, REML = TRUE, formula = math ~ 
+ 1 + (1 | schoolID))
> print("Show results with all relevant empty model options turned on")
[1] "Show results with all relevant empty model options turned on"
> lmerSummary(Model1, CI = TRUE, ICC1 = TRUE, LRT = TRUE, ICC2 = TRUE, 
+ randCI = TRUE, Neff = TRUE)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 47122.793  47143.433 -23558.397  47116.793 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    8.614   2.935
 Residual              39.148   6.257
Number of obs: 7185 , groups: schoolID, 160 

Fixed Effects (with 95% CI):
 Term          Est    SE      DF t.value p.value  Lower  Upper
 Intercept  12.637 0.244 156.647  51.707   <.001 12.154 13.120

ICC1 intraclass correlation:
    Adjusted ICC: 0.180 
  Unadjusted ICC: 0.180 

Design Effect and Effective Sample Size:
  Design Effect (using mean n = 44.906 per cluster): 8.919 
  Effective N: 805.622 

LRT for random effect(s):
               Npar     LL  -2LL   AIC   BIC     LRT DF p.value p.mix
<none>            3 -23558 47117 47123 47143                         
(1 | schoolID)    2 -24051 48103 48107 48121 986.124  1   <.001 <.001

ICC2 reliability of the cluster means (using mean n = 44.906 per cluster):
0.908 

Random Effect 95% CIs (Fixed Effect +/- 1.96*SD(random effect variance)):
 Term      FixedEffect Variance    SD Lower  Upper
 Intercept      12.637    8.614 2.935 6.884 18.390

> print("Model 2 using lme function from nlme package instead of lmer from lme4 package")
[1] "Model 2 using lme function from nlme package instead of lmer from lme4 package"
> print("nlme has more options for residual correlation and heterogeneous variance")
[1] "nlme has more options for residual correlation and heterogeneous variance"
> print("but doesn’t work as well with Satterthwaite DDF as lmerTest")
[1] "but doesn’t work as well with Satterthwaite DDF as lmerTest"
> Model2lme = lme(data = Example2, method = "REML", fixed = math ~ 
+ 1 + size100 + private, random = ~1 | schoolID)
> lmeSummary(Model2lme, CI = TRUE, effectSizes = TRUE)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 47092.313  47126.709 -23541.156  47082.313 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    6.618   2.573
 Residual              39.150   6.257
Number of obs: 7185 , groups: schoolID , 160 

Overall F-test for All Fixed Slopes:
 Model     NumDF   DenDF F.value p.value
 Model2lme     2 155.903  21.803   <.001

Fixed Effects (with 95% CI):
 Term          Est    SE      DF t.value p.value  Lower  Upper
 Intercept  11.182 0.321 155.973  34.827   <.001 10.548 11.816
 size100     0.061 0.039 155.903   1.574   0.117 -0.016  0.138
 private     3.151 0.490 152.073   6.436   <.001  2.184  4.119

Effect Sizes for Fixed Slopes:
 Term    Estimate    SE p.value Cohen.d Partial.r
 size100    0.061 0.039   0.117   0.252     0.125
 private    3.151 0.490   <.001   1.044     0.463

> print("Model 2: Add Main Effects of Level-2 School Size and Public vs Private")
[1] "Model 2: Add Main Effects of Level-2 School Size and Public vs Private"
> Model2 = lmer(data = Example2, REML = TRUE, formula = math ~ 
+ 1 + size100 + private + (1 | schoolID))
> lmerSummary(Model2, CI = TRUE, effectSizes = TRUE)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 47092.313  47126.711 -23541.156  47082.313 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    6.618   2.573
 Residual              39.150   6.257
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model2     2 154.470  21.803   <.001

Fixed Effects (with 95% CI):
 Term          Est    SE      DF t.value p.value  Lower  Upper
 Intercept  11.182 0.321 156.320  34.827   <.001 10.548 11.816
 size100     0.061 0.039 156.070   1.574   0.118 -0.016  0.138
 private     3.151 0.490 152.059   6.436   <.001  2.184  4.119

Effect Sizes for Fixed Slopes:
 Term    Estimate    SE p.value Cohen.d Partial.r
 size100    0.061 0.039   0.118   0.252     0.125
 private    3.151 0.490   <.001   1.044     0.463

> PseudoR2(baseModel = Model1, altModel1 = Model2)

Pseudo-R2
 Term       Model  Variance Pseudo.R2
 Intercept  Model1    8.614          
 Intercept  Model2    6.618     0.232
 Residual   Model1   39.148          
 Residual   Model2   39.150    -0.000

> print("Rights & Sterba R2 suite -- f2 is for random intercept")
[1] "Rights & Sterba R2 suite -- f2 is for random intercept"
> r2mlm::r2mlm(model = Model2, bargraph = FALSE)
$Decompositions
                   total within between
fixed, within   0.000000      0      NA
fixed, between  0.044375     NA 0.24307
slope variation 0.000000      0      NA
mean variation  0.138188     NA 0.75693
sigma2          0.817437      1      NA

$R2s
       total within between
f1  0.000000      0      NA
f2  0.044375     NA 0.24307
v   0.000000      0      NA
m   0.138188     NA 0.75693
f   0.044375     NA      NA
fv  0.044375      0      NA
fvm 0.182563     NA      NA

> TotalR2(altModel1 = Model2)

Total R2
 Model  Total.R2
 Model2    0.043

> print("Model 3: Add Interaction of Level-2 School Size and Public vs Private")
[1] "Model 3: Add Interaction of Level-2 School Size and Public vs Private"
> Model3 = lmer(data = Example2, REML = TRUE, formula = math ~ 
+ 1 + size100 * private + (1 | schoolID))
> lmerSummary(Model3, CI = TRUE)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 47095.398  47136.677 -23541.699  47083.398 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    6.569   2.563
 Residual              39.151   6.257
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model3     3 151.894  15.261   <.001

Fixed Effects (with 95% CI):
 Term                Est    SE      DF t.value p.value  Lower  Upper
 Intercept        11.288 0.329 155.264  34.299   <.001 10.637 11.938
 size100           0.031 0.045 157.645   0.685   0.495 -0.058  0.119
 private           3.252 0.494 150.004   6.589   <.001  2.277  4.227
 size100:private   0.124 0.090 149.283   1.377   0.171 -0.054  0.302

> print("Size slope for private schools = size100 + size100:private")
[1] "Size slope for private schools = size100 + size100:private"
> Model3con1D = contest1D(Model3, L = c(0, 1, 0, 1))
> con1DSummary(Model3con1D, CI = TRUE, effectSizes = TRUE)

Fixed Effect (with 95% CI):
   Est    SE      DF t.value p.value Lower Upper
 0.155 0.078 146.678   1.979   0.050 0.000 0.309

Effect Size:
   Est    SE p.value Cohen.d Partial.r
 0.155 0.078   0.050   0.327     0.161

> PseudoR2(baseModel = Model1, altModel1 = Model2, altModel2 = Model3)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model  Variance Pseudo.R2 Change.R2
 Intercept  Model1    8.614                    
 Intercept  Model2    6.618     0.232          
 Intercept  Model3    6.569     0.237     0.006
 Residual   Model1   39.148                    
 Residual   Model2   39.150    -0.000          
 Residual   Model3   39.151    -0.000    -0.000

> print("Rights & Sterba R2 and change -- f2 is for random intercept")
[1] "Rights & Sterba R2 and change -- f2 is for random intercept"
> r2mlm::r2mlm_comp(modelA = Model2, modelB = Model3, bargraph = FALSE)
$`Model A R2s`
       total within between
f1  0.000000      0      NA
f2  0.044375     NA 0.24307
v   0.000000      0      NA
m   0.138188     NA 0.75693
f   0.044375     NA      NA
fv  0.044375      0      NA
fvm 0.182563     NA      NA

$`Model B R2s`
       total within between
f1  0.000000      0      NA
f2  0.046869     NA 0.25497
v   0.000000      0      NA
m   0.136954     NA 0.74503
f   0.046869     NA      NA
fv  0.046869      0      NA
fvm 0.183823     NA      NA

$`R2 differences, Model B - Model A`
         total within   between
f1   0.0000000      0        NA
f2   0.0024941     NA  0.011901
v    0.0000000      0        NA
m   -0.0012336     NA -0.011901
f    0.0024941     NA        NA
fv   0.0024941      0        NA
fvm  0.0012605     NA        NA

> TotalR2(altModel1 = Model2, altModel2 = Model3)

Total-R2 and Change in Total-R2
 Model  Total.R2 Change.R2
 Model2    0.043          
 Model3    0.045     0.002

