Coverage Probability

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Coverage Probability

hubinho
Hello.

I'm allready this far. I have a function which is calculating the lower (l) and upper (u) limit for a confidence interval for the odds ratio.

For example for 5 simulated 2x2 tables the upper and lower limits are:

> u
[1] 2.496141 7.436524 8.209161 4.313587 3.318612
> l
[1] -0.9718608  1.1000713  1.5715373  0.1135158 -0.2700517

With (l[1]; u[1]) being the confidence interval for the odds ratio for the first simulated table and so on.

Now I want to compute the coverage probability. For that I've created a function which is return 1 if the odds ratio is in the interval and 0 if it isn't.

cover <- function(theta, u, l){
if(theta >= l && theta   <= u){z=1}
if(theta < l || theta   > u){z=0}; return(z)
}

This works but unfortunately not if I want to summarize the function and divide it with the sample size to get the coverage probability.

I tried it this way

for(for(x in 1:5) {a = (sum(cover(theta, u[x], l[x]))/5; return(a)}

Maybe someone can help me. Thank you

Reply | Threaded
Open this post in threaded view
|

Re: Coverage Probability

Jorge I Velez
Hi hubinho,

You are almost there.  Try this slightly modification of your function:

# theta, u and l are vectors of the same length
foo <- function(theta, u, l) mean(theta >= l & theta <= u, na.rm = TRUE)
foo(theta, u, l)

HTH,
Jorge.-


On Mon, Mar 19, 2012 at 12:55 PM, hubinho <> wrote:

> Hello.
>
> I'm allready this far. I have a function which is calculating the lower (l)
> and upper (u) limit for a confidence interval for the odds ratio.
>
> For example for 5 simulated 2x2 tables the upper and lower limits are:
>
> > u
> [1] 2.496141 7.436524 8.209161 4.313587 3.318612
> > l
> [1] -0.9718608  1.1000713  1.5715373  0.1135158 -0.2700517
>
> With (l[1]; u[1]) being the confidence interval for the odds ratio for the
> first simulated table and so on.
>
> Now I want to compute the coverage probability. For that I've created a
> function which is return 1 if the odds ratio is in the interval and 0 if it
> isn't.
>
> cover <- function(theta, u, l){
> if(theta >= l && theta   <= u){z=1}
> if(theta < l || theta   > u){z=0}; return(z)
> }
>
> This works but unfortunately not if I want to summarize the function and
> divide it with the sample size to get the coverage probability.
>
> I tried it this way
>
> for(for(x in 1:5) {a = (sum(cover(theta, u[x], l[x]))/5; return(a)}
>
> Maybe someone can help me. Thank you
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Coverage-Probability-tp4485511p4485511.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] 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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage Probability

hubinho
Thank you very much. This was, was i needed. Unfortunately I have one futher problem with this Code. I don't only need the coverage probability for one but for a range of different odds ratios. (for example [1;30]). I tried it with a loop but I get an error. I think again, that I'm almost there but having a little mistake. The complete code is:

#setting values

n1 <- 10
n2 <- 10
y <- 100
alpha <- 1
z<-1.96

# creating 2x2 table

for (i in 1:30)

{

theta <- i
x1 <- exp(alpha +theta)/ (1+  exp(alpha +theta))
x2 <- exp(alpha)/ (1+  exp(alpha))


n11 <- rbinom(y, 10, x1)
n12 <- n1 - n11
n21 <- rbinom(y, 10, 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))}


u <- gartu(z, n11[i],n22[i],n12[i],n21[i])
l <- gartl(z, n11[i],n22[i],n12[i],n21[i])

foo <- function(theta, u, l) mean(theta >= l & theta <= u, na.rm = TRUE)
 foo(theta, u, l)
}
Reply | Threaded
Open this post in threaded view
|

Re: Coverage Probability

Jorge I Velez
Hi hubinho,

You need to initialize the for() loop and then store the u and l values
properly:

# parameters
n1 <- 10
n2 <- 10
y <- 100
alpha <- 1
z<-1.96

# creating B 2x2 tables
B <- 50
u <- l <- vector('numeric', B)
for (i in 1:B){
    theta <- i
    x1 <- exp(alpha +theta)/ (1+  exp(alpha +theta))
    x2 <- exp(alpha)/ (1+  exp(alpha))

    n11 <- rbinom(y, 10, x1)
    n12 <- n1 - n11
    n21 <- rbinom(y, 10, 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))}

# store results
u[i] <- gartu(z, n11[i],n22[i],n12[i],n21[i])
l[i] <- gartl(z, n11[i],n22[i],n12[i],n21[i])
}

# coverage
theta <- 1:B
foo <- function(theta, u, l) mean(theta >= l & theta <= u, na.rm = TRUE)
foo(theta, u, l)
#  [1] 0.14

HTH,
Jorge.-


On Mon, Mar 19, 2012 at 2:25 PM, hubinho <> wrote:

> Thank you very much. This was, was i needed. Unfortunately I have one
> futher
> problem with this Code. I don't only need the coverage probability for one
> but for a range of different odds ratios. (for example [1;30]). I tried it
> with a loop but I get an error. I think again, that I'm almost there but
> having a little mistake. The complete code is:
>
> #setting values
>
> n1 <- 10
> n2 <- 10
> y <- 100
> alpha <- 1
> z<-1.96
>
> # creating 2x2 table
>
> for (i in 1:30)
>
> {
>
> theta <- i
> x1 <- exp(alpha +theta)/ (1+  exp(alpha +theta))
> x2 <- exp(alpha)/ (1+  exp(alpha))
>
>
> n11 <- rbinom(y, 10, x1)
> n12 <- n1 - n11
> n21 <- rbinom(y, 10, 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))}
>
>
> u <- gartu(z, n11[i],n22[i],n12[i],n21[i])
> l <- gartl(z, n11[i],n22[i],n12[i],n21[i])
>
> foo <- function(theta, u, l) mean(theta >= l & theta <= u, na.rm = TRUE)
>  foo(theta, u, l)
> }
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Coverage-Probability-tp4485511p4485865.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] 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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage Probability

hubinho
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.
Reply | Threaded
Open this post in threaded view
|

Re: Coverage Probability

Jorge I Velez
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.
>
> ______________________________________________
> [hidden email] 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]]

______________________________________________
[hidden email] 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.