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)
## id couple dob doe dox dosp fail group disab health sex
## 1 14464 93 1901-01-31 1981-03-11 1990-11-10 2000-01-01 1 1 0 2 1
## 2 2000 18 1896-06-27 1981-03-19 1981-10-02 1982-03-12 1 2 3 1 1
## 3 24486 187 1905-06-03 1981-04-13 1991-01-01 1982-03-03 0 2 0 2 1
## 4 24646 195 1895-11-20 1981-04-29 1985-08-17 2000-01-01 1 1 0 2 1
## 5 16810 119 1898-12-26 1981-03-18 1989-08-31 2000-01-01 1 2 0 2 1
## 6 6344 49 1904-10-31 1981-05-22 1989-02-19 1986-10-05 1 1 3 2 1
## age_entry att_age t_at_risk
## 1 80.10897 89.77659 9.6676158
## 2 84.72511 85.26448 0.5393714
## 3 75.86245 85.58208 9.7196364
## 4 85.43971 89.74099 4.3012813
## 5 82.22539 90.68010 8.4547147
## 6 76.55788 84.30621 7.7483299##Look at the five first couples
brv %>% select(couple, id, sex, doe, dosp, dox, fail) %>% filter(couple<=5) %>% arrange(couple, id)## 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6389 -1.3080 0.2636 0.8742 3.1495
##
## 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## # A tibble: 2 × 2
## sex meanAgeAtEntry
## <int> <dbl>
## 1 1 79.1
## 2 2 78.7(c) Breaking records into pre and post bereavement.
## Creating times relativ 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 y_after_sp_dth
## 1 1 197 1 1981-01-15 2000-01-01 1983-04-18 1 -18.9601358 -16.7068229
## 2 2 177 1 1981-04-27 1984-11-20 1983-01-24 1 -3.5675172 -1.8234585
## 3 2 270 2 1981-04-27 1983-01-24 1984-11-20 0 -1.7440587 0.0000000
## 4 2 270 2 1981-04-27 1983-01-24 1984-11-20 1 0.0000000 1.8234585
## 5 3 168 1 1981-01-20 1981-12-31 1981-08-03 1 -0.9445844 -0.4106889
## 6 3 384 2 1981-01-20 1981-08-03 1981-12-31 0 -0.5338955 0.0000000
## 7 3 384 2 1981-01-20 1981-08-03 1981-12-31 1 0.0000000 0.4106889
## 8 4 12 1 1981-01-20 1988-11-23 1991-01-01 0 -7.8414193 0.0000000
## 9 4 12 1 1981-01-20 1988-11-23 1991-01-01 0 0.0000000 2.1054649
## 10 4 300 2 1981-01-20 2000-01-01 1988-11-23 1 -18.9464462 -11.1050268
## 11 5 212 1 1981-05-07 1984-02-04 1984-02-04 1 -2.7461395 0.0000000
## t_sp_at_risk
## 1 2.2533129
## 2 1.7440587
## 3 1.7440587
## 4 1.8234585
## 5 0.5338955
## 6 0.5338955
## 7 0.4106889
## 8 7.8414193
## 9 2.1054649
## 10 7.8414193
## 11 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5442 -1.0377 -0.0657 0.8531 3.4313
##
## 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))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6325 -1.0016 0.1241 0.8594 3.1501
##
## 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))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4851 -0.9557 -0.5097 0.8569 3.4490
##
## 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6325 -0.9659 -0.2355 0.8583 3.4490
##
## 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8226 -0.7876 -0.5159 0.4976 3.3740
##
## 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7766 -0.7973 -0.4988 0.4702 3.3896
##
## 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7246 -0.7874 -0.4931 0.4611 3.3591
##
## 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