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)
. gllamm y x , i(id) link(logit) fam(bin) adapt
. 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))
