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]]

