Estimating sigma squared with ppmlhdfe

[Katinka Larsen]

Hi,

I am currently writing my master thesis where I am using the ppmlhdfe command to estimate a poisson model with multiple fixed effects. I want to find sigma squared, so I can do a quasi-LR test (quasi-likelihood ratio test), however I am not sure how to obtain it. Normally, if using a regular poisson model (without fixed effects), I would instead use the glm command and then get the (1/df) Pearson in the output. Do you know if I can do something similar with the ppmlhdfe model?


Thanks in advance,

Katinka

[Andrew Musau]

ppmlhdfe is from SSC (FAQ Advice #12). You can compute the Pearson statistic and then multiply it by (1/df).

$$X^{2} = \sum_{i=1}^{n} \frac{(y_i - \text{ exp}\{\text{X}_i\hat{\beta}\})^{2}}{\text{ exp}\{\text{X}_i\hat{\beta}\}}.$$

Here is an example:

webuse ships, clear
glm accident op_75_79 co_65_69 co_70_74 co_75_79 i.ship, family(poisson) nolog
ppmlhdfe accident op_75_79 co_65_69 co_70_74 co_75_79, absorb(ship) d
predict double yhat, xbd
gen double p= ((accident-exp(yhat))^2)/exp(yhat) if e(sample)
qui sum p
di "(1/df) Pearson ="  %7.1f  r(sum)*(1/e(df))

Res.:

. glm accident op_75_79 co_65_69 co_70_74 co_75_79 i.ship, family(poisson) nolog

Generalized linear models                         Number of obs   =         34
Optimization     : ML                             Residual df     =         25
                                                  Scale parameter =          1
Deviance         =  139.0852637                   (1/df) Deviance =   5.563411
Pearson          =  132.0602162                   (1/df) Pearson  =   5.282409

Variance function: V(u) = u                       [Poisson]
Link function    : g(u) = ln(u)                   [Log]

                                                  AIC             =   7.498581
Log likelihood   = -118.4758775                   BIC             =   50.92625

------------------------------------------------------------------------------
             |                 OIM
    accident |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    op_75_79 |   .2928003   .1127466     2.60   0.009      .071821    .5137796
    co_65_69 |   .5824489   .1480547     3.93   0.000     .2922671    .8726308
    co_70_74 |   .4627844    .151248     3.06   0.002     .1663437    .7592251
    co_75_79 |  -.1951267   .2135749    -0.91   0.361    -.6137258    .2234724
             |
        ship |
          2  |    1.79572   .1666196    10.78   0.000     1.469151    2.122288
          3  |  -1.252763   .3273268    -3.83   0.000    -1.894312   -.6112142
          4  |  -.9044563   .2874597    -3.15   0.002    -1.467867   -.3410457
          5  |  -.1462833   .2351762    -0.62   0.534    -.6072202    .3146537
             |
       _cons |   1.308451   .1972718     6.63   0.000     .9218049    1.695096
------------------------------------------------------------------------------

. ppmlhdfe accident op_75_79 co_65_69 co_70_74 co_75_79, absorb(ship) d
Iteration 1:   deviance = 1.6167e+02  eps = .         iters = 1    tol = 1.0e-04  min(eta) =  -2.14  P  
Iteration 2:   deviance = 1.3971e+02  eps = 1.57e-01  iters = 1    tol = 1.0e-04  min(eta) =  -2.56      
Iteration 3:   deviance = 1.3909e+02  eps = 4.46e-03  iters = 1    tol = 1.0e-04  min(eta) =  -2.69      
Iteration 4:   deviance = 1.3909e+02  eps = 1.15e-05  iters = 1    tol = 1.0e-04  min(eta) =  -2.70      
Iteration 5:   deviance = 1.3909e+02  eps = 2.66e-10  iters = 1    tol = 1.0e-05  min(eta) =  -2.70   S O
------------------------------------------------------------------------------------------------------------
(legend: p: exact partial-out   s: exact solver   h: step-halving   o: epsilon below tolerance)
Converged in 5 iterations and 5 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs      =         34
Absorbing 1 HDFE group                            Residual df     =         25
                                                  Wald chi2(4)    =       8.14
Deviance             =  139.0852637               Prob > chi2     =     0.0866
Log pseudolikelihood = -118.4758775               Pseudo R2       =     0.6674
------------------------------------------------------------------------------
             |               Robust
    accident |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    op_75_79 |   .2928003   .2736869     1.07   0.285    -.2436162    .8292168
    co_65_69 |   .5824489   .3063304     1.90   0.057    -.0179477    1.182846
    co_70_74 |   .4627844   .3591039     1.29   0.197    -.2410462    1.166615
    co_75_79 |  -.1951267   .3849763    -0.51   0.612    -.9496664     .559413
       _cons |    2.48605   .3500144     7.10   0.000     1.800034    3.172066
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
        ship |         5           0           5     |
-----------------------------------------------------+

.
. predict yhat, xbd

.
. gen double p= ((accident-exp(yhat))^2)/exp(yhat)
(6 missing values generated)

.
. qui sum p

.
. di "(1/df) Pearson ="  %7.1f  r(sum)*(1/e(df))
(1/df) Pearson =    5.3

.

[Katinka Larsen]

Thanks a lot! What does %7.1f mean in the last command?

[Andrew Musau]

The %f format

In %w.df, w is the total output width, including sign and decimal point, and d is the number of digits to appear to the right of the decimal point. The result is
right-justified.

The number 5.139 in %12.2f format displays as

----+----1--
5.14

So I am specifying that the result should be displayed with output width 7 and one digit to the right of the decimal point. See

help format