txt> print("Descriptive Statistics in Child-Level-1 Data")
[1] "Descriptive Statistics in Child-Level-1 Data"
txt> univSummary(data = Example6a, vars = c("achieve", "ses", "WPses", 
txt+ "WSses", "WPSses"))

 Variable    N   Mean    SD Variance    Min   Max Median
 achieve  1000  6.343 0.868    0.753  3.900 9.900  6.300
 ses      1000  4.098 1.398    1.954  1.000 6.000  4.000
 WPses    1000 -0.000 1.371    1.879 -3.684 2.529  0.200
 WSses    1000  0.000 1.374    1.888 -3.412 2.529  0.217
 WPSses   1000  0.000 0.793    0.629 -3.333 2.333  0.000

txt> corrSummary(data = Example6a, p.separate = TRUE, vars = c("achieve", 
txt+ "PMachieve", "SMachieve", "PSMachieve", "ses", "PMses", "SMses", 
txt+ "PSMses", "WPses", "WSses", "WPSses"))

Correlations:
 Variable   achieve PMachieve SMachieve PSMachieve    ses  PMses  SMses PSMses  WPses  WSses WPSses
 achieve      1.000     0.521     0.357      0.885  0.174 -0.013  0.103  0.155  0.180  0.158  0.082
 PMachieve    0.521     1.000     0.084      0.588 -0.005 -0.025 -0.020 -0.006 -0.000 -0.001  0.000
 SMachieve    0.357     0.084     1.000      0.403  0.053 -0.026  0.288  0.064  0.059 -0.000  0.000
 PSMachieve   0.885     0.588     0.403      1.000  0.144 -0.015  0.116  0.175  0.150  0.125  0.000
 ses          0.174    -0.005     0.053      0.144  1.000  0.197  0.184  0.824  0.980  0.983  0.567
 PMses       -0.013    -0.025    -0.026     -0.015  0.197  1.000  0.015  0.239 -0.000  0.197  0.000
 SMses        0.103    -0.020     0.288      0.116  0.184  0.015  1.000  0.224  0.185  0.000 -0.000
 PSMses       0.155    -0.006     0.064      0.175  0.824  0.239  0.224  1.000  0.792  0.796 -0.000
 WPses        0.180    -0.000     0.059      0.150  0.980 -0.000  0.185  0.792  1.000  0.963  0.578
 WSses        0.158    -0.001    -0.000      0.125  0.983  0.197  0.000  0.796  0.963  1.000  0.577
 WPSses       0.082     0.000     0.000      0.000  0.567  0.000 -0.000 -0.000  0.578  0.577  1.000

P-values:
 Variable   achieve PMachieve SMachieve PSMachieve   ses PMses SMses PSMses WPses WSses WPSses
 achieve                <.001     <.001      <.001 <.001 0.683 0.001  <.001 <.001 <.001  0.009
 PMachieve    <.001               0.008      <.001 0.877 0.432 0.520  0.851 1.000 0.971  1.000
 SMachieve    <.001     0.008                <.001 0.093 0.412 <.001  0.042 0.061 1.000  1.000
 PSMachieve   <.001     <.001     <.001            <.001 0.644 <.001  <.001 <.001 <.001  1.000
 ses          <.001     0.877     0.093      <.001       <.001 <.001  <.001 <.001 <.001  <.001
 PMses        0.683     0.432     0.412      0.644 <.001       0.626  <.001 1.000 <.001  1.000
 SMses        0.001     0.520     <.001      <.001 <.001 0.626        <.001 <.001 1.000  1.000
 PSMses       <.001     0.851     0.042      <.001 <.001 <.001 <.001        <.001 <.001  1.000
 WPses        <.001     1.000     0.061      <.001 <.001 1.000 <.001  <.001       <.001  <.001
 WSses        <.001     0.971     1.000      <.001 <.001 <.001 1.000  <.001 <.001        <.001
 WPSses       0.009     1.000     1.000      1.000 <.001 1.000 1.000  1.000 <.001 <.001       

txt> print("Number of Children per Unique Primary/Secondary Combination")
[1] "Number of Children per Unique Primary/Secondary Combination"
txt> freqSummary(data = Example6a[!duplicated(Example6a$UniqueID), 
txt+ ], vars = "NperUniqueID")

 Variable     Level   N Proportion
 NperUniqueID 1     394      0.604
 NperUniqueID 2     185      0.284
 NperUniqueID 3      60      0.092
 NperUniqueID 4      10      0.015
 NperUniqueID 5       2      0.003
 NperUniqueID 6       1      0.002

txt> print("Model 1a: Empty Means, Random Intercept for Primary School")
[1] "Model 1a: Empty Means, Random Intercept for Primary School"
txt> Model1a = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + (1 | PschoolID))
txt> lmerSummary(Model1a, ICC1 = TRUE, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2393.736  2408.460 -1193.868  2387.736 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 PschoolID  Intercept    0.176   0.419
 Residual                0.577   0.760
Number of obs: 1000 , groups: PschoolID, 50 

Fixed Effects:
 Term         Est    SE     DF t.value p.value
 Intercept  6.359 0.064 49.241  99.142   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.233 
  Unadjusted ICC: 0.233 

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                    3 -1194 2388 2394 2408                         
 Intercept | PschoolID     2 -1279 2558 2562 2572 170.632  1   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                  FixedEffect Variance    SD Lower Upper
 Intercept | PschoolID       6.359    0.176 0.419 5.538 7.180

txt> print("Model 1b: Empty Means, Crossed Random Intercepts for Primary and Secondary School")
[1] "Model 1b: Empty Means, Crossed Random Intercepts for Primary and Secondary School"
txt> Model1b = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + (1 | PschoolID) + (1 | SschoolID))
txt> lmerSummary(Model1b, ICC1 = TRUE, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2320.269  2339.900 -1156.134  2312.269 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 PschoolID  Intercept    0.175   0.418
 SschoolID  Intercept    0.079   0.281
 Residual                0.505   0.711
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Fixed Effects:
 Term         Est    SE     DF t.value p.value
 Intercept  6.342 0.081 66.879  78.468   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.334 
  Unadjusted ICC: 0.334 

Proportion of Variance by Source (adj = adjusted after fixed effects):
 Source        Variance Prop.Adj Prop.Unadj
 PschoolID        0.175    0.230      0.230
 SschoolID        0.079    0.104      0.104
 Residual         0.505    0.666      0.666
 Fixed Effects    0.000               0.000

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                    4 -1156 2312 2320 2340                         
 Intercept | PschoolID     3 -1250 2501 2507 2521 188.309  1   <.001 <.001
 Intercept | SschoolID     3 -1194 2388 2394 2408  75.468  1   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                  FixedEffect Variance    SD Lower Upper
 Intercept | PschoolID       6.342    0.175 0.418 5.523 7.161
 Intercept | SschoolID       6.342    0.079 0.281 5.792 6.892

txt> print("Model 1c: Model 1b + Random Intercept for Unique Primary/Secondary Combination")
[1] "Model 1c: Model 1b + Random Intercept for Unique Primary/Secondary Combination"
txt> Model1c = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + (1 | PschoolID) + (1 | SschoolID) + (1 | UniqueID))
txt> lmerSummary(Model1c, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2321.284  2345.822 -1155.642  2311.284 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 UniqueID   Intercept    0.027   0.164
 PschoolID  Intercept    0.174   0.417
 SschoolID  Intercept    0.078   0.279
 Residual                0.480   0.693
Number of obs: 1000 , groups: UniqueID, 652; PschoolID, 50; SschoolID, 33 

Fixed Effects:
 Term         Est    SE     DF t.value p.value
 Intercept  6.344 0.081 66.345  78.627   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.367 
  Unadjusted ICC: 0.367 

Proportion of Variance by Source (adj = adjusted after fixed effects):
 Source        Variance Prop.Adj Prop.Unadj
 UniqueID         0.027    0.036      0.036
 PschoolID        0.174    0.229      0.229
 SschoolID        0.078    0.103      0.103
 Residual         0.480    0.633      0.633
 Fixed Effects    0.000               0.000

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                    5 -1156 2311 2321 2346                         
 Intercept | PschoolID     4 -1232 2465 2473 2492 153.233  1   <.001 <.001
 Intercept | SschoolID     4 -1188 2376 2384 2404  64.961  1   <.001 <.001
 Intercept | UniqueID      4 -1156 2312 2320 2340   0.985  1   0.321 0.160

txt> print("Empty Means, Crossed Random Intercepts for Child SES Level-1 Predictor")
[1] "Empty Means, Crossed Random Intercepts for Child SES Level-1 Predictor"
txt> EmptySES = lmer(data = Example6a, REML = TRUE, formula = ses ~ 
txt+ 1 + (1 | PschoolID) + (1 | SschoolID) + (1 | UniqueID))
txt> lmerSummary(EmptySES, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting ses:
      AIC       BIC    logLik -2*logLik 
 3521.238  3545.777 -1755.619  3511.238 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 UniqueID   Intercept    0.032   0.179
 PschoolID  Intercept    0.000   0.000
 SschoolID  Intercept    0.000   0.000
 Residual                1.922   1.386
Number of obs: 1000 , groups: UniqueID, 652; PschoolID, 50; SschoolID, 33 

Fixed Effects:
 Term         Est    SE      DF t.value p.value
 Intercept  4.098 0.045 469.461  92.004   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.016 
  Unadjusted ICC: 0.016 

Proportion of Variance by Source (adj = adjusted after fixed effects):
 Source        Variance Prop.Adj Prop.Unadj
 UniqueID         0.032    0.016      0.016
 PschoolID        0.000    0.000      0.000
 SschoolID        0.000    0.000      0.000
 Residual         1.922    0.984      0.984
 Fixed Effects    0.000               0.000

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC    LRT DF p.value p.mix
 <none>                    5 -1756 3511 3521 3546                        
 Intercept | PschoolID     4 -1756 3511 3519 3539 -0.000  1   1.000 1.000
 Intercept | SschoolID     4 -1756 3511 3519 3539 -0.000  1   1.000 1.000
 Intercept | UniqueID      4 -1756 3511 3519 3539  0.096  1   0.756 0.378

txt> print("Model 2: Model 1b + Primary by Secondary School Denomination")
[1] "Model 2: Model 1b + Primary by Secondary School Denomination"
txt> Model2 = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + Pdenom + Sdenom + Pdenom:Sdenom + (1 | PschoolID) + (1 | 
txt+ 
txt+ SschoolID))
txt> lmerSummary(Model2)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2326.533  2360.887 -1156.267  2312.533 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 PschoolID  Intercept    0.172   0.415
 SschoolID  Intercept    0.073   0.269
 Residual                0.504   0.710
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Overall F-test for All Fixed Slopes:
 Model  NumDF   DenDF F.value p.value
 Model2     3 128.976   2.450   0.066

Fixed Effects:
 Term             Est    SE      DF t.value p.value
 Intercept      6.167 0.137  78.411  44.862   <.001
 Pdenom         0.073 0.144  75.523   0.506   0.614
 Sdenom         0.101 0.126  45.008   0.798   0.429
 Pdenom:Sdenom  0.161 0.101 949.193   1.599   0.110

txt> PseudoR2(baseModel = Model1b, altModel1 = Model2)

Pseudo-R2
 Term                   Model   Variance Pseudo.R2
 Intercept | PschoolID  Model1b    0.175          
 Intercept | PschoolID  Model2     0.172     0.015
 Intercept | SschoolID  Model1b    0.079          
 Intercept | SschoolID  Model2     0.073     0.079
 Residual               Model1b    0.505          
 Residual               Model2     0.504     0.002

txt> print("Model 3a: Add Child SES (Smushed, Centered at 4)")
[1] "Model 3a: Add Child SES (Smushed, Centered at 4)"
txt> Model3a = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + Pdenom + Sdenom + Pdenom:Sdenom + ses4 + (1 | PschoolID) + 
txt+ 
txt+ (1 | SschoolID))
txt> lmerSummary(Model3a)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2292.859  2332.121 -1138.429  2276.859 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 PschoolID  Intercept    0.173   0.415
 SschoolID  Intercept    0.068   0.260
 Residual                0.483   0.695
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model3a     4 162.926  12.661   <.001

Fixed Effects:
 Term             Est    SE      DF t.value p.value
 Intercept      6.171 0.135  79.189  45.568   <.001
 Pdenom         0.052 0.144  74.533   0.359   0.721
 Sdenom         0.073 0.123  45.283   0.592   0.557
 ses4           0.107 0.016 939.146   6.556   <.001
 Pdenom:Sdenom  0.199 0.099 947.930   2.016   0.044

txt> PseudoR2(baseModel = Model1b, altModel1 = Model2, altModel2 = Model3a)

Pseudo-R2 and Change in Pseudo-R2
 Term                   Model   Variance Pseudo.R2 Change.R2
 Intercept | PschoolID  Model1b    0.175                    
 Intercept | PschoolID  Model2     0.172     0.015          
 Intercept | PschoolID  Model3a    0.173     0.012    -0.003
 Intercept | SschoolID  Model1b    0.079                    
 Intercept | SschoolID  Model2     0.073     0.079          
 Intercept | SschoolID  Model3a    0.068     0.138     0.060
 Residual               Model1b    0.505                    
 Residual               Model2     0.504     0.002          
 Residual               Model3a    0.483     0.043     0.042

txt> print("Model 3b: Add Contextual Effects of Child SES (Each Centered at 4)")
[1] "Model 3b: Add Contextual Effects of Child SES (Each Centered at 4)"
txt> Model3b = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + Pdenom + Sdenom + Pdenom:Sdenom + ses4 + PMses4 + SMses4 + 
txt+ 
txt+ PSMses4 + (1 | PschoolID) + (1 | SschoolID))
txt> lmerSummary(Model3b)

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2304.899  2358.884 -1141.450  2282.899 

Model for the Variance:
 Groups     Name      Variance Std.Dev
 PschoolID  Intercept    0.177   0.421
 SschoolID  Intercept    0.066   0.258
 Residual                0.484   0.695
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Overall F-test for All Fixed Slopes:
 Model   NumDF  DenDF F.value p.value
 Model3b     7 95.393   7.468   <.001

Fixed Effects:
 Term              Est    SE      DF t.value p.value
 Intercept       6.166 0.138  77.484  44.592   <.001
 Pdenom          0.050 0.145  72.328   0.342   0.733
 Sdenom          0.054 0.123  42.604   0.442   0.661
 ses4            0.090 0.028 914.494   3.237   0.001
 PMses4         -0.036 0.230  48.395  -0.157   0.876
 SMses4          0.179 0.176  38.457   1.020   0.314
 PSMses4         0.023 0.034 915.094   0.680   0.497
 Pdenom:Sdenom   0.200 0.099 945.777   2.028   0.043

txt> PseudoR2(baseModel = Model1b, altModel1 = Model3a, altModel2 = Model3b)

Pseudo-R2 and Change in Pseudo-R2
 Term                   Model   Variance Pseudo.R2 Change.R2
 Intercept | PschoolID  Model1b    0.175                    
 Intercept | PschoolID  Model3a    0.173     0.012          
 Intercept | PschoolID  Model3b    0.177    -0.014    -0.026
 Intercept | SschoolID  Model1b    0.079                    
 Intercept | SschoolID  Model3a    0.068     0.138          
 Intercept | SschoolID  Model3b    0.066     0.155     0.017
 Residual               Model1b    0.505                    
 Residual               Model3a    0.483     0.043          
 Residual               Model3b    0.484     0.042    -0.001

txt> TotalR2(altModel1 = Model2, altModel2 = Model3a, altModel3 = Model3b)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model2     0.022          
 Model3a    0.053     0.031
 Model3b    0.058     0.005

txt> print("Model 4a: Add Random Within-Both-Centered Child SES Slope across Primary")
[1] "Model 4a: Add Random Within-Both-Centered Child SES Slope across Primary"
txt> Model4a = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + Pdenom + Sdenom + Pdenom:Sdenom + ses4 + PMses4 + SMses4 + 
txt+ 
txt+ PSMses4 + (1 + WPSses | PschoolID) + (1 | SschoolID))
txt> lmerSummary(Model4a, LRT = TRUE, randCI = TRUE, randCIfixed = "ses4", 
txt+ randCIrandom = "WPSses")

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2308.793  2372.594 -1141.396  2282.793 

Model for the Variance:
 Groups     Name      Variance Std.Dev  Corr
 PschoolID  Intercept    0.177   0.421      
            WPSses       0.000   0.010 1.000
 SschoolID  Intercept    0.066   0.258      
 Residual                0.484   0.695      
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Overall F-test for All Fixed Slopes:
 Model   NumDF  DenDF F.value p.value
 Model4a     7 95.281   7.492   <.001

Fixed Effects:
 Term              Est    SE      DF t.value p.value
 Intercept       6.166 0.138  77.500  44.610   <.001
 Pdenom          0.051 0.145  72.345   0.354   0.724
 Sdenom          0.054 0.123  42.626   0.439   0.663
 ses4            0.090 0.028 824.146   3.252   0.001
 PMses4         -0.036 0.230  48.395  -0.158   0.875
 SMses4          0.179 0.176  38.469   1.019   0.315
 PSMses4         0.023 0.034 877.592   0.665   0.506
 Pdenom:Sdenom   0.201 0.099 945.767   2.031   0.043

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC    LRT DF p.value p.mix
 <none>                   13 -1141 2283 2309 2373                        
 WPSses | PschoolID       11 -1141 2283 2305 2359  0.106  2   0.948 0.846
 Intercept | SschoolID    12 -1170 2341 2365 2423 57.724  1   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                             FixedEffect Variance    SD Lower Upper
 Intercept | PschoolID                  6.166    0.177 0.421 5.341 6.990
 Intercept | SschoolID                  6.166    0.066 0.258 5.660 6.671
 WPSses | PschoolID (fixed: ses4)       0.090    0.000 0.010 0.071 0.110

txt> print("Model 4b: Add Random Within-Both-Centered Child SES Slope across Secondary")
[1] "Model 4b: Add Random Within-Both-Centered Child SES Slope across Secondary"
txt> Model4b = lmer(data = Example6a, REML = TRUE, formula = achieve ~ 
txt+ 1 + Pdenom + Sdenom + Pdenom:Sdenom + ses4 + PMses4 + SMses4 + 
txt+ 
txt+ PSMses4 + (1 | PschoolID) + (1 + WPSses | SschoolID))
txt> lmerSummary(Model4b, LRT = TRUE, randCI = TRUE, randCIfixed = "ses4", 
txt+ randCIrandom = "WPSses")

Model Fit for Predicting achieve:
      AIC       BIC    logLik -2*logLik 
 2308.768  2372.568 -1141.384  2282.768 

Model for the Variance:
 Groups     Name      Variance Std.Dev  Corr
 PschoolID  Intercept    0.177   0.421      
 SschoolID  Intercept    0.067   0.258      
            WPSses       0.000   0.012 1.000
 Residual                0.484   0.695      
Number of obs: 1000 , groups: PschoolID, 50; SschoolID, 33 

Overall F-test for All Fixed Slopes:
 Model   NumDF  DenDF F.value p.value
 Model4b     7 94.986   7.451   <.001

Fixed Effects:
 Term              Est    SE      DF t.value p.value
 Intercept       6.164 0.138  77.594  44.612   <.001
 Pdenom          0.050 0.145  72.329   0.343   0.732
 Sdenom          0.059 0.123  42.763   0.481   0.633
 ses4            0.089 0.028 613.655   3.213   0.001
 PMses4         -0.035 0.230  48.387  -0.154   0.879
 SMses4          0.172 0.175  38.575   0.980   0.333
 PSMses4         0.024 0.034 762.649   0.689   0.491
 Pdenom:Sdenom   0.200 0.099 945.740   2.025   0.043

LRT for random effect(s):
 Random Effect Dropped  Npar    LL -2LL  AIC  BIC     LRT DF p.value p.mix
 <none>                   13 -1141 2283 2309 2373                         
 Intercept | PschoolID    12 -1238 2476 2500 2559 193.003  1   <.001 <.001
 WPSses | SschoolID       11 -1141 2283 2305 2359   0.132  2   0.936 0.827

Random Effect 95% CIs (Fixed Effect +/- 1.96*SQRT(random effect variance)):
 Term                             FixedEffect Variance    SD Lower Upper
 Intercept | PschoolID                  6.164    0.177 0.421 5.340 6.989
 Intercept | SschoolID                  6.164    0.067 0.258 5.659 6.670
 WPSses | SschoolID (fixed: ses4)       0.089    0.000 0.012 0.067 0.112

