*------------------------------------------------------------------------------*
* Fardousi et al. 2022
* PLOS Medicine
* Performance bonuses and the quality of primary health care delivered by 
* family health teams in Brazil: A difference-in differences analysis
*------------------------------------------------------------------------------*


*------------------------------------------------------------------------------*
* README
*------------------------------------------------------------------------------*
/* 
This Stata do file runs the analysis that was used to produce the results in 
the main text of Fardousi et al (2022). Some variables used in the analysis are 
not publicly available and are therefore not contained in the analytical 
dataset accompanying this do file. The restricted data are: the PMAQ scores for
each family health team (our measure of quality of care) and the responses 
from the survey of municipality managers about performance bonuses (our 
exposure variables). For these restricted data, we generate placeholder 
variables using a random number generator to allow the do file to run and for
illustrative purposes. Please note that the results contained within Fardousi 
et al (2022) can not be replicated without access to the restricted 
data. Requests for access to these data should directed at the Department for 
Family Health: telephone number +55 61 33159044 or email desf@saude.gov.br.
*/

*------------------------------------------------------------------------------*
* Preliminaries
*------------------------------------------------------------------------------*
version 16
clear all
set more off

global root "specify pathway"
global raw "$root\raw STATA"
global xls "$root\raw XLS"
global temp "$root\temp"
global clean "$root\clean"
global analysis "$root\analysis"
global out "$root\do-files\analysis\tables"

use "$analysis\pmaq_analytical_plos_med", clear
global f_controls "type1 staff_total"
global mun_control "mun_pmaqfunds gdp_pc hdi gini pop_total urban_perc  pop_share_u5 pop_share_o60 "

*------------------------------------------------------------------------------*
* Variable generation
*------------------------------------------------------------------------------*
svyset ibge
gen mun_n = 0
bysort ibge: replace mun_n = 1 if _n == 1
gen cnes_n = 0
bysort cnes: replace cnes_n = 1 if _n == 1


*------------------------------------------------------------------------------*
* Placeholder (random number generated) variables for restricted data
*------------------------------------------------------------------------------*
set seed 12345
g design_c3=uniform() < 0.5523 /*set at the proportion of the actual data */
set seed 12345
g bonus1=uniform() < 0.1949 if design_c3!=0 /*set at the proportion of the actual data */
set seed 12345
g bonus2=uniform() < 0.2967 if design_c3!=0 & bonus1!=1 /*set at the proportion of the actual data */
set seed 12345
g bonus3=uniform() < 0.2361 if design_c3!=0 & bonus1!=1 & bonus2!=1 /*set at the proportion of the actual data */
g bonus_amount=0 if design_c3==0
replace bonus_amount=1 if bonus1==1
replace bonus_amount=2 if bonus2==1
replace bonus_amount=3 if bonus3==1
drop bonus1 bonus2 bonus3
label define design 0 "Non-bonus municipalities" 1 "Bonus municipalities"
label values design_c3 design
label define amount 0 "No bonus" 1 "1 to 20%" 2 "21 to 50%" 3 "More than 50%"
label values bonus_amount amount

set seed 12345
g double score_c1=rnormal(60.68405,10.34623) if design_c3==0 /*set at mean and sd of the actual data */
set seed 12345
replace score_c1=rnormal(61.56753,8.831329) if bonus_amount==1 /*set at mean and sd of the actual data */
set seed 12345
replace score_c1=rnormal(61.07565,9.300989) if bonus_amount==2 /*set at mean and sd of the actual data */
set seed 12345
replace score_c1=rnormal(60.08015,10.09426) if bonus_amount==3 /*set at mean and sd of the actual data */
set seed 12345
replace score_c1=rnormal(61.94767,9.142843) if bonus_amount==. /*set at mean and sd of the actual data */
set seed 12345
g double score_c3=rnormal(59.4829,11.89947) if design_c3==0 /*set at mean and sd of the actual data */
set seed 12345
replace score_c3=rnormal(62.93029,12.9752) if bonus_amount==1 /*set at mean and sd of the actual data */
set seed 12345
replace score_c3=rnormal(66.51982,11.77422) if bonus_amount==2 /*set at mean and sd of the actual data */
set seed 12345
replace score_c3=rnormal(67.23538,10.01879) if bonus_amount==3 /*set at mean and sd of the actual data */
set seed 12345
replace score_c3=rnormal(63.35956,13.70382) if bonus_amount==. /*set at mean and sd of the actual data */
g pmaq_c3_c1=score_c3-score_c1
label var score_c1 "Placeholder PMAQ score in round 1"
label var score_c3 "Placeholder PMAQ score in round 3"
label var pmaq_c3_c1 "Change in PMAQ score (round 3 minus round 1)"


*------------------------------------------------------------------------------*
* Results
*------------------------------------------------------------------------------*
** Results behind Fig 1a
mean score_c1 score_c3 if sample_analytical==1, over(design_c3) vce(cluster ibge)
mean score_c1 score_c3 if sample_analytical==1, over(bonus_amount) vce(cluster ibge)

** Results behind Fig 1b
mean score_c1 score_c3 if sample_analytical==1 & sample_matched==1, over(design_c3) vce(cluster ibge)
mean score_c1 score_c3 if sample_analytical==1 & sample_matched==1, over(bonus_amount) vce(cluster ibge)

* Unadjusted matched diff-in-diff results reported in main text
mean pmaq_c3_c1 if sample_analytical==1 & sample_matched==1, over(design_c3) vce(cluster ibge)
mean pmaq_c3_c1 if sample_analytical==1 & sample_matched==1, over(bonus_amount) vce(cluster ibge)
reg pmaq_c3_c1 design_c3 if sample_analytical==1 & sample_matched==1, cluster(ibge)
reg pmaq_c3_c1 i.bonus_amount if sample_analytical==1 & sample_matched==1, cluster(ibge)

** Table 1
*Descriptive stats
tabstat score_c1 income_1000 if sample_analytical==1, by(design_c3) stat(mean sd n )
tabstat type1 staff_total if sample_analytical==1 & cnes_n==1, by(design_c3) stat(mean sd n )
tabstat $mun_control if sample_analytical==1 & mun_n==1, by(design_c3) stat(mean sd n )

tabstat score_c1 income_1000 if sample_analytical==1 & sample_matched==1, by(design_c3) stat(mean sd sum n )
tabstat mun_pmaqfunds $mun_control  if sample_analytical==1 & mun_n==1 & sample_matched==1, by(design_c3) stat(mean sd n )
tabstat type1 staff_total if sample_analytical==1 & cnes_n==1 & sample_matched==1, by(design_c3) stat(mean sd sum n )

*P values of difference
foreach var of varlist score_c1 income_1000 {
reg `var' design_c3 if sample_analytical==1, cluster(ibge)
}
foreach var of varlist staff_total {
reg `var' design_c3 if sample_analytical==1 & cnes_n==1, cluster(ibge)
}
svy: tab type1 design_c3 if sample_analytical==1 & cnes_n==1, col
foreach var of varlist $mun_control  {
reg `var' design_c3 if sample_analytical==1 & mun_n==1
}

foreach var of varlist score_c1 income_1000 {
reg `var' design_c3 if sample_analytical==1 & sample_matched==1, cluster(ibge)
}
foreach var of varlist staff_total {
reg `var' design_c3 if sample_analytical==1 & cnes_n==1 & sample_matched==1, cluster(ibge)
}
svy: tab type1 design_c3 if sample_analytical==1 & cnes_n==1 & sample_matched==1, col
foreach var of varlist $mun_control  {
reg `var' design_c3 if sample_analytical==1 & mun_n==1 & sample_matched==1
}

** Table 2
*Unmatched
reg pmaq_c3_c1 design_c3 i.income_q_c $f_controls $mun_control if sample_analytical==1, cluster(ibge)
outreg2 using table2.xls, replace label bdec(1) dec(1) pdec(4) stats(coef ci pval)  nocons sideway noas
reg pmaq_c3_c1 i.bonus_amount i.income_q_c $f_controls $mun_control if sample_analytical==1, cluster(ibge)
outreg2 using table2.xls, append label bdec(1) dec(1) pdec(4) stats(coef ci pval)  nocons sideway noas
mean score_c1 if sample_analytical==1 & bonus_amount!=.

*Matched
reg pmaq_c3_c1 design_c3  i.income_q_c_matched $f_controls $mun_control if sample_analytical==1 & sample_matched==1, cluster(ibge)
outreg2 using table2.xls, append label bdec(1) dec(1) pdec(4) stats(coef ci pval)  nocons sideway noas
mean score_c1 if sample_analytical==1 & sample_matched==1
reg pmaq_c3_c1 i.bonus_amount i.income_q_c_matched $f_controls $mun_control if sample_analytical==1 & sample_matched==1, cluster(ibge)
outreg2 using table2.xls, append label bdec(1) dec(1) pdec(4) stats(coef ci pval)  nocons sideway noas
mean score_c1 if sample_analytical==1 & sample_matched==1 & bonus_amount!=.

** Results behind Fig 2
*Unmatched, any bonus
reg pmaq_c3_c1 design_c3##i.income_q_c   $f_controls $mun_control if sample_analytical==1, cluster(ibge)
margins income_q_c, dydx(design_c3) atmeans
mean score_c1 if sample_analytical==1, over(income_q_c)
*Matched, any bonus
reg pmaq_c3_c1 design_c3##i.income_q_c_matched  $f_controls $mun_control if sample_analytical==1 & sample_matched==1, cluster(ibge)
margins income_q_c_matched, dydx(design_c3) atmeans
mean score_c1 if sample_analytical==1 & sample_matched==1, over(income_q_c_matched)
*Unmatched, bonus amount
reg pmaq_c3_c1 i.bonus_amount##i.income_q_c  $f_controls $mun_control if sample_analytical==1, cluster(ibge)
margins income_q_c, dydx(bonus_amount) atmeans
mean score_c1 if sample_analytical==1 & bonus_amount!=., over(income_q_c)
*Matched, bonus amount
reg pmaq_c3_c1 i.bonus_amount##i.income_q_c_matched  $f_controls $mun_control if sample_analytical==1 & sample_matched==1, cluster(ibge)
margins income_q_c_matched, dydx(bonus_amount) atmeans
mean score_c1 if sample_analytical==1 & bonus_amount!=. & sample_matched==1, over(income_q_c_matched)