We will illustrate the application of increment-decrement life tables to the analysis of contraceptive effectiveness using a classic dataset from Tietze (1967), Intra-Uterine Contraception: Recommended Procedures for Data Analysis. Studies in Family Planning Vol 1, No 18, Supplement, pp 1-6.
Women will enter the life table with the insertion of an IUD, which defines a segment of use. The segment ends with an event or “termination” or “continues” as of the analysis date. Some events (pregnancy) close the case and are called “closures”, but others (first expulsion) terminate the segment but not necessarily the case, as the device may be “reinserted”, in which case the woman re-enters the table with a new segment of use.
The data come in two tables, one showing first insertions per
calendar month, and another giving details of terminations and
reinsertions.
I did some
preliminary work to put the data in a more usable form. The
manipulations are simple but tedious. The end result is a file
tabulating insertions, reinsertions, terminations, closures and
continuations by ordinal month. The file is called
tietze.dta and is available in the datasets section in
Stata format.
The variables include the ordinal month, counts of
insertions and reinsertions (ins and reins),
terminations by accidental pregnancy (ap), first and later
expulsion (fe,le), removal for medical or
personal reasons or by investigator’s choice (rm,
rp, ric), and lost to follow up
(lfu), as well as counts of closures by the same reasons
(with names that add a c prefix, so cfe are
closures due to first expulsions.
We start by reading the file and listing some of the counts.
. use https://grodri.github.io/datasets/tietze, clear
. list month ins reins terms cont in 1/12, clean
month ins reins terms cont
1. 1 250 7 25 18
2. 2 . 8 9 14
3. 3 . 1 7 17
4. 4 . 1 7 16
5. 5 . 2 4 13
6. 6 . 9 11 10
7. 7 . 5 7 17
8. 8 . 2 3 9
9. 9 . 0 3 17
10. 10 . 1 2 12
11. 11 . 0 0 15
12. 12 . 7 7 14
.
> library(haven)
> library(dplyr)
> tietze <- read_dta("https://grodri.github.io/datasets/tietze.dta")
> select(tietze, month, ins, reins, terms, cont) |> slice(1:12)
# A tibble: 12 × 5
month ins reins terms cont
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1 250 7 25 18
2 2 NA 8 9 14
3 3 NA 1 7 17
4 4 NA 1 7 16
5 5 NA 2 4 13
6 6 NA 9 11 10
7 7 NA 5 7 17
8 8 NA 2 3 9
9 9 NA 0 3 17
10 10 NA 1 2 12
11 11 NA 0 0 15
12 12 NA 7 7 14
We have a total of 250 first insertions. In the first month there are 7 reinsertions, 25 terminations of various kinds, and 18 continuing segments.
The first task is to compute woman-months of exposure in each month. We will follow Tietze in assuming that events are uniformly distributed within each month of use. In the first month we have 250 segments at the start and 250 + 7 - 25 - 18 = 214 at the end, and exposure is simply the average, 232 woman-months.
. gen exit = ins[1] + sum(reins - terms - cont)
. gen enter = ins[1]
. replace enter = exit[_n - 1] in 2/L
(14 real changes made)
. gen expo = (enter + exit)/2
. list month enter reins terms cont exit expo in 1/12, clean
month enter reins terms cont exit expo
1. 1 250 7 25 18 214 232
2. 2 214 8 9 14 199 206.5
3. 3 199 1 7 17 176 187.5
4. 4 176 1 7 16 154 165
5. 5 154 2 4 13 139 146.5
6. 6 139 9 11 10 127 133
7. 7 127 5 7 17 108 117.5
8. 8 108 2 3 9 98 103
9. 9 98 0 3 17 78 88
10. 10 78 1 2 12 65 71.5
11. 11 65 0 0 15 50 57.5
12. 12 50 7 7 14 36 43
> last <- nrow(tietze) > tietze <- mutate(tietze, exit = ins[1] + cumsum(reins - terms - cont), + enter = c(ins[1], exit[-last]), + expo = (enter + exit)/2) > select(tietze, month, enter, reins, terms, cont, exit, expo) |> slice(1:12) # A tibble: 12 × 7 month enter reins terms cont exit expo <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 1 250 7 25 18 214 232 2 2 214 8 9 14 199 206. 3 3 199 1 7 17 176 188. 4 4 176 1 7 16 154 165 5 5 154 2 4 13 139 146. 6 6 139 9 11 10 127 133 7 7 127 5 7 17 108 118. 8 8 108 2 3 9 98 103 9 9 98 0 3 17 78 88 10 10 78 1 2 12 65 71.5 11 11 65 0 0 15 50 57.5 12 12 50 7 7 14 36 43
We will compute continuation probabilities using all “relevant” closures, defined as pregnancy, first and later expulsion, and removal for medical or personal reasons. We are used to computing rates and then converting to probabilities (see formula 3.2 in the textbook), but we can also compute the probability directly using as denominator an adjusted number at the start of each month obtained by adding half the relevant closures to the exposure.
It is easy to verify that the two approaches given exactly the same answer. For example in the first month we have 12 relevant closures. Dividing by the 232 woman-months of exposure we obtain a rate of 0.0517 which converts to a probability of 0.0504 (using formula 3.12 in the text: q = 2 m/(2 + m)). Or we can divide 12 by the adjusted denominator of 232 + 12/2 = 238 to obtain 0.0504 directly.
Once we have the probability of discontinuation in each month we obtain the probability of continuation as the complement, and the cumulative continuation probability (unfortunately often called cumulative continuation “rate”) as a cumulative product
. gen relclos = cap + cfe + cle + crm + crp
. gen nadj = expo + relclos/2
. gen q = relclos/nadj
. gen p = 1 - q
. gen cr = exp(sum(log(p)))
. list month nadj relclos q p cr in 1/12, clean
month nadj relclos q p cr
1. 1 238 12 .05042017 .94957983 .94957983
2. 2 208.5 4 .01918465 .98081535 .93136247
3. 3 190.5 6 .03149606 .96850394 .90202822
4. 4 167.5 5 .02985075 .97014925 .87510201
5. 5 147 1 .00680272 .99319728 .86914893
6. 6 133.5 1 .00749064 .99250936 .86263845
7. 7 118.5 2 .01687764 .98312236 .84807915
8. 8 104.5 3 .02870813 .97129187 .82373238
9. 9 89.5 3 .03351955 .96648045 .79612124
10. 10 72 1 .01388889 .98611111 .785064
11. 11 57.5 0 0 1 .785064
12. 12 43.5 1 .02298851 .97701149 .76701655
> tietze <- mutate(tietze, relclos = cap + cfe + cle + crm + crp, + nadj = expo + relclos/2, + q = relclos/nadj, + p = 1 - q, + cr = cumprod(p)) > select(tietze, month, nadj, relclos, q, p, cr) |> slice(1:12) # A tibble: 12 × 6 month nadj relclos q p cr <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 1 238 12 0.0504 0.950 0.950 2 2 208. 4 0.0192 0.981 0.931 3 3 190. 6 0.0315 0.969 0.902 4 4 168. 5 0.0299 0.970 0.875 5 5 147 1 0.00680 0.993 0.869 6 6 134. 1 0.00749 0.993 0.863 7 7 118. 2 0.0169 0.983 0.848 8 8 104. 3 0.0287 0.971 0.824 9 9 89.5 3 0.0335 0.966 0.796 10 10 72 1 0.0139 0.986 0.785 11 11 57.5 0 0 1 0.785 12 12 43.5 1 0.0230 0.977 0.767
So the probability of continuing use after 12 months is 76.7%.
We now compute so-called “net” rates of events, which correspond to
the multiple decrement life table and reflect the probabilities of
experiencing various events in the presence of the competing
risks. Here is the calculation for pregnancies (technically closures due
to accidental pregnancies or cap).
We obtain a conditional monthly probability of pregnancy dividing pregnancies by the adjusted denominator, and then unconditional probabilities multiplying by the continuation probability at the beginning of the month, which can be accumulated over duration of use.
. gen qap = cap/nadj
. gen lcr = cr[_n-1] // lagged continuation rate
(1 missing value generated)
. replace lcr = 1 in 1
(1 real change made)
. gen cnap = sum(lcr * qap)
. list month nadj cap qap cnap in 1/12, clean
month nadj cap qap cnap
1. 1 238 0 0 0
2. 2 208.5 0 0 0
3. 3 190.5 2 .01049869 .00977808
4. 4 167.5 3 .01791045 .02593381
5. 5 147 0 0 .02593381
6. 6 133.5 0 0 .02593381
7. 7 118.5 1 .00843882 .03321346
8. 8 104.5 0 0 .03321346
9. 9 89.5 0 0 .03321346
10. 10 72 1 .01388889 .0442707
11. 11 57.5 0 0 .0442707
12. 12 43.5 0 0 .0442707
> tietze <- mutate(tietze, qap = cap/nadj, + lcr = c(1, cr[-last]), + cnap = cumsum(lcr * qap)) > select(tietze, month, nadj, cap, qap, cnap) |> slice(1:12) # A tibble: 12 × 5 month nadj cap qap cnap <dbl> <dbl> <dbl> <dbl> <dbl> 1 1 238 0 0 0 2 2 208. 0 0 0 3 3 190. 2 0.0105 0.00978 4 4 168. 3 0.0179 0.0259 5 5 147 0 0 0.0259 6 6 134. 0 0 0.0259 7 7 118. 1 0.00844 0.0332 8 8 104. 0 0 0.0332 9 9 89.5 0 0 0.0332 10 10 72 1 0.0139 0.0443 11 11 57.5 0 0 0.0443 12 12 43.5 0 0 0.0443
The probability of accidental pregnancy in the first year of use is 4.4%. Try computing similar probabilities for the other relevant closures. You should get 4.3% for first expulsions, 3.6% for later expulsions, 6.7% for medical reasons and 4.3% for personal reasons. These add up to 23.3%, which is of course the complement of the probability of continuing use after 12 months.
These “rates” can also be computed for events other than closures using exactly the same procedure. Here is the probability of a first expulsion in the presence of the other risks:
. gen qfe = fe/nadj
. gen cnfe = sum(lcr * qfe)
. list month nadj fe qfe cnfe in 1/12, clean
month nadj fe qfe cnfe
1. 1 238 15 .06302521 .06302521
2. 2 208.5 4 .01918465 .08124257
3. 3 190.5 2 .01049869 .09102065
4. 4 167.5 0 0 .09102065
5. 5 147 2 .01360544 .1029268
6. 6 133.5 1 .00749064 .10943728
7. 7 118.5 3 .02531646 .13127623
8. 8 104.5 1 .00956938 .13939182
9. 9 89.5 0 0 .13939182
10. 10 72 0 0 .13939182
11. 11 57.5 0 0 .13939182
12. 12 43.5 1 .02298851 .15743927
> tietze <- mutate(tietze, qfe = fe/nadj, + cnfe = cumsum(lcr * qfe)) > select(tietze, month, nadj, fe, qfe, cnfe) |> slice(1:12) # A tibble: 12 × 5 month nadj fe qfe cnfe <dbl> <dbl> <dbl> <dbl> <dbl> 1 1 238 15 0.0630 0.0630 2 2 208. 4 0.0192 0.0812 3 3 190. 2 0.0105 0.0910 4 4 168. 0 0 0.0910 5 5 147 2 0.0136 0.103 6 6 134. 1 0.00749 0.109 7 7 118. 3 0.0253 0.131 8 8 104. 1 0.00957 0.139 9 9 89.5 0 0 0.139 10 10 72 0 0 0.139 11 11 57.5 0 0 0.139 12 12 43.5 1 0.0230 0.157
We get a 15.7% probability of first expulsion in the first year of use. (Recall that first expulsions are not “fatal”, as the IUD can be reinserted.) Try computing similar rates for later expulsions and removals for medical and personal reasons. You should get 6.1%, 8.3% and 4.3% respectively. For pregnancies the answer is the same as above because every termination is a closure.
We can also compute the equivalent of an associated single-decrement life table under the assumption of independence of risks. The results are often called “gross” rates. The calculation proceeds as before but we use a new adjusted denominator, which adds to the woman-months of exposure half the events of interest (as opposed to half of all relevant closures). Let us do the calculation for pregnancies
. gen nadjap = expo + 0.5 * cap // pregnancies are only risk
. gen gqap = cap/nadjap
. gen gcrap = exp(sum(log(1 - gqap))) // continuation
. gen gdrap = 1 - gcrap // discontinuation
. list month gqap gcrap gdrap in 1/12, clean
month gqap gcrap gdrap
1. 1 0 1 0
2. 2 0 1 0
3. 3 .01061008 .98938992 .01061008
4. 4 .01801802 .97156308 .02843692
5. 5 0 .97156308 .02843692
6. 6 0 .97156308 .02843692
7. 7 .00847458 .96332949 .03667051
8. 8 0 .96332949 .03667051
9. 9 0 .96332949 .03667051
10. 10 .01388889 .94994991 .05005009
11. 11 0 .94994991 .05005009
12. 12 0 .94994991 .05005009
> tietze <- mutate(tietze, nadjap = expo + cap/2, + gqap = cap/nadjap, + gcrap = cumprod(1 - gqap), + gdrap = 1 - gcrap) > select(tietze, month, gqap, gcrap, gdrap) |> slice(1:12) # A tibble: 12 × 4 month gqap gcrap gdrap <dbl> <dbl> <dbl> <dbl> 1 1 0 1 0 2 2 0 1 0 3 3 0.0106 0.989 0.0106 4 4 0.0180 0.972 0.0284 5 5 0 0.972 0.0284 6 6 0 0.972 0.0284 7 7 0.00847 0.963 0.0367 8 8 0 0.963 0.0367 9 9 0 0.963 0.0367 10 10 0.0139 0.950 0.0501 11 11 0 0.950 0.0501 12 12 0 0.950 0.0501
The probability of accidental pregnancy would be a bit higher, 5%, in the absence of the other risks. The calculation is based on the assumption of independence of risks, as the observed probabilities of pregnancy are applied to women who discontinued for other reasons.
Note that the demographic literature uses “net” for rates in the presence of other risks and “gross” for rates in the absence of other risks. The statistical literature reverses the usage.
A useful article in addition to Tietze’s classic is S. R. Karia, J. Toivonen and E. Arjas (1998), Analysis of Contraceptive Failure Data in Intrauterine Device Studies, Contraception, 58:361-374.