* Program to produce the sex ratio at birth for districts in the India NFHS-5 survey.
* Can be adapted for subpopulations in any DHS survey.

* Two definitions of the sex ratio are used.

* The time interval is the past five years, all births in the KR file, except those
*  in the month of interview.

* The program provides the lower and upper ends of 95% confidence intervals.  
* The intervals are wide.
* The estimates use svy, including subpop (within the state).
* Bayesian procedures would reduce the confidence intervals and move the estimates
*  toward the state value

* The results are saved into the workfile and then the workfile is reduced
*  to just the saved results.

* specify a workspace
cd e:\DHS\DHS_data\scratch

use "...IAKR7DFL.DTA" 
keep v001 v002 v003 v005 v008 v023 v024 v025 sdist b4

* Construct binary variable m_per_f
gen m_per_f=0
replace m_per_f=1 if b4==1

* Sex ratio defined as males per 100 females; do not construct until end
* m_per_100f

* Sex ratio defined as females per 1000 males; do not construct until end
* f_per_1000m

*************
* for testing; comment out the next line for a national run
keep if v024==24
*************

save IAKR7Dtemp.dta, replace

levelsof v024, local(lstates)
foreach ls of local lstates {
use IAKR7Dtemp.dta, clear
keep if v024==`ls'

* Trick: use this file to save the results
gen vstate=`ls'
gen vdist=.
gen vb=.
gen vL=.
gen vU=.
gen vcases=.

svyset v001 [pweight=v005], strata(v023) singleunit(centered)

* First do the state estimate
svy: logit m_per_f 
matrix T=r(table)
replace vb=T[1,1]  if _n==1
replace vL=T[5,1]  if _n==1
replace vU=T[6,1]  if _n==1
replace vcases=e(N) if _n==1

* Now loop through all the districts in this state
scalar sline=2
levelsof sdist, local(ldistricts)
  quietly foreach ld of local ldistricts {

  * Construct a variable for subpop to select the district
  gen select_dist=1 if sdist==`ld'

  svy, subpop(select_dist): logit m_per_f 
  matrix T=r(table)
  replace vdist=`ld' if _n==sline
  replace vb=T[1,1]  if _n==sline
  replace vL=T[5,1]  if _n==sline
  replace vU=T[6,1]  if _n==sline
  replace vcases=e(N) if _n==sline
  drop select_dist
  scalar sline=sline+1
  }

* Finished with a state; save the results for this state in a data file
drop if vb==.
keep vstate vdist vb vL vU vcases

rename v* *

* Re-attach the labels for state and district; must confirm the label names
label values state V024
label values dist SDIST

gen P_m_per_100f=100*exp(b)
gen L_m_per_100f=100*exp(L)
gen U_m_per_100f=100*exp(U)

gen P_f_per_1000m=1000*exp(-b)
gen L_f_per_1000m=1000*exp(-L)
gen U_f_per_1000m=1000*exp(-U)

save results_`ls'.dta, replace
}

format P_* L_* U_* %6.1f
list, table clean noobs