# 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

## (a)

Load the bereavement data and explore it.

    head(brv)

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

## (b)

Calculate the mortality rate per 1000 years for men and for women, and find the rate ratio comparing women (coded 2) with men (coded 1).

1. What dimension of time did we use as the timescale? Do you think this is a sensible choice?
2. Which gender has the highest mortality? Is this expected?
3. Could age be a potential confounder? Does age at entry differ between males and females? Later we will estimate the rate ratio while controlling for age.

    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)

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

    select(brv, sex, age_entry) %>%
        group_by(sex) %>%
        summarise(meanAgeAtEntry = mean(age_entry))

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

In these data a subject changes exposure status from not bereaved to bereaved when his or her spouse dies. The first stage of the analysis therefore is to partition each follow–up into a record describing the period of follow-up pre–bereavement and (for subjects who were bereaved during the study) the period 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)

## (d)

Now find the (crude) effect of bereavement.

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

## (e)

Since there is a strong possibility that the effect of bereavement is not the same for men as for women, use streg to estimate the effect of bereavement separately for men and women. Do this both by fitting separate models for males and females (e.g. streg brv if sex==1) as well as by using a single model with an interaction term (you may need to create dummy variables). Confirm that the estimates are identical for these two approaches.

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

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

    ## 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)
    eform(poisson22e3)

## (f) Controlling for age.

There is strong confounding by age. Expand the data by 5 year age–bands, and check that the rate is increasing with age. Find the effect of bereavement controlled for age. If you wish to study the distribution of age then it is useful to know that age at entry and exit. Translate your time scale to 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)

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

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

## (g)

Now estimate the effect of bereavement (controlled for age) separately for each sex.

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

## (h)

We have assumed that any effect of bereavement is both immediate and permanent. This is not realistic and we might wish to improve the analysis by further subdividing the post–bereavement follow–up. How might you do this? (you are not expected to actually do it)

## (i) Analysis using Cox regression.

We can also model these data using Cox regression. Provided we use the attained age as the time scale and split the data to obtain separate observations for the bereaved and non-bereaved person-time the following command will estimate the effect of bereavement adjusted for attained age.

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

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

That is, we do not have to split the data by attained age (although we can fit the model to data split by attained age and the results will be the same).

## (j)

Use the Cox model to estimate the effect of bereavement separately for males and females and compare the estimates to those obtained using Poisson regression.

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