*Rieger M, Wagner N., Polygyny and child health revisited. Proc Natl Acad Sci U S A. 2016 Mar 29;113(13):E1769-70. 
*rieger@iss.nl
*Comment on Lawson et al. No evidence that polygynous marriage is a harmful cultural practice in northern Tanzania. Proc Natl Acad Sci U S A. 2015 Nov 10;112(45):13827-32. 
*Data available from https://www.pnas.org/doi/10.1073/pnas.1507151112#supplementary-materials 
clear
import excel "...pnas.1507151112.sd01.xlsx", sheet("Sheet1") firstrow

tab SeasonofDataCollection, gen(sss)
tab HouseholdType, gen(poly)
tab ChildSex, gen(gender)

gen female=gender1

gen poly=poly2+poly4

gen poly_male=poly4
gen poly_female=poly2
gen mono_female=poly1

*Replicate Table S3 for HAZ
xtset VillageID
reg ChildHAZ poly ChildAgeinMonths ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears, cluster(VillageID) 
xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1  || VillageID:,  mle

*Create Figure 1
twoway (lpoly ChildHAZ ChildAgeinMonths if poly==1, legend(size(vsmall) lab(2 "Monogamous household" ) lab(1 "Polygynous household")) ///
ytitle("Height-for-age z-score") ylabel(-2(0.5)0.5) xmtick(0(5)60) xtitle("Age in months") xlabel(0(10)60) ///
xmtick(0(5)60) acol(blue) fcol(none) alw(thin) alpat(shortdash) col(blue) lpat(solid)) ///
(lpoly ChildHAZ ChildAgeinMonths if poly==0, acol(red) fcol(none) alw(thin) alpat(shortdash) col(red) lpat(solid))

*Age Heterogeneities
eststo haz_global: xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1  || VillageID:,  mle
eststo haz_old: xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1 || VillageID: if ChildAgeinMonths>=30,  mle
eststo haz_young: xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1 || VillageID: if ChildAgeinMonths<30,  mle

eststo haz_old_fe: xtreg ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1  if ChildAgeinMonths>=30,  fe
eststo haz_young_fe: xtreg ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1 if ChildAgeinMonths<30,  fe

eststo haz_old_d: xtmixed ChildHAZ poly_male poly_female mono_female ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1 || VillageID: if ChildAgeinMonths>=30,  mle
eststo haz_young_d: xtmixed ChildHAZ poly_male poly_female mono_female ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1 || VillageID: if ChildAgeinMonths<30,  mle

eststo haz_old_v: xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1    PercentageofHouseholdsSampled EstimatedMonthlyRainfall PercentageofHouseholdHeadswi DistanceofVillagefromDistric || VillageID: if ChildAgeinMonths>=30,  mle
eststo haz_young_v: xtmixed ChildHAZ poly ChildAgeinMonths  ChildAgeinMonthsSquared female AgeofHouseholdHeadinYears sss1   PercentageofHouseholdsSampled EstimatedMonthlyRainfall PercentageofHouseholdHeadswi DistanceofVillagefromDistric || VillageID: if ChildAgeinMonths<30,  mle

*Build Table 1
estout haz_global haz_old haz_young haz_old_v haz_young_v, cells(b(star fmt(2)) ci(fmt(2))) stats(N) starlevels(* 0.10 ** 0.05 *** 0.01)