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.
(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)
- The timescale is attained age, which would seem to be a reasonable choice.
- Males have the higher mortality which is to be expected.
- 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.1080471poisson22b <- 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##             exp(beta)     2.5 %    97.5 %
## (Intercept) 0.2058383 0.1446403 0.2929295
## sex         0.6559620 0.5125885 0.8394379## 
##  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##             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##             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##              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##             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.81209383poisson22f1 <- 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##               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.5373752poisson22f2 <- 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##               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##               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:
## 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## 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:
## 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