3

I'm currently working in a project, where I replicate a project that has been conducted with Stata. I, however, work with R.

The task is to estimate a multi level logit model. In Stata, the project was conducted with the gllamm package and its function.

My problem is now, that I'm unable to replicate what the gllamm function does in R. So far I used the glmer() function from the lme4 package.

Here's an example data set and these are the respective code snippets:

Stata:

gllamm y z1 z2 z3 z4 z5 ///
x1 x2 x3 x4 x5 x6 x7,  ///
i(rank) family(binomial) link(logit)

Here is the output:

enter image description here

And R:

library(haven)
library(lme4)

# Data
data <- read_dta("data.dta")

# Model with converence warning:
mod1 <- glmer(y ~ z1 + z2 + z3 + z4 + z5 + 
                x1 + x2 + x3 + x4 + x5 + x6 + x7 + 
                (1 | rank), 
             family = "binomial",
             data = dataw)

summary(mod1)

# Model with different maximizer:
mod2 <- glmer(y ~ z1 + z2 + z3 + z4 + z5 + 
                x1 + x2 + x3 + x4 + x5 + x6 + x7 + 
               (1 | rank), 
              family = "binomial",
              control=glmerControl(optimizer="nloptwrap", calc.derivs = FALSE),
             data = dataw)
summary(mod2)

Here is the summary(mod2) output:

enter image description here

Unfortunately, one can see that the results differ in respect to the size and certainty of the estimates.

  • Is there a possibility to get (approximately) the same results in R?
  • What does gllamm do differently than glmer?
  • Additional question: In a possible solution, is there an argument that would fit the pweight() argument of the gllamm function as well?
Penguin_Knight
  • 11,078
  • 29
  • 48
MNeumann
  • 33
  • 3

1 Answers1

2

You can try the GLMMadaptive package that fits the same model using the adaptive Gaussian quadrature rule by default, e.g.,

library("GLMMadaptive")
library("haven")

dataw <- read_dta("data.dta")

fm <- mixed_model(y ~ z1 + z2 + z3 + z4 + z5 + 
                      x1 + x2 + x3 + x4 + x5 + x6 + x7,
                  random = ~ 1 | rank, family = binomial(),
                  data = dataw, iter_EM = 0)

summary(fm)
Dimitris Rizopoulos
  • 17,519
  • 2
  • 16
  • 37