# Coverage Probability Classic List Threaded 6 messages Reply | Threaded
Open this post in threaded view
|

## Coverage Probability

 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  2.496141 7.436524 8.209161 4.313587 3.318612 > l  -0.9718608  1.1000713  1.5715373  0.1135158 -0.2700517 With (l; u) 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

 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 >  2.496141 7.436524 8.209161 4.313587 3.318612 > > l >  -0.9718608  1.1000713  1.5715373  0.1135158 -0.2700517 > > With (l; u) 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: Coverage Probability

 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

 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) #   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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: Coverage Probability

 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

 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.