[R] Confidence intervals for relative risk

[Terry Therneau]

Wolfgang,
  It is common to handle relative risk problems using Poisson regression.
  In your example you have 8 events out of 508 tries, and 0/500 in the second
data set.

> tdata <- data.frame(y=c(8,0), n=c(508,500), group=1:0)
> fit <- glm(y ~ group + offset(log(n)), data=tdata, family=poisson)

  Because of the zero, the standard beta/se(beta) confidence intervals don't
work.  In fact, the actual MLE for the ratio is infinity, which the above
glm decides is exp(33.85) = 5* 10^14 --- which is close enough to infinity
for me.  What you want is the lower limit of the interval, which you can
find with a profile likelihood.  That is, look at the curve of deviance vs
beta, draw a horizontal line 3.84 units down from the top (chisq on 1 df),
and where it itersects the curve is your confidence limit.

  For the above, since it is 2 groups with 2 coefficients, the deviance of
the full model happens to be 0.  Your approximation gave a lower limit of
.97 which is about exp(0), so I'll guess a solution between -2 and 2.

> xx <- seq(-2, 2, length=21)
> for (i in 1:21) {
	fit <- glm(y ~ offset(group*xx[i] + log(n)), poisson, tdata)
	print(c(xx[i], fit$deviance))
	}
	
[1]  -2.00000  33.80736
[1]  -1.80000  31.02994
[1]  -1.60000  28.33138
[1]  -1.40000  25.72327
[1]  -1.20000  23.21769
[1]  -1.00000  20.82691
[1]  -0.80000  18.56281
[1]  -0.60000  16.43629
[1]  -0.40000  14.45668
[1]  -0.20000  12.63108
[1]  0.00000 10.96387
[1] 0.20000 9.45639
[1] 0.400000 8.106805
[1] 0.600000 6.910274
[1] 0.800000 5.859303
[1] 1.000000 4.944278
[1] 1.200000 4.154088
[1] 1.400000 3.476757
[1] 1.60000 2.90003
[1] 1.80000 2.41186
[1] 2.000000 2.000785

 So the "exact" lower limit is somewhere between exp(1.4) and exp(1.2), which
is around 3.6.  One can easily refine this by tucking the fit into an
iterative search, or just resetting xx to a new range.

  The offset(log(n)) part of the model, BTW, is a well known trick in poisson
models.  It has to do with the fact that the likelihood is in terms of the
number of events, but we want coefficients in terms of rates, and E(y) = rate*n.
Adding an offset of xx[i] * group fits a model with the group effect fixed at
xx, but allowing the intercept to vary.  
  For more conceptual depth, you could look up 'rate regression' or 
'standardized mortality ratio' in the Encyclopedia of Biostatistics --- it is
in the latter computation that this approach is quite common.  The idea of
adding a fraction is found in Anscombe(1949), Transformations of Poisson, 
binomial and negative-binomial data.  Biometrika, vol 35, p246-254, but for
each single estimate not for the ratio.

	Terry Therneau
	Mayo Clinic