Biostatistics III in R

Exercise 6. Diet data: tabulating incidence rates and modelling with Poisson regression

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)
library(dplyr)    # for data manipulation

Load the diet data using time-on-study as the timescale.

head(diet)
summary(diet)

(a)

Tabulate CHD incidence rates per 1000 person-years for each category of hieng. Calculate (by hand) the ratio of the two incidence rates.

diet <- biostat3::diet
diet$y1k <- diet$y/1000

diet.ir6a <- survRate(Surv(y/1000,chd) ~ hieng, data=diet)
## or
diet %>%
    group_by(hieng) %>%
    summarise(Event = sum(chd), Time = sum(y1k), Rate = Event/Time,      # group sums
              CI_low = poisson.test(Event,Time)$conf.int[1],
              CI_high = poisson.test(Event,Time)$conf.int[2]) 

## IRR
with(diet.ir6a, poisson.test(event,tstop)) 
with(diet.ir6a, poisson.test(rev(event),rev(tstop)))

(b)

Fit a poisson model to find the incidence rate ratio for the high energy group compared to the low energy group and compare the estimate to the one you obtained in the previous question.

poisson6b <- glm( chd ~ hieng + offset( log( y1k ) ), family=poisson, data=diet)
summary(poisson6b)
eform(poisson6b)
eform(poisson6b, method="Profile")

## Waiting for profiling to be done...

Also write out the regression equation, defining your notation.

(c)

Grouping the values of total energy into just two groups does not tell us much about how the CHD rate changes with total energy. It is a useful exploratory device, but to look more closely we need to group the total energy into perhaps 3 or 4 groups.

In this example we shall use the cut points 1500, 2500, 3000, 4500. Compare with a normal distribution to check if these cutpoints seem reasonable.

hist6c <- hist(diet$energy, breaks=25, probability=TRUE, xlab="Energy (units)")
curve(dnorm(x, mean=mean(diet$energy), sd=sd(diet$energy)), col = "red", add=TRUE)

quantile(diet$energy, probs=c(0.01,0.05,0.1,0.25,0.5,0.75,0.90,0.95,0.99))
# For kurtosis and skewness, see package e1071

(d)

Create a new variable eng3 coded 1500 for values of energy in the range 1500–2499, 2500 for values in the range 2500–2999, and 3000 for values in the range 3000–4500.

diet$eng3 <- cut(diet$energy, breaks=c(1500,2500,3000,4500),labels=c("low","medium","high"), 
                 right = FALSE)
cbind(Freq=table(diet$eng3),
      Prop=table(diet$eng3)/nrow(diet))

(e)

Estimate the rates for different levels of eng3. Calculate (by hand) the ratio of rates in the second and third levels to the first level.

## rates and IRRs
diet.ir6e <- survRate(Surv(y/1000,chd) ~ eng3, data=diet)
print(diet.ir6e)

# calculate IRR and confidence intervals
with(diet.ir6e, rate[eng3=="medium"] / rate[eng3=="low"])
with(diet.ir6e[c(2,1),], { # compare second row with first row
  poisson.test(event, tstop)
})
with(diet.ir6e, rate[eng3=="high"] / rate[eng3=="low"])
with(diet.ir6e[c(3,1),], { # compare third row with first row
  poisson.test(event, tstop)
})

(f)

Create your own indicator variables for the three levels of eng3.

diet <- mutate(diet, 
               X1 = as.numeric(eng3 == "low"),
               X2 = as.numeric(eng3 == "medium"),
               X3 = as.numeric(eng3 == "high"))
# or
diet <- biostat3::addIndicators(diet, ~eng3+0) %>%
    mutate(X1 = eng3low, X2 = eng3medium, X3 = eng3high)
colSums(diet[c("X1","X2","X3")])

(g)

Check the indicator variables.

filter(diet, eng3=="low")    %>% select(c(energy,eng3,X1,X2,X3)) %>% head
filter(diet, eng3=="medium") %>% select(c(energy,eng3,X1,X2,X3)) %>% head
filter(diet, eng3=="high")   %>% select(c(energy,eng3,X1,X2,X3)) %>% head

(h)

Use poisson to compare the second and third levels with the first.

poisson6h <- glm( chd ~ X2 + X3 + offset( log( y1k ) ), family=poisson, data=diet )
summary(poisson6h)
eform(poisson6h)

Compare your estimates with those you obtained in part 6e. Write the linear predictor using pen and paper.

Also write out the regression equation, defining your notation.

(i)

Use a poisson model to compare the first and third levels with the second. Write the linear predictor using pen and paper.

poisson6i <- glm( chd ~ X1 + X3 + offset( log( y1k ) ), family=poisson, data=diet )
# or 
poisson6i <- glm( chd ~ I(eng3=="low") + I(eng3=="high") + offset( log( y1k ) ), family=poisson, data=diet )

summary( poisson6i )
eform( poisson6i )

(j)

Repeat the analysis comparing the second and third levels with the first but this time have R create the indicators automatically y using the categorical variable eng3.

poisson6j <- glm( chd ~ eng3 + offset( log( y1k ) ), family=poisson, data=diet )
summary( poisson6j )
eform( poisson6j )

(k)

Calculate the total number of events during follow-up, person-time at risk, and the crude incidence rate (per 1000 person-years)

summarise(diet, rate = sum(chd) / sum(y))