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,
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,
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,
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,
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,
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 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.
[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 || thetau){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.
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 || thetau){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.
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 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 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,
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.
[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.
