*cd [specify preferred folder path here]
*use [integrated longitudinal (wide) dataset here], clear

*Some formatting
foreach var of varlist fertpref_1 fertpref_2 fertpref_3 mcp_1 mcp_2 mcp_3 age_1 age_2 age_3 educattgen_1 educattgen_2 educattgen_3 marstat_1 marstat_2 marstat_3 score_1 score_2 score_3{
	replace `var'=. if `var'>=90
}

foreach var of varlist fertpref_1 fertpref_2 fertpref_3{
	replace `var'=. if `var'==3
}

*Create a variable denoting whether modern contraceptive use is persistent or intermittent at phases 2 & 3
gen mcpt=0 if mcp_2==0 & mcp_3==0
replace mcpt=1 if (mcp_2==1 & mcp_3==0) | (mcp_2==0 & mcp_3==1)
replace mcpt=2 if mcp_2==1 & mcp_3==1
label define mcpt 0 "No use" 1 "Intermittent use" 2 "Persistent use"
label values mcpt mcpt

*Create variables denoting whether fertility preferences are persistent or intermittent at phases 1 & 2
gen fertpreft=0 if fertpref_1==1 & fertpref_2==1
replace fertpreft=1 if (fertpref_1==2 & fertpref_2==1) | (fertpref_1==1 & fertpref_2==2)
replace fertpreft=2 if fertpref_1==2 & fertpref_2==2
label define fertpreft 0 "Persistent preference for more children" 1 "Intermittent preference to stop/delay" 2 "Persistent preference to stop/delay"
label values fertpreft fertpreft

gen fertpreft2=0 if fertpref_1==1 & fertpref_2==1
replace fertpreft2=1 if fertpref_1==2 & fertpref_2==1
replace fertpreft2=2 if fertpref_1==1 & fertpref_2==2
replace fertpreft2=3 if fertpref_1==2 & fertpref_2==2
label define fertpreft2 0 "Persistent preference for more children" 1 "Lost preference to stop/delay" 2 "Gained preference to stop/delay" 3 "Persistent preference to stop/delay"
label values fertpreft2 fertpreft2

*ANALYSES
*Question: do stated fertility preferences predict future adoption of modern contraceptives?
*Denominator: currently partnered women who are not pregnant and not currently using modern contraception

*In this initial analysis, we're only using data from phases 1 and 3, and can therefore use panel1_3_weight
svyset [pw=panel1_3_weight]
svy: logit mcp_3 i.country i.fertpref_1 i.educattgen_1 age_1 score_1 if mcp_1==0 & pregnant_1==0 & (marstat_1==21 | marstat_1==22), or
margins, at(fertpref_1==1) at(fertpref_1==2)
marginsplot, recast(scatter) xscale(range(0.5 2.5)) ///
	title("Predicted Probability of Using Modern Contraception (Phase 3)""by Fertility Preferences (Phase 1) (with 95% CI)") ///
	ytitle("Predicted Probability") xtitle("") ///
	xlabel(1 "Have Another Child" 2 "No More Children") scheme(s1mono)
graph export figure1-webinar.pdf, replace
	
*In these analyses, we're using data from all 3 phases and therefore rely on fullpanelweight
svyset [pw=fullpanelweight]
svy: logit mcp_3 i.country i.fertpreft2 i.educattgen_1 age_1 score_1 if mcp_1==0 & pregnant_1==0 & (marstat_1==21 | marstat_1==22), or
margins, at(fertpreft2=0) at(fertpreft2=1) at(fertpreft2=2) at(fertpreft2=3)
marginsplot, recast(scatter) xscale(range(0.5 4.5)) ///
	title("Predicted Probability of Using Modern Contraception (Phase 3)""by Fertility Preferences (Phases 1 and 2) (with 95% CI)") ///
	ytitle("Predicted Probability") xtitle("") ///
	xlabel(1 "PRO-PRO" 2 "CON-PRO" 3 "PRO-CON" 4 "CON-CON") scheme(s1mono)
graph export figure2-webinar.pdf, replace

svy: mlogit mcpt i.country i.fertpref_1 i.educattgen_1 age_1 score_1 if mcp_1==0 & pregnant_1==0 & (marstat_1==21 | marstat_1==22), rr
margins, at(fertpref_1==1) at(fertpref_1==2) predict(outcome(1)) predict(outcome(2))
marginsplot, recast(scatter) xscale(range(0.5 2.5)) ///
	title("Predicted Probability of Using Modern Contraception (Phases 2 and 3)""by Fertility Preferences (Phase 1) (with 95% CI)") ///
	ytitle("Predicted Probability") xtitle("") ///
	xlabel(1 "Have Another Child" 2 "No More Children") ///
	legend(order(1 "Intermittent Use" 2 "Persistent Use")) scheme(s1mono)
graph export figure3-webinar.pdf, replace