[R] lme4 print and summary errror
Hi all - (this is posted to r-help and R-SIG-MAC) OSX 10.3.7, R 2.0.1, lme4/Matrix/latticeExtra latest, fresh install of R. MASS loaded (or not). I am getting an error message for the print() and summary() commands with all lme models I try and run in lme4 (GLMM's work fine). Using the example from the lme help, summary and print produce the following errors, despite the model being fit, as indicated by VarCorr() and examination of str(fm). Any ideas? (I can't reproduce this on a windows95 install of 2.0.1, so I am guessing it may be a mac thing at the moment? This happens with the binary or the source installation of lme4.) Cheers andrew data(bdf) fm - lme(langPOST ~ IQ.ver.cen + avg.IQ.ver.cen, data = bdf, + random = ~ IQ.ver.cen | schoolNR) summary(fm) Error in verbose || attr(x, verbose) : invalid `y' type in `x || y' fm Linear mixed-effects model fit by Data: NULL Log-likelihood: NULL Fixed: list() NULL Length Class Mode 0 NULL NULL Number of Observations: Number of Groups: Error in 1:dd$Q : NA/NaN argument VarCorr(fm) Groups NameVariance Std.Dev. Corr schoolNR (Intercept) 8.07702 2.84201 IQ.ver.cen 0.20806 0.45614 -0.642 Residual 41.34942 6.43035 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] GLMM and crossed effects
Hi again. Perhaps a simple question this time I am analysing data with a dependent variable of insect counts, a fixed effect of site and two random effects, day, which is the same set of 10 days for each site, and then transect, which is nested within site (5 each). I am trying to fit the cross classified model using GLMM in lme4. I have, for potential use, created a second coding of transect with levels 1-5 for site 1 and 6-10 for site2. Likewise, if a groupedData object is necessary, there are als ts1 and ts2 dummy variables, as was necessary in the old lme. str(dat3) `data.frame': 100 obs. of 7 variables: $ site : Factor w/ 2 levels Here,There: 1 1 1 1 1 1 1 1 1 1 ... $ day : Factor w/ 10 levels 1,2,3,4,..: 1 1 1 1 1 2 2 2 2 2 ... $ trans : Factor w/ 5 levels 1,2,3,4,..: 1 2 3 4 5 1 2 3 4 5 ... $ count : int 77 109 81 124 115 84 90 85 130 106 ... $ trans2: Factor w/ 10 levels 1,2,3,4,..: 1 2 3 4 5 1 2 3 4 5 ... $ ts1 : Factor w/ 10 levels Here 1,Here 2,..: 1 2 3 4 5 1 2 3 4 5 ... $ ts2 : Factor w/ 10 levels Here 1,Here 2,..: 1 2 3 4 5 1 2 3 4 5 ... Might someone explain to me how I might reflect the fact that transects are different between sites, while days are not? #this does not work, though I thought it might be the best way to specify the model. GLMM(count~site,data=dat3,random=list(day=~1,trans=~1|site,family=poisso n) Error in GLMM(count ~ site, data = dat3, random = list(day = ~1, trans = ~1 | : subscript out of bounds In addition: Warning message: | not meaningful for factors in: Ops.factor(1, site) #This does... but also note the differences in the summary and VarCorr variance components... summary(GLMM(count~site,data=dat3,random=list(day=~1,trans2=~1),family= poisson)) Generalized Linear Mixed Model Family: poisson family with log link Fixed: count ~ site Data: dat3 AIC BIC logLik 103.1494 116.1753 -46.5747 Random effects: Groups NameVariance Std.Dev. trans2 (Intercept) 0.073011 0.27020 day(Intercept) 0.034373 0.18540 # of obs: 100, groups: trans2, 10; day, 10 Estimated scale (compare to 1) 0.6232135 Fixed effects: Estimate Std. Error z value Pr(|z|) (Intercept) 4.662800.13502 34.534 2e-16 *** siteThere -0.255720.17216 -1.485 0.1375 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Correlation of Fixed Effects: (Intr) siteThere -0.636 VarCorr(GLMM(count~site,data=dat3,random=list(day=~1,trans2=~1),family= poisson)) Groups NameVariance Std.Dev. trans2 (Intercept) 0.028936 0.17010 day (Intercept) 0.013623 0.11672 Residual 0.396322 0.62954 Many thanks andrew __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lme, glmmPQL, multiple random effects
, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 4.4834749486, 5.11362806892812, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214, 4.76190015454611, 5.59005352989214)), method = ML, weights = varFixed(~invwt)) 3: eval(expr, envir, enclos) 2: eval(mcall) 1: glmmPQL(Mate1 ~ Cross, data = gd, random = pdBlocked(list(pdIdent(~Female - 1), pdIdent(~Male - 1))), family = binomial) - Dr. Andrew Beckerman Department of Animal and Plant Sciences, University of Sheffield, Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK ph +44 (0)114 222 0026; fx +44 (0)114 222 0002 http://www.shef.ac.uk/beckslab -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] LME-glmmPQL formulation
Hi all - R2.0.1 on OSX;MASS library;nlme library I am trying to emulate the solution to a problem set that has normally been run in Genstat, using R. The problem that I am having at the moment is with the following glmm question (using glmmPQL from the MASS library): We have two different forest habitats (first rotation thicket, and high forest) which we want to survey for the presence of our study animal. We survey both habitats on each of 10 days, and within each habitat we have five transects. The sampling unit is the number of animals counted per transect. We therefore have two sources of random variation: · counts will vary between days due, say, to variation in weather · counts will vary between transects, within sites, for any number of (known and unknown) reasons There is no relationship between transects at the two sites: transect 1 in site 1 has no link to transect 1 in site 2, etc. The random term for transectis therefore nested within site, while the main effect of site, which is what we are interested in, is a fixed effect. summary(dat) site day trans count Here :50 Min. : 1.0 Min. :1 Min. : 48.00 There:50 1st Qu.: 3.0 1st Qu.:2 1st Qu.: 79.00 Median : 5.5 Median :3 Median : 95.00 Mean : 5.5 Mean :3 Mean : 95.85 3rd Qu.: 8.0 3rd Qu.:4 3rd Qu.:112.25 Max. :10.0 Max. :5 Max. :165.00 In Genstat, the (supposed) procedure is to fit a model with site as a fixed effect and then a random effects model of day+transect.site, where the transect.site indicates that there are 5 transects nested within each site. My first thought was the following: glmmPQL(count~site,data=dat,random=~day|site/transect, family=poisson) however, the random effects are not separated into day and site/transect. Instead, there is day|site and day|site %in% transect, which I realize makes sense in light of the model formulation. my second guess was glmmPQL(count~site,random=list(~day|site,~1|trans),family=poisson,data =dat2) which estimates a random effect on ~day|site and on ~1|trans%in%site. which seems more appropriate, but does not give the same answers as I have for the genstat; nor does it estimate the p-value for site. I guess my question is how to separate the two random effects so that there is a estimate for day and for transect/site. I would be happy to provide the data if anyone needs/wants IT. Cheers andrew - Dr. Andrew Beckerman Department of Animal and Plant Sciences, University of Sheffield, Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK ph +44 (0)114 222 0026; fx +44 (0)114 222 0002 http://www.shef.ac.uk/beckslab -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R sumo package suggestion
to install a large bundle of xemacs packages at one time (about a 120 modes including ESS). I think R should have a similar bundle. It would be so much easier than hunting/downloading/installing. Martin encouraged me to send this suggestion to r-help. In addition, he put together a few comments relating to the previous times that this, or a similar suggestion, has been brought up here. Martin wrote: If you search for install all CRAN packages on http://maths.newcastle.edu.au/~rking/R/ (the URL which is quickly found from the [Search] sidebar of http://www.R-project.org/) You find things like Greg Warnes 'Makefile' http://tolstoy.newcastle.edu.au/R/help/04/04/0723.html and http://tolstoy.newcastle.edu.au/R/help/04/04/0616.html which is from Tony and has the following small function: installNewCRANPackages - function() { ## (C) A.J. Rossini, 2002--2004 test2 - packageStatus()$avail[Status] install.packages(row.names(test2)[which(test2$Status==not installed)]) } -- Rodney Sparapani Medical College of Wisconsin Sr. Biostatistician Patient Care Outcomes Research [EMAIL PROTECTED] http://www.mcw.edu/pcor Was 'Name That Tune' rigged? WWLD -- What Would Lombardi Do __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html - Dr. Andrew Beckerman Department of Animal and Plant Sciences, University of Sheffield, Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK ph +44 (0)114 222 0026; fx +44 (0)114 222 0002 http://www.shef.ac.uk/beckslab -- __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bootstrap, lme, random effects
Hi there. OSX/R2.0 We are trying to implement a bootstrap of the coeffecients of a mixed effect model. In particular, we are interested in the intercept and slope of the random effects. Following from the basics for a linear model, we construct our lme models and a boot function: library(nlme) library(boot) data-read.csv(~/data.csv) bootcoef-function(data,index){ dat-data[index,] mod-lme(Frames~Man2+Manip+Strings+Date.+Cut.,random=~Man2|ID.,data=dat) fixef(mod) } boot.out-boot(data,bootcoef,99) boot.ci(boot.out) # produces information via boot.ci() that suggests this is not necessarily successful ORDINARY NONPARAMETRIC BOOTSTRAP Call: boot(data = data, statistic = bootcoef, R = 99) Bootstrap Statistics : original biasstd. error t1* 16.125904015 5.299827478 9.98818463 t2* 0.010901682 -0.004134585 0.01621935 t3* -0.038168126 -0.078833467 0.35778286 t4* 1.101486342 -0.021886290 0.45720400 t5* 0.005982241 -0.009140839 0.01563175 t6* 2.729537567 0.287663533 1.56150779 boot.ci(boot.out) BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 99 bootstrap replicates CALL : boot.ci(boot.out = boot.out) Intervals : Level Normal Basic Studentized 95% ( -8.75, 30.40 ) ( -8.40, 36.01 ) (-40.93, 35.15 ) Level PercentileBCa 95% (-3.75, 40.65 ) (-5.55, 28.65 ) Calculations and Intervals on Original Scale Some basic intervals may be unstable Some studentized intervals may be unstable Some percentile intervals may be unstable Warning : BCa Intervals used Extreme Quantiles Some BCa intervals may be unstable Warning messages: 1: NaNs produced in: sqrt(tv[, 2]) 2: Extreme Order Statistics used as Endpoints in: norm.inter(t, adj.alpha) As stated above, we are interested in the ranef(mod) components. including this instead of fixef(mod) results in the error: bootcoef-function(data,index){ + dat-data[index,] + mod- lme(Frames~Man2+Manip+Strings+Date.+Cut.+ID.*Man2+ID.*Manip,random=~Man2 |ID.,data=dat) + ranef(mod) + } boot.out-boot(data,bootcoef,99) Error: incorrect number of subscripts on matrix Suggesting that the setup of the ranef(mod) list is different (clearly) Any suggestions on any of this? I have a sneaking suspicion this is not a straightforward issue. Cheers andrew - Dr. Andrew Beckerman Department of Animal and Plant Sciences, University of Sheffield, Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK ph +44 (0)114 222 0026; fx +44 (0)114 222 0002 http://www.shef.ac.uk/beckslab -- __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lme with poly(x,2) terms
Hi there. Mac OSX 3.3.4 R 1.9.1 I am analysing a data set with the following model m4- lme(fixed=sr~time*poly(energy,2)*poly(dist,2),random=~time|pot,data=deh) where time is one of six months, pot is a jar in which the repeated measures of species number (sr) was made. energy and dist (disturbance) are fixed experimental treatments. We are trying to test the hypothesis that there is an interaction between energy and disturbance that varies through time, with the expectation that sr varies quadratically with energy and with disturbance. Our difficulty is interpreting the various outputs from the model, assuming it is specified correctly - sorry if this is more a stats question than a R mechanics question. summary(m1) and anova(m1) produce the tables below the . Q1) Am i correct to assume that the anova table is sequential? Q2) How does one interpret the fixed effects/coefficients table? Do the insignificant terms for poly(dist)2 all the way down (Up) to its main effect suggest that a quadratic function in dist is not significant? Q3) If we remove the quadratic term in dist and compare it to the model with poly(dist,2), the anova says the polynomial is significant anova(update(m2,~.,method=ML),update(m4,~.,method=ML)) Model df AIC BIClogLik Test L.Ratio p-value update(m2, ~., method = ML) 1 16 2781.683 2858.271 -1374.841 update(m4, ~., method = ML) 2 22 2771.380 2876.688 -1363.690 1 vs 2 22.303 0.0011 despite only the main effect of poly(dist,2) being significant in the terms. Is the best approach to use the anova test or the coefficients? How does one justify the insignificance of every term with poly(dist)2 in it? Many thanks in advance andrew - summary(m1) Linear mixed-effects model fit by REML Data: deh AIC BIClogLik 2687.974 2792.830 -1321.987 Random effects: Formula: ~time | pot Structure: General positive-definite, Log-Cholesky parametrization StdDevCorr (Intercept) 1.5503393 (Intr) time0.1858609 -0.862 Residual0.9234853 Fixed effects: sr ~ time * poly(energy, 2) * poly(dist, 2) Value Std.Error DF t-value p-value (Intercept)8.2424 0.14576 721 56.54737 0. time -1.1447 0.02376 721 -48.16926 0. poly(energy, 2)1 18.2052 4.34118 721 4.19361 0. poly(energy, 2)2 -43.8133 4.34213 721 -10.09028 0. poly(dist, 2)1-9.9600 4.34169 721 -2.29403 0.0221 poly(dist, 2)2 -10.6639 4.34198 721 -2.45599 0.0143 time:poly(energy, 2)1 1.7320 0.70705 721 2.44961 0.0145 time:poly(energy, 2)2 5.6245 0.70695 721 7.95608 0. time:poly(dist, 2)1 -0.6569 0.70701 721 -0.92908 0.3532 time:poly(dist, 2)20.0400 0.70697 721 0.05657 0.9549 poly(energy, 2)1:poly(dist, 2)1 356.6786 128.77967 721 2.76968 0.0058 poly(energy, 2)2:poly(dist, 2)1 -99.7288 128.60505 721 -0.77547 0.4383 poly(energy, 2)1:poly(dist, 2)2 -11.4295 129.65263 721 -0.08816 0.9298 poly(energy, 2)2:poly(dist, 2)2 149.5420 129.80979 721 1.15201 0.2497 time:poly(energy, 2)1:poly(dist, 2)1 -79.3803 20.96606 721 -3.78613 0.0002 time:poly(energy, 2)2:poly(dist, 2)1 59.4570 20.93577 721 2.83997 0.0046 time:poly(energy, 2)1:poly(dist, 2)2 -20.6131 21.10723 721 -0.97659 0.3291 time:poly(energy, 2)2:poly(dist, 2)2 -22.3304 21.13159 721 -1.05673 0.2910 anova(m4) numDF denDF F-value p-value (Intercept)1 721 888.6686 .0001 time 1 721 2321.2473 .0001 poly(energy, 2)2 721 77.1328 .0001 poly(dist, 2) 2 721 22.9940 .0001 time:poly(energy, 2) 2 721 34.6873 .0001 time:poly(dist, 2) 2 7210.4551 0.6345 poly(energy, 2):poly(dist, 2) 4 7212.5824 0.0361 time:poly(energy, 2):poly(dist, 2) 4 7216.1290 0.0001 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] converting apply output
Hi - platform powerpc-apple-darwin6.8 status major1 minor9.0 year 2004 I am trying to deal with the output of apply(). As indicated, when each call to 'FUN' returns a vector of length 'n', then 'apply' returns an array of dimension 'c(n, dim(X)[MARGIN])'. However, I would like this to be a list in the same format as is produced when 'FUN' return vectors of different lengths ('apply' returns a list of length 'dim(X)[MARGIN]'). e.g. tt1-c(0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0) m1-matrix(tt1,10,10) out-apply(m1,2,function(x) which(x==1)) produces an array, out [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,]897867666 7 [2,]9 108988 10 107 9 but I would like out as a list of 10 elements with two elements in each, e.g. [[1]] [1] 8 9 [[2]] [1] 9 10 etc. I have tried apply(out,2,function(x) list(x))), but the subsrcripting is not equal to the pattern when FUN returns a vectors of different length. Any help would be appreciated. Cheers andrew __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] converting apply output
Thanks eric... I figured this routine out as well. cvt-function(dat){ x-as.list(rep(0,dim(dat)[2])) for(i in 1:dim(dat)[2]){ x[[i]]-dat[,i] x}} # get ragged array of 1's dat-apply(mat,2,function(x) which(x==1)) # deal with this using cvt to creat list if(is.null(dim(dat))) dat2-dat else dat2-cvt(dat) On 23 Jun 2004, at 13:13, Eric Lecoutre wrote: Hi, Why not try with the data.frame structure, wich internally yet consists in a list: lapply(as.data.frame(m1),function(x) which(x==1)) $V1 [1] 8 9 $V2 [1] 9 10 [...] Eric At 12:53 23/06/2004, Andrew Beckerman wrote: Hi - platform powerpc-apple-darwin6.8 status major1 minor9.0 year 2004 I am trying to deal with the output of apply(). As indicated, when each call to 'FUN' returns a vector of length 'n', then 'apply' returns an array of dimension 'c(n, dim(X)[MARGIN])'. However, I would like this to be a list in the same format as is produced when 'FUN' return vectors of different lengths ('apply' returns a list of length 'dim(X)[MARGIN]'). e.g. tt1-c(0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0) m1-matrix(tt1,10,10) out-apply(m1,2,function(x) which(x==1)) produces an array, out [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,]897867666 7 [2,]9 108988 10 107 9 but I would like out as a list of 10 elements with two elements in each, e.g. [[1]] [1] 8 9 [[2]] [1] 9 10 etc. I have tried apply(out,2,function(x) list(x))), but the subsrcripting is not equal to the pattern when FUN returns a vectors of different length. Any help would be appreciated. Cheers andrew __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Eric Lecoutre UCL / Institut de Statistique Voie du Roman Pays, 20 1348 Louvain-la-Neuve Belgium tel: (+32)(0)10473050 [EMAIL PROTECTED] http://www.stat.ucl.ac.be/ISpersonnel/lecoutre If the statistics are boring, then you've got the wrong numbers. -Edward Tufte - Dr. Andrew Beckerman Department of Animal and Plant Sciences, University of Sheffield, Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK ph +44 (0)114 222 0026; fx +44 (0)114 222 0002 http://www.shef.ac.uk/beckslab -- __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html