st: two-parts model bootstrap using XTGEE

["Li, Elena (Zhonghe)" <zli@hsph.harvard.edu>]

Hi
We are fitting a two-parts model on health care cost using longitudinal data.
Part1 is modeling whether there is any cost, and
Part2 is modeling the cost when if there is a cost.

The goal is to bootstrap the difference in difference (post_pilot) in terms of the p0*p1 (the predicted probability of having cost * the predicted cost).
I am a SAS user and tried to put the following Stata program together.  The models converged the first time but when the replication starts the modeling does not converge:
"convergence not achieved
an error occurred when bootstrap executed tempa, posting missing values"

Can you please help?  Thank you.
Zhonghe Li

_____________________________________________________________
sort member qtr
by member qtr: drop if _n~=1
egen memberid=group(member)
egen practice=group(pcp)

gen newid=memberid
xtset newid qtr
generate sample=(1-missing(havecost, logcost, female, aet, uhc, ant,  Qtr_months, age1-age4, qtr, health,  post, pilot, post_pilot))
keep if sample

matrix define results = (1,2,3)


capture program drop tempa

program define tempa, rclass


xtgee havecost i.female aet uhc ant  Qtr_months age1-age4 qtr health  i.post i.pilot i.post_pilot,  family(binomial) link(logit) corr(ar 1) 

    
	gen truepost=post
	gen truepilot=pilot
	gen true_post_pilot=post_pilot

	replace post=0
	replace pilot=0
	replace post_pilot=0
	predict pre_control_p1

	replace pilot=1
	predict pre_pilot_p1

	replace post=1
	replace post_pilot=1
	predict post_pilot_p1

	replace pilot=0
	replace post_pilot=0
	predict post_control_p1
	
	replace post=truepost 
	replace pilot=truepilot 
	replace post_pilot=true_post_pilot
	
xtgee logcost female aet uhc ant  Qtr_months age1-age4 qtr health  post pilot post_pilot if logcost>0,  family(gamma) link(log) corr(ar 1) vce(robust) 

		
	replace post=0
	replace pilot=0
	replace post_pilot=0
	predict pre_control_p2, mu

	replace pilot=1
	predict pre_pilot_p2, mu

	replace post=1
	replace post_pilot=1
	predict post_pilot_p2, mu

	replace pilot=0
	replace post_pilot=0
	predict post_control_p2, mu

	drop post pilot post_pilot 
	rename truepost post
	rename truepilot pilot
	rename true_post_pilot post_pilot
	

gen phat_pre_control=pre_control_p1*pre_control_p2
gen phat_pre_pilot=pre_pilot_p1*pre_pilot_p2
gen phat_post_pilot=post_pilot_p1*post_pilot_p2
gen phat_post_control=post_control_p1*post_control_p2

gen diff_pilot=phat_post_pilot-phat_pre_pilot
gen diff_control=phat_post_control-phat_pre_control
gen diff_in_diff=diff_pilot-diff_control

summ diff_pilot
return scalar m_diff_pilot = r(mean)
local var1 = r(mean)
 
summ diff_control 
return scalar m_diff_control = r(mean)
local var2 = r(mean)
 
summ diff_in_diff 
return scalar m_diff_in_diff = r(mean)
local var3 = r(mean)

matrix define temp=(`var1',`var2',`var3')
matrix define results = results \ temp
              
matrix list results
 
quietly: drop diff*  pre_* post_control_p1 post_control_p2 post_pilot_p1 post_pilot_p2 phat* 
*quietly: scalar drop _all

xtset,clear

end



bootstrap mean_diff_pilot=r(m_diff_pilot) mean_diff_control=r(m_diff_control) mean_diff_in_diff=r(m_diff_in_diff), noisily strata(pilot) cluster(memberid) idcluster(newid) group(qtr) reps(100) seed(1234567) : tempa

svmat results, names(diff)
outsheet twopartdiff* using diff.csv, comma replace

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/