2
$\begingroup$

I'm running a fixed-effects Poisson regression and get different results in Stata and R. Unfortunately I cannot upload and share the data due to legal restrictions.

My code in R is (formula shortened for illustration):

library(pglm)    
pdf <- pdata.frame(data, index=c("id","timevar"))    
model <- pglm(y ~ x1 + x3 + x3_lag + x3_lag2 + season + x4 + x1*x3_lag + x1*x3_lag2+x1*x4,
              data = pdf,
              effect = "individual",
              model = "within",
              family = "poisson")

where season controls for seasonal effects, Y is a count variable and x3 is a treatment.

summary(model) yields this:

Maximum Likelihood estimation
Newton-Raphson maximisation, 4 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -47369.66 
19  free parameters
Estimates:
                   Estimate Std. error t value Pr(> t)
x1            -1.031e-04        Inf       0       1
x3            -1.196e-02        Inf       0       1
x3_lag        -4.783e-02        Inf       0       1
x3_lag2       -5.159e-02        Inf       0       1
season02      -7.038e-02        Inf       0       1
season03       9.323e-02        Inf       0       1
season04       1.257e-01        Inf       0       1
season05       1.427e-01        Inf       0       1
season06      -1.217e-01        Inf       0       1
season07      -1.566e+01        Inf       0       1
season08      -2.095e-01        Inf       0       1
season09      -1.886e-01        Inf       0       1
season10       4.488e-02        Inf       0       1
season11      -9.954e-02        Inf       0       1
season12       8.201e-02        Inf       0       1
x4             6.055e-01        Inf       0       1
x1:x3_lag      1.888e-05        Inf       0       1
x1:x3_lag2     3.529e-05        Inf       0       1
x1:x4          4.948e-04        Inf       0       1

In Stata14 I used:

xtset id timevar
xtpoisson x1 x3 x3_lag x3_lag2 season x4 x1#x3_lag x1#x3_lag2 x1#x4, fe

The results in Stata (see below) yield standard errors and where the standard errors are significant the coefficients are very much the same as in R, so I assume in both cases the same model was calculated. The issue is: why do I get standard errors in Stata but "Inf" in R?

Conditional fixed-effects Poisson regression    Number of obs     =    110,233
Group variable: id                              Number of groups  =     15,945

                                                Obs per group:
                                                              min =          2
                                                              avg =        6.9
                                                              max =         13

                                                Wald chi2(19)     =     816.49
Log likelihood  = -47369.663                    Prob > chi2       =     0.0000

--------------------------------------------------------------------------------------
           y         |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------------+----------------------------------------------------------------
                  x1 |  -.0001031   .0002119    -0.49   0.626    -.0005184    .0003121
                  x3 |  -.0119605   .0058094    -2.06   0.040    -.0233467   -.0005744
              x3_lag |  -.0478287   .0122489    -3.90   0.000    -.0718361   -.0238212
                     |
           x3_lag#x1 |   .0000189   .0000225     0.84   0.401    -.0000252     .000063
                     |
             x3_lag2 |  -.0515859   .0138096    -3.74   0.000    -.0786523   -.0245195
                     |
          x3_lag2#x1 |   .0000353   .0000245     1.44   0.149    -.0000127    .0000833
                     |
              season |
                 02  |  -.0703797   .0211975    -3.32   0.001    -.1119261   -.0288333
                 03  |    .093219   .0233858     3.99   0.000     .0473836    .1390544
                 04  |   .1256642   .0270927     4.64   0.000     .0725635    .1787649
                 05  |   .1426421   .0335767     4.25   0.000     .0768331    .2084512
                 06  |  -.1217462    .046838    -2.60   0.009     -.213547   -.0299454
                 07  |  -18.24382   741.5026    -0.02   0.980    -1471.562    1435.075
                 08  |  -.2095094   .0373193    -5.61   0.000    -.2826539   -.1363649
                 09  |  -.1886333   .0341701    -5.52   0.000    -.2556054   -.1216611
                 10  |    .044879   .0278778     1.61   0.107    -.0097604    .0995184
                 11  |  -.0995352   .0249665    -3.99   0.000    -.1484686   -.0506018
                 12  |   .0819983   .0214643     3.82   0.000      .039929    .1240675
                     |
                  x4 |   .6055111   .0754535     8.02   0.000      .457625    .7533972
                     |
               x4#x1 |   .0004948   .0001416     3.49   0.000     .0002173    .0007723
--------------------------------------------------------------------------------------

If you have any suggestions on how I could provide a replicable example (maybe with some sample panel data, any suggestions?), I will gladly do so.

EDIT

A smaller model produces identical results in R and Stata. It seems that with poor models R is more restrictive and does not show standard errors.

$\endgroup$

migrated from stackoverflow.com Jun 23 '16 at 14:16

This question came from our site for professional and enthusiast programmers.

  • 2
    $\begingroup$ I can't see this as a programming problem, as you have code and are not asking for programming advice. It's more of a statistical question of why a model fit is produced in one program but not another. Achieving a good fit in one program and not another is often a sign of a fragile model as much as, if not more than, a comment on the software. Although it often attracts flak, my suggestion is to try fitting a simpler model in R: the Stata results suggest which predictors may be omitted. As you realise the lack of a reproducible example constrains comment. $\endgroup$ – Nick Cox Jun 8 '16 at 16:49
  • $\begingroup$ Voting for migration to Cross Validated. $\endgroup$ – Nick Cox Jun 8 '16 at 16:50
  • $\begingroup$ You are right, simpler models do not produce this problem. $\endgroup$ – yoland Jun 8 '16 at 17:07
  • $\begingroup$ As I don't think I will get any help with my model fit without uploading the data, I will mark this post as answered - even though the exact reasons for the disparity remain unclear. Thank you! $\endgroup$ – yoland Jun 8 '16 at 17:11
  • 2
    $\begingroup$ The R code issued a warning ("Return Code 2") that should be heeded: it did not actually find a solution. Until you can develop a reproducible example, I doubt we can do any more about the problem here on CV than they could on SO. However, you could inspect the variance-covariance matrix produced by the model in R: that could be revealing. $\endgroup$ – whuber Jun 23 '16 at 14:28

Your Answer

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