# Biostatistics III in R

#### Code by Andreas Karlsson and Mark Clements

# Exercise 22. Estimating the effect of a time-varying exposure – the bereavement data

These data were used to study a possible effect of marital bereavement (loss of husband or wife) on all–cause mortality in the elderly. The dataset was extracted from a larger follow-up study of an elderly population and concerns subjects whose husbands or wives were alive at entry to the study. Thus all subjects enter as not bereaved but may become bereaved at some point during follow–up. The variable dosp records the date of death of each subject’s spouse and takes the value 1/1/2000 where this has not yet happened.

You may have to install the required packages the first time you use them. You can install a package by `install.packages("package_of_interest")` for each package you require.

    library(biostat3) # for Surv and survfit
    library(dplyr)    # for data manipulation
    library(rstpm2)

## (a)

    ## Look at the first five couples
    brv |> select(couple, id, sex, doe, dosp, dox, fail) |>
        filter(couple<=5) |> arrange(couple, id) |> kable("html")

| couple | id | sex | doe | dosp | dox | fail
--- | --- | --- | --- | --- | --- | --- | ---
197 | 1 | 48 | 1 | 1981-01-15 | 2000-01-01 | 1983-04-18 | 1
270 | 2 | 49 | 2 | 1981-04-27 | 1983-01-24 | 1984-11-20 | 1
177 | 2 | 263 | 1 | 1981-04-27 | 1984-11-20 | 1983-01-24 | 1
168 | 3 | 60 | 1 | 1981-01-20 | 1981-12-31 | 1981-08-03 | 1
384 | 3 | 63 | 2 | 1981-01-20 | 1981-08-03 | 1981-12-31 | 1
12 | 4 | 156 | 1 | 1981-01-20 | 1988-11-23 | 1991-01-01 | 0
300 | 4 | 220 | 2 | 1981-01-20 | 2000-01-01 | 1988-11-23 | 1
212 | 5 | 361 | 1 | 1981-05-07 | 1984-02-04 | 1984-02-04 | 1

## (b)

1. The timescale is attained age, which would seem to be a reasonable choice.
2. Males have the higher mortality which is to be expected.
3. Age could potentially be a confounder. Males are slightly older at diagnosis (although we haven’t studied pairwise differences).

    brv <- mutate(brv, age_entry = as.numeric(doe - dob) / 365.24, # Calc age at entry
               att_age = as.numeric(dox - dob) / 365.24,   # Calc attained age
               t_at_risk = att_age - age_entry)            # Calc time at risk

    ## crude rates
    survRate(Surv(age_entry, att_age, fail) ~ sex, data=brv)

    ##       sex    tstop event       rate      lower     upper
    ## sex=1   1 1340.521   181 0.13502210 0.11606748 0.1561895
    ## sex=2   2 1095.187    97 0.08856937 0.07182379 0.1080471

    poisson22b <- glm( fail ~ sex + offset( log( t_at_risk) ), family=poisson, data=brv)
    summary( poisson22b )

    ##
    ## Call:
    ## glm(formula = fail ~ sex + offset(log(t_at_risk)), family = poisson,
    ##     data = brv)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)  -1.5807     0.1800  -8.780  < 2e-16 ***
    ## sex          -0.4217     0.1258  -3.351 0.000806 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 526.28  on 398  degrees of freedom
    ## Residual deviance: 514.64  on 397  degrees of freedom
    ## AIC: 1074.6
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22b)

    ##             exp(beta)     2.5 %    97.5 %
    ## (Intercept) 0.2058383 0.1446403 0.2929295
    ## sex         0.6559620 0.5125885 0.8394379

    t.test(age_entry~sex, data=brv)

    ##
    ##  Welch Two Sample t-test
    ##
    ## data:  age_entry by sex
    ## t = 1.3355, df = 385.25, p-value = 0.1825
    ## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
    ## 95 percent confidence interval:
    ##  -0.1943401  1.0174920
    ## sample estimates:
    ## mean in group 1 mean in group 2
    ##        79.07153        78.65995

## (c) Breaking records into pre and post bereavement.

    ## Creating times relative to spouse death (year=0)
    brv2 <- mutate(brv,
                   id=NULL,
                   y_before_sp_dth =  as.numeric(doe -dosp) / 365.24,
                   y_after_sp_dth = as.numeric(dox - dosp) / 365.24)

    ## Splitting at spouse death (year=0)
    brvSplit <- survSplit(brv2, cut = 0, end="y_after_sp_dth", start="y_before_sp_dth",
                          id="id",event="fail")

    ## Calculating risk times
    brvSplit <- mutate(brvSplit,
                       t_sp_at_risk =   y_after_sp_dth - y_before_sp_dth,
                       brv = ifelse(y_after_sp_dth > 0, 1, 0))

    ## Look at the five first couples
    brvSplit |>
        select(couple, id, sex, doe, dosp, dox, fail, y_before_sp_dth,
               y_after_sp_dth, t_sp_at_risk) |>
        filter(couple<=5) |>
        arrange(couple, id)

    ##    couple  id sex        doe       dosp        dox fail y_before_sp_dth
    ## 1       1 197   1 1981-01-15 2000-01-01 1983-04-18    1     -18.9601358
    ## 2       2 177   1 1981-04-27 1984-11-20 1983-01-24    1      -3.5675172
    ## 3       2 270   2 1981-04-27 1983-01-24 1984-11-20    0      -1.7440587
    ## 4       2 270   2 1981-04-27 1983-01-24 1984-11-20    1       0.0000000
    ## 5       3 168   1 1981-01-20 1981-12-31 1981-08-03    1      -0.9445844
    ## 6       3 384   2 1981-01-20 1981-08-03 1981-12-31    0      -0.5338955
    ## 7       3 384   2 1981-01-20 1981-08-03 1981-12-31    1       0.0000000
    ## 8       4  12   1 1981-01-20 1988-11-23 1991-01-01    0      -7.8414193
    ## 9       4  12   1 1981-01-20 1988-11-23 1991-01-01    0       0.0000000
    ## 10      4 300   2 1981-01-20 2000-01-01 1988-11-23    1     -18.9464462
    ## 11      5 212   1 1981-05-07 1984-02-04 1984-02-04    1      -2.7461395
    ##    y_after_sp_dth t_sp_at_risk
    ## 1     -16.7068229    2.2533129
    ## 2      -1.8234585    1.7440587
    ## 3       0.0000000    1.7440587
    ## 4       1.8234585    1.8234585
    ## 5      -0.4106889    0.5338955
    ## 6       0.0000000    0.5338955
    ## 7       0.4106889    0.4106889
    ## 8       0.0000000    7.8414193
    ## 9       2.1054649    2.1054649
    ## 10    -11.1050268    7.8414193
    ## 11      0.0000000    2.7461395

## (d)

Now find the (crude) effect of bereavement.

    poisson22d <- glm(fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson, data=brvSplit)
    summary(poisson22d)

    ##
    ## Call:
    ## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
    ##     data = brvSplit)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept) -2.20379    0.07125 -30.932   <2e-16 ***
    ## brv          0.11970    0.13199   0.907    0.364
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 704.42  on 554  degrees of freedom
    ## Residual deviance: 703.61  on 553  degrees of freedom
    ## AIC: 1263.6
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22d)

    ##             exp(beta)      2.5 %   97.5 %
    ## (Intercept) 0.1103837 0.09599736 0.126926
    ## brv         1.1271537 0.87022603 1.459937

## (e)

    ## Poisson regression for sex==1
    poisson22e1 <- glm( fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson,
                       data=filter(brvSplit, sex==1))
    summary(poisson22e1)

    ##
    ## Call:
    ## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
    ##     data = filter(brvSplit, sex == 1))
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept) -2.00434    0.08248 -24.301   <2e-16 ***
    ## brv          0.01080    0.19030   0.057    0.955
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 341.95  on 294  degrees of freedom
    ## Residual deviance: 341.95  on 293  degrees of freedom
    ## AIC: 707.95
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22e1)

    ##             exp(beta)     2.5 %   97.5 %
    ## (Intercept) 0.1347495 0.1146361 0.158392
    ## brv         1.0108630 0.6961580 1.467834

    ## Poisson regression for sex==2
    poisson22e2 <- glm( fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson,
                       data=filter(brvSplit, sex==2))
    summary(poisson22e2)

    ##
    ## Call:
    ## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
    ##     data = filter(brvSplit, sex == 2))
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)  -2.6301     0.1414 -18.598   <2e-16 ***
    ## brv           0.4853     0.2032   2.389   0.0169 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 350.83  on 259  degrees of freedom
    ## Residual deviance: 345.21  on 258  degrees of freedom
    ## AIC: 543.21
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22e2)

    ##              exp(beta)      2.5 %     97.5 %
    ## (Intercept) 0.07206987 0.05462311 0.09508916
    ## brv         1.62461324 1.09098341 2.41925600

    ## Poisson regression, interaction with sex
    brvSplit2 <- mutate(brvSplit,
                        sex = as.factor(sex),
                        brv = as.factor(brv))
    poisson22e3 <- glm(fail ~ sex + brv:sex + offset(log(t_sp_at_risk)),
                       family=poisson, data=brvSplit2)
    summary(poisson22e3)

    ##
    ## Call:
    ## glm(formula = fail ~ sex + brv:sex + offset(log(t_sp_at_risk)),
    ##     family = poisson, data = brvSplit2)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept) -2.00434    0.08248 -24.301  < 2e-16 ***
    ## sex2        -0.62578    0.16371  -3.822 0.000132 ***
    ## sex1:brv1    0.01080    0.19030   0.057 0.954724
    ## sex2:brv1    0.48527    0.20316   2.389 0.016913 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 704.42  on 554  degrees of freedom
    ## Residual deviance: 687.16  on 551  degrees of freedom
    ## AIC: 1251.2
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22e3)

    ##             exp(beta)     2.5 %    97.5 %
    ## (Intercept) 0.1347495 0.1146361 0.1583920
    ## sex2        0.5348431 0.3880363 0.7371919
    ## sex1:brv1   1.0108630 0.6961580 1.4678335
    ## sex2:brv1   1.6246132 1.0909834 2.4192560

## (f) Controlling for age.

    ## Translate time scale from years from spouse death to ages
    brvSplit3 <- brvSplit2 |>
        mutate(age_sp_dth =  as.numeric(dosp - dob) / 365.24, # Age at spouse death
               age_start = age_sp_dth + y_before_sp_dth,      # Age at start of timeband
               age_end = age_sp_dth + y_after_sp_dth)         # Age at end of timeband

    age_cat <- seq(70,100,5) # Split at these ages
    brvSplit4 <- survSplit(brvSplit3, cut=age_cat, start="age_start", end="age_end",
                           event="fail", zero = 0)

    brvSplit4 <- mutate(brvSplit4,
                        t_at_risk = age_end- age_start, # Creating new time at risk
                        age = cut(age_end, age_cat))   # Creating age band category

    ## Calculate crude rates
    survRate(Surv(t_at_risk, fail) ~ age, data=brvSplit4)

    ##                   age       tstop event       rate      lower      upper
    ## age=(75,80]   (75,80]  703.612419    45 0.06395566 0.04664970 0.08557771
    ## age=(80,85]   (80,85] 1184.684043   123 0.10382515 0.08628885 0.12387811
    ## age=(85,90]   (85,90]  490.021356    95 0.19386910 0.15685168 0.23699492
    ## age=(90,95]   (90,95]   55.090352    12 0.21782399 0.11255283 0.38049467
    ## age=(95,100] (95,100]    2.299858     3 1.30442857 0.26900453 3.81209383

    poisson22f1 <- glm(fail ~ brv + age + offset(log(t_at_risk)), family=poisson,
                       data=brvSplit4)
    summary(poisson22f1)

    ##
    ## Call:
    ## glm(formula = fail ~ brv + age + offset(log(t_at_risk)), family = poisson,
    ##     data = brvSplit4)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)  -2.7359     0.1495 -18.299  < 2e-16 ***
    ## brv1         -0.1515     0.1371  -1.105  0.26928
    ## age(80,85]    0.5106     0.1757   2.907  0.00365 **
    ## age(85,90]    1.1627     0.1869   6.220 4.98e-10 ***
    ## age(90,95]    1.2847     0.3290   3.905 9.42e-05 ***
    ## age(95,100]   3.0431     0.5967   5.100 3.40e-07 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 1131.0  on 1035  degrees of freedom
    ## Residual deviance: 1074.4  on 1030  degrees of freedom
    ## AIC: 1642.4
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22f1)

    ##               exp(beta)      2.5 %     97.5 %
    ## (Intercept)  0.06483296 0.04836526  0.0869077
    ## brv1         0.85941225 0.65684485  1.1244503
    ## age(80,85]   1.66632969 1.18097507  2.3511543
    ## age(85,90]   3.19848093 2.21730933  4.6138264
    ## age(90,95]   3.61371336 1.89632619  6.8864335
    ## age(95,100] 20.97060543 6.51145075 67.5373752

    poisson22f2 <- glm(fail ~ brv +age + sex + offset(log(t_at_risk)),
                       family=poisson, data=brvSplit4)
    summary(poisson22f2)

    ##
    ## Call:
    ## glm(formula = fail ~ brv + age + sex + offset(log(t_at_risk)),
    ##     family = poisson, data = brvSplit4)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept) -2.57737    0.15364 -16.776  < 2e-16 ***
    ## brv1        -0.02676    0.14019  -0.191 0.848603
    ## age(80,85]   0.51641    0.17567   2.940 0.003286 **
    ## age(85,90]   1.15434    0.18626   6.197 5.74e-10 ***
    ## age(90,95]   1.29672    0.32900   3.941 8.10e-05 ***
    ## age(95,100]  3.32531    0.60227   5.521 3.37e-08 ***
    ## sex2        -0.49188    0.13054  -3.768 0.000165 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 1131.0  on 1035  degrees of freedom
    ## Residual deviance: 1059.6  on 1029  degrees of freedom
    ## AIC: 1629.6
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22f2)

    ##               exp(beta)      2.5 %     97.5 %
    ## (Intercept)  0.07597334 0.05621929  0.1026685
    ## brv1         0.97359228 0.73968454  1.2814678
    ## age(80,85]   1.67599651 1.18779672  2.3648527
    ## age(85,90]   3.17193748 2.20178944  4.5695502
    ## age(90,95]   3.65729031 1.91917581  6.9695399
    ## age(95,100] 27.80767040 8.54103516 90.5354583
    ## sex2         0.61147403 0.47343733  0.7897571

## (g)

    poisson22g <- glm(fail ~ age + sex + brv:sex + offset(log(t_at_risk)),
                      family=poisson, data=brvSplit4)
    summary(poisson22g)

    ##
    ## Call:
    ## glm(formula = fail ~ age + sex + brv:sex + offset(log(t_at_risk)),
    ##     family = poisson, data = brvSplit4)
    ##
    ## Coefficients:
    ##             Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)  -2.5398     0.1554 -16.343  < 2e-16 ***
    ## age(80,85]    0.5176     0.1757   2.946 0.003221 **
    ## age(85,90]    1.1410     0.1866   6.113 9.76e-10 ***
    ## age(90,95]    1.2962     0.3291   3.939 8.19e-05 ***
    ## age(95,100]   3.3586     0.6031   5.568 2.57e-08 ***
    ## sex2         -0.6221     0.1656  -3.756 0.000172 ***
    ## sex1:brv1    -0.1940     0.1925  -1.008 0.313628
    ## sex2:brv1     0.1823     0.2085   0.874 0.381981
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for poisson family taken to be 1)
    ##
    ##     Null deviance: 1131.0  on 1035  degrees of freedom
    ## Residual deviance: 1057.8  on 1028  degrees of freedom
    ## AIC: 1629.8
    ##
    ## Number of Fisher Scoring iterations: 6

    eform(poisson22g)

    ##               exp(beta)      2.5 %     97.5 %
    ## (Intercept)  0.07888104 0.05816888  0.1069682
    ## age(80,85]   1.67794302 1.18911789  2.3677154
    ## age(85,90]   3.12991468 2.17100707  4.5123602
    ## age(90,95]   3.65549666 1.91790712  6.9673112
    ## age(95,100] 28.74862893 8.81493354 93.7594891
    ## sex2         0.53681350 0.38801481  0.7426746
    ## sex1:brv1    0.82368705 0.56482095  1.2011954
    ## sex2:brv1    1.19991663 0.79745115  1.8055023

## (h)

We could split the post bereavement period into multiple categories (e.g., within one year and subsequent to one year following bereavement) and compare the risks between these categories.

## (i) Analysis using Cox regression.

Cox regression: effect of `brv` controlled for attained age:

    summary(coxph(Surv(age_start, age_end, fail) ~ brv,
                  data = brvSplit4))

    ## Call:
    ## coxph(formula = Surv(age_start, age_end, fail) ~ brv, data = brvSplit4)
    ##
    ##   n= 1036, number of events= 278
    ##
    ##         coef exp(coef) se(coef)      z Pr(>|z|)
    ## brv1 -0.2070    0.8131   0.1390 -1.488    0.137
    ##
    ##      exp(coef) exp(-coef) lower .95 upper .95
    ## brv1    0.8131       1.23    0.6191     1.068
    ##
    ## Concordance= 0.511  (se = 0.014 )
    ## Likelihood ratio test= 2.26  on 1 df,   p=0.1
    ## Wald test            = 2.22  on 1 df,   p=0.1
    ## Score (logrank) test = 2.22  on 1 df,   p=0.1

    summary(coxph(Surv(age_start, age_end, fail) ~ brv + sex,
                  data = brvSplit4))

    ## Call:
    ## coxph(formula = Surv(age_start, age_end, fail) ~ brv + sex, data = brvSplit4)
    ##
    ##   n= 1036, number of events= 278
    ##
    ##          coef exp(coef) se(coef)      z Pr(>|z|)
    ## brv1 -0.07842   0.92458  0.14245 -0.551 0.581971
    ## sex2 -0.47291   0.62318  0.13075 -3.617 0.000298 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##      exp(coef) exp(-coef) lower .95 upper .95
    ## brv1    0.9246      1.082    0.6993    1.2224
    ## sex2    0.6232      1.605    0.4823    0.8052
    ##
    ## Concordance= 0.56  (se = 0.018 )
    ## Likelihood ratio test= 15.85  on 2 df,   p=4e-04
    ## Wald test            = 15.21  on 2 df,   p=5e-04
    ## Score (logrank) test = 15.51  on 2 df,   p=4e-04

## (j)

Cox regression estimating the effect of `brv` for each sex, controlling for age:

    summary(coxph(Surv(age_start, age_end, fail) ~ sex + sex:brv,
                  data = brvSplit4))

    ## Call:
    ## coxph(formula = Surv(age_start, age_end, fail) ~ sex + sex:brv,
    ##     data = brvSplit4)
    ##
    ##   n= 1036, number of events= 278
    ##
    ##               coef exp(coef) se(coef)      z Pr(>|z|)
    ## sex2      -0.58156   0.55902  0.16556 -3.513 0.000444 ***
    ## sex1:brv1 -0.21678   0.80511  0.19303 -1.123 0.261415
    ## sex2:brv1  0.09791   1.10286  0.21191  0.462 0.644065
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##           exp(coef) exp(-coef) lower .95 upper .95
    ## sex2         0.5590     1.7888    0.4041    0.7733
    ## sex1:brv1    0.8051     1.2421    0.5515    1.1753
    ## sex2:brv1    1.1029     0.9067    0.7280    1.6707
    ##
    ## Concordance= 0.568  (se = 0.017 )
    ## Likelihood ratio test= 17.12  on 3 df,   p=7e-04
    ## Wald test            = 16.55  on 3 df,   p=9e-04
    ## Score (logrank) test = 16.92  on 3 df,   p=7e-04