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