We will look at the interval from second to third birth using Colombian WFS data, circa 1976. We start from an extract that has the dates of R’s birth, date of interview, the birth history, date of first union if any, and current and childhood type of place of residence.
. use https://grodri.github.io/datasets/cofertx, clear (COSR02 extract)
We will use intervals that start in the ten years before the interview and exclude twins. (The results are a bit different from my paper with Hobcraft, which used all birth intervals, but I can reproduce the earlier results by removing the period restriction.)
. keep if b022 >= v007-120 & b022 < v007 (4,144 observations deleted) . drop if b032==b022 // only 6 (6 observations deleted)
We construct time at risk starting in the middle of the month of birth of the second child and ending in the middle of the month when the third child is born or at the end of the month before the interview, whichever occurs first
. gen expo = b032 - b022 . replace expo = v007 - b022 -0.5 if v007 <= b032 (496 real changes made) . gen third = b032 < v007
If we stset the data we can take advantage of Stata’s
survival analysis facilities. For example it is very easy to plot a
Kaplan-Meier estimate of survival at parity 2:
. gen id=_n
. stset expo, fail(third) id(id)
Survival-time data settings
ID variable: id
Failure event: third!=0 & third<.
Observed time interval: (expo[_n-1], expo]
Exit on or before: failure
──────────────────────────────────────────────────────────────────────────
1,228 total observations
0 exclusions
──────────────────────────────────────────────────────────────────────────
1,228 observations remaining, representing
1,228 subjects
732 failures in single-failure-per-subject data
34,388 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 115.5
. sts graph
Failure _d: third
Analysis time _t: expo
ID variable: id
. graph export co3rdkm.png, replace
file co3rdkm.png saved as PNG format
We will take advantage of Stata’s facilities to split the exposure into 3 month segments. The cutpoints have the form 0 3.5 6.5 9.5 … 57.5 60.5 120.
. stsplit segment, at(0 3.5(3)60.5 120) (9,844 observations (episodes) created)
Stata’s built-in variables _t and _t0 have
the start and end of each segment and _d has the death
indicator. We use these to compute events and exposure. This is not a
bad time to save the data.
. replace expo = _t-_t0 (9,844 real changes made) . gen births = _d . save temp,replace file temp.dta saved
We then collapse to obtain total events and exposure by duration (segment) and childhood place of residence.
. drop if v703 > 3 // a few missing values (36 observations deleted) . collapse (sum) births (sum) expo , by(segment v703)
At this point we can easily calculate life tables.
I will compute the hazard dividing births by exposure and then use the constant hazard assumption to estimate the cumulative hazard and survival (or rather its complement, the birth function). You can get very similar results assuming that events occur half-way through each three-month interval. (Even more similar if we had used single months.)
. gen h = births/expo
. bysort v703 (segment): gen H = 3*sum(h)
. gen B = 1-exp(-H) // the birth function
. separate B, by(v703)
Variable Storage Display Value
name type format label Variable label
────────────────────────────────────────────────────────────────────────────────
B1 float %9.0g B, v703 == Rural
B2 float %9.0g B, v703 == Town
B3 float %9.0g B, v703 == City
. line B1 B2 B3 segment, lpat(solid longdash dash) ///
> title("The Third Child in Colombia, 1976") ///
> subtitle(Birth Functions by Childhood Residence) ///
> ytitle(Proportion having a 3rd birth) ///
> xlabel(12 24 36 48 60) xtitle(months since 2nd birth) ///
> legend(order(1 "Rural" 2 "Town" 3 "City") ///
> col(1) ring(0) pos(5))
. graph export co3rdbirth.png, replace
file co3rdbirth.png saved as PNG format
We see that the proportions moving to parity three are higher for women of rural origins than for others at every duration. Women who grew up in towns are a bit less likely to have a third child, and when they do they have it later. Women who grew up in cities are much less likely to have a third child in the long run, but they are also more likely than townfolk to have very short intervals, with a cross-over around 18 months.
To compute quintums we need to interpolate in the category 57.5-60.5. Linear interpolation would probably do, but with constant hazard a more exact calculation based on S(60) = S(57.5) exp{-2.5 h(57.5)} gives Q = 1-(1-B(57.5))exp{-2.5 h(57.5)}. Note that the way we calculated things, the survival and birth functions pertain to the end of each segment, which is why the code below looks at 54.5
. gen Q =1-(1-B[_n]) * exp(-2.5*h[_n+1]) if segment==54.5
(60 missing values generated)
. list v703 segment h B Q if segment > 54 & v703 < 9, sep(3)
┌──────────────────────────────────────────────────┐
│ v703 segment h B Q │
├──────────────────────────────────────────────────┤
19. │ Rural 54.5 .0207612 .8118494 .8230233 │
20. │ Rural 57.5 .0244898 .8251771 . │
21. │ Rural 60.5 .0086022 .829631 . │
├──────────────────────────────────────────────────┤
40. │ Town 54.5 .0363636 .7815775 .7930105 │
41. │ Town 57.5 .0215054 .7952242 . │
42. │ Town 60.5 .0059435 .7988432 . │
├──────────────────────────────────────────────────┤
61. │ City 54.5 .0137931 .6706035 .6766252 │
62. │ City 57.5 .0073801 .6778163 . │
63. │ City 60.5 .0223537 .6987138 . │
└──────────────────────────────────────────────────┘
So the quintums are 82.3, 79.3 and 67.7 for women who grew up in rural areas, towns and cities. (These are a bit lower than in my paper with Hobcraft, reflecting the fact that they pertain to more recent births.)
To obtain the trimean we use a similar interpolation strategy for computing the quartiles and then apply Tukey’s formula. I will illustrate with rural women. A quarter of the quintum of .823 is .206, so we need the durations by which 20.6, 41.2 and 61.7% of the women have moved to parity 3. For the first value we need to interpolate in category 9.5-12.5. Write 1 - 0.206 = (1-B(9.5)) exp{-2.5 h(9.5)(Q1-9.5)} and solve for Q1 to obtain Q1 = 9.5 -log((1-0.206)/(1-B(9.5)))/h(9.5). Because the calculations are repetitive I encapsulate them in a little program which takes as argument a category of v703, finds the row where we need to interpolate, and applies the above formula:
. capture program drop trimean
. program define trimean
1. args ctpr
2. quietly sum Q if v703==`ctpr'
3. local q = r(mean)/4
4. local T = 0
5. tempvar i
6. gen `i' = _n
7. forvalues k=1/3 {
8. local b = `k'*`q'
9. quietly sum `i' if v703==`ctpr' & B < `b'
10. local n = r(max)
11. local Q = seg[`n']-log((1-`b')/(1-B[`n']))/h[`n'+1]
12. display "Q`k' = `Q'"
13. local T = `T' + `Q'
14. if (`k' == 2) local T = `T' + `Q'
15. }
16. local T = `T'/4
17. display "Trimean: `T'"
18. end
We now run the computations for the three categories
. trimean 1 Q1 = 13.04612965753713 Q2 = 18.84358957514994 Q3 = 27.81360231842412 Trimean: 19.63672778156528 . trimean 2 Q1 = 14.75915079139585 Q2 = 21.72924126033396 Q3 = 30.01968865537609 Trimean: 22.05933049185996 . trimean 3 Q1 = 13.15086885691453 Q2 = 19.51578799192259 Q3 = 28.77935879083599 Trimean: 20.24045090789892
So the average birth intervals for those who go on to have a third child are 19.6, 22.1 and 20.2 months, respectively, for women who grew up in rural areas, towns and cities. (The intermediate value for the last group was also observed in the analysis using all intervals.)
The fact that the hazard functions are somewhat noisy is easily verified:
. separate h, by(v703)
Variable Storage Display Value
name type format label Variable label
────────────────────────────────────────────────────────────────────────────────
h1 float %9.0g h, v703 == Rural
h2 float %9.0g h, v703 == Town
h3 float %9.0g h, v703 == City
. gen d = segment+1.5
. line h1 h2 h3 d, lpat(solid longdash dash) ///
> title("The Third Child in Colombia, 1976") ///
> subtitle(Hazard Functions by Childhood Residence) ///
> xtitle(duration since 2nd birth) ///
> ytitle(hazard of having 3rd birth) ///
> legend(order(1 "Rural" 2 "Town" 3 "City") ring(0) col(1) pos(2))
. graph export co3rdhaz.png, replace
file co3rdhaz.png saved as PNG format
Fitting a proportional hazards model is quite easy using the Poisson trick. (This is covered in some detail in the GLM course.)
Let us construct dummy variables for duration and childhood residence and then run a regression for the first five years using the log of exposure time as an offset.
. gen town = v703==2
. gen city = v703==3
. forvalues d = 12(3)57 {
2. gen dur`d' = segment == `d'.5
3. }
. gen os = log(expo)
. gen touse = seg > 9 & seg < 60
. poisson births dur12-dur57 town city if touse, offset(os)
Iteration 0: log likelihood = -126.22561
Iteration 1: log likelihood = -125.00294
Iteration 2: log likelihood = -124.9991
Iteration 3: log likelihood = -124.9991
Poisson regression Number of obs = 51
LR chi2(18) = 120.26
Prob > chi2 = 0.0000
Log likelihood = -124.9991 Pseudo R2 = 0.3248
─────────────┬────────────────────────────────────────────────────────────────
births │ Coefficient Std. err. z P>|z| [95% conf. interval]
─────────────┼────────────────────────────────────────────────────────────────
dur12 │ .8387151 .1677175 5.00 0.000 .5099949 1.167435
dur15 │ 1.030044 .1672236 6.16 0.000 .7022913 1.357796
dur18 │ .9195063 .1758572 5.23 0.000 .5748326 1.26418
dur21 │ 1.00266 .1794974 5.59 0.000 .650851 1.354468
dur24 │ 1.131769 .1816092 6.23 0.000 .7758211 1.487716
dur27 │ .7913616 .2053167 3.85 0.000 .3889482 1.193775
dur30 │ .9017034 .209525 4.30 0.000 .4910419 1.312365
dur33 │ .3386295 .2624972 1.29 0.197 -.1758555 .8531145
dur36 │ .3342246 .2728632 1.22 0.221 -.2005775 .8690267
dur39 │ .3081459 .2853069 1.08 0.280 -.2510455 .8673372
dur42 │ .7543197 .2579441 2.92 0.003 .2487585 1.259881
dur45 │ -.1988996 .4022531 -0.49 0.621 -.9873012 .5895019
dur48 │ -.4680164 .4679286 -1.00 0.317 -1.38514 .4491068
dur51 │ .3255532 .3448924 0.94 0.345 -.3504236 1.00153
dur54 │ .3411703 .3606533 0.95 0.344 -.3656972 1.048038
dur57 │ .0717598 .4308883 0.17 0.868 -.7727656 .9162853
town │ -.2020243 .088543 -2.28 0.023 -.3755654 -.0284831
city │ -.3866947 .0977489 -3.96 0.000 -.578279 -.1951103
_cons │ -3.951717 .1419443 -27.84 0.000 -4.229923 -3.673512
os │ 1 (offset)
─────────────┴────────────────────────────────────────────────────────────────
. estat gof
Deviance goodness-of-fit = 51.44754
Prob > chi2(32) = 0.0161
Pearson goodness-of-fit = 45.39549
Prob > chi2(32) = 0.0587
The duration coefficients track an early increase in hazard up to about a year and a half followed by a steady decline thereafter. On average, at any given duration since second birth the risk of having a third child is 18.3% lower for women who grew up in towns and 32.1% lower for women who grew up in cities than for those who grew up in rural areas. The goodness of fit test indicates, however, that the hazards for the three categories of childhood residence are not proportional. This is consistent with the differences in timing between women with town and city origins noted earlier.