# Why does this mixed model produce discrepant output in SPSS and R?

I've been working with an example from the Heck et al (2013) textbook, which gives the following model:

$$Y_{ij} = \gamma_{00} + \gamma_{01}\text{ses_mean}_j + \gamma_{02}\text{pro4yrc}_j + \gamma_{03}\text{public}_j + \gamma_{10}\text{ses}_{ij} + u_{0j} + u_{1j}\text{ses}_{ij} + \epsilon_{ij}$$

I tried to run this in R using the following code

heckc3  <- read_sav(file ="https://github.com/user1205901/codeforposting/blob/main/ch3multilevel.sav?raw=true")

model1 <- lmer("math ~ 1 + ses + ses_mean + pro4yrc + public + (1 + ses|schcode)", data = heckc3)

summary(model1)


I found that the results did not match the SPSS output, though the output I got from SPSS did match the output in the book.

In my previous experience running uncomplicated models of this sort in SPSS and lme4 I had thought that parameter estimates tended to come out very similar or identical, but here there is a substantial difference.

Trying to investigate why this happened, I noticed something defective in the dataset such that within the Level 2 units 418 and 419 certain individuals have different values on the "public" (Level 2) predictor. However, when I run heckc3 <- subset(heckc3, schcode < 417) rerun the R analysis it still does not match and SPSS analysis with those Level 2 units deleted.

Heck, R. H., Thomas, S. L., & Tabata, L. N. (2012). Multilevel and longitudinal modeling with IBM SPSS. Routledge.

Okay, so I think that I have an answer here that is explaining the differences. The model in lme4 results in a singular fit. Looking at the model summary, this is happening because the correlation parameter between the random effects could not be estimated. When this parameter cannot be estimated, it will normally show up in the summary as 1.00 or -1.00. Here is -1.00 as shown in the printout of your model below.

As this parameter could not be estimated properly, I dropped the correlation parameter and fit the same model with an uncorrelated random slope below. The estimates from this model are very close to the estimates you posted from SPSS, and are definitely within the range of differences that I would expect from slightly different algorithms for maximum likelihood estimation.

What I believe is happening is that SPSS is automatically dropping this correlation parameter, or that SPSS defaults to fitting uncorrelated random slopes. This is somewhat supported by that not being reported in the table you posted, although I have not been able to confirm this in the SPSS documentation.

Here is the printout of the model fitted by your code on my machine.

summary(model1)
Linear mixed model fit by REML ['lmerMod']
Formula: math ~ 1 + ses + ses_mean + pro4yrc + public + (1 + ses | schcode)
Data: d

REML criterion at convergence: 48095.1

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.8525 -0.5590  0.1383  0.6541  5.7017

Random effects:
Groups   Name        Variance Std.Dev. Corr
schcode  (Intercept)  2.045   1.430
ses          1.317   1.147    -1.00
Residual             62.204   7.887
Number of obs: 6871, groups:  schcode, 419

Fixed effects:
Estimate Std. Error t value
(Intercept)  56.6683     0.4662 121.561
ses           3.1874     0.1684  18.926
ses_mean      2.4601     0.2796   8.798
pro4yrc       1.3146     0.4700   2.797
public       -0.3210     0.2609  -1.230

Correlation of Fixed Effects:
(Intr) ses    ses_mn pr4yrc
ses      -0.044
ses_mean  0.209 -0.521
pro4yrc  -0.875 -0.011 -0.246
public   -0.412 -0.002 -0.053  0.005
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')


Here is the code that fit the model with uncorrelated random slopes, and prints out the model summary. The numbers are almost exactly the same as SPSS:

model2 <- lmer(math ~ 1 + ses + ses_mean + pro4yrc + public + (ses||schcode),
data = d)

summary(model2)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ 1 + ses + ses_mean + pro4yrc + public + ((1 | schcode) +      (0 + ses | schcode))
Data: d

REML criterion at convergence: 48121.8

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.7910 -0.5538  0.1342  0.6587  5.7021

Random effects:
Groups    Name        Variance Std.Dev.
schcode   (Intercept)  2.112   1.453
schcode.1 ses          1.314   1.146
Residual              62.115   7.881
Number of obs: 6871, groups:  schcode, 419

Fixed effects:
Estimate Std. Error t value
(Intercept)  56.4698     0.4716 119.749
ses           3.1639     0.1689  18.734
ses_mean      2.6596     0.3136   8.480
pro4yrc       1.3602     0.4679   2.907
public       -0.1200     0.2744  -0.437

Correlation of Fixed Effects:
(Intr) ses    ses_mn pr4yrc
ses      -0.001
ses_mean  0.233 -0.486
pro4yrc  -0.868 -0.001 -0.255
public   -0.444 -0.002 -0.045  0.023