How do you calculate a paired risk ratio and its confidence interval?

Question

I have some data showing the number of times individuals make GP appointments at the week and at the weekend, over a period of about 5 years. The data also show the number of times each individual does not attend (DNA).

I'm trying to calculate a risk ratio for DNA (i.e. is the risk lower for appointments during the week?). One approach I am using is to look at the paired risk ratios, as this will control for individual-level time-invariant confounders (though not other confounders).

The data look like this:

+-----+------------+-----------+------------+-----------+
| ID  | Week: apts | Week: DNA | Wknd: apts | Wknd: DNA |
+-----+------------+-----------+------------+-----------+
|   1 |          1 |         0 |          3 |         1 |
|   2 |          3 |         2 |          2 |         2 |
|   3 |          2 |         1 |          4 |         0 |
| ... |        ... |       ... |        ... |       ... |
+-----+------------+-----------+------------+-----------+

I've generated some data in R and shown how I've attempted to calculate a paired risk ratio.

I calculated the paired risk ratio by logging the individual-level risk ratios, taking the mean, and then exponentiating.

I tried two approaches to calculating confidence intervals - firstly by bootstrapping and secondly using a meta analysis function.

I don't know if any of this is correct. How would you calculate the paired risk ratio and its confidence interval? I have not been able to find any published method. There are methods for paired binary or continuous data, but I can't find anything for this situation.

#----------------------
# generate example data
#----------------------

# we are interested in whether patients are more likely to miss (did not attend / DNA) their appointments at the weekend
# the data show that there is a lower proportion of DNAs during the week
# but this is confounded by an unobserved variable: whether the patient has dependents (e.g. children)
# those with dependents are overrepresented among weekend appointments, and also more likely to DNA

set.seed(13)
d_no_dependents <- data.frame(n.week = rpois(500, 4.5), dna.week = rpois(500, 1.5), n.wknd = rpois(500, 2), dna.wknd = rpois(500, 0.5))
d_dependents <- data.frame(n.week = rpois(250, 1.5), dna.week = rpois(250, 1.5), n.wknd = rpois(250, 8), dna.wknd = rpois(250, 8))
d <- rbind(d_no_dependents, d_dependents)

# adjust impossible cases
d[d$n.week == 0 | d$n.wknd == 0,] <- d[d$n.week == 0 | d$n.wknd == 0,] + 1
d$dna.week <- pmin(d$dna.week, d$n.week)
d$dna.wknd <- pmin(d$dna.wknd, d$n.wknd)

#-------------------------------
# crude risk ratio (unpaired RR)
#-------------------------------

library(epitools) # for function 'riskratio'
crude <- colSums(d)
crude <- matrix(crude, ncol = 2, byrow = T)
crude[,1] <- crude[,1] - crude[,2]
riskratio(crude, rev = 'rows')$measure # crude RR = 0.61 (95% CI 0.58-0.65)

#-----------------------------------------------
# 'real' answer (if we could observe dependents)
#-----------------------------------------------

#reshape data
dr <- d
dr$id <- seq_len(nrow(d))
dr$dependents <- rep(c('no', 'yes'), c(500, 250))
week <- data.frame(id = dr$id, n = dr$n.week, dna = dr$dna.week, dependents = dr$dependents, week = 'week')
wknd <- data.frame(id = dr$id, n = dr$n.wknd, dna = dr$dna.wknd, dependents = dr$dependents, week = 'wknd')
dr <- rbind(week, wknd)
dr$week <- factor(dr$week, c('wknd', 'week'))

# unadjusted model - should give same answer as crude risk ratio
m1 <- glm(dna ~ week + offset(log(n)), dr, family = 'poisson')
exp(cbind(m1$coef, confint(m1))) # RR = 0.61 (0.57-0.66)

# adjusted model
m2 <- glm(dna ~ week + dependents + offset(log(n)), dr, family = 'poisson')
exp(cbind(m2$coef, confint(m2))) # RR = 0.95 (0.88-1.04)

# or with mixed model clustering on patient
# confidence intervals are an approximation for speed
# results are the same

library(lme4)
m3 <- glmer(dna ~ week + dependents + offset(log(n)) + (1 | id), dr, family = 'poisson')
exp(cbind(fixef(m3), confint(m3, method = 'Wald')[-1,])) # RR = 0.95 (0.88-1.04)

#----------------------------------
# paired risk ratio: point estimate
#----------------------------------

dcc <- d
dcc[d$dna.week == 0 | d$dna.wknd == 0,] <- d[d$dna.week == 0 | d$dna.wknd == 0,] + 0.5 # continuity correction - add to all
dcc$RR <- (dcc$dna.week / dcc$n.week) / (dcc$dna.wknd / dcc$n.wknd)
exp(mean(log(dcc$RR))) # mean of individual risk ratios = 0.90

#-------------------------------
# bootstrap confidence intervals
#-------------------------------

N <- nrow(dcc)
B <- 10000 # number of resamples
ind <- sample(seq_len(N), B * N, replace = T) # bootstrap indices
boot.rr <- dcc$RR[ind]
boot.rr <- matrix(boot.rr, ncol = B)
boot.rr <- exp(colMeans(log(boot.rr)))
quantile(boot.rr, probs = c(0.025, 0.5, 0.975)) 
# RR = 0.90 (95% CI 0.84-0.97)

#-----------------------
# meta-analysis approach
#-----------------------

# different result to paired ratio, partly because it accounts for number of admissions per individual
# continuity correction is 0.5 by default

library(meta)
metabin(event.e = dcc$dna.week, n.e = dcc$n.week, event.c = dcc$dna.wknd, n.c = dcc$n.wknd)

# from random effect model: RR = 0.94 (0.89-0.98)

Answer 1

One simple way to proceed would be to do a binomial logistic regression to determine statistical significance and confidence intervals, and then convert the results from the log-odds scale of the logistic regression to probability scales if you wish to present results in terms of probability differences or risk ratios.

With the standard logit link function in binomial regression, you model the log-odds of the event (in your case, a DNA) as a function of predictors. You main predictor is part-of-week, in your example a 2-level categorical variable for which we can take "weekday" as the reference level. Then the intercept of the logistic regression is the log-odds of aDNA on a weekday, and the regression coefficient for part-of-week would be the difference in log-odds between weekend and weekday; equivalently, this regression coefficient is the log of the odds ratio for a DNA between weekend and weekday. If that regression coefficient is significantly different from 0 then the null hypothesis of no difference between weekend and weekday is not supported.

You would take the individuals into account with a random effect in a mixed logistic regression model. If the values available for each patient are already aggregated among the 5 years of data collection then you would be limited to a random effect for the intercept of the logistic regression model, allowing individuals to differ in the log-odds of a DNA on a weekday while modeling a single log-odds-ratio in the part-of-week coefficient.

Many people (including me) have difficulty in thinking about odds and odds ratios, but those values are easily transformed into probability values, differences in probabilities, or risk ratios to make your presentation easier for your audience to understand.

If you use the R glmer() function in the lme4 package for this analysis, you would re-organize your data into a long form, with one row representing data for a single individual and a single part-of-week, labeled by individual and part-of-week in columns. The corresponding appointment observations could be represented by two columns, one for the number of DNAs ("successes" if you model the probability of a DNA) and one for the number of non-DNAs; alternatively, you could provide the proportion of DNAs as the outcome and use the total number of appointments as a weight.

Answer 2

Having worked on this problem for a while, I believe the answer is to use conditional poisson regression.

# unadjusted model - gives "crude" risk ratio
m1 <- glm(dna ~ week + offset(log(n)), dr, family = 'poisson')
exp(cbind(m1$coef, confint(m1))) # RR = 0.61 (0.57-0.66)

# adjusted model - gives "correct" answer (using unobserved variable 'dependents')
m2 <- glm(dna ~ week + dependents + offset(log(n)), dr, family = 'poisson')
exp(cbind(m2$coef, confint(m2))) # RR = 0.95 (0.88-1.04)

# conditioning on individual - gives "correct" answer without using unobserved variable
library(gnm)
m3 <- gnm(dna ~ week + offset(log(n)), data = dr, eliminate = factor(id), family = poisson())
exp(c(m3$coefficients[1], confint(m3))) # RR = 0.94 (0.86-1.03)