  # Re: st: RRR with CI from logit model

 From Michael Ingre To statalist@hsphsun2.harvard.edu Subject Re: st: RRR with CI from logit model Date Wed, 3 Nov 2004 16:49:10 +0100

Dear list,

I would like to thank all of you for your help (none mentioned and none forgotten). Many suggestions and comments have been put forward and it has really helped me with this task. I think I have a solution and the paper I had in mind, seems to be possible to realise. I'm very interested in your comments and maybe they will help me defend my case against the reviewers.

I have tried three approaches. The first is a variant of the simulation approach suggested by Leonelo Bautista (http://www.stata.com/statalist/archive/2004-11/msg00071.html). The second is a variant of the first where the ratios have been log transformed and RRR and CIs are calculated from the exponent of the mean and std dev. The third approach was the one I proposed earlier where the ratio are log transformed before the wald-test with -nlcom-. At the bottom are some comparisons of RRRs with CIs with the three different approaches.

The main impression is that the two first give similar answers which supports the fact that the ratios can be log transformed and then treated as normally distributed. Also, -ladder- suggest that the ratios could be log transformed to have a distribution that does not deviate significantly from the normal distribution with 1000 observations.

The third approach with log transformation and -nlcom- has more narrow CIs with the point estimate right on target. Given the evidence above, that the ratios are log-normally distributed, I vote in favour for this approach.

What would be your comment if you reviewed a paper with these RRRs?

Thanks again for all the help.

Michael

// The code is only to illustrate. It has not been tested and might contain mistakes.
// level_a and level_b describes the levels of the covariates to calculate RRR.

1)

. set obs 10000
. local level_a = 9
. local level_b = 1
. lincom `level_a' * x + _cons
. gen log_odds_a = r(estimate)
. gen random_a = r(se) * invnorm(uniform())
. lincom `level_b' * x + _cons
. gen log_odds_b = r(estimate)
. gen random_b = r(se) * invnorm(uniform())
. gen random_a = log_odds_se * invnorm(uniform())
. gen random_b = log_odds_se * invnorm(uniform())
. gen risk_a = exp(log_odds_a + random_a) / (1 + exp(log_odds_a + random_a))
. gen risk_b = exp(log_odds_b + random_b) / (1 + exp(log_odds_b + random_b))
. gen ratio = risk_a/risk_b
. centile ratio , c(2.5 50 97.5)
. di "RRR :" r(c_2)
. di "lower CI: " r(c_1)
. di "upper CI: " r(c_3)

2)
. gen ratio_ln = ln(ratio)
. su ratio_ln
. di "RRR :" exp(r(mean))
. di "lower CI: " exp(r(mean)-1.96*r(sd))
. di "upper CI: " exp(r(mean)+1.96*r(sd))

3)
. nlcom ln( ///
. (exp(`level_a' * _b[x] + _b[_cons]) / ///
. (1 + exp(`level_a' * _b[x] + _b[_cons]))) / ///
. (exp(`level_b' * _b[x] + _b[_cons]) / ///
. (1 + exp(`level_b' * _b[x] + _b[_cons]))) ///
. )
. di "RRR :" exp(b[1,1])
. di "lower CI: " exp(b[1,1] - sqrt(V[1,1])*1.96)
. di "upper CI: " exp(b[1,1] + sqrt(V[1,1])*1.96)

Some comparisons:

RRR 95% CI
1) 19432 107 - 2704450
2) 18493 111 - 3072209
3) 22186 227 - 2170434

1) 53.45 0.21 - 11882.82
2) 52.96 0.22 - 12860.14
3) 50.70 1.35 - 1901.76

1) 0.226 0.006 - 8.137
2) 0.224 0.006 - 7.798
3) 0.215 0.129 - 0.359

1) 3.852 0.024 - 598.329
2) 3.792 0.024 - 599.047
3) 3.742 0.775 - 18.064

*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/