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]
Array


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

Code:
********************************************************************************
***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:
Code:
***********************
****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.

Array

Thank you for any advice you can provide!!