> print("Descriptive Statistics in School-Level-2 Data")
[1] "Descriptive Statistics in School-Level-2 Data"
> univSummary(data = Example3Means, vars = c("CMmath", "CMses", 
+ "CMnon"))

 Variable   N   Mean    SD Variance    Min    Max
 CMmath   160 12.621 3.118    9.720  4.240 19.719
 CMses    160 -0.006 0.414    0.171 -1.194  0.825
 CMnon    160  0.275 0.301    0.091  0.000  1.000

> print("Correlations in School-Level-2 Data")
[1] "Correlations in School-Level-2 Data"
> corrSummary(data = Example3Means, p.separate = FALSE, vars = c("CMmath", 
+ "CMses", "CMnon"))

Correlations (and p-values):
 Variable  CMmath   CMses   CMnon
 CMmath     1.000   0.785  -0.491
                  (<.001) (<.001)
 CMses      0.785   1.000  -0.483
          (<.001)         (<.001)
 CMnon     -0.491  -0.483   1.000
          (<.001) (<.001)        

> Example3$CMnon30 = Example3$CMnon - 0.3
> Example3$WCses = Example3$ses - Example3$CMses
> Example3$WCnon = Example3$nonwhite - Example3$CMnon
> Example3$WCmath = Example3$math - Example3$CMmath
> print("R Descriptive Statistics in Student-Level-1 Data")
[1] "R Descriptive Statistics in Student-Level-1 Data"
> univSummary(data = Example3, vars = c("math", "WCmath", "CMmath", 
+ "ses", "WCses", "CMses", "nonwhite", "WCnon", "CMnon"))

 Variable    N   Mean    SD Variance     Min    Max
 math     7185 12.748 6.878   47.310  -2.832 24.993
 WCmath   7185  0.000 6.187   38.275 -19.675 17.969
 CMmath   7185 12.748 3.006    9.035   4.240 19.719
 ses      7185  0.000 0.779    0.607  -3.758  2.692
 WCses    7185  0.000 0.661    0.436  -3.651  2.856
 CMses    7185  0.000 0.414    0.171  -1.194  0.825
 nonwhite 7185  0.275 0.446    0.199   0.000  1.000
 WCnon    7185  0.000 0.327    0.107  -0.982  0.982
 CMnon    7185  0.275 0.304    0.093   0.000  1.000

> print("Correlations in Student-Level-1 Data")
[1] "Correlations in Student-Level-1 Data"
> corrSummary(data = Example3, p.separate = TRUE, vars = c("math", 
+ "ses", "nonwhite", "WCmath", "WCses", "WCnon", "CMmath", 
+ "CMses", "CMnon"))

Correlations:
 Variable   math    ses nonwhite WCmath  WCses  WCnon CMmath  CMses  CMnon
 math      1.000  0.361   -0.268  0.899  0.210 -0.169  0.437  0.344 -0.212
 ses       0.361  1.000   -0.271  0.198  0.848 -0.141  0.417  0.531 -0.247
 nonwhite -0.268 -0.271    1.000 -0.137 -0.122  0.732 -0.331 -0.317  0.681
 WCmath    0.899  0.198   -0.137  1.000  0.234 -0.188 -0.000  0.000  0.000
 WCses     0.210  0.848   -0.122  0.234  1.000 -0.167 -0.000 -0.000  0.000
 WCnon    -0.169 -0.141    0.732 -0.188 -0.167  1.000 -0.000 -0.000  0.000
 CMmath    0.437  0.417   -0.331 -0.000 -0.000 -0.000  1.000  0.787 -0.485
 CMses     0.344  0.531   -0.317  0.000 -0.000 -0.000  0.787  1.000 -0.465
 CMnon    -0.212 -0.247    0.681  0.000  0.000  0.000 -0.485 -0.465  1.000

P-values:
 Variable  math   ses nonwhite WCmath WCses WCnon CMmath CMses CMnon
 math           <.001    <.001  <.001 <.001 <.001  <.001 <.001 <.001
 ses      <.001          <.001  <.001 <.001 <.001  <.001 <.001 <.001
 nonwhite <.001 <.001           <.001 <.001 <.001  <.001 <.001 <.001
 WCmath   <.001 <.001    <.001        <.001 <.001  1.000 1.000 1.000
 WCses    <.001 <.001    <.001  <.001       <.001  1.000 1.000 1.000
 WCnon    <.001 <.001    <.001  <.001 <.001        1.000 1.000 1.000
 CMmath   <.001 <.001    <.001  1.000 1.000 1.000        <.001 <.001
 CMses    <.001 <.001    <.001  1.000 1.000 1.000  <.001       <.001
 CMnon    <.001 <.001    <.001  1.000 1.000 1.000  <.001 <.001      

> print("Model 1: Empty Means, Random Intercept for Math")
[1] "Model 1: Empty Means, Random Intercept for Math"
> Model1 = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + (1 | schoolID))
> lmerSummary(Model1, CI = TRUE, ICC1 = TRUE, LRT = TRUE, randCI = 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 

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

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("Empty Means, Random Intercept for SES")
[1] "Empty Means, Random Intercept for SES"
> EmptySES = lmer(data = Example3, REML = TRUE, formula = ses ~ 
+ 1 + (1 | schoolID))
> lmerSummary(EmptySES, CI = TRUE, ICC1 = TRUE, LRT = TRUE, randCI = TRUE)

Model Fit for Predicting ses:
      AIC       BIC    logLik -2*logLik 
15052.016 15072.655 -7523.008 15046.016 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    0.161   0.401
 Residual               0.446   0.668
Number of obs: 7185 , groups: schoolID, 160 

Fixed Effects (with 95% CI):
 Term          Est    SE      DF t.value p.value  Lower Upper
 Intercept  -0.006 0.033 159.187  -0.173   0.863 -0.070 0.059

ICC1 intraclass correlation:
    Adjusted ICC: 0.265 
  Unadjusted ICC: 0.265 

LRT for random effect(s):
               Npar    LL  -2LL   AIC   BIC      LRT DF p.value p.mix
<none>            3 -7523 15046 15052 15073                          
(1 | schoolID)    2 -8407 16814 16818 16832 1768.398  1   <.001 <.001

Random Effect 95% CIs (Fixed Effect +/- 1.96*SD(random effect variance)):
 Term      FixedEffect Variance    SD  Lower Upper
 Intercept      -0.006    0.161 0.401 -0.792 0.780

> print("Empty Means, Random Intercept for Binary Nonwhite")
[1] "Empty Means, Random Intercept for Binary Nonwhite"
> print("Logit link, bernoulli level-1 distribution, 15 quadrature points")
[1] "Logit link, bernoulli level-1 distribution, 15 quadrature points"
> EmptyNon = glmer(data = Example3, family = binomial(link = "logit"), 
+ nAGQ = 15, nonwhite ~ 1 + (1 | schoolID))
> glmerSummary(EmptyNon, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting nonwhite:
      AIC       BIC    logLik -2*logLik 
 5486.867  5500.626 -2741.433  5482.867 

Model for the Variance:
 Groups    Name             Variance Std.Dev
 schoolID  Intercept           5.467   2.338
 Residual  (logit constant)  (3.290) (1.814)
Number of obs: 7185 , groups: schoolID, 160 

Fixed Effects:
 Term          Est    SE z.value p.value
 Intercept  -1.620 0.195  -8.306   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.624 
  Unadjusted ICC: 0.624 

LRTs for Random Effects (each removed one at a time):
               Model Npar    LL -2LL  AIC  BIC      LRT DF p.value p.mix
              <none>    2 -2741 5483 5487 5501                          
 1 in (1 | schoolID)    1 -4224 8448 8450 8457 2965.481  1   <.001 <.001

> print("Empty Means, Random Intercept for Binary Nonwhite")
[1] "Empty Means, Random Intercept for Binary Nonwhite"
> print("Logit link, bernoulli level-1 distribution, 1 quadrature point as default")
[1] "Logit link, bernoulli level-1 distribution, 1 quadrature point as default"
> EmptyNon = glmer(data = Example3, family = binomial(link = "logit"), 
+ nAGQ = 1, nonwhite ~ 1 + (1 | schoolID))
> glmerSummary(EmptyNon, ICC1 = TRUE, LRT = TRUE)

Model Fit for Predicting nonwhite:
      AIC       BIC    logLik -2*logLik 
 5492.077  5505.837 -2744.039  5488.077 

Model for the Variance:
 Groups    Name             Variance Std.Dev
 schoolID  Intercept           5.322   2.307
 Residual  (logit constant)  (3.290) (1.814)
Number of obs: 7185 , groups: schoolID, 160 

Fixed Effects:
 Term          Est    SE z.value p.value
 Intercept  -1.617 0.190  -8.496   <.001

ICC1 intraclass correlation:
    Adjusted ICC: 0.618 
  Unadjusted ICC: 0.618 

LRTs for Random Effects (each removed one at a time):
               Model Npar    LL -2LL  AIC  BIC      LRT DF p.value
              <none>    2 -2744 5488 5492 5506                    
 1 in (1 | schoolID)    1 -4224 8448 8450 8457 2960.271  1   <.001
 p.mix
      
 <.001

> print("Model 2a: Add L1 Cluster-MC Student SES")
[1] "Model 2a: Add L1 Cluster-MC Student SES"
> Model2a = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + WCses + (1 | schoolID))
> lmerSummary(Model2a)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46731.996  46759.515 -23361.998  46723.996 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    8.672   2.945
 Residual              37.010   6.084
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF    DenDF F.value p.value
 Model2a     1 7022.024 406.684   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  12.636 0.244  156.741  51.684   <.001
 WCses       2.191 0.109 7022.024  20.166   <.001

> PseudoR2(baseModel = Model1, altModel1 = Model2a)

Pseudo-R2
 Term       Model   Variance Pseudo.R2
 Intercept  Model1     8.614          
 Intercept  Model2a    8.672    -0.007
 Residual   Model1    39.148          
 Residual   Model2a   37.010     0.055

> TotalR2(altModel1 = Model2a)

Total R2
 Model   Total.R2
 Model2a    0.044

> print("Model 2b: Add L2 Cluster Mean of Student SES")
[1] "Model 2b: Add L2 Cluster Mean of Student SES"
> Model2b = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + WCses + CMses + (1 | schoolID))
> lmerSummary(Model2b)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46578.584  46612.983 -23284.292  46568.584 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    2.693   1.641
 Residual              37.019   6.084
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model2b     2 298.343 334.812   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  12.683 0.149  153.652  84.906   <.001
 WCses       2.191 0.109 7021.509  20.164   <.001
 CMses       5.866 0.362  153.367  16.218   <.001

> con1DSummary(contest1D(Model2b, L = c(0, -1, 1)), name = "L2 Contextual SES")

Fixed Effect: L2 Contextual SES
   Est    SE      DF t.value p.value
 3.675 0.378 182.286   9.731   <.001

> PseudoR2(baseModel = Model1, altModel1 = Model2a, altModel2 = Model2b)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1     8.614                    
 Intercept  Model2a    8.672    -0.007          
 Intercept  Model2b    2.693     0.687     0.694
 Residual   Model1    39.148                    
 Residual   Model2a   37.010     0.055          
 Residual   Model2b   37.019     0.054    -0.000

> TotalR2(altModel1 = Model2a, altModel2 = Model2b)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model2a    0.044          
 Model2b    0.162     0.118

> print("Model 3a: Use Smushed Main Effect of L1 Student SES Instead")
[1] "Model 3a: Use Smushed Main Effect of L1 Student SES Instead"
> Model3a = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + ses + (1 | schoolID))
> lmerSummary(Model3a)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46653.169  46680.688 -23322.585  46645.169 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    4.768   2.184
 Residual              37.034   6.086
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF    DenDF F.value p.value
 Model3a     1 6838.078 511.164   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  12.657 0.188  148.302  67.332   <.001
 ses         2.390 0.106 6838.078  22.609   <.001

> PseudoR2(baseModel = Model1, altModel1 = Model3a)

Pseudo-R2
 Term       Model   Variance Pseudo.R2
 Intercept  Model1     8.614          
 Intercept  Model3a    4.768     0.446
 Residual   Model1    39.148          
 Residual   Model3a   37.034     0.054

> TotalR2(altModel1 = Model3a)

Total R2
 Model   Total.R2
 Model3a    0.130

> print("Model 3b: L1 Grand-MC Student SES + L2 Cluster Mean SES")
[1] "Model 3b: L1 Grand-MC Student SES + L2 Cluster Mean SES"
> Model3b = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + ses + CMses + (1 | schoolID))
> lmerSummary(Model3b)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46578.584  46612.983 -23284.292  46568.584 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    2.693   1.641
 Residual              37.019   6.084
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model3b     2 362.097 334.812   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  12.683 0.149  153.652  84.906   <.001
 ses         2.191 0.109 7021.509  20.164   <.001
 CMses       3.675 0.378  182.286   9.731   <.001

> con1DSummary(contest1D(Model3b, L = c(0, 1, 1)), name = "L2 Between SES")

Fixed Effect: L2 Between SES
   Est    SE      DF t.value p.value
 5.866 0.362 153.367  16.218   <.001

> PseudoR2(baseModel = Model1, altModel1 = Model3a, altModel2 = Model3b)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1     8.614                    
 Intercept  Model3a    4.768     0.446          
 Intercept  Model3b    2.693     0.687     0.241
 Residual   Model1    39.148                    
 Residual   Model3a   37.034     0.054          
 Residual   Model3b   37.019     0.054     0.000

> TotalR2(altModel1 = Model3a, altModel2 = Model3b)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model3a    0.130          
 Model3b    0.162     0.032

> print("Model 4a: Add L1 Binary Student Nonwhite")
[1] "Model 4a: Add L1 Binary Student Nonwhite"
> Model4a = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + ses + CMses + nonwhite + (1 | schoolID))
> lmerSummary(Model4a)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46408.480  46449.759 -23198.240  46396.480 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    2.638   1.624
 Residual              36.139   6.012
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model4a     3 542.008 286.922   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  13.424 0.158  191.567  84.955   <.001
 ses         1.967 0.109 7063.565  18.102   <.001
 CMses       2.953 0.378  186.919   7.822   <.001
 nonwhite   -2.714 0.205 4400.818 -13.249   <.001

> PseudoR2(baseModel = Model1, altModel1 = Model3b, altModel2 = Model4a)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1     8.614                    
 Intercept  Model3b    2.693     0.687          
 Intercept  Model4a    2.638     0.694     0.006
 Residual   Model1    39.148                    
 Residual   Model3b   37.019     0.054          
 Residual   Model4a   36.139     0.077     0.022

> print("Model 4b: L1 Binary Student Nonwhite + L2 Cluster Mean Nonwhite")
[1] "Model 4b: L1 Binary Student Nonwhite + L2 Cluster Mean Nonwhite"
> Model4b = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + ses + CMses + nonwhite + CMnon30 + (1 | schoolID))
> lmerSummary(Model4b)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46404.554  46452.712 -23195.277  46390.554 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    2.559   1.600
 Residual              36.136   6.011
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model4b     4 399.685 218.054   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  13.511 0.161  219.691  83.906   <.001
 ses         1.952 0.109 7020.170  17.931   <.001
 CMses       3.375 0.417  172.761   8.090   <.001
 nonwhite   -2.896 0.220 7020.170 -13.152   <.001
 CMnon30     1.351 0.594  199.136   2.274   0.024

> print("Computing between = within + contextual")
[1] "Computing between = within + contextual"
> con1DTable(Model4b, L = list(`SES L2 Between` = c(0, 1, 1, 0, 
+ 0), `Nonwhite L2 Between` = c(0, 0, 0, 1, 1)))

Requested Fixed Effects:
 Name                   Est    SE      DF t.value p.value
 SES L2 Between       5.327 0.403 150.020  13.230   <.001
 Nonwhite L2 Between -1.544 0.552 148.270  -2.798   0.006

> PseudoR2(baseModel = Model1, altModel1 = Model4a, altModel2 = Model4b)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1     8.614                    
 Intercept  Model4a    2.638     0.694          
 Intercept  Model4b    2.559     0.703     0.009
 Residual   Model1    39.148                    
 Residual   Model4a   36.139     0.077          
 Residual   Model4b   36.136     0.077     0.000

> print("Model 4c: Use Cluster-MC Versions of L1 Predictors Instead")
[1] "Model 4c: Use Cluster-MC Versions of L1 Predictors Instead"
> Model4c = lmer(data = Example3, REML = TRUE, formula = math ~ 
+ 1 + WCses + CMses + WCnon + CMnon30 + (1 | schoolID))
> lmerSummary(Model4c)

Model Fit for Predicting math:
       AIC        BIC     logLik  -2*logLik 
 46404.554  46452.712 -23195.277  46390.554 

Model for the Variance:
 Groups    Name      Variance Std.Dev
 schoolID  Intercept    2.559   1.600
 Residual              36.136   6.011
Number of obs: 7185 , groups: schoolID, 160 

Overall F-test for All Fixed Slopes:
 Model   NumDF   DenDF F.value p.value
 Model4c     4 292.047 218.054   <.001

Fixed Effects:
 Term          Est    SE       DF t.value p.value
 Intercept  12.642 0.147  152.081  86.087   <.001
 WCses       1.952 0.109 7020.170  17.931   <.001
 CMses       5.327 0.403  150.020  13.230   <.001
 WCnon      -2.896 0.220 7020.170 -13.152   <.001
 CMnon30    -1.544 0.552  148.270  -2.798   0.006

> print("Computing contextual = between - within")
[1] "Computing contextual = between - within"
> con1DTable(Model4c, L = list(`SES L2 Contextual` = c(0, -1, 1, 
+ 0, 0), `Nonwhite L2 Contextual` = c(0, 0, 0, -1, 1), `Nonwhite L2 Between per 10%` = c(0, 
+ 0, 0, 0, 1/10), `Nonwhite L2 Contextual per 10%` = c(0, 0, 
+ 0, -1/10, 1/10)))

Requested Fixed Effects:
 Name                              Est    SE      DF t.value p.value
 SES L2 Contextual               3.375 0.417 172.761   8.090   <.001
 Nonwhite L2 Contextual          1.351 0.594 199.136   2.274   0.024
 Nonwhite L2 Between per 10%    -0.154 0.055 148.270  -2.798   0.006
 Nonwhite L2 Contextual per 10%  0.135 0.059 199.136   2.274   0.024

> PseudoR2(baseModel = Model1, altModel1 = Model2b, altModel2 = Model3b, 
+ altModel3 = Model4c)

Pseudo-R2 and Change in Pseudo-R2
 Term       Model   Variance Pseudo.R2 Change.R2
 Intercept  Model1     8.614                    
 Intercept  Model2b    2.693     0.687          
 Intercept  Model3b    2.693     0.687    -0.000
 Intercept  Model4c    2.559     0.703     0.016
 Residual   Model1    39.148                    
 Residual   Model2b   37.019     0.054          
 Residual   Model3b   37.019     0.054     0.000
 Residual   Model4c   36.136     0.077     0.023

> TotalR2(altModel1 = Model2b, altModel2 = Model3b, altModel3 = Model4c)

Total-R2 and Change in Total-R2
 Model   Total.R2 Change.R2
 Model2b    0.162          
 Model3b    0.162    -0.000
 Model4c    0.184     0.022

> packageVersion("lme4")
[1] ‘2.0.6’
> citation("lme4")
To cite lme4 in publications use:

  Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker
  (2015). Fitting Linear Mixed-Effects Models Using lme4.
  Journal of Statistical Software, 67(1), 1-48.
  doi:10.18637/jss.v067.i01.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {Fitting Linear Mixed-Effects Models Using {lme4}},
    author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
    journal = {Journal of Statistical Software},
    year = {2015},
    volume = {67},
    number = {1},
    pages = {1--48},
    doi = {10.18637/jss.v067.i01},
  }
> packageVersion("lmerTest")
[1] ‘3.2.1’
> citation("lmerTest")
To cite lmerTest in publications use:

  Kuznetsova A, Brockhoff PB, Christensen RHB (2017). “lmerTest
  Package: Tests in Linear Mixed Effects Models.” _Journal of
  Statistical Software_, *82*(13), 1-26.
  doi:10.18637/jss.v082.i13
  <https://doi.org/10.18637/jss.v082.i13>.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
    author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
    journal = {Journal of Statistical Software},
    year = {2017},
    volume = {82},
    number = {13},
    pages = {1--26},
    doi = {10.18637/jss.v082.i13},
  }
