Re: [R] Data framing question
Hi, You may also try: library(data.table) dt1 - data.table(dat1, key = identif) dt2 - copy(dt1) dt3 - copy(dt1) dt1[,c(roa1, eta1) := list(c(NA, roa[-.N]),c(NA, diff(eta))), by=identif] #or dt2[, `:=`(c(roa1, eta1), list(c(NA, roa[-.N]), c(NA, diff(eta, by = identif] # or dt3[, `:=`(roa1 = c(NA, roa[-.N]), eta1 = c(NA, diff(eta))), by = identif] identical(dt1, dt2) # [1] TRUE identical(dt1,dt3) #[1] TRUE identical(dat2, as.data.frame(dt1)) # [1] TRUE A.K. On Friday, June 6, 2014 10:47 PM, arun smartpink...@yahoo.com wrote: Hi, May be this helps: dat1 - read.table(text=identif roa eta 1 7 5 2 8 9 2 9 8 2 10 7 3 11 6 3 1 4 3 2 2 4 3 3 4 6 5,sep=,header=TRUE) dat2 -within(dat1, { eta1 - ave(eta, identif, FUN = function(x) c(NA, diff(x))) roa1 - ave(roa, identif, FUN = function(x) c(NA, x[-length(x)])) }) A.K. So here is my little problem. I don't know how to compute a new variable roa1, such that first data is split by identif then roa 1 will take the value of the line just before. so it will be something like roa1=(NA,NA,8,9,NA,11,1,NA,3) and the same for eta just that for this one first data is split by identif then deta will take the value of the difference between the actual value of eta with the previous something like (NA,NA,-1,-1,NA,-2,2,NA,-2) 21 mins · Like __ 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] finding average of entire rows to equal values in a vector
Hi, This is not clear. If this is a combination of rows using a specific formula as you showed, use ?combn dat - read.table(text=0.7, 0.3, 0.6, 0.9 1.0, 0.1, 0.4, 0.7 1.2, 0.8, 0.3, 0.1,sep=,,header=FALSE) indx - combn(seq(dim(dat)[1]), 2) vec1 - c(0.86, 0.46, 0.5, 0.63) indx1 - sapply(seq_len(ncol(indx)), function(i) { ind - indx[, i] val - unlist(round((2 * dat[ind[1], ] + dat[ind[2], ])/3, 2)) tol - 0.015 all(abs(val - vec1) = tol) }) dat[indx[, indx1], ] # V1 V2 V3 V4 #1 0.7 0.3 0.6 0.9 #3 1.2 0.8 0.3 0.1 A.K. Hi, I'm looking for a way to find the average of several rows, of a 4 x 7 matrix, that will eventually reach a fixed set of values in a vector while always averaging the entire row from the 4x7 to get there. So... in a matrix of 4 columns and 7 rows with values like 0.7, 0.3, 0.6, 0.9 1.0, 0.1, 0.4, 0.7 1.2, 0.8, 0.3, 0.1 ... .. . how can i write it to find the possible combinations of rows 1, 2, and 3 so that it matches the vector: 0.86, 0.46, 0.5, 0.63 with some small margin of error? in this case the answer is (2x(row 1)+(row 3))/3 best, -sunny0 __ 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] Package 'effects' plotting allEffects question
See ?Effect On Jun 6, 2014 9:38 PM, Nwinters nicholas.wint...@mail.mcgill.ca wrote: Is there a way to plot just the response function for one variable, instead of every variable on the same plot? my model has 13 variables and the plot produces a graph for each of these on the same plot, which is useless. -- View this message in context: http://r.789695.n4.nabble.com/Package-effects-plotting-allEffects-question-tp4691817.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] Help with factor levels and reference level
On Fri, 6 Jun 2014 11:16:11 AM Nwinters wrote: I have a variable coded in Stata as follows: ** *gen sat_pm25cat_=. replace sat_pm25cat_= 1 if (sat_pm25=4 sat_pm25=7.1 sat_pm25!=.) replace sat_pm25cat_= 2 if (sat_pm25=7.1 sat_pm25=10) replace sat_pm25cat_= 3 if (sat_pm25=10.1 sat_pm25=11.3) replace sat_pm25cat_= 4 if (sat_pm25=11.4 sat_pm25=12.1) replace sat_pm25cat_= 5 if (sat_pm25=12.2 sat_pm25=17.1) gen satpm25catR= A if sat_pm25cat_==1 replace satpm25catR= B if sat_pm25cat_==2 replace satpm25catR= C if sat_pm25cat_==3 replace satpm25catR= D if sat_pm25cat_==4 replace satpm25catR= E if sat_pm25cat_==5 *** my model for R is: ## *glm.PM25linB -glm(leuk ~ satpm25catR + sex + ageR, data=leuk, family=binomial, epsilon=1e-15, maxit=1000)* ## In the summary, satpm25catR is being reported as all levels: http://r.789695.n4.nabble.com/file/n4691823/Screen_Shot_2014-06-06_at_2.png *What I want is to make A the reference level, how do I do this??* Hi Nwinters, I get what you want with this example: leukdf- data.frame(leuk=sample(0:1,100,TRUE),sat_pm25=runif(100,0,17.1), sex=sample(c(M,F),100,TRUE),ageR=sample(20:75,100,TRUE)) leukdf$satpm25catR-factor(NA,levels=LETTERS[1:5]) leukdf$satpm25catR-factor(rep(NA,100),levels=LETTERS[1:5]) leukdf$satpm25catR[leukdf$sat_pm25 7.1]-A leukdf$satpm25catR[leukdf$sat_pm25 = 7.1 leukdf$sat_pm25 10.1]-B leukdf$satpm25catR[leukdf$sat_pm25 = 10.1 leukdf$sat_pm25 11.3]-C leukdf$satpm25catR[leukdf$sat_pm25 = 11.3 leukdf$sat_pm25 12.1]-D leukdf$satpm25catR[leukdf$sat_pm25 = 12.1 leukdf$sat_pm25 17.1]-E summary(glm(leuk ~ satpm25catR + sex + ageR, data=leukdf, family=binomial, epsilon=1e-15, maxit=1000)) Call: glm(formula = leuk ~ satpm25catR + sex + ageR, family = binomial, data = leukdf, epsilon = 1e-15, maxit = 1000) Deviance Residuals: Min 1Q Median 3Q Max -1.4813 -1.1798 0.7631 1.1347 1.5195 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 1.675650.87205 1.922 0.0547 satpm25catRB -0.522890.58578 -0.893 0.3721 satpm25catRC -0.799980.78405 -1.020 0.3076 satpm25catRD -0.364880.88162 -0.414 0.6790 satpm25catRE -0.653720.51461 -1.270 0.2040 sexM -0.540630.42073 -1.285 0.1988 ageR -0.020950.01455 -1.440 0.1500 (Dispersion parameter for binomial family taken to be 1) Null deviance: 138.59 on 99 degrees of freedom Residual deviance: 133.74 on 93 degrees of freedom AIC: 147.74 Number of Fisher Scoring iterations: 5 It may be a problem with the way you have calculated the categorical variable as David noted. However, if you haven't read a paper I had published a few years ago titled On the perils of categorizing responses, you might want to have a look. Jim __ 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] error in R program
REPLY TO ALL FOR THE R-HELP LIST!!! I apologize for the bluntness but you must realize that it is critical if you desire to get help on the mailing list beyond a single person your question must actually get to the mailing list. My expertise only goes so far and there is an infinitely larger community that could likely help you more quickly and often in a way that can increase performance. I am always happy to help when time allows but please make sure the r-help@r-project.org is a recipient as well. Regarding your code, after glancing at it again it is clear that you aren't actually submitting a covariance matrix to mvrnorm. The only place you specify phi1 is initially with NA's and then within the list init1. If you intend to use the list init1 you specify it within the function call: xi1-mvrnorm(1,c(0,0,0,0),init1$phi1, tol = 1e-6, empirical = FALSE, EISPACK = FALSE) Once again, please make sure you are not just sending these messages solely to me. Regards, Charles On Thu, Jun 5, 2014 at 9:04 AM, thanoon younis thanoon.youni...@gmail.com wrote: many thanks to you Dr. Charles Really i have a problem with simulation data in xi and now i have this erro r Error in mvrnorm(1, c(0, 0, 0), phi1, tol = 1e-06, empirical = FALSE, : incompatible arguments Regards On 5 June 2014 16:45, Charles Determan Jr deter...@umn.edu wrote: Hello again Thanoon, Once again, you should send these request not to me but to the r-help list. You are far more likely to get help from the greater R community than just me. Furthermore, it is not entirely clear where your error is. It is courteous to provide only the code that is run up to the error and commenting the location. I see you have a loop and all the more important to isolate where the error is taking place and commenting it out. My recommendations: 1. Set both t and i = 1 and try to run your loop code manually to locate the specific function providing the error. 2. Check your data inputs for the function. The error tells you that some data is missing or infinite. Perhaps a slight error in your data generation? 3. See if you can address the specific instance before running the full loops, it may be multiple instances or just one. I hope this provides some guidance, Charles On Thu, Jun 5, 2014 at 3:10 AM, thanoon younis thanoon.youni...@gmail.com wrote: Dear Dr. Charles i need your help to correct the winbugs code below to estimate parameters of SEM by using bayesian inference for two group. I have this errorError in eigen(Sigma, symmetric = TRUE, EISPACK = EISPACK) : infinite or missing values in 'x'. many thanks in advance R-CODE library(MASS) #Load the MASS package library(R2WinBUGS) #Load the R2WinBUGS package library(boa) #Load the boa package library(coda) #Load the coda package N1-200;N2=100; P1-10;P2-10 phi1-matrix(NA,nrow=4,ncol=4) #The covariance matrix of xi1 Ro1-matrix(NA,nrow=4,ncol=4) yo1-matrix(data=NA,nrow=200,ncol=10) phi2-matrix(NA,nrow=4,ncol=4) #The covariance matrix of xi2 Ro2-matrix(NA,nrow=4,ncol=4) yo2-matrix(data=NA,nrow=200,ncol=10) #Matrices save the Bayesian Estimates and Standard Errors Eu1-matrix(data=NA,nrow=100,ncol=10) SEu-matrix(data=NA,nrow=100,ncol=10) Elam1-matrix(data=NA,nrow=100,ncol=6) SElam1-matrix(data=NA,nrow=100,ncol=6) Egam1-matrix(data=NA,nrow=100,ncol=4) SEgam1-matrix(data=NA,nrow=100,ncol=4) Ephx1-matrix(data=NA,nrow=100,ncol=4) SEphx1-matrix(data=NA,nrow=100,ncol=4) Eb1-numeric(100); SEb-numeric(100) Esgd1-numeric(100); SEsgd-numeric(100) Eu2-matrix(data=NA,nrow=100,ncol=10) SEu2-matrix(data=NA,nrow=100,ncol=10) Elam2-matrix(data=NA,nrow=100,ncol=6) SElam2-matrix(data=NA,nrow=100,ncol=6) Egam2-matrix(data=NA,nrow=100,ncol=3) SEgam2-matrix(data=NA,nrow=100,ncol=3) Ephx2-matrix(data=NA,nrow=100,ncol=3) SEphx2-matrix(data=NA,nrow=100,ncol=3) Eb2-numeric(100); SEb2-numeric(100) Esgd2-numeric(100); SEsgd2-numeric(100) #Arrays save the HPD intervals mu.y1=array(NA, c(100,10,2)) lam1=array(NA, c(100,6,2)) gam1=array(NA, c(100,4,2)) sgd1=array(NA, c(100,2)) phx1=array(NA, c(100,4,2)) mu.y2=array(NA, c(100,10,2)) lam2=array(NA, c(100,6,2)) gam2=array(NA, c(100,3,2)) sgd2=array(NA, c(100,2)) phx2=array(NA, c(100,3,2)) DIC=numeric(100)#DIC values #Parameters to be estimated parameters-c (mu.y1,lam1,gam1,phi1,psi1,psd1,sgd1,phx1,mu.y2,lam2,gam2,phi2,psi2,psd2,sgd2,phx2) #Initial values for the MCMC in WinBUGS init1-list(uby1=rep(0.0,10),lam1=rep(0.0,10), gam1=c(1.0,1.0,1.0,1.0),psd1=1.0,phi1=matrix(data=c(1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0),ncol=4,byrow=TRUE),psi1=rep(0.0,10), xi1=matrix(data=rep(0.0,200),ncol=4),xi2=matrix(data=rep(0.0,200),ncol=4)) init2-list(mu.y1=rep(0.0,10),lam1=rep(0.0,10),
Re: [R] plot in package psych with function error.bars.by
Tham and Jim, As usual, my first response to this is when you find a problem with the psych package, write me (as author) as well as the R-help list. In addition, always include which version of psych you are running. That will help in the debugging. The current version 1.4.5 on CRAN draws “cats eyes” instead of error bars, unless you turn off that option. You might find that useful. Then, my comment to Tham, Yes, you have found a weakness (some would say a bug) in that I currently default the base plot character to be 15. I will correct this in the next release (which won’t be shipped to CRAN until mid to late July). In the interim, I will try to get a fix up on the personality-project.org/r repository by early next week. To Jim, I can not get your error at all. Tham was finding a problem with the basic example, which works, unless you try to specify the pch, which doesn’t work. Bill On Jun 3, 2014, at 4:10 PM, Jim Lemon j...@bitwrit.com.au wrote: On Mon, 2 Jun 2014 11:28:19 PM Tham Tran wrote: Hi, I have a problem with the function error.bars.by in package psych. This is the code for example of a graph: keys.list=list(Agree=c(-1,2:5),Conscientious=c(6:8,-9,-10),Extraversion=c(-1 1,-12,13:15),Neuroticism=c(16:20),Openness = c(21,-22,23,24,-25)) keys = make.keys(28,keys.list,item.labels=colnames(bfi)) scores = score.items(keys,bfi,min=1,max=6) require(psych) error.bars.by(scores$scores,round(bfi$age/10)*10,by.var=TRUE, main=BFI age trends,legend=3,labels=colnames(scores$scores), xlab=Age,ylab=Mean item score) I need to change the plotting character and line type of the graph according to the scores (Agree,Conscientious,Extraticism,Openness). I have tried with: error.bars.by(scores$scores,round(bfi$age/10)*10,by.var=TRUE, main=BFI age trends,legend=3,labels=colnames(scores$scores), pch=c(1,2,3,4), lty=1 ,xlab=Age,ylab=Mean item score) But there is a problem: Error in localWindow(xlim, ylim, log, asp, ...) : formal argument pch matched by multiple actual arguments Anyone can help me for this problem. Hi Tham, When I run your example, I get the following error: Error in seq.default(clim, -clim, 0.01) : 'from' cannot be NA, NaN or infinite In addition: Warning messages: 1: In qt(1 - alpha/2, group.stats[[g]]$n - 1) : NaNs produced 2: In dt(ln, n - 1) : NaNs produced 3: In qt(alpha/2, n - 1) : NaNs produced As I don't know what is happening here, I can't do anything about the pch problem. Jim __ 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. William Revellehttp://personality-project.org/revelle.html Professor http://personality-project.org Department of Psychology http://www.wcas.northwestern.edu/psych/ Northwestern Universityhttp://www.northwestern.edu/ Use R for psychology http://personality-project.org/r It is 5 minutes to midnighthttp://www.thebulletin.org __ 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] Package 'effects' plotting allEffects question
just use plot(Effect(myvar, mymodel)) On 6/6/2014 11:28 AM, Nwinters wrote: Is there a way to plot just the response function for one variable, instead of every variable on the same plot? my model has 13 variables and the plot produces a graph for each of these on the same plot, which is useless. -- View this message in context: http://r.789695.n4.nabble.com/Package-effects-plotting-allEffects-question-tp4691817.html Sent from the R help mailing list archive at Nabble.com. -- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. Chair, Quantitative Methods York University Voice: 416 736-2100 x66249 Fax: 416 736-5814 4700 Keele StreetWeb: http://www.datavis.ca Toronto, ONT M3J 1P3 CANADA __ 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] Interval for rnorm command?
Hello everyone! I want to create a data vector using the rnorm command. Let's say /x - rnorm(100, 1, 2)/. This gives me 100 values from the specific normal distribution. However, I want to generate a variable that has a certain range of values, e.g. values can't exceed the interval /[-3.7; 3.7]/. How can I implement such an interval? Thanks for your help and sorry if this is totally trivial, I am new to R. -- View this message in context: http://r.789695.n4.nabble.com/Interval-for-rnorm-command-tp4691856.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] Help on postscript
Hello everyone, I'm a beginner on R, and I recently met a paragraph of code that I don't quite understand, anyone could nicely offer me some detailed explanation on it? Many thanks, and I'll wait on-line. Here is the code: postscript(SPreturns_tplot.ps,width=6,height=5) # Figure 19.2 par(mfrow=c(1,1)) qqplot(SPreturn, qt(grid,df=2.9837),main=t-probability plot, df=2.9837,xlab=data,ylab=t-quantiles) abline(lm(qt(c(.25,.75),df=2.9837)~quantile(SPreturn,c(.25,.75 graphics.off() I suppose this to be a plotting process, however, after I run it in RStudio, there is no response like plotting, so I can't figure out what exactly is this code telling about. Thanks again! [[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] copula fitting
Hi guys, can i fit a copula to two marginal distributions with different sample size? like one has 2340 observations and other has 1912. thanks Mudit [[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] Discrete distributions in R
Hello there, I am a new R user. I would like to test eight discrete distributions as Barrigossi et al. 2001. Goodness-of-fit was indicated by chi-square test. They used a Fortran program for fitting the discrete frequency distributions: http://r.789695.n4.nabble.com/file/n4691859/discrete_distributions.png Please, can you guide me to obtain similar tests for other data sets? Barrigossi, JAF, et al. Spatial and probability distribution of Mexican bean beetle (Coleoptera: Coccinellidae) egg mass populations in dry bean. Environmental entomology 30.2 (2001): 244-253. Thank you. -- View this message in context: http://r.789695.n4.nabble.com/Discrete-distributions-in-R-tp4691859.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] rake() error message
I'm teaching myself how to use rake() in the R survey package, using Thomas Lumley's Complex Surveys book. Working with one of my data sets, I received an error that has me stumped. Here is the code I used: d473.raked - rake(sdes473, sample=list(~m01, ~m02, ~m03, ~m04, ~m05, ~m06c, ~m07c, ~m08, ~m09, ~m10c, ~m11, ~m12c, ~m13, ~m14c, ~m15c, ~m16c), population=list(pop.m01, pop.m02, pop.m03, pop.m04, pop.m05, pop.m06, pop.m07, pop.m08, pop.m09, pop.m10, pop.m11, pop.m12, pop.m13, pop.m14, pop.m15, pop.m16)) Here is the error I received: Error in array(dim = extent, dimnames = namelist) : vector is too large In addition: Warning messages: 1: In ngroup * (as.integer(index) - one) : NAs produced by integer overflow 2: In group + ngroup * (as.integer(index) - one) : NAs produced by integer overflow 3: In ngroup * nlevels(index) : NAs produced by integer overflow Any guidance on this would be much appreciated. Thanks, Michael Willmorth [[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] Help with R package
Dear R Community, I am in process of developing an R package which in turn depends on a package that is not available in CRAN but has to be downloaded from a web source (as follows). Could someone guide me on how to include the package listed in Depends of the Description file. I have the package installed in my computer: Source of Package: source(âhttp://bioconductor.org/biocLite.Râ) biocLite(âPackagenameâ) library(Packagename) Thank you in advance, Milan [[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] Help on postscript
On Sat, 7 Jun 2014 06:48:49 AM Yijia Wang wrote: Hello everyone, I'm a beginner on R, and I recently met a paragraph of code that I don't quite understand, anyone could nicely offer me some detailed explanation on it? Many thanks, and I'll wait on-line. Here is the code: postscript(SPreturns_tplot.ps,width=6,height=5) # Figure 19.2 par(mfrow=c(1,1)) qqplot(SPreturn, qt(grid,df=2.9837),main=t-probability plot, df=2.9837,xlab=data,ylab=t-quantiles) abline(lm(qt(c(.25,.75),df=2.9837)~quantile(SPreturn,c(.25,.75 graphics.off() I suppose this to be a plotting process, however, after I run it in RStudio, there is no response like plotting, so I can't figure out what exactly is this code telling about. Hi Yijia, The above should create a file named SPreturns_tplot.ps in the working directory of your R session (use getwd() to find it). You should be able to see the image using a Postscript display manager like Ghostview - gs in Linux, GSview in Windows. Jim __ 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] Interval for rnorm command?
On 08/06/14 02:23, zagreus wrote: Hello everyone! I want to create a data vector using the rnorm command. Let's say /x - rnorm(100, 1, 2)/. This gives me 100 values from the specific normal distribution. However, I want to generate a variable that has a certain range of values, e.g. values can't exceed the interval /[-3.7; 3.7]/. How can I implement such an interval? What you are asking for amounts to drawing a sample from a *truncated* normal distribution. Try RSiteSearch(truncated normal) to find appropriate functions to effect this procedure. cheers, Rolf Turner __ 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.