Calculate a permutation-based p-value for a risk ratio

Question

Summary

How do you calculate a two-sided p-value for a risk ratio/relative risk (obtained from a GEE logistic regression via predicting risks with and without a treatment) based on a permutation test, where the risk ratios are used as the test statistics, given that the resulting null distribution of permuted risk ratios is skewed and centred on 1?

Full details

Based on a randomised study design I am analysing the effect of a binary treatment/control condition on a binary outcome. I want to estimate the effects as risk (or proportion) differences and risk ratios/relative risks and compute the 95% confidence intervals and (as I must) the p-values for those effect estimates. Due to the longitudinal nature of the data I am using a suitable generalised estimating equation (logistic link and binomial errors) to get the point estimates, by computing the model predicted risks for each individual assuming they are exposed and then not exposed, and then either subtracting (for the risk difference) or dividing (for the risk ratio) those values, and then bootstrapping them to get the confidence intervals.

However, although I could use the bootstrap approach to calculate a "corresponding" p-value for the risk difference (by centring the bootstrap distribution and treating it as a null distribution), I don't believe this makes sense for the risk ratio given the skewed distribution of the bootstraps. Therefore, I was thinking about using a permutation approach where I do the standard permutation process: permute the group allocations and calculate the risk ratio as above, repeat this many times to obtain the null distribution of the risk ratios, and then calculate a two-sided p-value in the usual way (the proportion of absolute permuted values greater than the absolute estimate obtained from the original data).

This works fine for the risk differences as you get a nice symmetrical null distribution of risk differences centred around 0, but for the risk ratios this approach appears to face similar challenges as the resulting null distribution is unsurprisingly skewed and centred on 1. I have struggled to find any obvious solution so far, but possible solutions I have considered so far are:

1. Transform the permuted null distribution for the risk ratio, e.g. using a log transform, to make it symmetrical and centred around 0, then calculate the two-tailed p-value as usual for a permutation p-value. However, this solution would seem to depend on how well/not the transform worked, which is not ideal, and there are surely other problems I'm missing too.

2. Use an alternative test statistic based on the risk ratios, such as one that gives you a null distribution where you can just compute a one-sided p-value and avoid the issues above. However, I have not found any examples of ones that might work.

Thank you very much for any suggestions.

Answer 1

given that the resulting null distribution of permuted risk ratios is skewed and centred on 1?

You can express the risk ratio, the ratio of probabilities for the cases $x_0$ and $x_1$ given the other covariates $\mathbf{z}$, as a function of the parameter $\beta_x$ and $\mathbf{\boldsymbol{\beta}_z}$.

$$RR = \frac{1+\text{exp}( - \beta_x x_1 - \mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z})}{1+\text{exp}(- \beta_x x_0 -\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z})}$$

To compute a p-value you could just as well compute the p-value of the hypothesis $\beta_x = 0$. The hypothesis $RR = 1$ is equivalent to the hypothesis $\beta_x = 0$.

But possibly you want to compute confidence intervals. In this case it is a bit more difficult...

...If the other $\mathbf{\boldsymbol{\beta}_z}$ would have been fixed, then you could convert the confidence intervals for $\beta_x$ directly to intervals for $RR$ because the function of $RR$ is a monotonous function of $\beta_x$. But, here you have a function of $\beta_x$ and also of $\mathbf{\boldsymbol{\beta}_z}$ and the true value of $\mathbf{\boldsymbol{\beta}_z}$ is not known, as they are estimates.

We could visualize this with a 2d plot of the value $\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z}$ and $\beta_x$ and on top of that isolines of the value of $RR$. For a given value of $\beta_x \neq 0$ the $RR$ will also depend on the value of $\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z}$. We can consider this $\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z}$ as a nuisance parameter and compute p-values for alternative hypotheses $RR = a \neq 1$ by considering

• All points $\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z}$ and $\beta_x$ on the iso-line $RR = a$.
• the sample distribution of the estimates of $\mathbf{\boldsymbol{\beta}_z} \cdot \mathbf{z}$ and $\beta_x$, which will be correlated and are approximately normal distributed.

and compute the worst case p-value. Or use some average value, or some Bayesian value.

Alternatively one could estimate a standard error for the distribution of $RR$ by using the Delta method. Then use the Wald interval by using the estimated $RR$ plusminus some constant times the estimated standard error.

set.seed(1)

n = 300
x = rbinom(n,1,0.5)
zv = rnorm(n)
z = cbind(rep(1,n) , zv)

betax = 0.1
betaz = c(1,1)

logistic = function(x) {
   (1+exp(-x))^-1
}

### gives beta as function of RR and other parameters z
l_inv = function(RR,z) {
   log(RR/(1-exp(z)*(RR-1)))
}

experiment = function() {
   p = logistic(x * betax + z %*% betaz)
   y = rbinom(n,1,p)
   mod = glm(y ~ 0 + x + z, family = binomial)
   b = coefficients(mod)
   bx = b[1]
   bz = colMeans(z) %*% b[2:3]
   return (c(bx,bz))
}

b = replicate(10^3, experiment())
plot(t(b), xlab = "beta_x", ylab = "beta_z * z", pch = 21, col = rgb(0,0,0,0.1), bg = rgb(0,0,0,0.1), 
main = "risk ratio as function of beta_x and beta_z*z \n scatter plot of simulations")


for (RR in seq(0.84,1.12,0.04)) {
   t = seq(0.3,1.7,0.01)
   lines(l_inv(RR,t),t)
   text(l_inv(RR,1.7),1.7,RR,pos=3)
}