inst/doc/q22.R

In biostat3: Utility Functions, Datasets and Extended Examples for Survival Analysis

## Purpose: To do the solution for Biostat III exercises in R
## Author: Andreas Karlsson, 2015-03-02
## Revised: Mark Clements, 2017-11-03
###############################################################################

## Install needed packages only need to be done once
## install.packages("foreign")
## install.packages("survival")
## install.packages("dplyr")


###############################################################################
## Exercise 22
###############################################################################
## @knitr loadDependencies
library(biostat3) # for Surv and survfit
library(dplyr)    # for data manipulation
library(rstpm2)

## @knitr loadPreprocess


## @knitr 22.a

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

## @knitr 22.b
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)

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

## @knitr 22.c

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

## @knitr 22.d
poisson22d <- glm(fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson, data=brvSplit)
summary(poisson22d)
eform(poisson22d)

## @knitr 22.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)
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)

## @knitr 22.f
## 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)

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

## @knitr 22.i
summary(coxph(Surv(age_start, age_end, fail) ~ brv,
              data = brvSplit4))


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

## @knitr 22.j
summary(coxph(Surv(age_start, age_end, fail) ~ sex + sex:brv,
              data = brvSplit4))