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.

  • $\begingroup$ You are using GEE to deal with correlated outcomes. Do you also consider this correlation in the permutation test or can it get lost due to scrambling of the data? $\endgroup$ Commented Jun 27, 2022 at 10:35
  • $\begingroup$ Yes I forgot to mention that. I permute the allocation of individuals to the treatment/control arm in the dataset to mirror the structure of the dataset, so their repeated outcome measurements (each across 7 time periods) remain untouched. $\endgroup$ Commented Jun 28, 2022 at 19:22

1 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.



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) {

### gives beta as function of RR and other parameters z
l_inv = function(RR,z) {

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)
  • $\begingroup$ Much appreciated but I was looking for an approach that included a reference/references as it had been done before and I could more easily defend for the bounty. Also, I can compute a (covariate-adjusted) empirical null distribution via permutations for my estimated parameter, so I still find it hard to believe there is not a reasonably straight-forward and proven way to then convert that skewed and "1 centered" distribution into a p-value! $\endgroup$ Commented Jul 4, 2022 at 15:31
  • $\begingroup$ @JupiterM104 did you want p-values, only for the hypothesis $RR = 1$, or more generally a way to compute p-values for other hypotheses as well? For this first case my solution is easy, you do not need to take into account the skewed distribution of the RR, you can just as well work with the distribution of the underlying $\beta$ instead. I am not sure whether there is a reference. I don't believe there is article stuff about this, but maybe it occurs somewhere in a textbook. $\endgroup$ Commented Jul 4, 2022 at 19:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.