Software Crosswalk: One Model in Every Program, and Troubleshooting

The chapters fit their models in R and close with a crosswalk box that gives the same models in other programs. This appendix goes one step further for two models: the random intercept and slope model with cross-level interactions for the High School and Beyond data (model M3 of Chapter 6) and the linear growth model for the Riesby data (model G1 of Chapter 19). It shows each model in every program, compares the estimates that the programs produce, explains why output that looks different is often the same result, and ends with a guide to the warnings and errors that multilevel software prints. The R fits, with lme4 (Bates et al., 2015), lmerTest (Kuznetsova et al., 2017), pbkrtest (Halekoh & Højsgaard, 2014), nlme (Pinheiro & Bates, 2000), glmmTMB (Brooks et al., 2017), and brms (Bürkner, 2017), and a Python fit with statsmodels (Seabold & Perktold, 2010) were run for this appendix. The syntax for the other programs was not run, because those programs were not available where the book was produced; it was written from the documentation of recent versions (SAS/STAT for SAS 9.4, IBM SPSS Statistics 29, Stata 18, Mplus 8, and HLM 8) and should be checked against the output of the version at hand. Book-length treatments of multilevel models in particular programs include Singer (1998) and Littell et al. (2006) for SAS, West et al. (2014) for several programs, and the Mplus User’s Guide (Muthén & Muthén, 1998--2017).

E.1The Programs and Their Defaults#

Most apparent disagreements between programs come from defaults rather than from different mathematics. Table E.1 lists, for each program, the default estimator, the default method for tests of fixed effects, the default structure of the covariance matrix of the random effects, and the sample size used in the BIC. A model fitted in two programs with the same estimator and the same covariance structure should reproduce the same estimates to several decimals; if it does not, check the data, the centering, and the coding of the predictors, and then whether each fit reached the maximum of the likelihood (Section E.6).

Three defaults cause most of the trouble. First, SAS and SPSS give random effects a diagonal (variance components) covariance matrix unless an unstructured matrix is requested, whereas the R packages estimate the covariance of a random intercept and slope by default, Stata assumes independence, and Mplus fixes their covariance at zero unless WITH frees it; a model with a random slope therefore differs across programs unless the structure is stated. Second, Stata, glmmTMB, Mplus, Julia, and MLwiN estimate by ML by default, whereas lme4, nlme, SAS, SPSS, HLM, JASP, jamovi, and statsmodels use REML; the variance components differ systematically between the two (Chapter 9). Third, the programs test fixed effects in different ways, from large-sample $z$ tests to $t$ tests with Satterthwaite or Kenward–Roger degrees of freedom (Chapter 10), so the same estimate can come with different $p$ values.

E.2One Model in Every Program: The HSB Cross-Level Model#

Model M3 of Chapter 6 predicts mathematics achievement from group-centered SES, the school mean of SES, the Catholic sector, and the two cross-level interactions, with a random intercept and a random slope of SES that are allowed to covary:

The data are the merged HSB file of Appendix C, which contains the group-centered SES (ses_cwc) and the sample school mean (meanses), so that every program uses the same variables (7,185 pupils in 160 schools). The syntax below requests REML and an unstructured covariance matrix wherever the program allows it.

R. The four packages use the same model formula, with minor differences in syntax:

library(lme4); library(lmerTest)          # lmerTest adds Satterthwaite df
m1 <- lmer(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
           data = hsb, REML = TRUE)
summary(m1)                                # summary(m1, ddf = "Kenward-Roger")
library(nlme)
m2 <- lme(mathach ~ (meanses + catholic) * ses_cwc, random = ~ ses_cwc | school,
          data = hsb, method = "REML")
library(glmmTMB)
m3 <- glmmTMB(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
              data = hsb, REML = TRUE)
library(brms)                              # Bayesian, default priors (Chapter 13)
m4 <- brm(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
          data = hsb, seed = 20270001)

SAS. PROC MIXED with TYPE=UN and Satterthwaite degrees of freedom:

proc mixed data=hsb method=reml covtest;
  class school;
  model mathach = meanses catholic ses_cwc meanses*ses_cwc catholic*ses_cwc
        / solution ddfm=satterthwaite;
  random intercept ses_cwc / subject=school type=un;
run;

SPSS. MIXED with COVTYPE(UN); the variables after WITH are covariates:

MIXED mathach WITH meanses catholic ses_cwc
  /FIXED=meanses catholic ses_cwc meanses*ses_cwc catholic*ses_cwc
  /METHOD=REML
  /PRINT=SOLUTION TESTCOV
  /RANDOM=INTERCEPT ses_cwc | SUBJECT(school) COVTYPE(UN).

Stata. mixed with factor-variable notation, REML, and Satterthwaite degrees of freedom:

mixed mathach c.meanses c.ses_cwc i.catholic c.meanses#c.ses_cwc ///
      i.catholic#c.ses_cwc || school: ses_cwc, ///
      covariance(unstructured) reml dfmethod(satterthwaite)

Mplus. Mplus has no REML; with ESTIMATOR = ML its estimates correspond to the ML fits of the R packages (REML = FALSE), and the default MLR changes only the standard errors, which then differ from the model-based standard errors of the other programs. The variable names are those of the Mplus file in the export archive of Appendix C:

DATA:     FILE = hsb.dat;
VARIABLE: NAMES = school minority female ses mathach size catholic pracad
                  disclim himinty meanses_ n_j meanses ses_cwc;
          USEVARIABLES = mathach ses_cwc meanses catholic;
          WITHIN = ses_cwc;  BETWEEN = meanses catholic;
          CLUSTER = school;  MISSING = ALL (-999);
ANALYSIS: TYPE = TWOLEVEL RANDOM;  ESTIMATOR = ML;
MODEL:    %WITHIN%
          s | mathach ON ses_cwc;
          %BETWEEN%
          mathach s ON meanses catholic;
          mathach WITH s;

HLM. HLM 8 is driven by menus. The multivariate data matrix (MDM) is built from a pupil file and a school file sorted by school, for example the SPSS files hsb_pupils.sav and hsb_schools.sav of the export archive. The model is specified equation by equation:

with SES entered group-mean centered and $u_{1j}$ activated. The school file of the archive contains the school mean SES as distributed with the data (meanses_rb), which differs slightly from the sample mean used here (Appendix C); to reproduce Table E.2 exactly, the sample mean has to be added to the school file.

JASP and jamovi. In JASP, the Linear Mixed Models analysis of the Mixed Models module takes mathach as the dependent variable, meanses, catholic, and ses_cwc as fixed-effects variables with the two interactions added in the model terms, and school as the random-effects grouping factor with a random slope of ses_cwc and the Correlations option checked. In jamovi, the Linear Mixed Model analysis of the GAMLj module takes the same variables, with the random effects entered as Intercept | school and ses_cwc | school and correlated. Both fit the model with lme4 (JASP through afex), but their defaults change the parameterization: JASP codes nominal predictors with sum-to-zero contrasts and includes random slopes for every fixed effect that varies within schools, and GAMLj centers covariates at their grand means. To reproduce the lme4 column of Table E.2, enter catholic as a 0/1 scale variable (or choose treatment coding), keep only the random intercept and the ses_cwc slope, and in GAMLj set the scaling of covariates to none; when only the coding or centering differs, the random part and the other coefficients agree, but the intercept, the SES slope, and (under sum coding) the Catholic terms do not.

Python and Julia. In Python, statsmodels fits the model with a formula and a grouping variable; in Julia, MixedModels uses a formula in lme4 style and estimates by ML unless REML is requested:

import statsmodels.formula.api as smf
m = smf.mixedlm("mathach ~ (meanses + catholic) * ses_cwc", hsb,
                groups=hsb["school"], re_formula="~ses_cwc")
fit = m.fit(reml=True)

using MixedModels, CSV, DataFrames
hsb = CSV.read("hsb.csv", DataFrame)
m = fit(MixedModel, @formula(mathach ~ (meanses + catholic) * ses_cwc +
        (1 + ses_cwc | school)), hsb; REML = true)

The estimates compared. Table E.2 shows the REML estimates of lme4, nlme, glmmTMB, and statsmodels, and the posterior means of brms. The four likelihood-based fits agree: the largest difference between them in any estimate is about 0.0004, and the REML log-likelihood is the same in all four ($-2\log L_R$ = 46,503.66). The ML fits of lme4, nlme, and statsmodels agree to the same degree (largest difference 0.0005). The standard errors of glmmTMB and statsmodels are slightly larger for some coefficients, most for the terms involving SES (for the interaction of mean SES and SES, 0.2996 against 0.2989), because these two programs derive them from the observed information matrix of all parameters jointly, which includes the uncertainty of the covariance parameters, whereas lme4 and nlme compute them from $({\mathbf{X}}^{\prime }{\hat{\mathbf{V}}}^{-1}\mathbf{X}{)}^{-1}$ with the covariance parameters treated as known; when that uncertainty is ignored, glmmTMB also gives 0.2989.

The Bayesian fit agrees closely for the fixed effects, whose posterior standard deviations are close to the REML standard errors, but not for the covariance parameters: the posterior of the slope variance is wide and strongly skewed (median 0.064, 95% credible interval [0.000, 0.475]), so its mean, 0.113, lies above the REML estimate, 0.101, and well above its own median, and the covariance is pulled toward zero (0.123 against 0.192), because the data carry little information about so small a variance (Chapter 13). Every $\hat{R}$ was below 1.01, and the smallest bulk or tail effective sample size over all parameters was 775.

The table also checks the syntax for the other programs: SAS, SPSS, and Stata with REML and an unstructured covariance matrix, and HLM with the sample school mean, should reproduce its first column to the precision they print, and Mplus with ESTIMATOR = ML the ML estimates (for example, a slope variance of 0.065 instead of 0.101).

E.3The Riesby Growth Model#

Model G1 of Chapter 19 is a linear growth model for the Hamilton depression ratings of the Riesby data, with a random intercept and a random slope of week that covary, estimated by ML from the 375 observed ratings of 66 patients. In long format, every program specifies it like the HSB model, with the occasion as level 1 and the patient as level 2:

lmer(hamd ~ week + (1 + week | id), data = riesby, REML = FALSE)

proc mixed data=riesby method=ml covtest;
  class id;
  model hamd = week / solution ddfm=satterthwaite;
  random intercept week / subject=id type=un;
run;

mixed hamd week || id: week, covariance(unstructured)

In SPSS, /METHOD=ML and /RANDOM=INTERCEPT week | SUBJECT(id) COVTYPE(UN); in Mplus, the same two-level setup as for HSB with week as the within-level variable, or a latent growth model in wide format (Chapter 22); in HLM, a two-level model with WEEK uncentered at level 1. Table E.3 shows that the four packages agree to within 0.002.

E.4Reconciling Output Across Programs#

When the estimates of two programs differ, five questions usually find the reason.

Same estimator? REML and ML estimates of the variance components differ, most for the smallest ones: the SES slope variance of the HSB model is 0.101 under REML and 0.065 under ML. The fixed effects differ much less. By default, Mplus and Stata estimate by ML, so their default output must be compared with ML fits; the MLR default of Mplus also changes the standard errors.

Same parameters? Programs report the random part in different forms: lme4, glmmTMB, and Julia print variances, standard deviations, and correlations; nlme prints standard deviations and correlations; SAS, SPSS, Mplus, statsmodels, and Stata (unless the option stddeviations is given) print variances and covariances, and HLM prints both forms. Squaring the standard deviations and multiplying the correlation by both standard deviations converts one form into the other. Programs also label the elements of the covariance matrix differently: the covariance ${\tau }_{01}$ appears as UN(2,1) in SAS, UN (2,1) in SPSS, cov(ses_cwc,_cons) in Stata, MATHACH WITH S at the between level in Mplus, and as the off-diagonal element of the $\tau$ matrix in HLM.

Same tests? The fixed effects are the same, but their tests differ. For the Catholic-by-SES interaction of the HSB model, the Satterthwaite approximation gives 143.2 degrees of freedom and the Kenward–Roger approximation 142.0, whereas nlme assigns the interaction the level-1 degrees of freedom (7,022) because it contains a level-1 variable, ignoring that the SES slope varies across schools, and HLM, which gives the coefficients of a random slope $J-q-1$ degrees of freedom, assigns it 157; for the sector effect itself, nlme uses 157 and Satterthwaite 149.6. With 160 schools the conclusions are the same, but with few clusters the choice matters (Chapter 10).

Same fit statistics? A REML log-likelihood cannot be compared with an ML log-likelihood, and REML likelihoods of models with different fixed effects cannot be compared at all (Chapters 9 to 11). lme4 reports the REML criterion ($-2\log L_R$ = 46,503.66 for the HSB model; extract it with REMLcrit(), not deviance(fit, REML = TRUE)) for REML fits and the deviance ($-2\log L$ = 46,496.43) for ML fits. The information criteria also differ in their penalty: for the ML fit of the HSB model, AIC is 46,516.43, BIC with the number of pupils is 46,585.23, and BIC with the number of schools, as SAS computes it, is 46,547.18; Stata gives the latter with estat ic, n(160), and in other programs it can be computed from the log-likelihood. Only differences between models fitted in the same program with the same definitions are meaningful.

Same data? Programs delete incomplete cases differently when variables are missing, and centering must be done identically; for the Riesby data, all programs must use the same 375 observed ratings. The HLM example above is a reminder that a variable with the same name, here the school mean of SES, can be computed in more than one way.

E.5Syntax Crosswalk for Common Model Features#

Table E.4 gives the syntax for the most common features of the models in this book in four programs. The chapters’ crosswalk boxes show complete examples; the table is meant as a quick reference.

E.6Troubleshooting#

Warnings and errors from multilevel software fall into a few groups, and each group has its usual causes and remedies (Table E.5). Table E.6 lists 15 messages that the examples of this appendix reproduce in R and Python, with the exact wording printed by the versions used in this book. Four of the examples use a small simulated data set of 10 clusters of five observations without any cluster variance, which is the typical situation in which a random slope cannot be estimated: lme4 reports a boundary (singular) fit, glmmTMB a non-positive-definite Hessian and singular convergence, and nlme stops with an error; the fourth example, with two predictors that are multiples of each other, adds a message about rank deficiency to the singular fit.

Two further cases from the HSB model show that a warning is not the only sign of trouble, and that its absence is no guarantee. Estimated by ML, the model gave glmmTMB a non-positive-definite Hessian and a slope variance of 0.018 (logLik() returns NA for this fit; its objective gives a log-likelihood of $-$23,248.24, 0.03 below the maximum), whereas lme4, nlme, and statsmodels agreed on 0.065; the likelihood is flat in this small variance, and the warning correctly flags that the optimizer stopped short of the maximum. Estimated by REML with the alternative optimizer of nlme (lmeControl(opt = "optim")), the model stopped at a slope variance of 0.0003 with a REML log-likelihood 0.56 below the maximum, without any message. Comparing the log-likelihood across optimizers or programs, which takes one line of code, finds both problems, provided that the fits maximize the same likelihood (the same estimator and, for generalized models, the same approximation). statsmodels with the L-BFGS optimizer alone also stopped early, but it said so (Table E.6); its default sequence of optimizers reached the maximum.