Going back to the case-control study

Here is the case-control table we made at the beginning.

However, we can compare the odds of smoking among cases and controls.

odds_case <- 90 / 10
odds_control <- 60 / 40
odds_ratio_smoking <- odds_case / odds_control
odds_ratio_smoking

## [1] 6

# Cross-products give the same OR.
(90 * 40) / (10 * 60)

## [1] 6

The odds of smoking among cases are 6 times those among controls. With appropriate case-control sampling, this estimates the disease odds ratio for smokers compared to non-smokers. This does not mean 6 times the risk.

What happens if we recruit 1,000 controls with the same smoking proportion?

odds_control_more <- 600 / 400
odds_case / odds_control_more

## [1] 6

The OR stays the same! This calculation does not remove confounding or selection bias; the way we recruit and measure the subjects still matters.

Question

The maternal/infant HIV transmission example has 13 transmissions among 180 infants in the AZT group and 40 among 183 in the placebo group.

Calculate RD, RR, and OR for AZT compared to placebo. Try it before running the next chunk.

p_azt <- 13 / 180
p_placebo <- 40 / 183

rd_azt <- p_azt - p_placebo
rr_azt <- p_azt / p_placebo
or_azt <- (p_azt / (1 - p_azt)) /
        (p_placebo / (1 - p_placebo))

data.frame(RD = rd_azt,
           RR = rr_azt,
           OR = or_azt)

The risk is about 14.6 percentage points lower in the AZT group. The RR is about 0.33, and the OR is about 0.28. Be careful when choosing between the words risk and odds.

Question

Calculate log(2), log(1 / 2), and log(1). Can you explain the three answers? Here we use positive, finite ratios; log(0) is -Inf.

Different follow-up times

The three-subject PBC example uses follow-up times of 7, 2, and 3 years on the study clock. Only the first subject died during follow-up.

dataset_follow_up <- data.frame(
        subject = 1:3,
        time_years = c(7, 2, 3),
        death = c(1, 0, 0)
)
dataset_follow_up

sum(dataset_follow_up$death) / sum(dataset_follow_up$time_years)

## [1] 0.08333333

The rate is 1 death per 12 person-years, or about 0.083 deaths per person-year. The two censored subjects still contribute their observed follow-up time.

mean(dataset_follow_up$time_years)

## [1] 4

mean(dataset_follow_up$death)

## [1] 0.3333333

These are the mean follow-up time and the observed fraction who died. The first is not the mean time to death. The second ignores the different lengths of follow-up. Censored subjects are event-free up to their last observation; we do not know what happened after that.

Changing time units

The Nepal infant mortality example has 644 deaths and 1,627,725 person-days of follow-up.

deaths_nepal <- 644
person_days_nepal <- 1627725

ir_nepal_day <- deaths_nepal / person_days_nepal
ir_nepal_year <- ir_nepal_day * 365

data.frame(deaths_per_person_day = ir_nepal_day,
           deaths_per_person_year = ir_nepal_year,
           deaths_per_500_person_years = ir_nepal_year * 500)

Here, 365 days is our conversion for one year. Expressing the rate per person-year does not turn the six-month study into a one-year study, or give us a one-year probability of death.

Question

Express the Nepal rate per 1,000 person-days. Do you get the same answer if you convert back from the rate per person-year?

PBC trial

Now, use the PBC totals in the lecture. Keep the direction as DPCA compared to placebo.

dataset_pbc_summary <- data.frame(
        treatment = c("DPCA", "Placebo"),
        deaths = c(65, 60),
        person_years = c(872.5, 842.5)
)
dataset_pbc_summary$incidence_rate <- dataset_pbc_summary$deaths /
        dataset_pbc_summary$person_years
dataset_pbc_summary

irr_dpca_to_placebo <- dataset_pbc_summary$incidence_rate[1] /
        dataset_pbc_summary$incidence_rate[2]
irr_dpca_to_placebo

## [1] 1.046084

From these totals, the IRR is about 1.046: the observed death rate in the DPCA group is about 4.6% higher. This is a descriptive comparison; we have not measured its uncertainty yet.

The hazard ratio mentioned in the ART paper compares instantaneous event rates. It is not generally the same quantity as the overall IRR.

Smoking cessation data

These are the 12 subjects from the lecture. The event is quitting smoking. event = 1 means that the subject quit, and event = 0 means that the observation was censored.

dataset_smoking <- data.frame(
        time = c(2, 3, 6, 8, 9, 10, 15, 16, 18, 24, 27, 30),
        event = c(1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0)
)
dataset_smoking

Here, “survival” means remaining event-free, or not having quit smoking yet. It does not always mean staying alive.

Calculate by hand first

At day 2, all 12 subjects are at risk and one quits. At day 6, only 10 subjects are at risk, because one quit at day 2 and one was censored at day 3.

s_day2 <- (12 - 1) / 12
s_day6 <- s_day2 * (10 - 1) / 10
s_day8 <- s_day6 * (9 - 1) / 9

c(day2 = s_day2, day6 = s_day6, day8 = s_day8)

##      day2      day6      day8 
## 0.9166667 0.8250000 0.7333333

At each event time, multiply the previous estimate by the fraction of the subjects at risk who remain event-free. Censoring changes later risk sets, but it does not make the curve drop at the censoring time.

Now, use R

Surv() puts time and event status together. survfit() estimates the Kaplan-Meier curve. The ~ 1 means that we want one curve for the whole sample.

km_smoking <- survival::survfit(
        survival::Surv(time, event) ~ 1,
        data = dataset_smoking
)

km_summary <- summary(km_smoking)
data.frame(time = km_summary$time,
           at_risk = km_summary$n.risk,
           events = km_summary$n.event,
           survival = km_summary$surv)

Check the first three rows against our calculations above.

plot(km_smoking,
     conf.int = FALSE,
     mark.time = TRUE,
     lwd = 2,
     xlab = "Time after workshop (days)",
     ylab = "Estimated proportion not yet quit",
     main = "Smoking cessation: Kaplan-Meier curve",
     xlim = c(0, 30),
     ylim = c(0, 1))

The curve is a step function. The marks show censoring times. We use both event and censored observations, assuming censoring is unrelated to the future event time within the group being analyzed.

Reading the curve

Let’s look at days 8 and 18.

summary(km_smoking, times = c(8, 18))$surv

## [1] 0.7333333 0.3771429

# First time the estimated event-free proportion reaches 0.5 or lower.
km_smoking$time[which(km_smoking$surv <= 0.5)[1]]

## [1] 18

The estimated median time to quitting is 18 days. At day 16 the estimate is still slightly above 0.5, even though it rounds to 0.50 in the lecture table. Use the unrounded values when finding the median.

The other way to show the same curve

1 - S(t) shows the estimated proportion who have had the event by time t. In this example, it is the estimated proportion who have quit.

plot(km_smoking,
     fun = "event",
     conf.int = FALSE,
     mark.time = TRUE,
     lwd = 2,
     xlab = "Time after workshop (days)",
     ylab = "Estimated proportion who have quit",
     main = "Smoking cessation: 1 - S(t)",
     xlim = c(0, 30),
     ylim = c(0, 1))

Comparing two groups

The survival package includes individual records from the PBC trial. Let’s select the 312 randomized subjects and draw one curve per treatment. Here, status = 2 means death; transplantation and other censored observations are treated as censoring for this example.

dataset_pbc <- subset(survival::pbc, !is.na(trt))
dataset_pbc$time_years <- dataset_pbc$time / 365
dataset_pbc$death <- as.integer(dataset_pbc$status == 2)
dataset_pbc$treatment <- factor(dataset_pbc$trt,
                                levels = c(1, 2),
                                labels = c("DPCA", "Placebo"))

table(dataset_pbc$treatment, dataset_pbc$death)

##          
##            0  1
##   DPCA    93 65
##   Placebo 94 60

Replace ~ 1 with ~ treatment to estimate a curve for each group.

km_pbc <- survival::survfit(
        survival::Surv(time_years, death) ~ treatment,
        data = dataset_pbc
)

plot(km_pbc,
     conf.int = FALSE,
     mark.time = TRUE,
     col = c("blue", "red"),
     lwd = 2,
     xlab = "Follow-up time (years)",
     ylab = "Estimated survival probability",
     main = "PBC trial: Kaplan-Meier curves",
     ylim = c(0, 1))

legend("bottomleft",
       legend = levels(dataset_pbc$treatment),
       col = c("blue", "red"),
       lwd = 2,
       bty = "n")

These curves assume that censored subjects, including those receiving a transplant, would have the same future death experience as subjects remaining under observation in their treatment group.

Question

Why does the smoking curve not drop at day 3? Why can the curve stay above zero even when few subjects are still under observation? Try adding fun = "event" to the PBC plot and change the y-axis label to match.