Disabled Development:
Where are Disabled Women amid the Sustainable Development Agenda in Costa Rica?

************************************************************
*********************** LOADING DATA *******************
************************************************************

clear all
use "C:\Sync\Data\mics_00007.dta" , clear

/*This is an IPUMS MICS data extract, courtesy of ipums.org,
All samples, all rounds, life satisfaction vars, discrim vars,
demographics, technical vars, March 2023 
*/

************************************************************
***********************  WEIGHTING   *********************
************************************************************

svyset [pweight=weightwm]
*help svy_estimation

************************************************************
***********************  PROGRAMS   *********************
************************************************************
* I used Stata SE version 16.1.
* The following programs are needed to run this do-file
ssc install elabel, replace //by Daniel Klein 
ssc install _GWTMEAN, replace //by David Kantor
ssc install grstyle, replace //by Ben Jann
ssc install mylabels, replace //by Nick Cox
ssc install nicelabels, replace //by Nick Cox
ssc install egenmore, replace //by Nick Cox
ssc install gologit2, replace //by Richard Williams
ssc install estout, replace //by Ben Jann

************************************************************
***********************   CLEANING   *********************
************************************************************

// recoding don't know (7), no response (8) and NIU (9) missing
recode resultwm area glasses hearaid diff* urban discrimethnic ///
everbirth discrimsexorien discrimage discrimrelig discrimdis ///
discrimother discrimgender healthinsur pregnant lit (7/9=.)

// labeling all functional disability vars thru a foreach loop
*STUDENTS - SEE BLOG POST ON LOCAL MACROS IN STATA DO FILES
local diffvars "see hear walk remcon com care"
local difflabels 0 "no difficulty" 1 "some difficulty" ///
			2 "a lot of difficulty" ///
			3 "cannot do this function at all"

foreach var of local diffvars {
    recode diff`var' (1=0) (2=1) (3=2) (4=3)
    label define diff`var' `difflabels'
    label values diff`var' diff`var'
}

//recoding no response 98 and NIU 99 as missing
recode marst everbirth birth2yr agewm edlevelwm (98/99=.)

// recoding no response (998) and NIU (999) as missing 
recode husbandage (998/999=.)

// recoding no response (8 or 98) and NIU (9 or 99) as missing
recode lsladder (98/99=.)
rename lsladder ladder
recode ls* (7/9=.)
rename ladder lsladder

/*The Cantril ladder ordinal outcome levels are 0-10, but very
negatively skewed here, with only 10% of scores being in 
levels 0-5, which are thus combined in the lsbin variable.
*/
gen lsbin = lsladder
recode lsbin (0/5=5)
label define lsbin 5 "ladder 0-5" 6 "ladder 6" 7 "ladder 7" ///
 8 "ladder 8" 9 "ladder 9" 10 "ladder 10"
label values lsbin lsbin
tab lsbin lsladder //verify this transformation

//recoding "NIU or missing" (0) as missing 
recode windex5 (0=.)

//combining wealth quintiles to simplify this predictor for model fit
//consistent with the distribution
gen wealth = windex5
recode wealth (1/2=0) (3/5=1)
label define wealth 0 "lower wealth" 1 "higher wealth"
label values wealth wealth 
tab wealth windex5 //verify the transformation 

//reverse-coding overall happiness so higher numbers mean more happiness 
generate lshappyorig = lshappy //backup of the orig for checking
recode lshappy (5=1) (4=2) (3=3) (2=4) (1=5)
label define lshappy 1 "very unhappy" 2 "a little unhappy" ///
                     3 "neither unhappy nor happy" 4 "somewhat happy" ///
                     5 "very happy"
label values lshappy lshappy 
tab lshappyorig lshappy //verify the transformation

*reverse coding lslastyr so that higher numbers mean more satisfaction 
generate lslastyrorig = lslastyr //backup of the original for checking 
recode lslastyr (3=1) (1=3)
label define lslastyr 1 "worsened" 2 "more or less the same" 3 "improved"
label values lslastyr lslastyr 
tab lslastyrorig lslastyr //verify the transformation 

*reverse coding lsnextyr so that higher numbers mean more satisfaction 
/* IMPORTANT NOTE: This item is not ordinal with "depends on God" 
 in it; I still reverse coded to flip the "worse" and "better" to be
 low and high respectively */  
generate lsnextyrorig = lsnextyr //backup of the orig for checking
recode lsnextyr (3=1) (2=3) (4=2) (1=4)
label define lsnextyr 1 "worse" 2 "depends on God" ///
3 "more or less the same" 4 "better"
label values lsnextyr lsnextyr 
tab lsnextyrorig lsnextyr //verify the transformation 

/* relabeling MICS round 6 countries (only) w/their ISO3 code for use
 in graphs. Note that I labeled Costa Rica with its full name */
elabel list (country)
label define COUNTRY 12 "DZ", replace
label define COUNTRY 32 "AR", add 
label define COUNTRY 50 "BD", add 
label define COUNTRY 112 "BY", add 
label define COUNTRY 140 "CF", add 
label define COUNTRY 148 "TD", add
label define COUNTRY 180 "CD", add
label define COUNTRY 188 "Costa Rica", add
label define COUNTRY 192 "CU", add
label define COUNTRY 214 "DO", add
label define COUNTRY 268 "GE", add
label define COUNTRY 270 "GM", add
label define COUNTRY 288 "GH", add
label define COUNTRY 296 "KI", add
label define COUNTRY 328 "GY", add
label define COUNTRY 340 "HN", add
label define COUNTRY 368 "IQ", add
label define COUNTRY 417 "KG", add
label define COUNTRY 418 "LA", add
label define COUNTRY 426 "LS", add
label define COUNTRY 450 "MG", add
label define COUNTRY 454 "MW", add
label define COUNTRY 496 "MN", add
label define COUNTRY 499 "ME", add
label define COUNTRY 524 "NP", add
label define COUNTRY 586 "PK", add
label define COUNTRY 624 "GW", add
label define COUNTRY 678 "ST", add
label define COUNTRY 688 "RS", add
label define COUNTRY 694 "SL", add
label define COUNTRY 716 "ZW", add
label define COUNTRY 740 "SR", add
label define COUNTRY 764 "TH", add
label define COUNTRY 768 "TG", add
label define COUNTRY 776 "TO", add
label define COUNTRY 788 "TN", add
label define COUNTRY 795 "TM", add
label define COUNTRY 796 "TC", add
label define COUNTRY 798 "TV", add
label define COUNTRY 807 "MK", add
label define COUNTRY 882 "WS", add
label define COUNTRY 990 "XK", add

************************************************************
**********************   INDICATORS   *********************
************************************************************

//creating a 0/1 indicator for MICS round 6 
gen round_6 = round
label variable round_6 "mics6"
recode round_6 (2/5=0) (6=1)
label define round_6 0 "not mics6" 1 "mics6"
label values round_6 round_6
tab round_6 round //verify this transformation

//creating the disability indicator of focus in the present study 
/* From this construction, the indicator "disabled" includes
 respondents with 1+ of "some functional difficulty in 1+ of the 
six domains measured by the survey: sight, 
hearing, walking, communicating, remembering/concentrating, or self-care.
*/

*STUDENTS - SEE BLOG POST ON  LOCAL MACROS IN STATA DO FILES
local diffvaris ///
"diffsee diffhear diffwalk diffremcon diffcare diffcom"
egen diffsum = rowtotal(`diffvaris') 
label variable diffsum "total disability score" 
recode diffsum (0=0) (.=.) (else=1) , gen(disabled)
label variable disabled "disability indicator"
label define disabled 0 "not disabled" 1 "disabled"
label values disabled disabled
tab diffsum disabled //verify the transformation 

//creating a discrimination indicator
/* From this construction, the indicator "discriminated" includes respondents 
who responded "yes" to 1+ of the 7 surveyed types of discrimination:
gender, sexual orientation, age, religion, disability, ethnicity/immigration, 
or other reason.
*/

*STUDENTS - SEE BLOG POST ON LOCAL MACROS IN STATA DO FILES
local discrimvars discrimethnic discrimgender discrimsexorien ///
discrimrelig discrimage discrimdis discrimother
egen discrimsum = rowtotal(`discrimvars') 
label variable discrimsum "total discrimination score" 
recode discrimsum (0=0) (.=.) (else=1) , gen(discriminated)
label variable discriminated "discriminated indicator"
label define discriminated 0 "not discriminated" 1 "discriminated"
label values discriminated discriminated
tab discrimsum discriminated //verify the transformation 

/*The diff* vars' distributions are very heavy with zeroes, meaning most 
report "no difficulty" in most functional domains. Thus the other levels 
require combining because there are relatively few responses in higher 
levels. Statistical detail is not lost by this, and power is gained.
*/

*Creating binary indicators for each diif* var via loop
*STUDENTS - SEE BLOG POST ON LOCAL MACROS IN STATA DO FILES
local diffvaris "diffcare diffcom diffremcon diffwalk diffhear diffsee"

foreach var of local diffvaris {
    generate `var'yn = `var'
    label variable `var'yn "difficulty `var' Y/N"
    recode `var'yn (1/3=1)
    label define `var'yn 0 "no difficulty" 1 "yes difficulty"
    label values `var'yn `var'yn
    tab `var'yn `var'
}

************************************************************
********************   WEIGHTED MEANS   ****************
************************************************************

//Note: dif* functioning vars and discrim vars are only in MICS 6

*mean women's Cantril ladder score by country in Round 6 = lsladderM
egen lsladderM = wtmean(lsladder) if round_6==1, ///
by( country ) weight( weightwm )
label variable lsladderM "mean women's national Cantril ladder"
label values lsladderM lsladderM

*mean overall women's happiness score by country in Round 6 = lshappyM
egen lshappyM = wtmean(lshappy) if round_6==1, ///
by( country ) weight( weightwm )
label variable lshappyM "mean women's national overall happiness"
label values lshappyM lshappyM

*mean overall women's disability total by country in Round 6 = disabilityM
egen disabilityM = wtmean(diffsum) if round_6==1, ///
by( country ) weight( weightwm )
label variable disabilityM "mean women's national total disabilities"
label values disabilityM disabilityM

*national proportion of disabled women in Round 6 = disabledM
egen disabledM = wtmean(disabled) if round_6==1, ///
by( country ) weight( weightwm )
label variable disabledM "mean national disabled women proportion"
label values disabledM disabledM

*national proportion of women exper. discrim. in Round 6 = discriminatedM
egen discriminatedM = wtmean(discriminated) if round_6==1, ///
by( country ) weight( weightwm )
label variable discriminatedM "mean natl women's discrimination proportion"
label values discriminatedM discriminatedM

*mean women's LS versus last year 
egen lslastyrM = wtmean(lslastyr) if round_6==1, ///
by( country ) weight( weightwm )
label variable lslastyrM "mean women's national ls vs last yr"
label values lslastyrM lslastyrM

*mean women's LS expected one year from now 
egen lsnextyrM = wtmean(lsnextyr) if round_6==1, ///
by( country ) weight( weightwm )
label variable lsnextyrM "mean women's national ls expected next yr"
label values lsnextyrM lsnextyrM

*Note: The rest of the ls* vars are not in mics6 

************************************************************
********************   STYLING GRAPHS   ******************
************************************************************

grstyle clear //good practice, resets any prev grstyle in mem
grstyle init
grstyle set horizontal //sets y axis tick labels horizon/readable
grstyle anglestyle vertical_tick horizontal 
grstyle set plain, nogrid
grstyle set legend 4, inside //clock position - choose best layout option
grstyle set graphsize 10in 16in //h x w
grstyle set margin ".25 .25 .5 .5", in: graph  // L R B T graph region
grstyle set margin .5in: twoway  // all four equal plot region
grstyle set symbolsize small
grstyle set size 32pt: heading
grstyle set size 30pt: subheading axis_title
grstyle color background white
grstyle color plotregion none
grstyle symbolsize p small
grstyle gsize axis_title_gap small //adds space bet ticks&axis titles
grstyle linewidth major_grid vthin

************************************************************
******************   DESCRIPTIVE GRAPHS   ***************
************************************************************

/*Countries are labeled in these graph by their alpha-2 codes 
per ISO 3166 international standards at the time of this pub.
*/

//this a manual process, tedious but necessary, to avoid overlapping
gen clockpos=12
replace clockpos=1 if country==796 | country==288  | ///
country==586 | country==426 //TC, GH, PK, LS
replace clockpos=2 if country==524 | country==678 //NP,ST
replace clockpos=3 if country==499 | country==716 //ME, ZW  
replace clockpos=4 if country==12 //DZ 
replace clockpos=7 if country==454 | country==798 //MW, TV 
replace clockpos=9 if country==180 //CD
replace clockpos=10 if country==140 | country==496 | ///
country==882 | country==624 //CF, MN, WS, GW
replace clockpos=11 if country==768 | country==188 | ///
country==368 | country==807 //TG, CR, IQ, MK

*ladderM by disabledM, twoway graph for Women by Country
*STUDENTS - SEE BLOG POST ON LOCAL MACROS IN STATA DO FILES
local fig1common ///
"mlabel(country) mlabgap(1.0) msize(small) mlabvpos(clockpos)"

mylabels 4(2)10, myscale(@) local(ylab2)
mylabels 0(10)50, myscale(@/100) local(ylab1) suffix(" %")
*help mylabels 

twoway (scatter lsladderM disabledM if round==6 & country!=188, ///
`fig1common' ms(o) mlabcolor(navy) mcolor(navy) mlabsize(medsmall)) ///
(scatter lsladderM disabledM if round==6 & country==188, ///
`fig1common' ms(D) mlabcolor(pink) mlabgap(1.5) mlabsize(large) ///
mcolor(pink) msize(large) ), ///
ylabel(`ylab2') xlabel(`ylab1') legend(off) ///
ytitle("Mean Life Satisfaction" "(Cantril Ladder 0-10 Scale)") ///
xtitle("Mean Disability Prevalence" ///
"(Indicator of 1+ Functional Impairments)") ///
title("Life Satisfaction by Disability Prevalance for Women by Country" ///
, span linegap(1.2) ) ///
subtitle("IPUMS MICS Data 2017-2019, Round 6" ///
"(n = 513,744, in 29 national samples)", span) ///
|| lfit lsladderM disabledM if round==6 & country!=188, ///
lwidth(thin) lcolor(navy) lpattern("-#-#.#-#..#") ///
name(paafigure1, replace)
graph export "C:\Sync\Data\paafigure1.jpg", ///
as(jpg) name("paafigure1") width(4000) quality(100) replace
*The svg will be huge. Uncomment to run if desired. Edit file path.
//graph export "C:\Sync\Data\paafigure1.svg", as(svg) replace


************************************************************
******************   DESCRIPTIVE TABLES   ****************
************************************************************

* sorting weighted mean for women's disability prevalance, ranking
preserve
collapse (mean) disabledM = disabled [pw=weightwm], ///
by(round_6 country)
gsort round_6 -disabledM
bysort round_6: list country disabledM if round_6==1
restore 
/*disabledM per country: 
37 countries measured in MICS round 6 
top five highest proportions of disability: 
#1: CR (.48898) ~ 49%
#2: CF (.46132) ~ 46%
#3: GE (.43178) ~ 43%
#4: ST (.43107) ~ 43%
#5: HN (.42639) ~ 43%
*/

* sorting weighted mean for women's LS, ranking 
preserve
collapse (mean) lsladderM = lsladder [pw=weightwm], ///
by(round_6 country)
gsort round_6 -lsladderM
bysort round_6: list country lsladderM if round_6==1
restore 
/* lsladderM per country:
30 countries measured, all mics round 6 only
top five highest ladder means by country:
#1: TO 8.637  Tonga 
#2: ME 8.342  Montenegro 
#3: WS 8.300  Samoa 
#4: CR: 8.091 Costa Rica 
#5: RS 7.815  Serbia 
*/

program savecompressed
    version 11.2
    syntax [ , REPLACE ]
    if (c(changed) & ("`replace'"!="replace")) error 4
    local width = c(width)
    compress
    if (c(width)<`width') save , replace
end

*To save cleaned & recoded data:
label data "IPUMS MICS All Samples Rounds PAA Cleaned"
saveold "C:\Sync\Data\ipums-mics-clean", replace //edit to yours


************************************************************
*******************   COSTA RICA ONLY   ******************
************************************************************

clear all
use "C:\Sync\Data\ipums-mics-clean.dta" , clear 

//country==188 for Costa Rica 
keep if country==188 
recode country (188=1) (else=0)
label define country 0 "not costa Rica" 1 "Costa Rica"
label values country country
tab country 

//this drops all variables with no data values in Costa Rica 

foreach var of varlist _all {
     capture assert mi(`var')
     if !_rc {
        drop `var'
     }
 }

*To save cleaned & recoded data for CR only:
label data "IPUMS MICS CR only"
saveold "C:\Sync\Data\ipums-mics-CR-clean", replace  //edit to yours

clear all
use "C:\Sync\Data\ipums-mics-CR-clean" , clear //edit to yours

/*Separately, I did unicode encoding for value labels display:
cd c:/sync/data/       //change directory to current dir 
unicode analyze ipums-mics-CR-clean.dta   //shows labels w/issues
unicode encoding set "Latin1"    //sets default enc to iso-Latin1 
unicode translate "ipums-mics-CR-clean.dta", transutf8 //translates
*/

************************************************************
**********************   CLEANING CR  ********************
************************************************************

* recoding missing geo1_cr
recode geo1_cr (188998=.) //missing fm incomplete questionnaire
* Create a new var "region" and assign it the values of "geo1_cr"
generate region = geo1_cr 
label variable region "Costa Rica Province"
label define region 188001 "San José" 188002 "Alajuela" ///
                    188003 "Cartago" 188004 "Heredia" ///
                    188005 "Guanacaste" 188006 "Puntarenas" ///
                    188007 "Limón"
label values region region 
tab region geo1_cr //verify this transformation 

* recoding missing ethnic_cr
recode ethnic_cr (98=.) //missing
generate ethnicity = ethnic_cr
label variable ethnicity "ethnicity of HH" //head of household
label define ethnicity 1 "Indigenous" ///
				   2 "Black/Afro-Costa Rican" ///
				  90 "Other Group Not-Spec."
label values ethnicity ethnicity
tab ethnicity ethnic_cr //verify this transformation 

*because 90% are in "Other" category, must combine v.small groups 
generate ethnicitybi = ethnicity 
label variable ethnicitybi "Black/Indigenous Indicator"
recode ethnicitybi (1/2=1) (90=0)
label define ethnicitybi 0 "Other Group Not-Spec." ///
1 "Black/Afro-CR or Indigenous"
label values ethnicitybi ethnicitybi 
tab ethnicitybi ethnicity //verify the transformation

gen married = marst 
lab var married "married indicator"
recode married (10=1) (20/30=0)
lab def married 0 "not married" 1 "married"
lab val married married 
tab married marst //verify the transformation 

************************************************************
**********************  WEIGHTING CR  *******************
************************************************************

svyset [pweight=weightwm]

************************************************************
********************   DESCRIPTIVES CR  ******************
************************************************************

*mean women's Cantril ladder score by province = lsladderCRR
egen lsladderCRR = wtmean(lsladder) if round_6==1, ///
by( region ) weight( weightwm )
label variable lsladderCRR "mean regional Cantril ladder"
label values lsladderCRR lsladderCRR

*mean women's overall happiness score by province - lshappyCRR
egen lshappyCRR = wtmean(lshappy) if round_6==1, ///
by( region ) weight( weightwm )
label variable lshappyCRR "mean regional overall happiness"
label values lshappyCRR lshappyCRR

*mean women's overall disability total by province = disabilityCRR
egen disabilityCRR = wtmean(diffsum) if round_6==1, ///
by( region ) weight( weightwm )
label variable disabilityCRR "mean regional total disabilities"
label values disabilityCRR disabilityCRR

*regional proportion of disabled women = disabledCRR
egen disabledCRR = wtmean(disabled) if round_6==1, ///
by( region ) weight( weightwm )
label variable disabledCRR "mean regional disabled proportion"
label values disabledCRR disabledCRR

*regional proportion of women exper. discrim = discriminatedCRR
egen discriminatedCRR = wtmean(discriminated) if round_6==1, ///
by( region ) weight( weightwm )
label variable discriminatedCRR "mean regional discrimination proportion"
label values discriminatedCRR discriminatedCRR

************************************************************
*******************   STYLING GRAPHS CR  ***************
************************************************************

grstyle clear //resets any previous grstyle
grstyle init
grstyle set horizontal //sets y axis tick labels horizontal/readable yay!
grstyle anglestyle vertical_tick horizontal 
grstyle set plain, nogrid
grstyle set legend 4, inside //clock position - choose best layout option
grstyle set graphsize 10in 12in //h x w
grstyle set margin ".25 .25 .5 .5", in: graph  // L R B T graph region
grstyle set margin .5in: twoway  // all four equal for plot region
grstyle set symbolsize medium
grstyle set margin medium: heading
grstyle set margin medium: subheading
grstyle set size 32pt: heading
grstyle set size 24pt: subheading axis_title
grstyle color background white
grstyle color plotregion none
grstyle symbolsize p small
grstyle gsize axis_title_gap small //adds small space bet ticks, axis titles
grstyle linewidth major_grid vthin

************************************************************
*****************   DESCRIPTIVE GRAPHS CR  *************
************************************************************

*Figure 2

/*SE would be wrong here using pweight as aweight but the weighted
 mean point estimates are correct, and all we are plotting. This graph is
 not formatted but depicts life satisfaction gaps by province */
graph dot (mean) lsladder [aweight = weightwm], ///
over(disabled) by(region) ///
blabel(bar, format(%4.2f)) exclude0 yscale(range(7(1)9)) ///
marker(1, ms(O) msize(*1.2) mcolor(navy) mlabpos(12)) ///
name(paafigure2, replace)
graph export "C:\Sync\Data\paafigure2.jpg", ///
as(jpg) width(4000) quality(100) replace

*Figure 3

gen clockposn=12
replace clockposn=1 if region==188005 | region==188007 
replace clockposn=2 if region==188001 | region==188003 
replace clockposn=3 if region==188002 

mylabels 40(5)60, myscale(@/100) local(ylab1) suffix(" %") 
mylabels 7.9(.1)8.4, myscale(@) local(ylab3) format(%2.1f)

* lsladderCRR by disabledCRR  twoway graph women by province  
twoway (scatter lsladderCRR disabledCRR if round==6 & ///
region!=188007, ///
mlabel(region) ms(o) mlabcolor(navy) mcolor(navy) mlabgap(1.0) ///
mlabsize(medium) msize(medsmall) mlabvpos(clockposn)) ///
(scatter lsladderCRR disabledCRR if round==6 & region==188007, ///
mlabel(region) ms(D) mlabcolor(pink) mcolor(pink) mlabgap(1.5) ///
mlabsize(huge) msize(large) mlabvpos(clockposn)), ///
ylabel(`ylab3') xlabel(`ylab1') legend(off) ///
ytitle("Mean Life Satisfaction" "(Cantril Ladder 0-10 Scale)")  ///
xtitle("Mean Disability Prevalence" ///
"(Indicator of 1+ Functional Impairments)") ///
title("Life Satisfaction and Disability for Women by Province" ///
, span linegap(1.2) ) ///
subtitle("IPUMS MICS Data 2018, Round 6 (n = 8,193)", span) ///
|| lfit lsladderCRR disabledCRR if round==6, lwidth(thin) ///
lcolor(gray) lpattern("-#-#.#-#..#") ///
name(paafigure3, replace)
graph export "C:\Sync\Data\paafigure3.jpg", ///
as(jpg) width(4000) quality(100) replace

************************************************************
**********************   ANALYSIS CR  *********************
************************************************************

* chi square test for independence
tabulate region disabled, chi2
//these two categorical variables (region and disabled) are related 

*one-way ANOVA 
anova lsladder region 
//sig differences between LS means of samples among levels of region  

drop if round_6==0 //done with the past

svy:ologit lsbin i.wealth i.disabled##i.region i.healthinsur ///
i.discriminated i.married, coeflegend // ,coeflegend to get names
eststo lsladder1

*Using gologit2 to fit partial proportional odds model for Wald tests
gsvy: gologit2 lsbin i.wealth i.disabled##i.region i.healthinsur ///
i.discriminated i.married, nolabel auto

/*An Adjusted Wald test on the final probability-weighted model 
was run using the user-writted Stata program gologit2. This was 
an insignificant Wald test which suggests that the final model does 
not violate the proportional odds, or parallel lines assumption of 
ordered logistic regression models, F(60, 8128) = 1.19, p = .148, n.s.
*/

*regression table output
esttab lsladder1, b(%5.3f) ci label nonumbers drop (cut*) ///
title("Life Satisfaction Ordered Logit Model Coefficients") ///
mtitles("lsladder1") nobaselevels

*coefplot - a graphical depiction of a regression table, which I prefer

// Defining coefficient labels
local coeflabels ///
1.disabled#188002.region = "Disabled in Alajuela" ///
1.disabled#188003.region = "Disabled in Cartago" ///
1.disabled#188004.region = "Disabled in Heredia" ///
1.disabled#188005.region = "Disabled in Guanacaste" ///
1.disabled#188006.region = "Disabled in Puntarenas" ///
1.disabled#188007.region = "Disabled in Limón"

//Defining headings 
local headings 1.disabled = "{bf:Disability} {it:(base = no)}" ///
1.disabled#188002.region = "{bf:Disability x Province}" ///
188002.region = "{bf:Provinces} {it:(base = San José)}" ///
1.wealth = "{bf:Wealth} {it:(base = lower)}"

// Defining variables to not graph here
local dropvars "_cons 1.healthinsur 1.discriminated 1.married"

// Defining titles
local title "Life Satisfaction Predictors"
local subtitle ///
`"Controls: marital status; discrimination; health insurance"'
local caption ///
`"{it:Note:} * p < .05; ** p < .01; *** p < .001; Caps are 95% CIs."'

coefplot ///
(lsladder1, label(1) mcolor(navy) mlabcolor(navy) msymbol(O) ///
ciopts(color(navy) recast(rcap))), ///
drop(_cons | 1.healthinsur | 1.discriminated | 1.married ) ///
xlabel(, notick) ///
grid(between glpattern(dash) glwidth(*0.5) glcolor(gs8)) ///
ylabel(,labsize(small) notick) ///
xline(0, lwidth(thin) lcolor(gray) lpattern(solid)) ///
order (1.disabled *.region *.disabled#*.region) ///
headings(`headings', ///
labcolor(navy) labsize(medsmall)) /// 
msize(medium) ///
mlabel(cond(@pval<.001, "***", ///
cond(@pval<.01, "**", ///
cond(@pval<.05, "*", ///
cond(@pval>.05, "", ///
string(@pval,"%9.3f")))))) ///
mlabposition(3) mlabsize(medium) /// 
coeflabels(`coeflabels') ///
title(`title', span) ///
subtitle(`subtitle', span ) ///
caption(`caption', span size(*.8)) ///
name(paafigure4, replace)

// Saving the graph to rasterized and vector file formats 
graph export "C:\Sync\Data\paafigure4.jpg", ///
as(jpg) name("paafigure4") width(4000) quality(100) replace
graph export "C:\Sync\Data\paafigure4.svg", as(svg) replace


************************************************************
********************  CROSS VALIDATION  ********************
************************************************************
*This was done earlier, but dropped down in this do-file 

*Gallup World Poll is used here for validating life satisfaction 

*OPENING DATA
clear all 
use "C:\Sync\Data\GWP_021723.dta" , clear

*WEIGHTING
svyset [pweight=wgt]

*CLEANING ONLY WHAT IS NEEDED
*some nondestuctive naming
generate income = INCOME_2
generate lifetoday = WP16

//lifetoday is the Cantril ladder item in Gallup World Poll
//98=DK, 99=refused
//add values labels so I know what my new cats mean
recode lifetoday (98=.) (99=.)
label define lifetoday ///
0 "worst possible" 10 "best possible"
label values lifetoday lifetoday 
tab1 WP16 lifetoday //verify that this is recoded right

*healthprob
//orig coding: 1=yes 2=no 3=DK 4=ref
generate healthprob = WP23
label variable healthprob "major health problems"
recode healthprob (2/4=0) (1=1)
label define healthprob 0 "no health prob" 1 "yes health prob"
label values healthprob healthprob
tab1 WP23 healthprob //verify this transformation

/*Note that for Gallup, this variable indicates having "major 
health problems" which differs from the present study, in which 
the disability indicator is for anyone with "some difficulty" in 
1+ functional domains. I am not using GWP to cross-validate 
disability data for this reason. */

*disability
gen disability = healthprob
label variable disability "disability indicator"
label define disability 0 "no disability" 1 "yes disability"
label values disability disability
tab1 disability healthprob //verify this transformation

*year
generate year = YEAR_WAVE if !missing(YEAR_WAVE)
label variable year "year"
label values year year
tab1 YEAR_WAVE year //verify this transformation

//country
generate country = WP5 if !missing(WP5)
label variable country "country"
label values country WP5
label define country 82 "CR", replace  // Costa Rica


/*********************************************************
                    DESCRIPTIVE MEASURES
**********************************************************/

/* Must use the respondent-level weight variable (wgt) if 
aggregating for national averages*/

*meanlifetoday //Cantril ladder today question (WP16)
egen meanlifetoday = wtmean(lifetoday), ///
by (country year) weight(wgt)
label variable meanlifetoday "annual mean Cantril ladder"
label values meanlifetoday meanlifetoday 

*meandis //binary disability indicator as above
egen meandis = wtmean(disability), ///
by (country year) weight(wgt)
label variable meandis "annual mean disability proportion"
label values meandis meandis 

*meancrlifetoday - regional annual mean in Costa Rica 
egen meancrlifetoday = wtmean(lifetoday), ///
by (REGION_CR year) weight(wgt)
label variable meancrlifetoday "mean Cantril ladder by region"
label values meancrlifetoday meancrlifetoday 

/* Overall, Gallup distributions affirm the present study's 
life satisfaction distribution, as aggregated by country and, 
as expected, Gallup's disability prevalance is lower than the 
present study, as it is using a different severity scale of functional disability */

exit