# Biostatistics III in R

#### Code by Annika Tillander, Andreas Karlsson and Mark Clements

# Exercise 1: Life tables and Kaplan-Meier estimates of survival

## (a) Hand calculation: Life table and Kaplan-Meier estimates of survival

The results are contained in an Excel file and are also shown in the R output below.

## (b) Using R to validate the hand calculations done in part 1 (a)

First, load the `biostat3` library:

    library(biostat3)

Following are the life table estimates. Note that in the lectures, when we estimated all-cause survival, there were 8 deaths in the first interval. One of these died of a cause other than cancer so in the cause-specific survival analysis we see that there are 7 ‘deaths’ and 1 censoring (Stata uses the term ‘lost’ for lost to follow-up) in the first interval.

    print(lifetab2(Surv(floor(surv_yy), status == "Dead: cancer")~1, colon_sample, breaks=0:10), digits=2)

    ##        nsubs nlost nrisk nevent surv   pdf hazard se.surv se.pdf se.hazard
    ## 0-1       35     1  34.5      7 1.00 0.203   0.23   0.000  0.068     0.085
    ## 1-2       27     3  25.5      1 0.80 0.031   0.04   0.068  0.031     0.040
    ## 2-3       23     4  21.0      5 0.77 0.182   0.27   0.073  0.073     0.120
    ## 3-4       14     1  13.5      2 0.58 0.086   0.16   0.090  0.058     0.113
    ## 4-5       11     1  10.5      0 0.50 0.000   0.00   0.095    NaN       NaN
    ## 5-6       10     0  10.0      0 0.50 0.000   0.00   0.095    NaN       NaN
    ## 6-7       10     3   8.5      0 0.50 0.000   0.00   0.095    NaN       NaN
    ## 7-8        7     1   6.5      0 0.50 0.000   0.00   0.095    NaN       NaN
    ## 8-9        6     4   4.0      1 0.50 0.124   0.29   0.095  0.110     0.283
    ## 9-10       1     1   0.5      0 0.37 0.000   0.00   0.129    NaN       NaN
    ## 10-Inf     0     0   0.0      0 0.37    NA     NA   0.129     NA        NA

Following is a table of Kaplan-Meier estimates. Although it’s not clear from the table, the person censored (lost) at time 2 was at risk when the other person dies at time 2. On the following is a graph of the survival function.

    mfit <- survfit(Surv(surv_mm, status == "Dead: cancer") ~ 1, data = colon_sample) # make Kaplan-Meier estimates
    summary(mfit)                                                  # print Kaplan-Meier table

    ## Call: survfit(formula = Surv(surv_mm, status == "Dead: cancer") ~ 1,
    ##     data = colon_sample)
    ##
    ##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
    ##     2     35       1    0.971  0.0282        0.918        1.000
    ##     3     33       1    0.942  0.0398        0.867        1.000
    ##     5     32       1    0.913  0.0482        0.823        1.000
    ##     7     31       1    0.883  0.0549        0.782        0.998
    ##     8     30       1    0.854  0.0605        0.743        0.981
    ##     9     29       1    0.824  0.0652        0.706        0.962
    ##    11     28       1    0.795  0.0692        0.670        0.943
    ##    22     24       1    0.762  0.0738        0.630        0.921
    ##    27     22       1    0.727  0.0781        0.589        0.898
    ##    28     20       1    0.691  0.0823        0.547        0.872
    ##    32     19       2    0.618  0.0882        0.467        0.818
    ##    33     16       1    0.579  0.0908        0.426        0.788
    ##    43     13       1    0.535  0.0941        0.379        0.755
    ##    46     12       1    0.490  0.0962        0.334        0.720
    ##   102      4       1    0.368  0.1284        0.185        0.729

    plot(mfit,                                                     # plot Kaplan-Meier curve
         ylab="S(t)",
         xlab="Time since diagnosis in months",
         main = "Kaplan−Meier estimates of cause−specific survival")