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 |
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. |
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) } |
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. |
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. |
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. |
Free forum by Nabble | Edit this page |