Germán Rodríguez
Demographic Methods Princeton University

Life Table Analysis of Contraceptive Use

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. Both tables are available in the datasets section of the course website as TietzeInsertions.dat and TietzeTerminations.dat. 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 http://data.princeton.edu/eco572/datasets/Tietze

. 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(foreign)

> tietze <- read.dta("http://data.princeton.edu/eco572/datasets/Tietze.dta")

> select(tietze, month, ins, reins, terms, cont) %>% slice(1:12)
   month ins reins terms cont
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.

Woman-Months of Use

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)
   month enter reins terms cont exit  expo
1      1   250     7    25   18  214 232.0
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.0
5      5   154     2     4   13  139 146.5
6      6   139     9    11   10  127 133.0
7      7   127     5     7   17  108 117.5
8      8   108     2     3    9   98 103.0
9      9    98     0     3   17   78  88.0
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.0

Continuation "Rates"

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   .0504202   .9495798   .9495798  
  2.       2   208.5         4   .0191847   .9808154   .9313625  
  3.       3   190.5         6   .0314961    .968504   .9020283  
  4.       4   167.5         5   .0298507   .9701493    .875102  
  5.       5     147         1   .0068027   .9931973    .869149  
  6.       6   133.5         1   .0074906   .9925094   .8626385  
  7.       7   118.5         2   .0168776   .9831223   .8480791  
  8.       8   104.5         3   .0287081   .9712918   .8237324  
  9.       9    89.5         3   .0335196   .9664804   .7961212  
 10.      10      72         1   .0138889   .9861111    .785064  
 11.      11    57.5         0          0          1    .785064  
 12.      12    43.5         1   .0229885   .9770115   .7670165  
> 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)
   month  nadj relclos           q         p        cr
1      1 238.0      12 0.050420168 0.9495798 0.9495798
2      2 208.5       4 0.019184652 0.9808153 0.9313625
3      3 190.5       6 0.031496063 0.9685039 0.9020282
4      4 167.5       5 0.029850746 0.9701493 0.8751020
5      5 147.0       1 0.006802721 0.9931973 0.8691489
6      6 133.5       1 0.007490637 0.9925094 0.8626385
7      7 118.5       2 0.016877637 0.9831224 0.8480792
8      8 104.5       3 0.028708134 0.9712919 0.8237324
9      9  89.5       3 0.033519553 0.9664804 0.7961212
10    10  72.0       1 0.013888889 0.9861111 0.7850640
11    11  57.5       0 0.000000000 1.0000000 0.7850640
12    12  43.5       1 0.022988506 0.9770115 0.7670166

So the probability of continuing use after 12 months is 76.7%.

Net Rates of Events

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   .0104987   .0097781  
  4.       4   167.5     3   .0179104   .0259338  
  5.       5     147     0          0   .0259338  
  6.       6   133.5     0          0   .0259338  
  7.       7   118.5     1   .0084388   .0332135  
  8.       8   104.5     0          0   .0332135  
  9.       9    89.5     0          0   .0332135  
 10.      10      72     1   .0138889   .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)
   month  nadj cap         qap        cnap
1      1 238.0   0 0.000000000 0.000000000
2      2 208.5   0 0.000000000 0.000000000
3      3 190.5   2 0.010498688 0.009778084
4      4 167.5   3 0.017910448 0.025933813
5      5 147.0   0 0.000000000 0.025933813
6      6 133.5   0 0.000000000 0.025933813
7      7 118.5   1 0.008438819 0.033213462
8      8 104.5   0 0.000000000 0.033213462
9      9  89.5   0 0.000000000 0.033213462
10    10  72.0   1 0.013888889 0.044270702
11    11  57.5   0 0.000000000 0.044270702
12    12  43.5   0 0.000000000 0.044270702

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   .0630252   .0630252  
  2.       2   208.5    4   .0191847   .0812426  
  3.       3   190.5    2   .0104987   .0910207  
  4.       4   167.5    0          0   .0910207  
  5.       5     147    2   .0136054   .1029268  
  6.       6   133.5    1   .0074906   .1094373  
  7.       7   118.5    3   .0253165   .1312762  
  8.       8   104.5    1   .0095694   .1393918  
  9.       9    89.5    0          0   .1393918  
 10.      10      72    0          0   .1393918  
 11.      11    57.5    0          0   .1393918  
 12.      12    43.5    1   .0229885   .1574393  
> tietze <- mutate(tietze, qfe = fe/nadj,   
+ cnfe = cumsum(lcr * qfe))

> select(tietze, month, nadj, fe, qfe, cnfe) %>% slice(1:12)
   month  nadj fe         qfe       cnfe
1      1 238.0 15 0.063025210 0.06302521
2      2 208.5  4 0.019184652 0.08124257
3      3 190.5  2 0.010498688 0.09102065
4      4 167.5  0 0.000000000 0.09102065
5      5 147.0  2 0.013605442 0.10292680
6      6 133.5  1 0.007490637 0.10943728
7      7 118.5  3 0.025316456 0.13127623
8      8 104.5  1 0.009569378 0.13939182
9      9  89.5  0 0.000000000 0.13939182
10    10  72.0  0 0.000000000 0.13939182
11    11  57.5  0 0.000000000 0.13939182
12    12  43.5  1 0.022988506 0.15743927

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.

Gross Rates of Events

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   .0106101   .9893899   .0106101  
  4.       4    .018018   .9715631   .0284369  
  5.       5          0   .9715631   .0284369  
  6.       6          0   .9715631   .0284369  
  7.       7   .0084746   .9633295   .0366705  
  8.       8          0   .9633295   .0366705  
  9.       9          0   .9633295   .0366705  
 10.      10   .0138889   .9499499   .0500501  
 11.      11          0   .9499499   .0500501  
 12.      12          0   .9499499   .0500501  
> 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)
   month        gqap     gcrap      gdrap
1      1 0.000000000 1.0000000 0.00000000
2      2 0.000000000 1.0000000 0.00000000
3      3 0.010610080 0.9893899 0.01061008
4      4 0.018018018 0.9715631 0.02843692
5      5 0.000000000 0.9715631 0.02843692
6      6 0.000000000 0.9715631 0.02843692
7      7 0.008474576 0.9633295 0.03667051
8      8 0.000000000 0.9633295 0.03667051
9      9 0.000000000 0.9633295 0.03667051
10    10 0.013888889 0.9499499 0.05005009
11    11 0.000000000 0.9499499 0.05005009
12    12 0.000000000 0.9499499 0.05005009

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.