Equivalent to glmer (R) in stata? for a GLM mixed models (Gamma distribution)

[Rodrigo Badilla]

Hi all,

I am trying to get regression parameters from for a simple experiment for time response with a mixed model (person as random effect), I get a lot of heterocedasty and normality residuals problems, for this reason now I focus my effort in a analysis with GLM mixed model (family: Gamma(link=log).

I was trying to get this regression with R with glmer (lme4 library). I would try to replicate my results with stata (14), I tried with meglm but my results are very differents...

I am using a wrong command? exist a equivalent to glmer (R) in stata?

Thanks in advance.

[Scott Merryman]

I am not familiar with lme4:glmer, but it appears using the nAGQ = 0 produces similar results. The example below uses E.F. Haghish's -rcall- (see https://github.com/haghish/rcall)

. webuse pig,clear
(Longitudinal analysis of pig weights)

. meglm weight week || id:, family(gamma) nolog

Mixed-effects GLM                               Number of obs     =        432
Family:                   gamma
Link:                       log
Group variable:              id                 Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =    9540.82
Log likelihood = -1194.0189                     Prob > chi2       =     0.0000
------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |   .1319252   .0013506    97.68   0.000     .1292781    .1345724
       _cons |   3.203014   .0131258   244.03   0.000     3.177288     3.22874
-------------+----------------------------------------------------------------
       /logs |  -2.648873   .0360505                      -2.71953   -2.578215
-------------+----------------------------------------------------------------
id           |
   var(_cons)|   .0055246   .0012417                      .0035562    .0085825
------------------------------------------------------------------------------
LR test vs. gamma model: chibar2(01) = 207.77         Prob >= chibar2 = 0.0000

. rcall clear
(R memory cleared)

. rcall: library(lme4)

. rcall: weight <-st.var(weight)

. rcall: week <- st.var(week)

. rcall: id <- st.var(id)

. rcall: m<-glmer(weight ~ week + ( 1 | id) ,family=Gamma(link="log"), nAGQ=0 )

. rcall: summary(m)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: Gamma  ( log )
Formula: weight ~ week + (1 | id)

     AIC      BIC   logLik deviance df.resid 
  2332.9   2349.2  -1162.5   2324.9      428 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7677 -0.5597  0.1470  0.6909  2.1890 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.002364 0.04862 
 Residual             0.005548 0.07448 
Number of obs: 432, groups:  id, 48

Fixed effects:
            Estimate Std. Error t value Pr(>|z|)    
(Intercept) 3.202922   0.010500  305.04   <2e-16 ***
week        0.131914   0.001388   95.04   <2e-16 ***

[Rodrigo Badilla]

Thanks Scott por you reply, I fix my problem both method gave me the same results...