1
$\begingroup$

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.

SPSS Estimates of Fixed Effects

SPSS Estimates of Covariance Parameters

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.

$\endgroup$
0

1 Answer 1

4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.