---
title: "Biostatistics III in R"
author:
- Code by Andreas Karlsson and Mark Clements
format:
  html:
    minimal: true
filters:
  - webr
execute:
  echo: true
  message: false
  cache: false
  fig-width: 7
  fig-height: 6
---

# 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) ## 

Load the bereavement data and explore it.

```{webr-r}
#| autorun: true
#| results: "asis"
library(biostat3)
library(knitr)
library(dplyr)
brv |> head() |> kable("html")
table(table(brv$couple)) |> kable("html")
```

Showing data for the first five couples:

```{webr-r}
#| results: "asis"
brv |> select(couple, id, sex, doe, dosp, dox, fail) |>
    filter(couple<=5) |> arrange(couple, id) |> kable("html")
```

## (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).

  i. What dimension of time did we use as the timescale? Do you think this is a sensible choice?
  ii. Which gender has the highest mortality? Is this expected?
  iii. 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.

```{webr-r}
#| autorun: true
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
```

```{webr-r}
## 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)

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

## (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.

```{webr-r}
#| autorun: true
## Splitting at spouse death (year=0)
brvSplit <- mutate(brv,
                   id=NULL,
                   y_before_sp_dth =  as.numeric(doe -dosp) / 365.24,
                   y_after_sp_dth = as.numeric(dox - dosp) / 365.24) |>
    survSplit(cut = 0, end="y_after_sp_dth", start="y_before_sp_dth",
              id="id", event="fail") |>
    mutate(t_sp_at_risk = y_after_sp_dth - y_before_sp_dth,
           brv = ifelse(y_after_sp_dth > 0, 1, 0),
           sex = as.factor(sex),
           brv = as.factor(brv),
           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
```

Look at the first five couples:

```{webr-r}
#| results: "asis"
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) |> kable("html")
```

## (d) ## 
Now find the (crude) effect of bereavement.


```{webr-r}
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 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.

```{webr-r}
poisson22e1 <- glm(fail ~ brv + offset(log(t_sp_at_risk)), family=poisson,
                   data=brvSplit, subset=sex==1)
summary(poisson22e1)
eform(poisson22e1)
```

```{webr-r}
poisson22e2 <- glm(fail ~ brv + offset(log(t_sp_at_risk)), family=poisson,
                   data=brvSplit, subset=sex==2)
summary(poisson22e1)
eform(poisson22e1)
```


```{webr-r}
poisson22e3 <- glm(fail ~ sex + brv:sex + offset(log(t_sp_at_risk)), family=poisson,
                   data=brvSplit)
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.

```{webr-r}
#| autorun: true
age_cat <- seq(70,100,5) # Split at these ages
brvSplit2 <- survSplit(brvSplit, cut=age_cat, start="age_start", end="age_end",
                       event="fail", zero = 0) |>
    mutate(t_at_risk = age_end- age_start, # Creating new time at risk
           age = cut(age_end, age_cat))   # Creating age band category
```

```{webr-r}
## Calculate crude rates
survRate(Surv(t_at_risk, fail) ~ age, data=brvSplit2)
## model for brv and age
poisson22f1 <- glm(fail ~ brv + age + offset(log(t_at_risk)), family=poisson,
                   data=brvSplit2)
summary(poisson22f1)
eform(poisson22f1)
## model for brv, age and sex
poisson22f2 <- glm(fail ~ brv + age + sex + offset(log(t_at_risk)), family=poisson,
                   data=brvSplit2)
summary(poisson22f2)
eform(poisson22f2)
```

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

```{webr-r}
poisson22g <- glm( fail ~ age + sex + brv:sex + offset( log(t_at_risk) ),
                  family=poisson, data=brvSplit2)
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.

```{webr-r}
summary(coxph(Surv(age_start, age_end, fail) ~ brv,
              data = brvSplit2))

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

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.

```{webr-r}
summary(coxph(Surv(age_start, age_end, fail) ~ sex + sex:brv,
              data = brvSplit2))
```
