newey-west correction for count data?

[Gabi3]

Hi,

This is my first time posting a question, so I'm a bit nervous 🙂 

I am a clinician not a statistician, so please excuse any naivete. I am modelling an interrupted time series analysis with 2 breakpoints, using GLM and an impact model with 2 slope changes and 2 level changes, no other covariates included. The dependent variable is a count, Poisson distributed. I found autocorrelation at lag 3 only. I'm using SAS 9.4.

I read about using the Newey-West correction for autocorrelation in time series, however I wasn't sure if that was appropriate for a count distribution? All the SAS codes I've seen were applied to continuous DV, so I essentially adapted those codes to my count model - but i'm not sure this is right...more specifically I'm asking:

- is newey west correction appropriate for a count model?

- if it is, did I code it correctly? (I've essentially taken the newey-west correction i've seen used in a proc model but for a continuous DV and applied it to mine)

There are no errors that come up when i run this or anything else weird but the reason I'm asking is that my mentor is doing the same analysis using STATA and is getting VERY different standard errors. So we're trying to figure out if I made a mistake in my code. 

Thanks very much!

Gabriela

Explanation of variables:

slopech1_2 = slope change from period 1 to 2

jump1_2 = level change from period 1 to 2 (same for periods 2 to 3)

time = month1 to 50 of 50 months (monthly measurements)

Exposure = exposure variable for each month

y= observed variable of interest (rate)

proc model data=Trends;

label b0="intercept" b1="time" b2="slopech1_2" b3="jump1_2" b4="slopech2_3" b5="jump2_3";

parms b0 = 0.5 b1 = 0.1 b2=0.1 b3=0.1 b4=0.1 b5=0.1;

yhat=exp(b0+b1*time+b2*slopech1_2+b3*jump1_2+b4*slopech2_3+b5*jump2_3)*Exposure;

y=yhat;

lk=exp(-yhat)*(yhat**y)/fact(y);

ll=-log(lk);

errormodel clots~general(ll);

fit clots/gmm kernel=(bart,4,0) vardef=n;

run;

Attached is a sample of the dataset. 

[SteveDenham]

I don't know if it is a mistake, but it seems unusual that these two lines are in there:

lk=exp(-yhat)*(yhat**y)/fact(y);

ll=-log(lk);

That would imply that the log likelihood is just

lk=(yhat)*(yhat**y)/fact(y);

and that doesn't seem right somehow.

SteveDenham

[Rick_SAS]

Although Steve's computations are a bit off

(I believe the correct computation is ll = -yhat + y*log(yhat) - LFACT(y)?)

his main point is valid: check the log-likelihood function.  Furthermore, for numerical stability, it is best to avoid forming the full likelihood and then taking the log. Instead, compute the log-likelihood directly, as shown in this article:

Two simple ways to construct a log-likelihood function in SAS

Some tips on specifying lo-likelihoods are available at 

Two ways to compute maximum likelihood estimates in SAS

This article talks about PROC NLMIXED, but the tips apply to other SAS procedures, too.

[SteveDenham]

Thanks, @Rick_SAS .  I missed that there was no closing parenthesis, so that the exp(yhat) is the only term that transforms back to the original scale in this representation.

So, this code calculates the log-likelihood as -yhat+y*log(yhat)-LFACT(y), and thus is off by a factor of y/2 after removing the identical parts of the calculation (check me here, Rick).

SteveDenham

[Rick_SAS]

I don't know the correct LL, but the OP's program has the line

y=yhat;

which seems very strange. The OP's comments indicate that y should be an observed variable (the count).

[Gabi3]

Hi,

Thanks guys!! I'm sorry to ask this, but I am a clinician so not really
mathematically savvy. Could you re-write the code with the proper syntax
for me? It would take me quite a bit of time to figure it out and would be
of great help (since i'm sure its super easy for you).

Also, any thoughts on the newey correction? Both in regards to the
applicability to the count model AND the syntax?

Thanks very much, appreciate it,

Gabi

[Ksharp]

You might post it at Forecasting Forum , since PROC MODEL is under SAS/ETS .

[Gabi3]

Thanks, good idea.
Is there a quick way of doing this, or should I copy / paste the entire thing into a different forum?

thanks
Gabi

[Ksharp]

Yeah. Post your question at this forum. Some expert about SAS/ETS might give you a hand.

https://communities.sas.com/t5/SAS-Forecasting-and-Econometrics/bd-p/forecasting_econometrics