/*Data quality measures preparation*/
clear all
quietly save "N:\quality_data.dta", replace emptyok
foreach waveprefix  in a b c d e f g h i j k l m n{	
	use "N:\HILDA - Wave 14 - Unconfidentialised\STATA 140u\Rperson_`waveprefix'150u.dta", clear
	renpfix `waveprefix'
	gen wave = index("abcdefghijklmnopqrstuvwxyz","`waveprefix'")
	gen no_SCQ = gh1 == -8 if gh1!=.
	gen follow_up_interview = hhiwin==2 if hhiwin!=.
	foreach var of varlist _all{
		local label : variable label `var'
		if strmatch("`label'", "*DV:*") | strmatch("`label'", "*SCQ*") | strmatch("`label'", "HF*") | strmatch("`label'", "PQ*") capture drop `var'
	}
	capture drop io*
	capture drop hhpmid
	gen Q_relevant = 0
	gen Q_bad_answer = 0
	foreach var of varlist _all{
		capture replace Q_relevant = Q_relevant + 1 if `var'>=0 & `var'!=.
		capture replace Q_bad_answer = Q_bad_answer + 1 if (`var'==-3|`var'==-4|`var'==-6|`var'==-7) // Don’t know / Refused or not answered / Value implausible / Unable to determine value
	}
	gen perc_prob_answ = (Q_bad_answer / Q_relevant) * 100
	keep xwaveid wave perc_prob_answ no_SCQ follow_up_interview Q_bad_answer Q_relevant
	append using "N:\HILDA - Wave 14 - Unconfidentialised\STATA 140u\quality_data.dta"
	quietly save "N:\HILDA - Wave 14 - Unconfidentialised\STATA 140u\quality_data.dta", replace
	disp in green "Wave `waveprefix' done!"
}
table wave, c(mean perc_prob_answ mean no_SCQ mean follow_up_interview ) f(%4.3f)
clear all

/*Full HILDA data preparation*/
ssc install hildasetup
hildasetup, saving("N:\") data("N:\HILDA - Wave 14 - Unconfidentialised\STATA 140u") ///
var(hgivw ghsf6d ghpf ghrp ghbp ghgh ghvt ghsf ghre ghmh hemirh pdk10s henec heslu ghmcs hhpxid anatsi anbcob hgage hgsex mrcurr edhigh esbrd hhsos hhra helth hifdip hifdin hhstate hhadult hh0_4 hh5_9 hh10_14 hhad10 hhpqlen hhpqmo hhintpq ahhura io*) filename(HILDA_MH_Survey_quality) rel(14) unconf

/*Preparations (I)*/
clear all
set more off
use "N:\HILDA_MH_Survey_quality.dta", clear
merge 1:1 xwaveid wave using "N:\HILDA - Wave 14 - Unconfidentialised\STATA 140u\quality_data.dta"
drop _merge
destring xwaveid, replace

/*Attrition*/
gen not_interviewed = hgivw>=5
gen attrition = .
bysort xwaveid (wave): replace attrition = 0 if not_interviewed==0 & not_interviewed[_n+1]==0
bysort xwaveid (wave): replace attrition = 1 if not_interviewed==0 & not_interviewed[_n+1]!=0
replace attrition = . if wave==14

/*Preparations (II)*/
drop if edhigh==-10
recode anatsi -10/1=0 2/4=1, gen (atsi)
recode hemirh -1=0
recode henec -1=0
recode heslu -1=0
mvdecode _all, mv(-10/-1)

/*Mental health*/  
rename ghmh SF36_MH
rename pdk10s K10
rename hemirh mental_illness
rename henec nerv_emot_cond
rename heslu diff_learn_underst

/*Other outcome variables*/
recode iosusp 1=0 2/3=1, gen(suspicious)
rename ioundst understanding
rename iocoop cooperation
recode understanding 1/2=0 3/5=1
recode cooperation 1/2=0 3/5=1
recode perc_prob_answ 100/1000=. 
gen refusals = perc_prob_answ>2 if perc_prob_answ!=.
rename no_SCQ scq_no

/*Control variables*/
replace anbcob=4 if atsi==1
tab anbcob, gen (xxx)
rename (xxx1 xxx2 xxx3 xxx4) (AB ESB NESB ATSI)
gen reading_language_probs = (iopread+iopeng+ioplang)
rename hhadult n_adults_in_HH
gen n_children_in_HH = hh0_4+hh5_9+hh10_14
rename hhstate state
tab state, gen(x)
rename x1 NSW
rename x2 VIC
rename x3 QLD
rename x4 SA
rename x5 WA
rename x6 TAS
rename x7 NT
rename x8 ACT
replace hifdip = hifdip/0.707 if wave==1
replace hifdip = hifdip/0.728 if wave==2
replace hifdip = hifdip/0.746 if wave==3
replace hifdip = hifdip/0.765 if wave==4
replace hifdip = hifdip/0.786 if wave==5
replace hifdip = hifdip/0.812 if wave==6
replace hifdip = hifdip/0.836 if wave==7
replace hifdip = hifdip/0.867 if wave==8
replace hifdip = hifdip/0.885 if wave==9
replace hifdip = hifdip/0.909 if wave==10
replace hifdip = hifdip/0.936 if wave==11
replace hifdip = hifdip/0.957 if wave==12
replace hifdip = hifdip/0.983 if wave==13
replace hifdip = hifdip/1 if wave==14
replace hifdip = hifdin if hifdin!=0 & hifdin!=.
replace hifdip = 1 if hifdip>=0 & hifdip<=1
rename hifdip income
gen income_in_10000s = income/10000
rename helth condition
recode condition 2=0
gen female=hgsex==2 & hgsex!=.
rename hgage age
gen age_squared = age^2
gen partnered=.
replace partnered =1 if mrcurr==1|mrcurr==2
replace partnered =0 if mrcurr>2 & mrcurr!=.
recode edhigh 1/3=1 4/7=2 8=3 9=4 10=., gen(education)
tab education, gen(xxx)
rename xxx1 degree
rename xxx2 VET
rename xxx3 year_12
rename xxx4 under_year_12
gen employed = esbrd == 1 if esbrd!=.
gen ue = esbrd == 2 if esbrd!=.
gen nilf = esbrd == 3 if esbrd!=. 
rename esbrd employment_status
gen majorurban = hhsos==0
gen MC = hhra==0
gen IR = hhra==1
gen OR_R_VR = hhra>1 & hhra<9
recode hhad10 1/2=1 3/4=2 5/6=3 7/8=4 9/10=5, gen(deprivation)
lab def deprivation 1 "1. Lowest" 2 "2" 3 "3" 4 "4" 5 "5. Highest"
lab val deprivation deprivation 
tab deprivation, gen(deprivation)
bysort xwaveid (wave): gen times_interviewed = _n
gen year = wave+2000
rename hhintpq interviewer_ID
bysort interviewer_ID wave: gen interv_workload = _N
bysort interviewer_ID xwaveid (wave): gen first_matching = _n==1

/*Housekeeping*/
ssc install runmlwin
capture rename (influenced suspicious understanding cooperation) (inf susp und coop)
capture rename (mental_illness nerv_emot diff_learn) (mental nerv learn)
global mental_health "SF36_MH K10 mental learn nerv"
global controls 	"female age age_squared partnered n_adults_in_HH n_children_in_HH ESB NESB ATSI VET year_12 under_year_12 nilf ue income_in_10000s IR OR_R_VR deprivation2 deprivation3 deprivation4 deprivation5 VIC QLD SA WA TAS NT ACT times_interviewed first_matching interv_workload year"
global outcomes 	"susp und coop"
global outcomes2	"attr ref scq"
global MLwiN_path "C:\Program Files (x86)\MLwiN v2.32\i386\MLwiN.exe"
keep xwaveid wave interviewer_ID $mental_health $controls $outcomes $outcomes2 AB degree employed MC deprivation1 NSW
quietly compress

/*Analysis: Simple logistic regression models and predicted probabilities*/
quietly{
	foreach y in susp und coop attr ref{
		foreach mh in SF36_MH K10{
			eststo b_`y'_`mh': quietly logit `y' `mh'
			foreach m in 10 25 50 75 90{
				quietly margins, at((p`m') `mh')
				quietly matrix table = r(table)
				local `y'_`mh'_`m' = table[1,1]
			}
			matrix input `y'_`mh' = (``y'_`mh'_10',``y'_`mh'_25',``y'_`mh'_50',``y'_`mh'_75',``y'_`mh'_90')
		}
	}
	foreach y in susp und coop attr scq ref{
		foreach mh in mental nerv learn{
			eststo b_`y'_`mh': quietly logit `y' `mh'
			foreach m in 0 1{
				quietly margins, at(`mh'=`m')
				quietly matrix table = r(table)
				local `y'_`mh'_`m' = table[1,1]
			}
			matrix input `y'_`mh' = (``y'_`mh'_0',``y'_`mh'_1')
		}
	}
	matrix final = und_SF36_MH\und_K10\coop_SF36_MH\coop_K10\susp_SF36_MH\susp_K10\attr_SF36_MH\attr_K10\ref_SF36_MH\ref_K10
	matrix colnames final = "p10" "p25" "p50" "p75" "p90"
	matrix rownames final = und_MHI und_K10 coop_MH coop_K10 susp_MH susp_K10 attr_SF36_MH attr_K10 ref_SF36_MH ref_K10
	matrix final2 = und_mental\und_learn\und_nerv\coop_mental\coop_learn\coop_nerv\susp_mental\susp_learn\susp_nerv\attr_mental\attr_learn\attr_nerv\scq_mental\scq_learn\scq_nerv\ref_mental\ref_learn\ref_nerv
	matrix colnames final2 = "p=0" "p=1"
	matrix rownames final2 = und_mental und_learn und_nerv coop_mental coop_learn coop_nerv susp_mental susp_learn susp_nerv attr_mental attr_learn attr_nerv scq_mental scq_learn scq_nerv ref_mental ref_learn ref_nerv
}
matrix list final, ti(Predicted probabilities) noh f(%4.3f)
matrix list final2, ti(Predicted probabilities) noh f(%4.3f)
*** user-written command esttab must be installed
esttab b_susp_SF36_MH b_susp_K10 b_susp_mental b_susp_learn b_susp_nerv using "N:\Simple logits.rtf", replace ti(susp) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not
esttab b_und_SF36_MH b_und_K10 b_und_mental b_und_learn b_und_nerv using "N:\Simple logits.rtf", append ti(und) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not
esttab b_coop_SF36_MH b_coop_K10 b_coop_mental b_coop_learn b_coop_nerv using "N:\Simple logits.rtf", append ti(coop) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not
esttab b_attr_SF36_MH b_attr_K10 b_attr_mental b_attr_learn b_attr_nerv using "N:\Simple logits.rtf", append ti(attr) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not
esttab b_ref_SF36_MH b_ref_K10 b_ref_mental b_ref_learn b_ref_nerv using "N:\Simple logits.rtf", append ti(ref) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not
esttab b_scq_mental b_scq_learn b_scq_nerv using "N:\Simple logits.rtf", append ti(scq) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not

/*Analysis: 3-level cross-classified random-effect logistic regression models*/
capture gen cons = 1
foreach key_X in $mental_health{
	preserve
	foreach i in `key_X' $controls $outcomes{
		display "`i'"
		drop if `i'==.
	}
	capture gen obs_id = _n
	sort interviewer_ID xwaveid obs_id
	recast float SF36_MH, force
	***
	runmlwin susp cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store in_sus_`key_X'
	runmlwin und cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store in_und_`key_X'
	runmlwin coop cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store in_coo_`key_X'
	runmlwin attr cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store in_attr_`key_X'
	runmlwin ref cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store in_ref_`key_X'
	if "`key_X'" != "K10" & "`key_X'" != "SF36_MH"{
		runmlwin scq cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) nopause or discrete(distribution(binomial) link(logit) denom(cons))
		est store in_scq_`key_X'
	}
	***
	runmlwin susp cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_sus_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store fu_sus_`key_X'
	runmlwin und cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_und_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store fu_und_`key_X'
	runmlwin coop cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_coo_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store fu_coo_`key_X'
	runmlwin attr cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_attr_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store fu_attr_`key_X'
	runmlwin ref cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_ref_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
	est store fu_ref_`key_X'
	if "`key_X'" != "K10" & "`key_X'" != "SF36_MH"{
		runmlwin scq cons `key_X' $controls, level3(interviewer_ID: cons) level2(xwaveid: cons) level1(obs_id: ) mcmc(cc) initsmodel(in_scq_`key_X') nopause or discrete(distribution(binomial) link(logit) denom(cons))
		est store fu_scq_`key_X'
	}
	***
	capture log close
	log using "N:\Regression results - Mental health_`key_X'.log", replace
	est replay fu*`key_X'
	capture log close
	restore
}

esttab fu_sus_SF36_MH fu_sus_K10 fu_sus_mental fu_sus_learn fu_sus_nerv using "N:\Full logits.rtf", replace ti(susp) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not order(SF36_MH K10 mental learn nerv) mti(SF36 K10 men learn ner)
esttab fu_und_SF36_MH fu_und_K10 fu_und_mental fu_und_learn fu_und_nerv using "N:\Full logits.rtf", append ti(und) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not order(SF36_MH K10 mental learn nerv) mti(SF36 K10 men learn ner)
esttab fu_coo_SF36_MH fu_coo_K10 fu_coo_mental fu_coo_learn fu_coo_nerv using "N:\Full logits.rtf", append ti(coop) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not order(SF36_MH K10 mental learn nerv) mti(SF36 K10 men learn ner)
esttab fu_attr_SF36_MH fu_attr_K10 fu_attr_mental fu_attr_learn fu_attr_nerv using "N:\Full logits.rtf", append ti(attr) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not order(SF36_MH K10 mental learn nerv) mti(SF36 K10 men learn ner)
esttab fu_ref_SF36_MH fu_ref_K10 fu_ref_mental fu_ref_learn fu_ref_nerv using "N:\Full logits.rtf", append ti(ref) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not order(SF36_MH K10 mental learn nerv) mti(SF36 K10 men learn ner)
esttab fu_scq_mental fu_scq_learn fu_scq_nerv using "N:\Full logits.rtf", append ti(scq) stats(N) star((*) 0.1 * 0.05 ** 0.01 *** 0.001) b(%12.3f) nogaps nonum eform not

/*Analysis: Descriptive statistics*/
tabstat $mental_health $outcomes $outcomes2 female age partnered n_adults_in_HH n_children_in_HH AB ESB NESB ATSI degree VET year_12 under_year_12 employed nilf ue income_in_10000s MC IR OR_R_VR deprivation1 deprivation2 deprivation3 deprivation4 deprivation5 NSW VIC QLD SA WA TAS NT ACT times_interviewed first_matching interv_workload year, f(%10.2f) st(n mean p50 sd min max) c(st)