Clustering standard error for hierarchical data

[jgreenberg321]

I am working with a large healthcare dataset from 15 states to conduct a difference-in-difference analysis. I have questions related to appropriate treatment of standard errors and interpretation of genmod results. I am attempting to develop hierarchical fixed effect regression models predicting the impact of a policy intervention on treatment frequency (procedures per 100,000 people), adjusting for patient-level confounders. However, I want to ensure that the standard errors appropriately reflect clustering at the state level. The dataset is not a sampling or survey data, but rather complete capture from each state. Reading through the different posts on this site, I have seen some people advocate using the absorb statement in proc glm, but I am unsure if this produces accurate standard errors. Does anyone know if that is the case?

Alternatively, others advocate using proc genmod, with example variables and code below:

outcome=treatment_freq

intervention=policy (1=intervention; 0=no intervention)

time-period (pre vs. post-intervention)=period (1=post-intervention, 0=pre-intervention)

confounders=conf1, conf2, etc.

State variable=state

proc genmod data=have;

class policy (ref='0') period (ref='0') conf1 (ref='0') conf2 (ref='0') state;

model treatment_freq= conf1 conf2 policy period policy*period;

repeated subject= state/type=ind;

run;

This produces the output shown in the attachment. Am I correct in interpreting the policy*period parameter estimate of 0.0434 as reflecting a 4.34% increase in treatment frequency associated with the intervention, with the standard errors accounting for clustering at the state level? 

Thank you!!

[SteveDenham]

There are two considerations here that make your interpretation a bit off, in my opinion.  First, all of the parameters are for the model on the logit scale, so the 0.0434 doesn't directly convert to a 4.34% change due to the intervention.  Second, the 0.0434 represents the additional change due to the policy intervention at the post-treatment period as compared to the pre-treatment period, above the 0.0137 change due to just post-treatment vs pre-treatment.

You probably ought to consider using either LSMEANS with the diff and ilink options.  This would give the difference between the groups on the original scale, with the confounding variables set at the reference value. To get differences at various combinations of the confounding variables, you would likely need to change your model to include interactions.  The point here is that the transformation back to the original scale is non-linear, so the 0.0434 will have different effects on the predicted proportions depending on the specific combination of confounders.  A plot of the lsmeans on the original scale as a function of the values of the confounders would probably make this clearer than my words.

SteveDenham.

[jgreenberg321]

Thank you very much for the reply. I just wanted to clarify a couple of points related to your reply:

1) My understanding was that, with the code shown, proc genmod would fit an OLS model and not a logistic regression model (which I thought required the link=logit statement). Your point about the interpretation of parameter estimates on the logit scale for a logistic regression is well taken, but my understanding is that proc genmod would produce a linear rather than logistic regression model with the specifications I showed. Is that not correct?

2) I can see how the LSmeans statement would help me obtain mean outcome values for different levels of the interaction, standardized for the other confounders in the model. However, my understanding is that for a difference-in-difference design (as I've attempted to implement), it is the coefficient of the interaction term that best captures the treatment effect. Otherwise, I would think I would need to manually compute a "difference in difference" between intervention and control groups in the pre and post-periods. Please let me know if this sounds off base, assuming I have reasonably conveyed what I have been trying to.

Thank you again for your time.

[SteveDenham]

I apologize.  Ordinarily, I don't use GENMOD for analysis of endpoints with normal errors.  As a result, your first point is right on the mark, and my concern about non-linearity was completely out of line. As a result, the interaction term you refer to is a very good measure of the additional effect and your conclusion of 4.3% makes sense.

So, with an endpoint with normal errors, why GENMOD for a hierarchical model, rather than PROC MIXED?  I suspect that your subjects are nested within states. I suppose if someone came to me with this design, I would not immediately think of a difference in differences approach, but rather a repeated measures analysis.  My approach would be to restructure the data set, with 'pre' and 'post' as levels of time, and their values put into a single response variable.  The code would look something like:

proc mixed data=have;

class policy (ref='0') period (ref='0') conf1 (ref='0') conf2 (ref='0') state subjid;

model treatment_freq= conf1 conf2 policy period policy*period/solution;

random state/subject=subjid;

repeated period/type=un subject=subjid

run;

This would model the covariance structure of period.  A difference in differences assumes complete independence after differencing, while this approach allows that assumption to be relaxed, and a heterogeneous variance per time period, with a non-zero covariance.  I am still pondering the confounders, and whether they should interact with period or not, but that is a subject matter decision.  In the end, the policy*period interaction coefficient will give the expected change in response above the period effect.

Thanks for pointing out where I was making my mistake.  

SteveDenham

[jgreenberg321]

Thank you for the quick and detailed reply. These are interesting and very helpful points.

Per your major point in #2: I do not have a strong rationale for using proc genmod compared to proc mixed. My rationale is partly that it seems that in the policy/difference-in-difference literature, people often use fixed effects rather than mixed models. I assume this is because fixed effects models require fewer assumptions, and I am interested in primarily "accounting for" rather than modeling the cluster correlations.

My knowledge of mixed model development is somewhat limited, so I just wanted to check that the hypothetical code you gave would provide standard errors that account for the clustering both within states and study period? 

Finally, per your point about the confounder-period interactions: I can understand the rationale here and have seen it done in some but not all studies. Practically, I don't expect this to have a large impact on my analysis but something to explore further as a secondary consideration. Thank you again for your time and thoughtfulness! 

[SteveDenham]

As period appears in the MODEL statement, it is a fixed effect, and the REPEATED statement allows for modeling the covariance between the residuals as period changes, marginalized over the random effects (state).  Actually fewer assumptions than are needed for a difference in difference analysis, although the assumptions of normality and independence are still in place.  The point here is that the only random effect is 'state', and if the literature is not treating this as a random effect, there are likely a lot of policies being touted that are not truly significant when implemented across a broader population than the described subjects in the states included in the design.  This is just one of the causes of the "non-reproducibility" crisis, especially in medical trials.

SteveDenham

[jgreenberg321]

Thank you. That all makes a lot of sense. One minor clarification: in the code you provided, "subjid" would just be referring to a unique id assigned to each subject (to denote the level of individual observation)? Thank you again!

[SteveDenham]

That is correct. If you happened to have duplicate id's you would almost certainly have to augment them so that they are distinguishable from one another.

SteveDenham

[jgreenberg321]

That makes sense. Thank you for clarifying!

Jacob

[jgreenberg321]

When I tried implementing this code, I reached the following error, likely reflecting the fact I have hundreds of thousands of individual observations. Is that a limitation of this approach?

WARNING: Class levels for subjid are not printed because of excessive size.

Thank you and sorry for the hassle!

Jacob

[SteveDenham]

Add NOCLPRINT to the PROC MIXED statement.  That will suppress the printing of all the class levels, and is really handy in situations like yours.  It will turn off the WARNING: as well, and make the output a reasonable size.

SteveDenham

[jgreenberg321]

Thank you. That seems to have taken me from the previous error message to the following. Perhaps it is just a limitation of my system's memory to run the model with so many observations? Thank you!


ERROR: The SAS System stopped processing this step because of insufficient memory.

[SteveDenham]

yeah, I hate that error.  It turns up a fair amount when working with mixed models.  If you have the ability to edit your sasv9.cfg file, find  -MEMSIZE.  It is probably set to 2G.  Edit that to MAX. The next time you invoke SAS you should have access to as much memory as your system will allow (depending on what other programs you have running).  And if that doesn't solve the problem, consider proportionally sampling subjects from each state to get an analyzable number of records.  You could then "model average/bootstrap" across multiple realizations of the sampled data to get a parameter estimate and bootstrapped confidence limits on the interaction term.

You have run into one of the real sticking points for mixed effects models as opposed to wholly fixed effects models.  I feel like a used car salesman for not addressing that drawback.

SteveDenham

[jgreenberg321]

Ha! Not at all. Thank you for explaining. That is all extremely helpful to understand!