The concrete problem is that I am refereeing
a paper where a confidence interval is
presented for the risk ratio and I do not find
it credible. I show below my attempts to
do this in R. The example is slightly changed
from the authors'.
I can obtain a confidence interval for
the odds ratio from fisher.test of
course
=== fisher.test example ==
> outcome <- matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE)
> fisher.test(outcome)
Fisher's Exact Test for Count Data
data: outcome
p-value = 0.00761
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.694792 Inf
sample estimates:
odds ratio
Inf
=== end example ==
but in epidemiology authors often
prefer to present risk ratios.
Using the facility on CRAN to search
the site I find packages epitools and Epi
which both offer confidence intervals
for the risk ratio
=== Epi example ==
> library(Epi)
> twoby2(outcome[c(2,1),c(2,1)])
2 by 2 table analysis:
------------------------------------------------------
Outcome : Col 1
Comparing : Row 1 vs. Row 2
Col 1 Col 2 P(Col 1) 95% conf. interval
Row 1 8 500 0.0157 0.0079 0.0312
Row 2 0 500 0.0000 0.0000 NaN
95% conf. interval
Relative Risk: Inf NaN Inf
Sample Odds Ratio: Inf NaN Inf
Conditional MLE Odds Ratio: Inf 1.6948 Inf
Probability difference: 0.0157 0.0027 0.0337
Exact P-value: 0.0076
Asymptotic P-value: NaN
------------------------------------------------------
=== end example ==
So Epi gives me a lower limit of NaN but the same confidence
interval and p-value as fisher.test
=== epitools example ==
> library(epitools)
> riskratio(outcome)
$data
Outcome
Predictor Disease1 Disease2 Total
Exposed1 500 0 500
Exposed2 500 8 508
Total 1000 8 1008
$measure
risk ratio with 95% C.I.
Predictor estimate lower upper
Exposed1 1 NA NA
Exposed2 Inf NaN Inf
$p.value
two-sided
Predictor midp.exact fisher.exact chi.square
Exposed1 NA NA NA
Exposed2 0.00404821 0.007610478 0.004843385
$correction
[1] FALSE
attr(,"method")
[1] "Unconditional MLE & normal approximation (Wald) CI"
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(xx, correct
correction)
=== end example ==
And epitools also gives a lower limit
of NaN.
=== end all examples ==
I would prefer not to have to tell the authors of the
paper I am refereeing that
I think they are wrong unless I can help them with what they
should have done.
Is there another package I should have tried?
Is there some other way of doing this?
Am I doing something fundamentally wrong-headed?
Michael Dewey
http://www.aghmed.fsnet.co.uk
Viechtbauer Wolfgang (STAT)
2006-Nov-10 15:54 UTC
[R] Confidence interval for relative risk
Hello,
A common way to calculate the CI for a relative risk is the following. Given a
2x2 table of the form:
a b
c d
the log of the risk ratio is given by:
lrr = log[ a/(a+b) / ( c/(c+d) ) ]
which is asymptotically normal with variance:
vlrr = 1/a - 1/(a+b) + 1/c - 1/(c+d).
So an approximate 95% CI for the risk ratio is given by:
exp[ lrr - 1.96*sqrt(vlrr) ], exp[ lrr + 1.96*sqrt(vlrr) ].
A common convention is to add 1/2 to each cell when there are zeros.
So, for the table:
Col 1 Col 2
Row 1 8 500
Row 2 0 500
lrr = log[ 8.5/509 / ( 0.5/501 ) ] = 2.817
vllr = 1/8.5 - 1/509 + 1/0.5 - 1/501 = 2.1137
exp[ 2.817-1.96*sqrt(2.1137) ] = .97
exp[ 2.817+1.96*sqrt(2.1137) ] = 289.04
Maybe that is what the authors did.
Best,
--
Wolfgang Viechtbauer
?Department of Methodology and Statistics
?University of Maastricht, The Netherlands
?http://www.wvbauer.com/
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-
> bounces at stat.math.ethz.ch] On Behalf Of Michael Dewey
> Sent: Friday, November 10, 2006 15:43
> To: r-help at stat.math.ethz.ch
> Subject: [R] Confidence interval for relative risk
>
> The concrete problem is that I am refereeing
> a paper where a confidence interval is
> presented for the risk ratio and I do not find
> it credible. I show below my attempts to
> do this in R. The example is slightly changed
> from the authors'.
>
> I can obtain a confidence interval for
> the odds ratio from fisher.test of
> course
>
> === fisher.test example ==>
> > outcome <- matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE)
> > fisher.test(outcome)
>
> Fisher's Exact Test for Count Data
>
> data: outcome
> p-value = 0.00761
> alternative hypothesis: true odds ratio is not equal to 1
> 95 percent confidence interval:
> 1.694792 Inf
> sample estimates:
> odds ratio
> Inf
>
> === end example ==>
> but in epidemiology authors often
> prefer to present risk ratios.
>
> Using the facility on CRAN to search
> the site I find packages epitools and Epi
> which both offer confidence intervals
> for the risk ratio
>
> === Epi example ==>
> > library(Epi)
> > twoby2(outcome[c(2,1),c(2,1)])
> 2 by 2 table analysis:
> ------------------------------------------------------
> Outcome : Col 1
> Comparing : Row 1 vs. Row 2
>
> Col 1 Col 2 P(Col 1) 95% conf. interval
> Row 1 8 500 0.0157 0.0079 0.0312
> Row 2 0 500 0.0000 0.0000 NaN
>
> 95% conf. interval
> Relative Risk: Inf NaN Inf
> Sample Odds Ratio: Inf NaN Inf
> Conditional MLE Odds Ratio: Inf 1.6948 Inf
> Probability difference: 0.0157 0.0027 0.0337
>
> Exact P-value: 0.0076
> Asymptotic P-value: NaN
> ------------------------------------------------------
>
> === end example ==>
> So Epi gives me a lower limit of NaN but the same confidence
> interval and p-value as fisher.test
>
> === epitools example ==>
> > library(epitools)
> > riskratio(outcome)
> $data
> Outcome
> Predictor Disease1 Disease2 Total
> Exposed1 500 0 500
> Exposed2 500 8 508
> Total 1000 8 1008
>
> $measure
> risk ratio with 95% C.I.
> Predictor estimate lower upper
> Exposed1 1 NA NA
> Exposed2 Inf NaN Inf
>
> $p.value
> two-sided
> Predictor midp.exact fisher.exact chi.square
> Exposed1 NA NA NA
> Exposed2 0.00404821 0.007610478 0.004843385
>
> $correction
> [1] FALSE
>
> attr(,"method")
> [1] "Unconditional MLE & normal approximation (Wald) CI"
> Warning message:
> Chi-squared approximation may be incorrect in: chisq.test(xx, correct >
correction)
>
> === end example ==>
> And epitools also gives a lower limit
> of NaN.
>
> === end all examples ==>
> I would prefer not to have to tell the authors of the
> paper I am refereeing that
> I think they are wrong unless I can help them with what they
> should have done.
>
> Is there another package I should have tried?
>
> Is there some other way of doing this?
>
> Am I doing something fundamentally wrong-headed?
>
>
>
> Michael Dewey
> http://www.aghmed.fsnet.co.uk
Michael Dewey <info at aghmed.fsnet.co.uk> wrote>>> Subject: [R] Confidence interval for relative risk >>> >>> The concrete problem is that I am refereeing >>> a paper where a confidence interval is >>> presented for the risk ratio and I do not find >>> it credible. I show below my attempts to >>> do this in R. The example is slightly changed >>> from the authors'. >>> >>> I can obtain a confidence interval for >>> the odds ratio from fisher.test of >>> course >>> >>> === fisher.test example ==>>> >>> > outcome <- matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE) >>> > fisher.test(outcome) >>> >>> Fisher's Exact Test for Count Data >>> >>> data: outcome >>> p-value = 0.00761 >>> alternative hypothesis: true odds ratio is not equal to 1 >>> 95 percent confidence interval: >>> 1.694792 Inf >>> sample estimates: >>> odds ratio >>> Inf >>> >>> === end example ==>>>Since RR = p1 + 1 p1=a/(a+b), p2=c/(c+d) -------------- p2 1 ----- + --- 1-p2 OR you can backsolve for the RR CI.> x <- function(p1, p2, oddsr) p1 + 1/(p2/(1-p2) + 1/oddsr) > 1/x(8/508, 0, 1/1.694792)[1] 1.650734 Another approach is the Cornfield method -- which is a profile likelihood based CI using the Gibbs Chi-square (LRTS), but IIRC originally used the Pearson chi-square. I don't think there are R implementations. David Duffy. -- | David Duffy (MBBS PhD) ,-_|\ | email: davidD at qimr.edu.au ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v
Michael Dewey
2006-Nov-22 16:31 UTC
[R] Summary, was Re: Confidence interval for relative risk
At 14:43 10/11/2006, Michael Dewey wrote: After considerable help from list members and some digging of my own I have prepared a summary of the findings which I have posted (see link below). Broadly there were four suggestions 1 - Wald-type intervals, 2 - transforming the odds ratio confidence interval, 3 - method based on score test, 4 - method based on likelihood. Method 3 was communicated to me off-list ==========================I haven't followed all of the thread either but has someone already pointed out to you confidence intervals that use the score method? For example Agresti (Categorical Data Analysis 2nd edition, p. 77-78) note that 'although computationally more complex, these methods often perform better'. However, they also note that 'currently they are not available in standard software'. But then again, R is not standard software: the code (riskscoreci) can be found from here: http://www.stat.ufl.edu/~aa/cda/R/two_sample/R2/index.html best regards, Jukka Jokinen ===============================================and so I reproduce it here. Almost a candidate for the fortunes package there. The other three can be found from the archive under the same subject although not all in the same thread. Methods 3 and 4 seem to have more going for them as far as I can judge. Thanks to David Duffy, Spencer Graves, Jukka Jokinen, Terry Therneau, and Wolfgan Viechtbauer. The views and calculations presented here and in the summary are my own and any errors are my responsibility not theirs. The summary document is available from here http://www.zen103156.zen.co.uk/rr.pdf Original post follows. ===================================>The concrete problem is that I am refereeing>a paper where a confidence interval is >presented for the risk ratio and I do not find >it credible. I show below my attempts to >do this in R. The example is slightly changed >from the authors'. > >I can obtain a confidence interval for >the odds ratio from fisher.test of >course > >=== fisher.test example ==> > > outcome <- matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE) > > fisher.test(outcome) > > Fisher's Exact Test for Count Data > >data: outcome >p-value = 0.00761 >alternative hypothesis: true odds ratio is not equal to 1 >95 percent confidence interval: > 1.694792 Inf >sample estimates: >odds ratio > Inf > >=== end example ==> >but in epidemiology authors often >prefer to present risk ratios. > >Using the facility on CRAN to search >the site I find packages epitools and Epi >which both offer confidence intervals >for the risk ratio > >=== Epi example ==> > > library(Epi) > > twoby2(outcome[c(2,1),c(2,1)]) >2 by 2 table analysis: >------------------------------------------------------ >Outcome : Col 1 >Comparing : Row 1 vs. Row 2 > > Col 1 Col 2 P(Col 1) 95% conf. interval >Row 1 8 500 0.0157 0.0079 0.0312 >Row 2 0 500 0.0000 0.0000 NaN > > 95% conf. interval > Relative Risk: Inf NaN Inf > Sample Odds Ratio: Inf NaN Inf >Conditional MLE Odds Ratio: Inf 1.6948 Inf > Probability difference: 0.0157 0.0027 0.0337 > > Exact P-value: 0.0076 > Asymptotic P-value: NaN >------------------------------------------------------ > >=== end example ==> >So Epi gives me a lower limit of NaN but the same confidence >interval and p-value as fisher.test > >=== epitools example ==> > > library(epitools) > > riskratio(outcome) >$data > Outcome >Predictor Disease1 Disease2 Total > Exposed1 500 0 500 > Exposed2 500 8 508 > Total 1000 8 1008 > >$measure > risk ratio with 95% C.I. >Predictor estimate lower upper > Exposed1 1 NA NA > Exposed2 Inf NaN Inf > >$p.value > two-sided >Predictor midp.exact fisher.exact chi.square > Exposed1 NA NA NA > Exposed2 0.00404821 0.007610478 0.004843385 > >$correction >[1] FALSE > >attr(,"method") >[1] "Unconditional MLE & normal approximation (Wald) CI" >Warning message: >Chi-squared approximation may be incorrect in: chisq.test(xx, correct >correction) > >=== end example ==> >And epitools also gives a lower limit >of NaN. > >=== end all examples ==> >I would prefer not to have to tell the authors of the >paper I am refereeing that >I think they are wrong unless I can help them with what they >should have done. > >Is there another package I should have tried? > >Is there some other way of doing this? > >Am I doing something fundamentally wrong-headed? > > >Michael Dewey http://www.aghmed.fsnet.co.uk