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

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)

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)

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)

• I had the same problem and I found a recent paper on confidence intervals for RR based on matched pairs - they also have a link to a ShinyApp: DelRocco, N., Wang, Y., Wu, D. et al. New Confidence Intervals for Relative Risk of Two Correlated Proportions. Stat Biosci 15, 1–30 (2023). doi.org/10.1007/s12561-022-09345-7
– Jens
Commented Sep 6, 2023 at 14:49
• This is the link to the ShinyApp mentioned in DelRocco et al. 2022: dongyuanwu.shinyapps.io/PairedRR
– Jens
Commented Sep 6, 2023 at 14:50

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.

• Thanks very much for this helpful answer. I think we can fit the mixed model using the reshaped data I've called dr. Are you thinking of a model along the lines of: glmer(dna ~ week + (1|id) + offset(log(n)), data = dr, family = 'poisson'), i.e. random intercepts for the individuals? This gives RR = 0.68. I've been trying to work out if the fixed effect of week is equivalent to the paired risk ratio - partly because it's different to both the 'correct' RR, 0.95, calculated from the simulated data (where we know the unobserved confounder) and also to the average individual-level RRs.
– Dan
Commented Feb 27, 2020 at 13:50
• @Dan I was proposing a logistic regression, for example glmer(cbind(DNAcounts, nonDNAcounts) ~ part_of_week + (1|id), data = dr, family = 'binomial') with the default logit link. It could be that your unobserved confounder has to to with a random slope effect related to the log(odds-ratio) rather than just the random intercept for the log(odds) at the reference level. That would depend on how you modeled your confounder. You might need multiple observations on the same individuals to estimate random slopes.
– EdM
Commented Feb 27, 2020 at 15:08

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)