Problem with graphing interrupted time series with random intercept -melogit-

[Dana Fletcher]

Dear Stata Gurus,

I am having a problem with correctly graphing output of an interrupted time series conducted using -melogit- with random intercept, and I would greatly appreciate some guidance. I’ve tried what feels like hundreds of configurations, and yet I just can't quite get it right.

This analysis is run on Stata v 15.1 in Windows 10 and is an evaluation of the association between a reimbursement change and utilization. Patients are clustered by hospital. There is a 4 month lead-in period before the payment change, during which patients were excluded.


Problems: When I try to plot both fixed and posterior means for the random effects (which is what I actually want to show), the pre-period lines for my predicted and counterfactual estimates do not overlap AND the y-intercept of the plots fall well outside the range of the unadjusted data.

Note that when I plot only fixed effects--as a benchmark for troubleshooting the graph, but not as my main desired output--my counterfactual and predicted plots in the pre-period lie on top of one another, as I would expect, but the plot still falls well below the unadjusted data.

My questions are:

1. What is Stata actually doing with this code?
2. What (and how) should I be asking it to plot the predicted values and counterfactuals correctly?
3. Is there a data characteristic (rather than operator naivete), that would explain why data would behave this way?

Details:
Elapsed_month is the time variable, 0 to 40, and represents the slope of the pre-period
Policy is the intervention, 0 or 1, and represents the level shift at the time the intervention was implemented
Policy#elapsed_month represents the change in the slope between the pre-period and the post period
Cohort4num is the binary dependent variable representing utilization of the treatment as 0,1

See code and resulting graphic below.

The -melogit- model and excerpt of the output:

[CODE]

melogit cohort4num c.elapsed_month i.policy c.elapsed_month#i.policy cr_onsite `patient'
`hospital' `market' `season', || new_hosp_num:, vce(robust)
[CODE]

Below is one way I've tried to graph the index functions:

********************************************************************************
***USING AVERAGES OF ALL COVARIATES***
*saves original covariate values, calculates mean, an puts that in the covariate field for each x

foreach x of varlist   cr_onsite `patient' `hospital' `market'  `season' {
gen `x'_save=`x'
egen `x'_mean=mean(`x') 
replace `x'=`x'_mean
}
predict average_pred7, mu

*******************************************************************************
**ESTIMATE RANDOM EFFECTS MEANS FOR COVARIATE MEANS
predict r7e, reffects
egen r7e_mean=mean(r7e)
summ r7e

***save so you can use the actual variable later ****/
gen elapsed_month_save=elapsed_month

*******************************************************************************
****ESTIMATE INDEX FUNCTIONS WITHOUT POSTERIOR MEANS RANDOM EFFECTS
forval t=0/40 {
replace elapsed_month = `t'

predict xb7_`t', xb 
}

gen mu7_xb =1 / (1+exp(-1*(xb7_1))) if elapsed_month_save == 0
****mu7_xb is the predicted probability at the mean of  covariates at time=0


****However, we need to compute at each value of time:
forval t=0/40 {
replace mu7_xb =1 /(1+exp(-1*(xb7_`t')))   if elapsed_month_save==`t'

}
*Here mu7_xb is the predicted probability for each time period based on x_mean*b^hat.


********************************************************************************
***TO INCLUDE RANDOM EFFECTS MEANS
***Recalculate mu7_xb with the re:

gen mu7_xb_re=1 /(1+exp(-1*(xb7_1+ r7e_mean)))   if elapsed_month_save==0

**** mu7_xb _re is the predicted probability at the mean of  covariates and mean of random effects  at time=0

****Recompute at each value of time:
forvalues t=0/40{
replace mu7_xb_re =1 /(1+exp(-1*(xb7_`t'+ r7e_mean )))   if elapsed_month_save==`t'

}

Next, is the same approach used to predict the counterfactual:

***********************
****COUNTERFACTUAL*****

*reset elapsed_month back to original
foreach x in elapsed_month {
replace `x'=`x'_save
drop elapsed_month_save
}
********************************************************************************    
*calculate counterfactual without effect of policy or policyXtime    
foreach x in  policy elapsed_month  {
gen `x'_save=`x'
replace `x'=0
}

predict  average_pred8, mu

*******************************************************************************
** PREDICT RANDOM EFFECTS MEANS [FOR COUNTERFACTUAL]
predict r8e, reffects
egen r8e_mean=mean(r8e)


*********************************************************************************
****PREDICT INDEX FUNCTIONS WITHOUT POSTERIOR MEANS RANDOM EFFECTS
forval t=0/40 {
replace elapsed_month = `t'
 
predict xb8_`t', xb
}

***The above is the index function without the random effects. 

gen mu8_xb =1 / (1+exp(-1*(xb8_1))) if elapsed_month_save==0

****mu8 is the predicted probability [OF THE COUNTERFACTUAL] at the mean of  covariates at time=0
****However, you need to compute at each value of time:
forval t=0/40 {
replace mu8_xb =1 /(1+exp(-1*(xb8_`t')))    if elapsed_month_save==`t'
}

*Here mu8_xb is the predicted probability [FOR THE COUNTERFACTUAL] for each time period based on x_mean*b^hat.

*******************************************************************************
***TO INCLUDE RANDOM EFFECTS MEANS [FOR THE COUNTERFACTUAL]
***Recalculate mu8_xb with the re:

gen mu8_xb_re=1 /(1+exp(-1*(xb8_1+ r8e_mean)))   if elapsed_month_save==0

**** mu8_xb _re is the predicted probability [OF THE COUNTERFACTUAL] at the mean of  covariates and mean of random effects  at time=0

****Recompute at each value of time:

forvalues t=0/40{
replace mu8_xb_re =1 /(1+exp(-1*(xb8_`t'+ r8e_mean )))   if elapsed_month_save==`t'

}

Which produces this graph, where the counterfactual (orange) and predicted index (blue) plots do not overlap in the pre-period when both the fixed and posterior means of the random effects are included, and the y-intercept lies much lower than the data.

Thank you for any advice you can provide!!

[Dana Fletcher]

Hi All,
I would like to revise my question. I revisited this post from 2014 (Interaction Effects and Models from OP Victoria Limon) because I suspect at least part of my answer resides here.

The part that struck me newly when I re-read this was, "The group (newid) level variance shows up in the graph as this kind of noise. If you want a graph of the more systematic relationship, specify the -xb- option of -predict-, and then you can transform it to the probability metric with the invlogit() function before graphing it."

I have been under the impression that the most appropriate presentation of this model was to include the random effects, but am I interpreting this statement correctly in that to represent the relationships between the pre-period trends, the policy implementation, and the post-period trends, rather than the noise or the effects of the groups (hospitals), it is acceptable to present the -xb- predictions?



When I change the code to simply read:

**USING AVERAGES OF ALL COVARIATES***
*saves original covariate values, calculates mean, an puts that in the covariate field for each x
foreach x of varlist   cr_onsite `patient' `hospital' `market'  `season' {
gen `x'_save=`x'
egen `x'_mean=mean(`x')
replace `x'=`x'_mean
}

*********************************************************
*PREDICT XB THEN INVLOGIT
predict xb7_xb, xb
gen phat7 = invlogit(xb7_xb)
twoway (line  phat7 elapsed_month)


*******************************************************
****COUNTERFACTUAL*****
*******************************************************    
*calculates counterfactual without effect of policy or policyXtime    
foreach x in  policy policyXmontime2  {
gen `x'_save=`x'
replace `x'=0
}

******************************************************
*PREDICT XB THEN INVLOGIT
predict xb8_xb, xb
gen phat8 = invlogit(xb8_xb)
twoway (line  phat8 elapsed_month)



twoway (line  phat7 elapsed_month, lcolor(red)) ||(line phat8 elapsed_month, lcolor (blue))

Then I get a very similar graph as to when I calculated only fixed effects (or if I model it using -logit- with clustered SE.)

phat7 is the index function, phat8 is the counterfactual, the lead-in period is not marked on this version of the graph, but is from month 19 to month 24.




I still don't know how to interpret why the predicted values fall so far outside the range of the unadjusted data points. For reference, the scatter plot is the number of eligible patients who participated in treatment / the total number of eligible patients in each month.

Thank you!!!

P.S. and my apologies for the obnoxiously large graphics--I haven't quite figured out how to size them yet.

[Dana Fletcher]

Hello All,
I think I have come to a resolution of how to present / interpret these predicted values, and that relates to presenting a more complete picture of the results by including and explaining both derivations of the graphs. Thanks to the StataList platform for helping me to formalize my question and answer it myself! Having said that, if anyone has any thoughts on presentation of this type of chart, interpretation of random effects in a time series, I'd love to hear it.