replicating Stata's xtmixed with lmer in R

Hi all. I am trying to replicate a Stata paper in R and am having an issue where the results of a multi-level mixed-effects linear regression are very similar, but slightly off (especially the standard errors). I am wondering if anyone knows why this is happening. I am using the exact same dataset in both specifications (which can be found here). The Stata code I am trying to replicate is:

use "/Users/mariaburzillo/Desktop/GOV1006/final_project/racial_polarization_winners.dta", clear

xtmixed biggestsplit H_citytract_multi_i diversityinterp pctasianpopinterp pctblkpopinterp pctlatinopopinterp medincinterp pctrentersinterp pctcollegegradinterp biracial nonpartisan primary logpop  i. year south midwest west if winner==1||geo_id2:

Which gives the following output:

The reprex for my attempt to replicate this in R is:

library(haven)
library(tidyverse)
library(lme4)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
library(reprex)

# load data

rp <- read_dta("/Users/mariaburzillo/Desktop/GOV1006/final_project/racial_polarization_winners.dta")


# create factor variable for year

rp$year.f <- as.factor(rp$year)

# apply condition that winner == 1 as in Stata regression

rp_sub <- rp %>%
  filter(winner == 1)

# multi-level mixed effects regression with random effects for geo_id2

m1 <- lmer(biggestsplit ~ H_citytract_multi_i + diversityinterp + pctasianpopinterp + pctblkpopinterp + pctlatinopopinterp + medincinterp + pctrentersinterp +  pctcollegegradinterp + biracial + nonpartisan + primary + logpop + year.f + south + midwest + west + (1 | geo_id2), data = rp_sub)
#> Warning: Some predictor variables are on very different scales: consider
#> rescaling
#> boundary (singular) fit: see ?isSingular

# summarize

summary(m1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: 
#> biggestsplit ~ H_citytract_multi_i + diversityinterp + pctasianpopinterp +  
#>     pctblkpopinterp + pctlatinopopinterp + medincinterp + pctrentersinterp +  
#>     pctcollegegradinterp + biracial + nonpartisan + primary +  
#>     logpop + year.f + south + midwest + west + (1 | geo_id2)
#>    Data: rp_sub
#> 
#> REML criterion at convergence: -5.2
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -2.0918 -0.4281  0.0000  0.5283  1.9850 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  geo_id2  (Intercept) 0.00000  0.000   
#>  Residual             0.02312  0.152   
#> Number of obs: 91, groups:  geo_id2, 25
#> 
#> Fixed effects:
#>                        Estimate Std. Error t value
#> (Intercept)          -2.419e-01  7.126e-01  -0.339
#> H_citytract_multi_i   9.320e-01  4.937e-01   1.888
#> diversityinterp       3.845e-01  4.530e-01   0.849
#> pctasianpopinterp    -1.149e-01  6.605e-01  -0.174
#> pctblkpopinterp      -4.318e-01  3.373e-01  -1.280
#> pctlatinopopinterp   -1.905e-01  3.213e-01  -0.593
#> medincinterp         -4.225e-06  8.304e-06  -0.509
#> pctrentersinterp     -5.796e-01  5.282e-01  -1.097
#> pctcollegegradinterp  3.277e-01  8.903e-01   0.368
#> biracial              2.104e-01  4.647e-02   4.527
#> nonpartisan          -8.956e-02  8.312e-02  -1.077
#> primary              -9.161e-02  3.998e-02  -2.291
#> logpop                3.466e-02  6.904e-02   0.502
#> year.f1991            7.068e-02  7.788e-02   0.908
#> year.f1993            6.612e-02  1.159e-01   0.570
#> year.f1994            1.153e-01  2.075e-01   0.556
#> year.f1995            4.397e-02  9.424e-02   0.467
#> year.f1996            1.255e-01  1.846e-01   0.680
#> year.f1997            5.762e-02  1.097e-01   0.526
#> year.f1998            2.933e-02  1.963e-01   0.149
#> year.f1999           -2.390e-02  1.027e-01  -0.233
#> year.f2000            2.958e-01  1.871e-01   1.581
#> year.f2001            8.088e-02  1.192e-01   0.679
#> year.f2002            6.029e-02  1.964e-01   0.307
#> year.f2003           -3.516e-02  1.233e-01  -0.285
#> year.f2004           -6.759e-02  1.459e-01  -0.463
#> year.f2005            1.175e-02  1.399e-01   0.084
#> year.f2006            5.294e-02  2.908e-01   0.182
#> year.f2007           -9.040e-02  1.831e-01  -0.494
#> year.f2009           -2.597e-02  2.317e-01  -0.112
#> south                 2.331e-01  1.194e-01   1.953
#> midwest               1.534e-01  1.328e-01   1.156
#> west                  8.498e-02  1.148e-01   0.740
#> 
#> Correlation matrix not shown by default, as p = 33 > 12.
#> Use print(x, correlation=TRUE)  or
#>     vcov(x)        if you need it
#> fit warnings:
#> Some predictor variables are on very different scales: consider rescaling
#> convergence code: 0
#> boundary (singular) fit: see ?isSingular

In response to a request for an example from a built-in dataset, etc., here is an alternative reprex using the EmplUK dataset from the plm package:

library(tidyverse)
library(lme4)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
library(reprex)
library(plm)
#> 
#> Attaching package: 'plm'
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, lag, lead

# load data from plm package

data("EmplUK", package="plm")

# multi-level mixed effects regression with random effects for firm


m4 <- lmer(wage ~ as.factor(year) + sector + emp + capital + 
             output + (1 | firm), data = EmplUK)

# summarize

summary(m4)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: wage ~ as.factor(year) + sector + emp + capital + output + (1 |  
#>     firm)
#>    Data: EmplUK
#> 
#> REML criterion at convergence: 4841.1
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -6.5171 -0.5067 -0.0280  0.4502  7.9279 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  firm     (Intercept) 21.489   4.636   
#>  Residual              3.793   1.948   
#> Number of obs: 1031, groups:  firm, 140
#> 
#> Fixed effects:
#>                     Estimate Std. Error t value
#> (Intercept)         26.57338    1.55552  17.083
#> as.factor(year)1977 -1.99478    0.27950  -7.137
#> as.factor(year)1978 -2.63032    0.27979  -9.401
#> as.factor(year)1979 -2.51523    0.27983  -8.988
#> as.factor(year)1980 -2.33948    0.28846  -8.110
#> as.factor(year)1981 -1.37792    0.32335  -4.261
#> as.factor(year)1982 -0.33180    0.33480  -0.991
#> as.factor(year)1983  0.49106    0.37275   1.317
#> as.factor(year)1984  0.33148    0.43801   0.757
#> sector              -0.78632    0.14900  -5.277
#> emp                 -0.09098    0.02393  -3.801
#> capital              0.17433    0.05542   3.145
#> output               0.03097    0.01132   2.736
#> 
#> Correlation matrix not shown by default, as p = 13 > 12.
#> Use print(x, correlation=TRUE)  or
#>     vcov(x)        if you need it


# write to csv for use in stata

write.csv(EmplUK, file = "empUK.csv" )

The corresponding stata command would be:

xtmixed wage i.year sector emp capital output || firm:

With this example, the outputs are again similar, but they are slightly off. I am trying to see if I can exactly reproduce the Stata code in R.

Thank you for the help! Happy to answer any more questions.

Here's a stern reply with kind-hearted intentions:

Screenshots of a reprex are not the thing itself. The whole point of a reprex is that it can be generated, then cut-and-pasted into the dialog box for the question, cut-and-pasted out and then pasted into an R session to see what's happening.

I'm not sure if I can help with this or not, but I'm pretty sure that no one is going to if they have to print out the image and type in the code.

Thanks! Can you share enough of rp to run the code or a built-in data set that will throw the same error?

Update: I tried my original code with

and this solved the problem with my dataset! I am not sure why it didn't solve the problem in the other example I posted, but it may have been another issue within that example. All in all, thank you!!!

1 Like

Great. Please mark the solution for the benefit of those to follow.

1 Like

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

Hi. Thank you for the feedback, I don't know why I did that! I just edited the post to include the actual reprex. Thank you.

1 Like

I used the tab delimited file and received the following

suppressPackageStartupMessages(library(lme4)) 
suppressPackageStartupMessages(library(plm)) 

# load data from plm package

data("EmplUK", package="plm")

# multi-level mixed effects regression with random effects for firm

m4 <- lmer(wage ~ as.factor(year) + sector + emp + capital + 
             output + (1 | firm), data = EmplUK)

# summarize

summary(m4)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: wage ~ as.factor(year) + sector + emp + capital + output + (1 |  
#>     firm)
#>    Data: EmplUK
#> 
#> REML criterion at convergence: 4841.1
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -6.5171 -0.5067 -0.0280  0.4502  7.9279 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  firm     (Intercept) 21.489   4.636   
#>  Residual              3.793   1.948   
#> Number of obs: 1031, groups:  firm, 140
#> 
#> Fixed effects:
#>                     Estimate Std. Error t value
#> (Intercept)         26.57338    1.55552  17.083
#> as.factor(year)1977 -1.99478    0.27950  -7.137
#> as.factor(year)1978 -2.63032    0.27979  -9.401
#> as.factor(year)1979 -2.51523    0.27983  -8.988
#> as.factor(year)1980 -2.33948    0.28846  -8.110
#> as.factor(year)1981 -1.37792    0.32335  -4.261
#> as.factor(year)1982 -0.33180    0.33480  -0.991
#> as.factor(year)1983  0.49106    0.37275   1.317
#> as.factor(year)1984  0.33148    0.43801   0.757
#> sector              -0.78632    0.14900  -5.277
#> emp                 -0.09098    0.02393  -3.801
#> capital              0.17433    0.05542   3.145
#> output               0.03097    0.01132   2.736
#> 
#> Correlation matrix not shown by default, as p = 13 > 12.
#> Use print(x, correlation=TRUE)  or
#>     vcov(x)        if you need it

Created on 2020-04-02 by the reprex package (v0.3.0)

This doesn't seem at all comparable to the Stata screen shots.

Pretty much the same result when calling lmer with REML = FALSE

suppressPackageStartupMessages(library(lme4)) 
suppressPackageStartupMessages(library(plm)) 

# load data from plm package

data("EmplUK", package="plm")

# multi-level mixed effects regression with random effects for firm

m4 <- lmer(wage ~ as.factor(year) + sector + emp + capital + 
             output + (1 | firm), data = EmplUK, REML = FALSE)

# summarize

summary(m4)
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> Formula: wage ~ as.factor(year) + sector + emp + capital + output + (1 |  
#>     firm)
#>    Data: EmplUK
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>   4842.4   4916.5  -2406.2   4812.4     1016 
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -6.5547 -0.5098 -0.0281  0.4528  7.9738 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  firm     (Intercept) 21.054   4.588   
#>  Residual              3.749   1.936   
#> Number of obs: 1031, groups:  firm, 140
#> 
#> Fixed effects:
#>                     Estimate Std. Error t value
#> (Intercept)         26.57625    1.54452  17.207
#> as.factor(year)1977 -1.99492    0.27790  -7.179
#> as.factor(year)1978 -2.63046    0.27818  -9.456
#> as.factor(year)1979 -2.51538    0.27822  -9.041
#> as.factor(year)1980 -2.33981    0.28680  -8.158
#> as.factor(year)1981 -1.37846    0.32148  -4.288
#> as.factor(year)1982 -0.33238    0.33287  -0.999
#> as.factor(year)1983  0.49026    0.37060   1.323
#> as.factor(year)1984  0.33059    0.43548   0.759
#> sector              -0.78633    0.14750  -5.331
#> emp                 -0.09100    0.02376  -3.830
#> capital              0.17440    0.05506   3.168
#> output               0.03094    0.01125   2.750
#> 
#> Correlation matrix not shown by default, as p = 13 > 12.
#> Use print(x, correlation=TRUE)  or
#>     vcov(x)        if you need it

Created on 2020-04-02 by the reprex package (v0.3.0)

All I can suggest is looking at the other arguments to lmer to see if any of them would show the results on a comparable basis.

Hi again, I posted the link to the data at the top of the post (or here). I am not sure what the ettiquette is, so if all data must be within the reprex/built in, I just added a similar example with built-in data below my main post. Thanks!

Thank you for the help!

1 Like