Re: [R] Coverage Probability
Muhammad Kashif uaf.edu.pk> writes: > > Dears > Can anyone help me to solve the issue. > > By using" boot" and "boot.ci" package in R we can construct bootstrap confidence intervals. How we > calculate the coverage probability of these intervals. Calculating coverage probability for any but the simplest cases requires simulations. You need to simulate a particular scenario; for each simulation, use your estimation and confidence-interval calculation procedure; and then score the particular simulation as 1 (estimated CI included the true parameter value) or 0 (true parameter value fell outside the estimated CI). Across a large number of simulations, the proportion of ones is an estimate of the coverage. This is discussed more here: http://ms.mcmaster.ca/~bolker/emdbook/chap5A.pdf __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
[R] Coverage Probability
Dears Can anyone help me to solve the issue. By using" boot" and "boot.ci" package in R we can construct bootstrap confidence intervals. How we calculate the coverage probability of these intervals. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, Some suggestion about the arguments of the function defined below. Since theta is calculated with the value of lambda1 and lambda2, there is no need to include theta in the argument. Or, your function can be defined as function(lambda1, lambda2, significance.level) cover <- function(theta, lambda1, lambda2, significance.level) { s1 <- rpois(1,lambda1) s2 <- rpois(1,lambda2) theta <- lambda2/(lambda1+lambda2) s <- s1+s2 z <- qnorm(1-0.05/2) k <- z^2 pi <- s2/s root <- ((s2/s)*(1-s2/s)+k/(4*s))^(1/2) low <- (s2+k/2)/(s+k)-((z*sqrt(s))/(s+k))*root hig <- (s2+k/2)/(s+k)+((z*sqrt(s))/(s+k))*root if (theta >= low & theta <= hig){1} else {0} } -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703248.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, Given the function cover, it's very likely that you will get 0 for both s1 and s1 with small value of lambda1 and lambda2. In that case the sum s will be 0. With s being 0, you will have issue with the expression in pi <- s2/s and root <- ((s2/s)*(1-s2/s)+k/(4*s))^(1/2). You need to take care of the case that s is 0 before proceeding calculating pi and root. cover <- function(theta, lambda1, lambda2, significance.level) { s1 <- rpois(1,lambda1) s2 <- rpois(1,lambda2) theta <- lambda2/(lambda1+lambda2) s <- s1+s2 z <- qnorm(1-0.05/2) k <- z^2 pi <- s2/s root <- ((s2/s)*(1-s2/s)+k/(4*s))^(1/2) low <- (s2+k/2)/(s+k)-((z*sqrt(s))/(s+k))*root hig <- (s2+k/2)/(s+k)+((z*sqrt(s))/(s+k))*root if (theta >= low & theta <= hig){1} else {0} } -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703238.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, In your function cover, lambda1 and lambda2 are used but not in the argument of the function. I suppose that you need to have lambda1 and lambda2 in the argument of the function cover, like function(lambda1, lambda2, n, significance.level). Give it a try. cover <- function(lambda, n, significance.level) { s1 <- rpois(1,lambda1) s2 <- rpois(1,lambda2) theta <- lambda2/(lambda1+lambda2) s <- s1+s2 z <- qnorm(1-0.05/2) k <- z^2 pi <- s2/s lower <- (pi+(k/(2*s))-z*sqrt((pi*(1-pi)+(k/4*s))/s))/(1+k/s) upper <- (pi+(k/(2*s))+z*sqrt((pi*(1-pi)+(k/4*s))/s))/(1+k/s) if (theta >= lower & theta <= upper){1} else {0} } -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4703230.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, After some thought, I found the treatment of sample mean equal 0 was not appropriate. I modified the function likelihood.ratio.test.Poisson. resulting.matrix now has 0.0512 as the average of type I error. function(lambda, sample.size, significance.level) { reject <- 0 sample.mean <- mean(rpois(sample.size, lambda)) if (sample.mean == 0) { test.statistics <- 2 * sample.size * lambda } else { test.statistics <- 2 * sample.size * (lambda - sample.mean + sample.mean * log(sample.mean / lambda)) } if (test.statistics >= qchisq(1 - significance.level, 1)) {reject <- 1} else {reject <- 0} return(reject) } > for (i in 1:500){ + resulting.matrix[i,1] <- 0.01 * i + resulting.matrix[i,2] <- mean(sapply(1:100,function(x) likelihood.ratio.test.Poisson(0.01*i,10,0.05))) + } > mean(resulting.matrix[,2]) [1] 0.05102 -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702909.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, My first question is what is your n when you say fixed n. I assume the lambda is the mean of the poisson distribution that you want to take sample from. Another question is about the sample size. It does not make too much sense to make a sample of size 1. Let's assume that you want to fix the sample size to be 100 and change lambda from 0.1 to 5 with an increment of 0.1. For each lambda, plan to run, say, 1000 times. Then the following will be my approach. Recall that the function cover returns 1 when lambda is in the confidence interval and 0 otherwise. resulting_matrix is created with size 50 x 2 with 0 populated. The matrix is to store lambda and the proportion of samples with lambda inside the confidence interval calculated from samples. With the resulting matrix, one can see that lambdas are in the first column with values of 0.1 to 5 with increment of 0.1. The corresponding proportions are in the second column. All of the proportions are from 0.917 to 0.969 as the last line shows. Hope this helps. > cover <- function(lambda, sample.size, significance.level) { + x <- rpois(sample.size,lambda) + estimate <- mean(x) + lower <- estimate - qnorm(1 - significance.level/2) * sqrt(estimate/sample.size) + upper <- estimate + qnorm(1 - significance.level/2) * sqrt(estimate/sample.size) + if (lambda > lower & lambda < upper){1}else{0} + } > resulting.matrix <- matrix(0, nrow=50,ncol=2) > for (i in 1:50) + { + resulting.matrix[i,1] <- 0.1 * i + resulting.matrix[i,2] <- mean(sapply(1:1000,function(x) cover(0.1*i,100,0.05))) + } > resulting.matrix [,1] [,2] [1,] 0.1 0.917 [2,] 0.2 0.949 [3,] 0.3 0.928 [4,] 0.4 0.939 [5,] 0.5 0.943 [6,] 0.6 0.949 [7,] 0.7 0.942 [8,] 0.8 0.939 [9,] 0.9 0.945 [10,] 1.0 0.943 [11,] 1.1 0.962 [12,] 1.2 0.933 [13,] 1.3 0.947 [14,] 1.4 0.951 [15,] 1.5 0.946 [16,] 1.6 0.939 [17,] 1.7 0.946 [18,] 1.8 0.953 [19,] 1.9 0.964 [20,] 2.0 0.943 [21,] 2.1 0.937 [22,] 2.2 0.944 [23,] 2.3 0.945 [24,] 2.4 0.950 [25,] 2.5 0.954 [26,] 2.6 0.946 [27,] 2.7 0.945 [28,] 2.8 0.949 [29,] 2.9 0.956 [30,] 3.0 0.953 [31,] 3.1 0.941 [32,] 3.2 0.949 [33,] 3.3 0.943 [34,] 3.4 0.956 [35,] 3.5 0.950 [36,] 3.6 0.944 [37,] 3.7 0.952 [38,] 3.8 0.958 [39,] 3.9 0.938 [40,] 4.0 0.944 [41,] 4.1 0.950 [42,] 4.2 0.945 [43,] 4.3 0.948 [44,] 4.4 0.962 [45,] 4.5 0.969 [46,] 4.6 0.956 [47,] 4.7 0.950 [48,] 4.8 0.955 [49,] 4.9 0.946 [50,] 5.0 0.945 > range(resulting.matrix[,2]) [1] 0.917 0.969 -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702551.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] Coverage probability for a Poisson parameter
Hi, Roughly reading the code, I find this statement "phat <- x / m" is probably incorrect since this will give you the set of 100 observed x values /100. I redefine the function cover with three inputs: lambda for the parameter of the poisson distribution, sample.size and significance.level. The output is 1 or 0, depending on whether lambda is inside the confidence interval or not. With 5% level of significance I expect to get 95% of the time the parameter will be included in the confidence interval. The two runs below shows 96% and 94.8%, pretty close. > cover <- function(lambda, sample.size, significance.level) + { + x <- rpois(sample.size,lambda) + estimate <- mean(x) + lower <- estimate - qnorm(1 - significance.level/2) * sqrt(estimate/sample.size) + upper <- estimate + qnorm(1 - significance.level/2) * sqrt(estimate/sample.size) + if (lambda > lower & lambda < upper){1}else{0} + } > mean(sapply(1:100, function(x)cover(2.5,100,0.05))) [1] 0.96 > mean(sapply(1:1000, function(x)cover(2.5,10,0.05))) [1] 0.948 -- View this message in context: http://r.789695.n4.nabble.com/Coverage-probability-for-a-Poisson-parameter-tp4702535p4702537.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
Re: [R] 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. > > __ > 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.
Re: [R] 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. -- 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.
Re: [R] 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) # [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. > > __ > 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.
Re: [R] 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) } -- 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. __ 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.
Re: [R] 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 > [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. > > __ > 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.
[R] 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 [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. __ 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.
[R] Coverage probability for the normal distribution in "plot.spec.coherency" function.
Hello to the members of the list, I am using the 'spectrum' function from 'stats' package, to calculate the squared coherency between two time series. The function 'plot.spec.coherency' provides information for the coverage probability for the normal distribution. It seems that the calculation is based on Enochson and Goodman, 1965, "Gaussian approximations to the distribution of sample coherence". AFFDL-TR-65-57, Whright-Patterson Air Force base. Could someone tell me if this reference was really used or another one? In advance, thank you. Pascal __ 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.