Hi hubinho,
This starts to look as homework to me so this will be my last try in
helping you.
The general strategy would be along the lines of (1) write a function that
does what you want for a value of theta and (2) sapply() that function to
the vector of theta values you would like to evaluate:
# function
# -- B is the number of tables
foo2 <- function(theta, n1, n2, B = 1000, alpha = 1, z = 1.96){
# 2x2 tables
x1 <- exp(alpha +theta)/ (1+ exp(alpha +theta))
x2 <- exp(alpha)/ (1+ exp(alpha))
n11 <- rbinom(B, n1, x1)
n12 <- n1 - n11
n21 <- rbinom(B, n2, x2)
n22 <- n2 - n21
# upper and lower limit gart interval
gartu <-function(z,d,e, f, g) log(((d+.5)*(g+.5))/((e+.5)*(f+.5)))+
z*sqrt(1/(d+.5)+1/(e+.5)+1/(f+.5)+1/(g+.5))
gartl <-function(z,d,e, f, g) log(((d+.5)*(g+.5))/((e+.5)*(f+.5)))-
z*sqrt(1/(d+.5)+1/(e+.5)+1/(f+.5)+1/(g+.5))
# calculations and results
u <- gartu(z, n11, n22, n12, n21)
l <- gartl(z, n11, n22, n12, n21)
theta >= l & theta <= u # TRUE if theta is in (l, u)
}
# example
# -- B is the number of tables
res <- foo2(theta = 1, n1 = 10, n2 = 10, B = 1000)
res
# coverage
mean(res)
# different values of theta
Theta <- 1:30
colMeans(sapply(Theta, foo2, n1 = 10, n2 = 10, B = 1000))
HTH,
Jorge.-
On Mon, Mar 19, 2012 at 4:24 PM, hubinho <> wrote:
> Thank you very much again. But in this case I get the coverage probability
> as
> an average over all values for the odds ratio.
>
> I need a coverage probability for every value for the odds ratio.
>
> So the coverage probability for odds ratio = 1, than for odds ratio = 2 and
> so on.
>
> Sorry to bother you again but I have some problems with loops.
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Coverage-Probability-tp4485511p4486264.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
