[R] analysis of life tables

[Göran Broström]

On Sun, Aug 15, 2004 at 02:28:03PM +0200, Christoph Scherber wrote:
> Dear all,
> 
> How can I analyze a life table (e.g. for a cohort of insects) in R?
> 
> I have 20 insects in 200 cages with two different treatments, whose 
> survival is followed over time, such that, e.g., in one treatment, the 
> number of animals surviving is c(20,18,16,12,10,8,4,0), while in the 
> other treatment the survival is c(20,20,18,18,16,15,15,14) at 8 
> subsequent time intervals.

One way of doing it would be to create an individual-based data frame
from your two life tables, which would give you 

> dat
   enter exit event treatment
1      0    1     1         1
2      0    1     1         1
3      0    2     1         1
4      0    2     1         1
5      0    3     1         1
6      0    3     1         1
7      0    3     1         1
8      0    3     1         1
9      0    4     1         1
10     0    4     1         1
11     0    5     1         1
12     0    5     1         1
13     0    6     1         1
14     0    6     1         1
15     0    6     1         1
16     0    6     1         1
17     0    7     1         1
18     0    7     1         1
19     0    7     1         1
20     0    7     1         1
21     0    2     1         2
22     0    2     1         2
23     0    4     1         2
24     0    4     1         2
25     0    5     1         2
26     0    7     1         2
27     0    7     0         2
28     0    7     0         2
29     0    7     0         2
30     0    7     0         2
31     0    7     0         2
32     0    7     0         2
33     0    7     0         2
34     0    7     0         2
35     0    7     0         2
36     0    7     0         2
37     0    7     0         2
38     0    7     0         2
39     0    7     0         2
40     0    7     0         2

where 'exit' represents the ordered time intervals. Then you can choose 
between a discrete-time (recommended) and a continuous-time Cox regression.
In package 'eha' that corresponds to the functions 'mlreg' and 'coxreg',
respectively. Output from mlreg:
----------------------------------------------------------------------
> mlreg(Surv(enter, exit, event) ~ treatment, data = dat)
Call:
mlreg(formula = Surv(enter, exit, event) ~ treatment, data = dat)

Covariate           Mean       Coef  Rel.Risk       Wald p
treatment
               1    0.419     0         1           (reference)
               2    0.581    -2.121     0.120        0.000

Events                    26
Total time at risk           210
Max. log. likelihood      -64.734
LR test statistic         21.9
Degrees of freedom        1
Overall p-value           2.82840e-06
-----------------------------------------------------------------------

This corresponds essentially to a binomial regression approach, where you 
regard the number of deaths in each time interval for each treatment as 
the outcome of a binomial experiment with n = riskset size at the beginning
of the interval.

An alternative would be Poisson regression with riskset size as an offset.
-- 
 Göran Broström                    tel: +46 90 786 5223
 Department of Statistics          fax: +46 90 786 6614
 Umeå University                   http://www.stat.umu.se/egna/gb/
 SE-90187 Umeå, Sweden             e-mail: gb at stat.umu.se