Re: [R] Visualize Sparse Matrix.
Hi, Thanks for your help. I used the SparseM package http://www.econ.uiuc.edu/~roger/research/sparse/SparseM.pdf <http://www.econ.uiuc.edu/~roger/research/sparse/SparseM.pdf> First of all, I create a class for sparse matrices stored in Compressed Sparse Row (CSR) with as.matrix.csr(matrix). After that, I plot the non-zero entries of a matrix of class matrix.csr with image(m.csr) This is my code: library(SparseM) data <- read.csv(pathCSV, header = FALSE, sep = ",") numcol <- ncol(data) dMatrix <- matrix(unlist(data), ncol = numcol, byrow = TRUE) dMatrix.csr <- as.matrix.csr(dMatrix) image(dMatrix.csr, col=c("white","blue")) After clustering, I will have the same matrix but each row (vector) has a tag to represent a cluster id. So, how could I plot my matrix to show a different color for cluster id? This is an example of my results: 213 0 0 0 0.213 0.3423 345 0 0 0.32 0 0 84 0 0.4 0 0.540 84 0.860 0 0 0 213 0 0.98 0 0 0.45 345 0 0.57 0 0 0.4 Cheers. > On Jun 10, 2016, at 18:51, Amos Elberg wrote: > > Sparse matrix visualization is a feature of my largeVis package: > https://github.com/elbamos/largeVis/tree/0.1.6 > <https://github.com/elbamos/largeVis/tree/0.1.6> > > > > On Thu, Jun 9, 2016 at 6:27 PM, FRANCISCO XAVIER SUMBA TORAL > mailto:xavier.sumb...@ucuenca.ec>> wrote: > Hi, > > First of all, sorry for my question it could be so basic for a common user in > R, but I am starting with this new environment. > > I have done a clustering job and I would like to visualize my vectors. I have > a matrix of TF-IDF weights of 4602 x 1817. I store the values in a CSV file. > How can I visualize my vectors in a 2D-space? > > After that, I execute a clustering algorithm and I got a label for each > cluster. How can I visualize my vectors resulting base on a color or figure > for each cluster? > > This is the code that I am having trying to accomplish my graphs: > > data <- read.csv(pathFile,header = FALSE, sep = ",”) > dMatrix <- matrix(unlist(data), ncol = 4602, byrow = TRUE) # Use a matrix to > use melt. > # Graph my data > ggplot(melt(dMatrix), aes(Var1,Var2, fill=value)) + geom_raster() + > scale_fill_gradient2(low='red', high=‘black', mid=‘white') + theme_bw() + > xlab("x1") + ylab("x2") > > > Cheers. > [[alternative HTML version deleted]] > > __ > R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To > UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > <https://stat.ethz.ch/mailman/listinfo/r-help> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > <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 -- 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] Visualize Sparse Matrix.
Hi Jim, Thanks for your answer. I try your code example, but it is basically the same that I had it. I want to visualise my matrix something like this image: With the graphics that I already have is difficult to visualise my data. I am getting this results: 1) With my first code, I got this: 2) With Jim’s code. I got this: Ho can I make my graphs more observable as in the first figure? My graphs shows points as if my screen was dirty. Cheers. > On Jun 10, 2016, at 04:39, Jim Lemon wrote: > > Hi Francisco, > I tried this just to see if it would work. It did, after a while. > > wtmat<-matrix(rnorm(4602*1817),nrow=4602) > library(plotrix) > x11(width=5,height=13) > color2D.matplot(wtmat,c(1,1,0),c(0,1,0),0,border=FALSE) > > Jim > > On Fri, Jun 10, 2016 at 8:27 AM, FRANCISCO XAVIER SUMBA TORAL > wrote: >> Hi, >> >> First of all, sorry for my question it could be so basic for a common user >> in R, but I am starting with this new environment. >> >> I have done a clustering job and I would like to visualize my vectors. I >> have a matrix of TF-IDF weights of 4602 x 1817. I store the values in a CSV >> file. How can I visualize my vectors in a 2D-space? >> >> After that, I execute a clustering algorithm and I got a label for each >> cluster. How can I visualize my vectors resulting base on a color or figure >> for each cluster? >> >> This is the code that I am having trying to accomplish my graphs: >> >> data <- read.csv(pathFile,header = FALSE, sep = ",”) >> dMatrix <- matrix(unlist(data), ncol = 4602, byrow = TRUE) # Use a matrix to >> use melt. >> # Graph my data >> ggplot(melt(dMatrix), aes(Var1,Var2, fill=value)) + geom_raster() + >> scale_fill_gradient2(low='red', high=‘black', mid=‘white') + theme_bw() + >> xlab("x1") + ylab("x2") >> >> >> Cheers. >>[[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. __ 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] Visualize Sparse Matrix.
Hi, First of all, sorry for my question it could be so basic for a common user in R, but I am starting with this new environment. I have done a clustering job and I would like to visualize my vectors. I have a matrix of TF-IDF weights of 4602 x 1817. I store the values in a CSV file. How can I visualize my vectors in a 2D-space? After that, I execute a clustering algorithm and I got a label for each cluster. How can I visualize my vectors resulting base on a color or figure for each cluster? This is the code that I am having trying to accomplish my graphs: data <- read.csv(pathFile,header = FALSE, sep = ",”) dMatrix <- matrix(unlist(data), ncol = 4602, byrow = TRUE) # Use a matrix to use melt. # Graph my data ggplot(melt(dMatrix), aes(Var1,Var2, fill=value)) + geom_raster() + scale_fill_gradient2(low='red', high=‘black', mid=‘white') + theme_bw() + xlab("x1") + ylab("x2") Cheers. [[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.
[R] Good pointers for understanding the R language implementation
Dear All, I'm currently working on a project with the purpose of remotely executing R code, which requires me to have to work with the code of R itself. I've searched the Internet for good information that will help me understand how R is implemented but what I've got so far isn't detailed enough. I've looked specifically at CRAN's manuals on the official website but they only address this issue briefly. I've also looked at other contents online but so far nothing has turned up that has the level of detail that I need to properly understand the inner workings of R. For example, I need to understand how exactly an expression is parsed and evaluated, because I will need to intervene in the process to decide whether to execute it remotely or not. Does anyone know of good pointers that would help me understand this? Thanks for any help! Best regards, Francisco [[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.
[R] Question about CHAID
Hallo there, I would like to work with CHAID, but the newest version of R does have it. So I thought I could use an older version of R which accepts or has the library CHAID. Could you tell me which version it is and where to download it? Thanks a lot in advance. Francisco __ 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] nlsystemfit help
Hello, For my master thesis I have fitted an individual tree diameter growth model and a survival probability model separately using R, but I was told that simultaneous estimation of these two models would minimize overall errors and provide a variance-covariance matrix as a whole. In that respect, can you please tell me if I can do it with nlssystemfit using SUR (seemingly unrelated regression) method? If not, do you know how can I do it in R? My equations are: - tree diameter growth model d_richards_k<-d2~A*(1-exp(-(k0+k1*puro+k2*GL1/100+k3*1/G1+k4*Tmedmax/100+k5* Perc_G_ec1))*(1-(d1/A)^(1-m)))^(1/(1-m)) nls_d_richards_k<-nlsLM(d_richards_k,start=list(A=100,k0=0.6,k1=0,k2=0,k3=0, k4=0,k5=0,m=0.6),control=nls.lm.control(maxiter=500)) - survival probability model mortal_final<-glm(sobrev~Verao+alto_fuste:puro+IC1arv+G1+Perc_G_ec1,family=b inomial) Thank you! Best regards, Francisco Goes [[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] I need help computing PRESS statistics (qpcR package) of a nlsLM model (minpack.lm package)
Dear all, I have used nlsLM in minpack.lm package for fitting a model using 22358 observations. Now I would like to compute PRESS residuals and have been trying it using qpcR package but it does not seem to work (I have fitted my data with a lm model and used qpcR PRESS and it worked fine). The message I get is: Error in as.data.frame.default(x[[i]], optional = TRUE) : cannot coerce class ""formula"" to a data.frame And also, with a different and smaller database, I got this message: .Error in PRESS.res[i] <- DATA[i, RESP.pos] - y.hat : replacement has length zero Can anyone please help me? Best regards, Francisco Goes [[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] (no subject)
Hey! http://expertiza-tv.ru/_1b3a_your.time.came_e5c2_.html?hicjpocix148039 [[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] Linear combination of time series for approximating another one
Dear R community, I would like to approximate a time serie as linear combination of a set of many time series, minimizing the number of the time series involved in the linear combination. Can you recommend to me any method for this? Any paper about this? I will appreciate much your suggestions. Thanks in advance [[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] Linear combination of time series for approximating another one.
Dear R community, I would like to approximate a time serie as linear combination of a set of many time series, minimizing the number of the time series involved in the linear combination. Can you recommend to me any method for this? Any paper about this? I will appreciate much your suggestions. Thanks in advance [[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] From Java to R OOP
Thanks MW! Your comments were very helpful. I also see how the NextMethod() works now. On Mar 25, 2013, at 5:03 AM, R. Michael Weylandt wrote: > On Mon, Mar 25, 2013 at 6:51 AM, Francisco J. Bido wrote: >> Hi, I'm new to OOP in R so please forgive the naiveness of some of the >> questions. Here are a couple of them. It would be great if you can >> contrast to OOP in Java. > > Java is not the end-all of OOP (in fact S is a good bit older than > Java) and you might find that the Lisp or Dylan object systems are a > better analogy. (I'm only going by hearsay on Dylan; never used it > myself) You might also quickly breeze through: > https://github.com/hadley/devtools/wiki/S3 > >> >> 1. R's S4 appears to centered around a dispatch mechanism which in my >> understanding is just a way to implement polymorphism. Now, here's the snag, >> I thought polymorphism was an aspect of OOP not by itself the definition of >> OOP. What am I missing here? Is any language that implements polymorphism >> automatically OO? >> > > If you accept the immutability of objects, then arguably yes, I > suppose polymorphism gives you a great deal of it. The remaining > weaknesses are generally addressed by the S4 object system. > > (Not immutability of bindings like Haskell, but the fact that x <- y > <- 1:5; y[3] <- 10 won't change x. In theory this is done by creating > a new y with the modified 3rd element and binding the name y to that; > not entirely thus in practice for performance reasons ) > >> 2. Can someone provide a simple example of how NextMethod() works? I read >> some things about but I can't make any sense out of it. >> It's supposed to facilitate inheritance but how? Why is it needed, what >> happens if it's ignored? An example would be useful. Is there a Java >> equivalent of NextMethod()? > > Grepping through R's source, it seems that the print system uses a > fair amount of NextMethod for the AsIs and noquote print methods. You > might take a look at those: also, section 7 of > http://cran.r-project.org/doc/manuals/r-devel/R-exts.html > > MW > > >> >> Many Thanks! >> >> __ >> 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@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] From Java to R OOP
Hi, I'm new to OOP in R so please forgive the naiveness of some of the questions. Here are a couple of them. It would be great if you can contrast to OOP in Java. 1. R's S4 appears to centered around a dispatch mechanism which in my understanding is just a way to implement polymorphism. Now, here's the snag, I thought polymorphism was an aspect of OOP not by itself the definition of OOP. What am I missing here? Is any language that implements polymorphism automatically OO? 2. Can someone provide a simple example of how NextMethod() works? I read some things about but I can't make any sense out of it. It's supposed to facilitate inheritance but how? Why is it needed, what happens if it's ignored? An example would be useful. Is there a Java equivalent of NextMethod()? Many Thanks! __ 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] Fwd: How to conditionally remove dataframe rows?
Hi, I have a data frame with two columns. I need to remove duplicated rows in first column, but I need to do it conditionally to values of the second column. Example: Point_counts Psi_Sp 1A 0 2A 1 3B 1 4B 2 5B 0 6C 1 7D 1 8D 2 I need to turn this data frame in one without duplicated rows at point-counts (one visit per point) but maintain the ones with maximum value at Psi_Sp, e.g. remove row 1 and maintain 2 or remove rows 3 and 5 and maintain 4. At the end I want a data frame like the one below: Point_counts Psi_Sp 1 A 1 2 B 2 3 C 0 4 D 2 How can I do it? I found several ways to edit data frames, but unfortunately I cound not use none of them. I appreciate Francisco [[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] Multistate occupancy models using Jags
Hi everyone, I am trying to run the example of multistate occupancy model from the book Bayesian Population Analysis using WinBUGS (Marc Kéry and Michael Schaub): This example is available on http://www.vogelwarte.ch/code-for-running-bpa-using-jags.html When I try to run the first line of this section of the model : # Initial values zst <- apply(y, 1, max, na.rm = TRUE) zst[zst == "-Inf"] <- 1 inits <- function(){list(z = zst)} I receive this warning message from R: 1: In FUN(newX[, i], ...) : no non-missing argument to max; returning -Inf (Original message in Portuguese: Mensagens de aviso perdidas: 1: In FUN(newX[, i], ...) : nenhum argumento não faltante para max; retornando -Inf ) Does someone know what is happening? I couldn't find anything related to this model and error. Thanks Francisco [[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] RWeka - weighted instances
Is there any way to train Weka classifiers using weighted instances in RWeka? I have done it in Weka in Java, but I don't see any way to pass the weights to the classifier in the RWeka documentation. Thanks __ 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] Problems with getURL (RCurl) to obtain list files of an ftp directory
Dear all, I have a problem with the command 'getURL' from the RCurl package, which I have been using to obtain a ftp directory list from the MOD16 (ET, DSI) products, and then to download them. (part of the script by Tomislav Hengl, spatial-analyst). Instead of the list of files (from ftp), I am getting the complete html code. Anyone knows why this might happen? This are the steps i have been doing: > MOD16A2.doy<- ' ftp://ftp.ntsg.umt.edu/pub/MODIS/Mirror/MOD16/MOD16A2.105_MERRAGMAO/' > items <- strsplit(getURL(MOD16A2.doy, .opts=curlOptions(ftplistonly=TRUE)), "\n")[[1]] >items #results [1] "http://www.w3.org/TR/html4/loose.dtd\";>\n\n\n\nFTP Directory: ftp://ftp.ntsg.umt.edu/pub/MODIS/Mirror/MOD16/MOD16A2.105_MERRAGMAO/\n\n<!--BODY{background-color:#ff;font-family:verdana,sans-serif}-->\n\n\nFTP Directory: ftp://ftp.ntsg.umt.edu/pub/MODIS/Mirror/MOD16/MOD16A2.105_MERRAGMAO/\n\nhttp://localhost:3128/squid-internal-static/icons/anthony-dirup.gif\"; ALT=\"[DIRUP]\"> Parent Directory \nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> GEOTIFF_0.05degree . . . . . . . Jun 3 18:00\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> GEOTIFF_0.5degree. . . . . . . . Jun 3 18:01\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2000. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2001. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2002. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2003. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2004. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2005. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2006. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2007. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2008. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2009. . . . . . . . . . . . . . Dec 23 2010\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2010. . . . . . . . . . . . . . Feb 20 2011\nhttp://localhost:3128/squid-internal-static/icons/anthony-dir.gif\"; ALT=\"[DIR] \"> Y2011. . . . . . . . . . . . . . Mar 12 2012\n\n\n\nGenerated Wed, 10 Oct 2012 13:43:53 GMT by localhost (squid/2.7.STABLE9)\n\n" The curious is that the command getURL was working well until I don't know what happened. And using the same command in Windows works fine. The sessionInfo() have given me the next: R version 2.14.1 (2011-12-22) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] MODIS_0.5-8 maptools_0.8-16 lattice_0.20-0 foreign_0.8-48 date_1.2-32 [6] RCurl_1.95-0.1 bitops_1.0-4.1 rgdal_0.7-19raster_2.0-12 sp_0.9-99 loaded via a namespace (and not attached): [1] grid_2.14.1 tools_2.14.1 Kind regard for all Francisco Zambrano Bigiarini INIA Quilamapu, Chillán, *Chile* [[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] bootstrapped CI for nonlinear models using nlsBoot from nlstools
Hi all I´m trying to get confidence intervals for parameters from nls modeling. I fitted a nls model to the following variables: > x [1] 2 1 1 5 4 6 13 11 13 101 101 101 > y [1] 1.281055090 1.563609934 0.001570796 2.291579783 0.841891853 [6] 6.553951324 14.243274230 14.519899320 15.066473610 21.728809880 [11] 18.553054450 23.722637370 The model fitted was: model<-nls(y~SSgompertz(x,a,b,c)) and it worked OK, with the following results: Formula: y ~ SSgompertz(x, a, b, c) Parameters: Estimate Std. Error t value Pr(>|t|) a 21.2426 0.9689 21.924 4.03e-09 *** b 5.3330 1.4722 3.622 0.00555 ** c 0.8045 0.0274 29.364 3.01e-10 *** Then, trying to get confidence intervals for the parameters using the nlsBoot function in the nlstools package I got this error: > modelboot<-nlsBoot(model) Error en data2[, var1] <- fitted1 + sample(scale(resid1, scale = FALSE), : objeto de tipo 'environment' no es subconjunto I´ve tried with another response variable and other self starting function (SSlogis) and got the same error. Any suggestions? Francisco -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] error using boxcox.nls during non linear estimation
Hi all I´ve got a problem using boxcox.nls function in nlrwr packagge. I´m fitting several non linear models to these data: > x [1] 2 1 1 5 4 6 13 11 13 101 101 101 > y [1] 1.281055090 1.563609934 0.001570796 2.291579783 0.841891853 [6] 6.553951324 14.243274230 14.519899320 15.066473610 21.728809880 [11] 18.553054450 23.722637370 I used nls function with self starters, for example SSgompertz: model<-nls(y~SSgompertz(x,a,b)) The model runs OK, with the following summary: Formula: y ~ SSgompertz(x, a, b, c) Parameters: Estimate Std. Error t value Pr(>|t|) a 21.2426 0.9689 21.924 4.03e-09 *** b 5.3330 1.4722 3.622 0.00555 ** c 0.8045 0.0274 29.364 3.01e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.686 on 9 degrees of freedom Number of iterations to convergence: 5 Achieved convergence tolerance: 1.508e-06 But then, trying to use boxcox.nls(model) gives the following error: > boxcox.nls(model) Error in eval(object$data)[, as.character(formula(object)[[2]])[2]] : object type 'environment' is not a subset (this last part traduced from spanish) I´ve tryed other functions, even without self starters but it does not work. Then it seems to be there something in the data, but I don´t know what. Any suggestions? Francisco ------ Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] scale_y_logit not present in ggplot2 0.9?
I have updated from ggplot2 0.8.9 on Windows to ggplot 0.9.0 (and then 0.9.1) and I can't find scale_y_logit() anymore ... On Mac, I can't see scale_y_logit() in ggplot 0.9.0 either. Am I missing anything? Thanks __ 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] confidence intervals for nls or nls2 model
Thanks! Now it is clear. Francisco On Wed, 16 May 2012 07:32:56 -0400, Gabor Grothendieck wrote > On Tue, May 15, 2012 at 11:20 PM, Gabor Grothendieck > wrote: > > On Tue, May 15, 2012 at 8:08 PM, Francisco Mora Ardila > > wrote: > >> Hi all > >> > >> I have fitted a model usinf nls function to these data: > >> > >>> x > >> [1] 1 0 0 4 3 5 12 10 12 100 100 100 > >> > >>> y > >> [1] 1.281055090 1.563609934 0.001570796 2.291579783 0.841891853 > >> [6] 6.553951324 14.243274230 14.519899320 15.066473610 21.728809880 > >> [11] 18.553054450 23.722637370 > >> > >> The model fitted is: > >> > >> modellogis<-nls(y~SSlogis(x,a,b,c)) > >> > >> It runs OK. Then I calculate confidence intervals for the actual data > >> using: > >> > >> dataci<-predict(as.lm(modellogis), interval = "confidence") > >> > >> BUt I don´t get smooth curves when plotting it, so I want to get other > >> "confidence > >> vectors" based on a new x vector by defining a new data to do predictions: > >> > >> x0 <- seq(0,15,1) > >> dataci<-predict(as.lm(modellogis), newdata=data.frame(x=x0), interval = "confidence") > >> > >> BUt it does not work: I get the same initial confidence interval > >> > >> Any ideas on how to get tconfidence and prediction intervals using new X > >> data on a > >> previous model? > >> > > > > as.lm is a linear model between the response variable and the gradient > > of the nonlinear model and as we see below x is not part of that > > linear model so x can't be in newdata when predicting from the tangent > > model. We can only make predictions at the original x points. For > > other x's we could use Interpolation. See ?approx (?spline can also > > work in smooth cases but in the example provided the function has a > > kink and that won't work well with splines.) > > > >> as.lm(modellogis)$model > > y a b c (offset) > > 1 1.281055090 0.06601796 -4.411829e-01 1.168928e+00 1.397153 > > 2 1.563609934 0.04798815 -3.268846e-01 9.766080e-01 1.015584 > > 3 0.001570796 0.04798815 -3.268846e-01 9.766080e-01 1.015584 > > 4 2.291579783 0.16311227 -9.767241e-01 1.597189e+00 3.451981 > > 5 0.841891853 0.12203013 -7.665928e-01 1.512752e+00 2.582551 > > 6 6.553951324 0.21464369 -1.206154e+00 1.564573e+00 4.542552 > > 7 14.243274230 0.74450055 -1.361047e+00 -1.455630e+00 15.756031 > > 8 14.519899320 0.59707858 -1.721353e+00 -6.770205e-01 12.636107 > > 9 15.066473610 0.74450055 -1.361047e+00 -1.455630e+00 15.756031 > > 10 21.728809880 1. -2.943955e-13 -9.073765e-12 21.163223 > > 11 18.553054450 1. -2.943955e-13 -9.073765e-12 21.163223 > > 12 23.722637370 1. -2.943955e-13 -9.073765e-12 21.163223 > > > > I have added a FAQ to the home page since this isn't the first time > this question has come up: > > http://nls2.googlecode.com#FAQs > > -- > Statistics & Software Consulting > GKX Group, GKX Associates Inc. > tel: 1-877-GKX-GROUP > email: ggrothendieck at gmail.com -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] confidence intervals for nls or nls2 model
On Tue, 15 May 2012 20:33:02 -0400, David Winsemius wrote > On May 15, 2012, at 8:08 PM, Francisco Mora Ardila wrote: > > > Hi all > > > > I have fitted a model usinf nls function to these data: > > > >> x > > [1] 1 0 0 4 3 5 12 10 12 100 100 100 > > > >> y > > [1] 1.281055090 1.563609934 0.001570796 2.291579783 0.841891853 > > [6] 6.553951324 14.243274230 14.519899320 15.066473610 21.728809880 > > [11] 18.553054450 23.722637370 > > > > The model fitted is: > > > > modellogis<-nls(y~SSlogis(x,a,b,c)) > > > > It runs OK. Then I calculate confidence intervals for the actual > > data using: > > > > dataci<-predict(as.lm(modellogis), interval = "confidence") > > > > BUt I don´t get smooth curves when plotting it, so I want to get > > other "confidence > > vectors" based on a new x vector by defining a new data to do > > predictions: > > > > x0 <- seq(0,15,1) > > dataci<-predict(as.lm(modellogis), newdata=data.frame(x=x0), > > interval = "confidence") > > > > BUt it does not work: > > I am really getting tired of seeing the phrase "doesn't work". And on > my machine it throws an error saying there is no as.lm function> This > seemes to give "smooth results: Sorry for that "doesn't work". Also, I forgot to say function as.lm() is in tha package nls2, which also fit nls models. > > dataci<-predict(as.lm(modellogis), newdata=data.frame(x=x0), interval > = "confidence") > > And notice that this message is in ?predict.nls > > interval > "At present this argument is ignored. Yeah, that´s why I´m using nls2 package > > > I get the same initial confidence interval > > > > Any ideas on how to get tconfidence and prediction intervals using > > new X data on a > > previous model? > > > > Thanks > > > > Francisco > > > > > > > > -- > > Francisco Mora Ardila > > Estudiante de Doctorado > > Centro de Investigaciones en Ecosistemas > > Universidad Nacional Autónoma de México > > > > __ > > 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. > > David Winsemius, MD > West Hartford, CT -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] confidence intervals for nls or nls2 model
Hi all I have fitted a model usinf nls function to these data: > x [1] 1 0 0 4 3 5 12 10 12 100 100 100 > y [1] 1.281055090 1.563609934 0.001570796 2.291579783 0.841891853 [6] 6.553951324 14.243274230 14.519899320 15.066473610 21.728809880 [11] 18.553054450 23.722637370 The model fitted is: modellogis<-nls(y~SSlogis(x,a,b,c)) It runs OK. Then I calculate confidence intervals for the actual data using: dataci<-predict(as.lm(modellogis), interval = "confidence") BUt I don´t get smooth curves when plotting it, so I want to get other "confidence vectors" based on a new x vector by defining a new data to do predictions: x0 <- seq(0,15,1) dataci<-predict(as.lm(modellogis), newdata=data.frame(x=x0), interval = "confidence") BUt it does not work: I get the same initial confidence interval Any ideas on how to get tconfidence and prediction intervals using new X data on a previous model? Thanks Francisco -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] Cook's distance for lme?
Hi Is there any function that I can use to calculate Cook's distance for an lme? Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk<mailto:beatriz.defranci...@sams.ac.uk> http://www.smi.ac.uk/beatriz-de-franciso The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. [[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] How can I brake a label in two lines when using expression()?
I making an xyplot and the y label is too long and needs to be in two rows, but when I brake it there is a huge gap between the last text string and the expression, and I can't get rid of it. Any ideas? Data: structure(list(Temp = c(8L, 8L, 8L, 8L, 8L, 8L, 12L, 12L, 12L, 12L, 12L, 12L), CO2 = c(380L, 380L, 380L, 750L, 750L, 750L, 380L, 380L, 380L, 750L, 750L, 750L), Treat = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("12-380", "12-750", "8-380", "8-750"), class = "factor"), Week = c(1L, 3L, 8L, 1L, 3L, 8L, 1L, 3L, 8L, 1L, 3L, 8L), Mean.Rate = c(2.125909389, 1.905870003, 1.417687602, 3.110439984, 2.31043989, 1.849232493, 2.546747098, 3.290235064, 3.000717599, 2.694901409, 3.852590547, 2.964084249 ), lower = c(1.846641409, 1.44072624, 1.185304427, 2.56408099, 2.02644683, 1.606374443, 2.253928482, 2.759177284, 2.49014747, 2.168437604, 3.075977559, 2.438453415), upper = c(2.405177369, 2.371013766, 1.650070777, 3.656798978, 2.59443295, 2.092090543, 2.839565714, 3.821292844, 3.511287728, 3.221365214, 4.629203535, 3.489715083), fTemp = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("8", "12"), class = "factor"), fCO2 = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("380", "750"), class = "factor"), fTreat = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("8-380", "8-750", "12-380", "12-750"), class = c("ordered", "factor" )), fWeek = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", "3", "8"), class = "factor")), .Names = c("Temp", "CO2", "Treat", "Week", "Mean.Rate", "lower", "upper", "fTemp", "fCO2", "fTreat", "fWeek"), row.names = c(NA, -12L), class = "data.frame") xyplot(cbind(Mean.Rate,lower,upper)~fWeek|fTreat, resp.week.mean.rate, as.table=TRUE, xlab="Week", ylab=expression("Mean Reapisration Rate (umol."*L^-1*".g (AFDM)"^-1*")"), scales=list(alternating=FALSE, tick.number=10, tck=c(-1,0)), layout=c(4,1), ylim=1:5, auto.key=list(title="Treatment", lines=TRUE, cex.title=1, columns=2), panel=function(x, y,...){ panel.errbars(x,y,make.grid="none",ewidth=0.2,type="p",...) panel.loess(x[resp.week.mean.rate$Treat=="8-380"],y[resp.week.mean.rate$Treat=="8-380"],span = 5, degree = 1,lwd=2,...) panel.loess(x[resp.week.mean.rate$Treat=="8-750"],y[resp.week.mean.rate$Treat=="8-750"],span = 5, degree = 1,lwd=2,...); panel.loess(x[resp.week.mean.rate$Treat=="12-380"],y[resp.week.mean.rate$Treat=="12-380"],span = 5, degree = 1,lwd=2,...); panel.loess(x[resp.week.mean.rate$Treat=="12-750"],y[resp.week.mean.rate$Treat=="12-750"],span = 5, degree = 1,lwd=2,...) } ) Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk<mailto:beatriz.defranci...@sams.ac.uk> http://www.smi.ac.uk/beatriz-de-franciso The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. [[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] braking a label in two lines when using expression()
I making an xyplot and the y label is too long and needs to be in two rows, but when I brake it there is a huge gap between the last text string and the expression, and I can't get rid of it. Any ideas? Data: structure(list(Temp = c(8L, 8L, 8L, 8L, 8L, 8L, 12L, 12L, 12L, 12L, 12L, 12L), CO2 = c(380L, 380L, 380L, 750L, 750L, 750L, 380L, 380L, 380L, 750L, 750L, 750L), Treat = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("12-380", "12-750", "8-380", "8-750"), class = "factor"), Week = c(1L, 3L, 8L, 1L, 3L, 8L, 1L, 3L, 8L, 1L, 3L, 8L), Mean.Rate = c(2.125909389, 1.905870003, 1.417687602, 3.110439984, 2.31043989, 1.849232493, 2.546747098, 3.290235064, 3.000717599, 2.694901409, 3.852590547, 2.964084249 ), lower = c(1.846641409, 1.44072624, 1.185304427, 2.56408099, 2.02644683, 1.606374443, 2.253928482, 2.759177284, 2.49014747, 2.168437604, 3.075977559, 2.438453415), upper = c(2.405177369, 2.371013766, 1.650070777, 3.656798978, 2.59443295, 2.092090543, 2.839565714, 3.821292844, 3.511287728, 3.221365214, 4.629203535, 3.489715083), fTemp = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("8", "12"), class = "factor"), fCO2 = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("380", "750"), class = "factor"), fTreat = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("8-380", "8-750", "12-380", "12-750"), class = c("ordered", "factor" )), fWeek = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", "3", "8"), class = "factor")), .Names = c("Temp", "CO2", "Treat", "Week", "Mean.Rate", "lower", "upper", "fTemp", "fCO2", "fTreat", "fWeek"), row.names = c(NA, -12L), class = "data.frame") xyplot(cbind(Mean.Rate,lower,upper)~fWeek|fTreat, resp.week.mean.rate, as.table=TRUE, xlab="Week", ylab=expression("Mean Reapisration Rate (umol."*L^-1*".g (AFDM)"^-1*")"), scales=list(alternating=FALSE, tick.number=10, tck=c(-1,0)), layout=c(4,1), ylim=1:5, auto.key=list(title="Treatment", lines=TRUE, cex.title=1, columns=2), panel=function(x, y,...){ panel.errbars(x,y,make.grid="none",ewidth=0.2,type="p",...) panel.loess(x[resp.week.mean.rate$Treat=="8-380"],y[resp.week.mean.rate$Treat=="8-380"],span = 5, degree = 1,lwd=2,...) panel.loess(x[resp.week.mean.rate$Treat=="8-750"],y[resp.week.mean.rate$Treat=="8-750"],span = 5, degree = 1,lwd=2,...); panel.loess(x[resp.week.mean.rate$Treat=="12-380"],y[resp.week.mean.rate$Treat=="12-380"],span = 5, degree = 1,lwd=2,...); panel.loess(x[resp.week.mean.rate$Treat=="12-750"],y[resp.week.mean.rate$Treat=="12-750"],span = 5, degree = 1,lwd=2,...) } ) Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk<mailto:beatriz.defranci...@sams.ac.uk> http://www.smi.ac.uk/beatriz-de-franciso The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. [[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] error bars for a barchart
Walmes, Thank you so much!!! I am still trying to understand all of your code but it works. I have changed it a bit so that I get upper and lower limits for the error bar, and that the origin starts at 0 so the negative values are plotted correctly. barchart(Change~fTreat,groups=Process,change, stderr=change$stderr, ylab="Pocertage change", xlab="Treatment", #ylim=-115:50, scales=list(alternating=FALSE, tick.number=7, tck=c(-1,0)), prepanel=function(y, stderr, subscripts=subscripts, ...){ uy <- as.numeric(y+stderr[subscripts]) ly <- as.numeric(y-stderr[subscripts]) list(ylim=range(y,uy,ly, finite=TRUE)) }, panel= function(x, y, subscripts, groups, stderr, box.ratio, ...){ panel.barchart(x, y, subscripts=subscripts, groups=groups, box.ratio=box.ratio,origin=0, ...) panel.abline(h=0,col="black",...) d <- 1/(nlevels(groups)+nlevels(groups)/box.ratio) g <- (as.numeric(groups[subscripts])-1); g <- (g-median(g))*d panel.arrows(as.numeric(x)+g,y-stderr[subscripts], as.numeric(x)+g, y+stderr[subscripts], code=3,angle=90, length=0.025) } ) I am very new to creating function and would be great if you could explain what the d and g elemens actually do? this is just for me to understand and later maybe make my own functions. I am assuming that g centers the error bars? but d? Regards Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk http://www.smi.ac.uk/beatriz-de-franciso From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of ilai [ke...@math.montana.edu] Sent: 02 May 2012 04:14 To: Walmes Zeviani Cc: r-help@r-project.org Subject: Re: [R] error bars for a barchart Thank you for your example. I only skimmed it, but since both solutions use nlevels and box.ratio it is no surprise we end up at the same place (although I do think your g-median is nicer than my 3/4). Thing is, I wouldn't call either of these "simple"... would be nice if one could just query the "new" centers, but I don't know if there is a way without hacking panel.barchart itself ? Cheers On Tue, May 1, 2012 at 1:34 PM, Walmes Zeviani wrote: > I have a repoducibe example here > > http://ridiculas.wordpress.com/2011/11/23/media-e-desvio-padrao-de-muitas-variaveis-separado-por-grupos/ > > Sorry for it be in Portuguese. > > Walmes. > > == > Walmes Marques Zeviani > LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) > Departamento de Estatística - Universidade Federal do Paraná > fone: (+55) 41 3361 3573 > VoIP: (3361 3600) 1053 1173 > e-mail: wal...@ufpr.br > twitter: @walmeszeviani > homepage: http://www.leg.ufpr.br/~walmes > linux user number: 531218 > == __ 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. The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. __ 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] error bars for a barchart
Hi I have the following barchart to which I want to add error bars. library(lattice) barchart(Change~fTreat,groups=Process,change, auto.key=list(points=FALSE,rectangles=TRUE), panel=function(x, y,...){ panel.barchart(x,y,origin = 0,...); panel.abline(h=0,col="black",...); } ) I have tried using the panel.errbars from the memisc package which works great for xyplots, but when I add it to my code it does not respect the groups. library(memisc) barchart(cbind(Change,lower,upper)~fTreat,groups=Process,change, ylab="Pocertage change", ylim=-115:50, scales=list(alternating=FALSE, tick.number=7, tck=c(-1,0)), panel=function(x, y,groups,...){ panel.barchart(x,y=change$Change,groups=change$Process,origin = 0,...); panel.abline(h=0,col="black",...); panel.errbars(x,y,make.grid="none",ewidth=0.2,type="n",...) } ) Any ideas of how to add error bars to my plot either using the panel.errbars or any other function? The data: structure(list(Treat = structure(c(3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L), .Label = c("12-380", "12-750", "8-380", "8-750"), class = "factor"), Process = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("Resp", "Cal"), class = c("ordered", "factor")), Change = c(-33.05, -34.74, 20.94, 18.06, 6.85, -28.57, -8.1, -78.72), upper = c(-13.22896628, -28.61149669, 31.29930461, 27.30173776, 39.73271282, 9.458372948, 13.11035572, -47.03745704), lower = c(-52.86120694, -40.87446411, 10.57421563, 8.822042178, -26.03144161, -66.60447035, -29.30563327, -110.3973761), fTreat = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("8-380", "8-750", "12-380", "12-750"), class = c("ordered", "factor"))), .Names = c("Treat", "Process", "Change", "upper", "lower", "fTreat"), row.names = c(NA, -8L), class = "data.frame") Regards Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk<mailto:beatriz.defranci...@sams.ac.uk> http://www.smi.ac.uk/beatriz-de-franciso The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. [[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] error bars for a barchart
Hi I have the following barchart to which I want to add error bars. library(lattice) barchart(Change~fTreat,groups=Process,change, auto.key=list(points=FALSE,rectangles=TRUE), panel=function(x, y,...){ panel.barchart(x,y,origin = 0,...); panel.abline(h=0,col="black",...); } ) I have tried using the panel.errbars from the memisc package which works great for xyplots, but when I add it to my code it does not respect the groups. library(memisc) barchart(cbind(Change,lower,upper)~fTreat,groups=Process,change, ylab="Pocertage change", ylim=-115:50, scales=list(alternating=FALSE, tick.number=7, tck=c(-1,0)), panel=function(x, y,groups,...){ panel.barchart(x,y=change$Change,groups=change$Process,origin = 0,...); panel.abline(h=0,col="black",...); panel.errbars(x,y,make.grid="none",ewidth=0.2,type="n",...) } ) Any ideas of how to add error bars to my plot either using the panel.errbars or any other function? The data: structure(list(Treat = structure(c(3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L), .Label = c("12-380", "12-750", "8-380", "8-750"), class = "factor"), Process = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("Resp", "Cal"), class = c("ordered", "factor")), Change = c(-33.05, -34.74, 20.94, 18.06, 6.85, -28.57, -8.1, -78.72), upper = c(-13.22896628, -28.61149669, 31.29930461, 27.30173776, 39.73271282, 9.458372948, 13.11035572, -47.03745704), lower = c(-52.86120694, -40.87446411, 10.57421563, 8.822042178, -26.03144161, -66.60447035, -29.30563327, -110.3973761), fTreat = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("8-380", "8-750", "12-380", "12-750"), class = c("ordered", "factor"))), .Names = c("Treat", "Process", "Change", "upper", "lower", "fTreat"), row.names = c(NA, -8L), class = "data.frame") Regards Beatriz de Francisco Mora PhD Student The Scottish Association for Marine Science Scottish Marine Institute Oban PA37 1QA Tel: 06131 559000 (switchboard) Fax: 01631559001 E. beatriz.defranci...@sams.ac.uk<mailto:beatriz.defranci...@sams.ac.uk> http://www.smi.ac.uk/beatriz-de-franciso The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. [[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] error using nls with logistic derivative
It seems to work! Thanks and sorry for the personal message Francisco On Tue, 17 Apr 2012 09:09:24 +0200, peter dalgaard wrote > On Apr 17, 2012, at 06:23 , Francisco Mora Ardila wrote: > > > Hi > > > > I´m trying to fit a nonlinear model to a derivative of the logistic function > > > > y = a/(1+exp((b-x)/c)) (this is the parametrization for the SSlogis > > function with nls) > > > > The derivative calculated with D function is: > > > >> logis<- expression(a/(1+exp((b-x)/c))) > >> D(logis, "x") > > a * (exp((b - x)/c) * (1/c))/(1 + exp((b - x)/c))^2 > > > > So I enter this expression in the nls function: > > > > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2, > > start=list(a = 21.16322, b = 8.83669, c = 2.957765), > > ) > > > > The data is: > > > >> Y > > [1] 5.5199668 1.5234525 3.3557000 6.7211704 7.4237955 1.9703127 > > [7] 4.3939336 -1.4380091 3.2650180 3.5760906 0.2947972 1.0569417 > >> X > > [1] 1 0 0 4 3 5 12 10 12 100 100 100 > > > > The problem is that I got the next error: > > > > Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2, : > > step factor 0.000488281 reduced below 'minFactor' of 0.000976563 > > > > I trien to change the minFactor using the control argument inside nls > > > > control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048 > > > > but got a new error message: > > > > > > Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2, : > > step factor 0.000244141 reduced below 'minFactor' of 0.000488281 > > > > So it seems that as I modify minFactor, the step factor reduces also and I > > can never > > reach a solution. > > > > Does anybody Know what am I doing wrong? Is there a problem with the > > formula? How can I > > solve it? I tried some suggestions on R-help related topics but did not > > work. > > (Please don't send private messages. I don't do free consulting outside of > the > mailing lists.) > > With nls(), there's really no substitute for good data and good starting > values. Did you actually try plotting those data? At best, they are extremely > noisy. How did you obtain the starting values? They seem remarkably accurate > for something that fits data so poorly. > > What is happening is that the algorithm shoots off into a region of the > parameter space where it can't make any progress: > > > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,+ > > start=list(a = 21.16322, b = 8.83669, c = 2.957765),+ trace=TRUE) > 151.098 : 21.163220 8.836690 2.957765 > 127.1149 : 103.49326 11.43274 20.29663 > 111.344 : 833.02386 -45.86244 140.32985 > 111.3375 : 978.97214 -76.20571 159.90833 > 111.3374 : 1097.1376 -101.6771 174.2037 > 111.3228 : 1179.8451 -119.7416 183.3794 > > Notice that the "b" parameter which is supposed to be the maximum point > starts > off well to the right of the position of the largest Y values, then shoots > into large negative values. > > With a little eyeballing, I can do better: > > > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2, > + start=list(a = 40, b = 3.5, c = 10),trace=T) > 139.5173 : 40.0 3.5 10.0 > 97.28614 : 26.513424 4.052639 2.267709 > 81.53127 : 32.384910 3.473411 2.307504 > 65.11387 : 39.53542 3.01097 2.07678 > 50.66328 : 36.223529 2.590102 1.133965 > 50.50921 : 35.984466 2.565698 1.067731 > 50.50162 : 36.149993 2.573420 1.081673 > 50.50145 : 36.129962 2.572195 1.079504 > 50.50144 : 36.133656 2.572402 1.079862 > 50.50144 : 36.133061 2.572368 1.079803 > > However, the fit isn't exactly stellar. Try this: > > X <- c(1, 0, 0, 4, 3, 5, 12, 10, 12, 100, 100, 100) > Y <- c(5.5199668, 1.5234525, 3.3557, 6.7211704, 7.4237955, 1.9703127, >4.3939336, -1.4380091, 3.265018, 3.5760906, 0.2947972, 1.0569417) > plot(X,Y) > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2, > start=list(a = 40, b = 3.5, c = 10),trace=T) > x0 <- seq(0,100,,1000) > lines(x0,predict(ratelogis,newdata=data.frame(X=x0))) > lines(x0,evalq(a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,envir=list(a = > 21.16322, b = 8.83669, c = 2.957765, X=x0)), lty="dashed") > > -- > Peter Dalgaard, Professor, > Center for Statistics, Copenhagen Business School > Solbjerg Plads 3, 2000 Frederiksberg, Denmark > Phone: (+45)38153501 > Email: pd@cbs.dk Priv: pda...@gmail.com -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] error using nls with logistic derivative
Hi I´m trying to fit a nonlinear model to a derivative of the logistic function y = a/(1+exp((b-x)/c)) (this is the parametrization for the SSlogis function with nls) The derivative calculated with D function is: > logis<- expression(a/(1+exp((b-x)/c))) > D(logis, "x") a * (exp((b - x)/c) * (1/c))/(1 + exp((b - x)/c))^2 So I enter this expression in the nls function: ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2, start=list(a = 21.16322, b = 8.83669, c = 2.957765), ) The data is: > Y [1] 5.5199668 1.5234525 3.3557000 6.7211704 7.4237955 1.9703127 [7] 4.3939336 -1.4380091 3.2650180 3.5760906 0.2947972 1.0569417 > X [1] 1 0 0 4 3 5 12 10 12 100 100 100 The problem is that I got the next error: Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2, : step factor 0.000488281 reduced below 'minFactor' of 0.000976563 I trien to change the minFactor using the control argument inside nls control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048 but got a new error message: Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2, : step factor 0.000244141 reduced below 'minFactor' of 0.000488281 So it seems that as I modify minFactor, the step factor reduces also and I can never reach a solution. Does anybody Know what am I doing wrong? Is there a problem with the formula? How can I solve it? I tried some suggestions on R-help related topics but did not work. Thanks Francisco -- Francisco Mora Ardila Estudiante de Doctorado Centro de Investigaciones en Ecosistemas Universidad Nacional Autónoma de México __ 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] phangorn and calculation of a rate matrix
Hi, I'm trying to calculate a ratematrix for a RNA aligment (U instead of T) in order to use it as a ratematrix in Profidst (a phylogenetic program which takes into account both the primary sequence and the secondary structure of the RNA, in my case rRNA). The sequence-structure aligment has been made in 4SALE (a java app) and saved as one-letter encoded (using a 12 letters alphabet, "a","c","d","e","g","h","i","k","l","n","q","r" instead of the conventional nucleotide codes). My intention is to calculate the ratematrix (which is a 12x12 matrix) for this special aligment with ape and phangorn, however I've repeatdly fail to do it. The aligment is in a phyDat object containing 30 sequences and 12 states (by using "user-defined" character states consisting on the 12 letters indicated above). I follow the steps described in the phangorn-specials vignette but the ratematrix (under a GTR substitution model) is not calculated. May the problem be that phangorn only accepts "a","g","c" and "t" as valid states for calculating the matrix? And in case phangorn could calculate the matrix, how could I do it? I have a very basic knowledge of R so please I would greatly appreciate (if possible) a step-by-step explanation. Thanks very much for the help. -- View this message in context: http://r.789695.n4.nabble.com/phangorn-and-calculation-of-a-rate-matrix-tp4550495p4550495.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] defining non linear predictors from nls in gam?
Hi all I´m trying to analize the role of time since abandonement (continuous variable) and biophysical environmental conditions on the recovery of the vegetation trough succession. First, I used non-linear least squares with nls function to model the effect of time on vegetation attributes. I tried several self-starting sigmoidal functions as data seems to conform to this type of models, and then chose the best model based on the minimum RSE and AIC. An example of the formula used is: model1<-nls(vegetattrib~SSlogis(time, a, b, c)) Then, I wanted to add to the model the effect of some biophysical attribute (ie, soil type). But three problems arise: 1) I can´t include factors in the nls, 2) even if it were a continuous predictor, nls model try to fit it as a nonlinear predictor and I don ´t have an apriori reason to think that kind of relation exist, 3) nls doesn’t give r2´s (the statistical reason for this is described by Douglas Bates and can be found in https://stat.ethz.ch/pipermail/r-help/2000-August/007778.html). But the point is that for me is interesting to have an idea about how well the model describe the patter in data, as the r2 does. Is correct to calculate the % of deviance instead? Or something else? So I decided to use a gam approach, were I can create an additive model with time and soil type. But gam creates a smoothing function for the relationship between time and vegetattr. The question is: Can I establish in gam the form of the relationship between time and vegetattr as, for example, a logistic relationship with the parameters estimated with the self starting nls function? I´ve revised the book from S.Wood about GAM´s in R, but haven´t find something like that. Any suggestions about how to model (and test) the effect of time as a nonlinear predictor plus other variables (preferably as linear predictors)? Thanks in advance Francisco -- oikos.unam.mx __ 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] subset select="variable with a list of names"
Hello, I would like to make a function which extracts a subset, from a dataset, with only the columns that I want (specifying their names). For example, having this matrix: > mydata<-matrix(c(22,1,3,2001,24,5,7,2002,26,7,8,2002,28,5,7,2003), byrow=TRUE, ncol=4, dimnames=list(c(1,2,3,4), c("age","day","month","year"))) > mydata age day month year 1 22 1 3 2001 2 24 5 7 2002 3 26 7 8 2002 4 28 5 7 2003 I would like to create a function like: x<-function(names) {subset(mydata, select=names) } So I can choose every time which columns select, i.e. when I call: x("age,day") it would returns: age day 1 22 1 2 24 5 3 26 7 4 28 5 Obviously it is not working, and I don't know how to do to fix it. Do you have any suggestion? Thank you very much __ 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] find a row identical to another
Hello, I have a dataset with many rows, starting from a row that I choose I would like to find the other rows in the dataset which are identical to this row (with the same values per each column) and assign them to a variable. How could I do? Thank you __ 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] drop columns whose rows are all 0
Hello, I have a dataset with 40 variables, some of them are always 0 (each row). I would like to make a subset containing only the columns which values are not all 0, but I don't know how to do it. I tried: for(cut_column in 1:40) { if(sum(dataset[,cut_column])!=0) { columns_useful<-c(columns_useful,dataset[cut_column]) } } sorted_dataset<-subset(dataset, select=columns_useful) But it doesn't work. Thank you Francisco __ 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] quantile type 1 perhaps?
Hello, I need to analyse some data coming from a questionnaire which have for each item a likert scale 1-5. I need to find the lowest scores in the distribution, and for this purpose I thought to use the quantile() function to identify the participants belonging to the 5% with lowest scores (who have a score < than quantile, am I right?). The problem is: which type of quantile should I calculate? I suppose type 1, since it is a discontinuous distribution, but I don't understand very well the differences between type 2 and 3. Thank you, Francisco __ 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] Word Wrap
Hello, I have to write a big sentence with cat() and I would like that R automatically adds a new line when it is needed (when the text arrives at the end of the window), the same as Windows Notepad does (Word Wrap). How could I do? Thank you F. __ 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] read.table as integer
Hello, I have a csv file with many variables, both characters and integers. I would like to load it on R and do some operations on integer variables, the problem is that R loads the entire dataset considering all variables as characters, instead I would like that R makes the distinction between the two types, because there are too many variables to do: x1<-as.integer(x1) x2<-as.integer(x2) x3<-as.integer(x3) ... I tried to specify read.table(... stringsAsFactors=FALSE) but it doesn't work. Thanks, Best Regards __ 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] testing significance of axis loadings from multivariate dudi.mix
Hi all I´m trying to tests the significance of loadings from a ordination of 46 variables (caategorical, ordinal and nominal). I used dudi.mix from ade4 for the ordination. A years ago Jari Oksanen wrote this script implementing Peres-Neto et al. 2003 (Ecology) bootstraping method: netoboot <- function (x, permutations=1000, ...) { pcnull <- princomp(x, cor = TRUE, ...) res <- pcnull$loadings out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { pc <- princomp(x[sample(N, replace=TRUE), ], cor = TRUE ...) pred <- predict(pc, newdata = x) r <- cor(pcnull$scores, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- pc$loadings[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } I tried to aply it to the case of dudi.mix instead of princomp this way: netoboot1<-function (x, permutations=1000,...) { dudinull <- dudi.mix(x, scannf = FALSE, nf = 3) res <- dudinull$c1 out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { dudi <- dudi.mix(x[sample(N, replace=TRUE), ], scannf = FALSE, nf = 3) pred <- predict(dudi, newdata = x) r <- cor(dudinull$li, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- dudi$c1[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } But a problem arised with the predict function: it doesn´t seem to work with an object from dudi.mix and I dont understand why. Can somebody tell me why? Any suggestions to modify the script or to use other method? Thanks in advance. Francisco Francisco Mora Ardila Laboratorio de Biodiversidad y Funcionamiento del Ecosistema Centro de Investigaciones en Ecosistemas UNAM-Campus Morelia Tel 3222777 ext. 42621 Morelia , MIchoacán, México. -- Open WebMail Project (http://openwebmail.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] Select some, but not all, variables stepwise
Hi Andre I don´t know if it will work, but I´ve tried the MuMIn package, were you can evaluate all possible models (usin for example AIC) at one time. Maybe you can focus on comparing those models which retain the explanators you want. Best wishes Francisco On Fri, 4 Nov 2011 13:06:09 +, Andre Easom wrote > Hi, > > I would like to fit a linear model where some but not all explanators are > chosen stepwise - ie I definitely want to include some terms, but others only > if they are deemed significant (by AIC or whatever other approach is > available) > . For example if I wanted to definitely include x1 and x2, but only include > z1 and z2 if they are significant, something like this: > > df <- data.frame(y=c(4,2,6,7,3,9,5,7,6,2), x1=c(2,3,4,0,5,8,8,1,1,2), > x2=c(0,0, > 0,0,1,1,0,0,0,1), z1=c(0,1,0,0,0,1,1,0,1,1), z2=c(1,1,1,0,0,1,1,1,1,0)) model > <- lm(y ~ x1 + x2 + stepwise(z1 + z2)) > > Any help would be appreciated. > > Cheers, > Andre > ** > This email and any attachments are confidential, protect...{{dropped:22}} > > __ > 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. Francisco Mora Ardila Laboratorio de Biodiversidad y Funcionamiento del Ecosistema Centro de Investigaciones en Ecosistemas UNAM-Campus Morelia Tel 3222777 ext. 42621 Morelia , MIchoacán, México. -- Open WebMail Project (http://openwebmail.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.
[R] Inputs from vectors in cor.test
Hello, I have two numeric vectors in R and used cor.test function with them. Is it possible to know the input of a particular row of the vectors to the total correlation value and significance? Of course I just could take out that row from the vectors and run the test again to see how the results vary but I wonder if that is the right way to do it. Thanks for your time. Regards, ---sram __ 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] Per sites contributions from a correlation test
Hello, I have two numeric vectors in R and used cor.test function with them. Is it possible in R to know how much contributed a particular row of the vectors to the total correlation value and significance? Of course I just could take out that row from the vectors and run the test again to see how the results vary but I wonder if that is the right way to do it. Thanks for your time. Regards, ---sram __ 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] A very big matrix inside a function
Dear listers: As a part of a intermediate process, I need to use and modify a very big matrix (some 3x3) inside the body of a function. If the matrix is defined in the function, R shows a error message "Cannot allocate vector of size 6.7 Gb". But if I define the matrix before the function is called (as soon as the dimension can be calculated) R can allocate it. The problem is that although my function can use the matrix, cannot modify it. The schema is the following: # OPTION 1 # Function definition Myfunction<-function(parameters, ...) { ... Verybigmatrix<-matrix(0, n, n) # Matrix definition inside the function ... } # Main process ... Value.returned<-Myfunction(parameters, ...) # R issues an error message " Cannot allocate vector of size 6.7 Gb" ... # OPTION 2 # Function definition... Myfunction<-function(parameters, ...) { ... Verybigmatrix[index1, index2]<<-newvalue # Issues an error message " Cannot allocate vector of size 6.7 Gb" # It assignation is made with `<-` error is issued too ... } # Main process ... Verybigmatrix<-matrix(0, n, n) # Definition before Myfunction is called doesn't issue an error ... Value.returned<-Myfunction(parameters, ...) # R sees Verybigmatrix but if I try to modify it issues an error message " Cannot allocate vector of size 6.7 Gb" ... Characteristics of the machine I'm working on: Operating system: Redhat Enterprise Linux 6 Kernel and CPU: Linux 2.6.32-131.6.1.el6.x86_64 on x86_64 Processor information: Intel(R) Xeon(TM) CPU 3.60GHz, 4 cores Memory: 8 GB Thanks in advance for your help. Francisco __ 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] Bar Chart
Hi, It's difficult to know what is going wrong from what you say below (please include some reproducible code in the future as indicated in the posting guide). If you want to produce a bar chart of the numbers with the corresponding names as labels for these numbers, you can do something like this: x <- rnorm(10) y <- letters[1:10] # use the first 10 letters of the alphabet as the labels barplot(x, names.arg = y) HTH, Francisco On Wed, Mar 23, 2011 at 4:26 PM, blutack wrote: > How do you do a bar chart of 2 vectors? > I have one vector which has 10 numbers, and another which has 10 names. > The numbers are the frequency of the corresponding name, but when I do a > bar > chart it says that there is no height. Thanks. > > -- > View this message in context: > http://r.789695.n4.nabble.com/Bar-Chart-tp3399924p3399924.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] create data set from selection of rows
Hi, What you are after is: datasubset <- dataset[ dataset[,3] == "text3", ] Equivalently, you can use: datasubset <- subset(dataset, subset = dataset[,3] == "text3") HTH, Francisco On Tue, Mar 15, 2011 at 3:09 PM, e-letter wrote: > Readers, > > For a data set: > > text1,23,text2,45 > text1,23,text3,78 > text1,23,text3,56 > text1,23,text2,45 > > The following command was entered: > > datasubset<-data.frame(dataset[,3]=="text3") > > The result of > > datasubset > > is > > TRUE > TRUE > > The required result is > > text1,23,text3,78 > text1,23,text3,56 > > What is the correct command to use please? > > Thanks in advance. > > __ > 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] Selecting ranges of dates from a dataframe
Benjamin, A more elegant "R-style" solution would be to use one of R's "apply"/aggregation routines, of which there are many. For example, the "by" function can split a data.frame by some factor/categorical variable(s), and then apply a function to each "slice". The result can then be pieced back together. See below for an example in which this factor is simply a parallel vector of pure dates: # extract pure date component of time and date dates <- format(serv$datum, "%Y-%m-%d") # write auxilliary function to aggregate a "slice" of the data.frame # x will be a "slice" of data from a single day aggregateDf <- function(x) { # return a one-row data.frame data.frame(datum = format(x$datum[1], "%Y-%m-%d"), write = sum(x$write), read = sum(x$read) ) } # now process each "slice" of the serv data.frame using "by" splitVals <- by(serv, dates, aggregateDf ) # bind back into a single data.frame values <- do.call(rbind, splitVals) The difference in execution speed is pretty negligible on my machine, so it's a more concise solution but I don't know if it is much faster. HTH, Francisco On Thu, Mar 10, 2011 at 1:23 PM, Benjamin Stier < benjamin.st...@ub.uni-tuebingen.de> wrote: > Hello list! > > I have a data.frame which looks like this: > > serv > datum op.read op.write read write > 1 2011-01-29 10:00:00 00 0 0 > 2 2011-01-29 10:00:01 00 0 0 > 3 2011-01-29 10:00:02 00 0 0 > 4 2011-01-29 10:00:03 04 0 647168 > 5 2011-01-29 10:00:04 00 0 0 > 6 2011-01-29 10:00:05 0 14 0 1960837 > 7 2011-01-29 10:00:06 00 0 0 > ... > 115 2011-01-30 10:00:54 00 0 0 > 116 2011-01-30 10:00:55 00 0 0 > 117 2011-01-30 10:00:56 00 0 0 > 118 2011-01-30 10:00:57 540 29184 0 > 119 2011-01-30 10:00:58 2040 122880 0 > 120 2011-01-30 10:00:59 00 0 0 > ... > > I want to compare read/write from each day. I already have a solution, but > it > is pretty slow. > > # read the data > serv <- read.delim("cut.inp") > > # Reformat the dates from the file > serv$datum <- strptime(serv$datum, "%Y-%m-%d %H:%M:%S") > > # select all single days > dates.serv <- unique(strptime(serv$datum, format="%Y-%m-%d")) > > # create a data.frame > values <- data.frame(row.names=1, datum=numeric(0), write=numeric(0), > read=numeric(0)) > for(i in as.character(dates.serv)) { ># build up a values for a day-range >searchstart <- as.POSIXlt(paste(i, "00:00:00", sep=" ")) >searchend <- as.POSIXlt(paste(i, "23:59:59", sep=" ")) ># select all values from a specific day >day <- serv[(serv$datum >= searchstart & serv$datum <= searchend),] >write <- as.numeric(sum(as.numeric(day$write))) >read <- as.numeric(sum(as.numeric(day$read))) ># add to the data.frame >values <- rbind(values, data.frame(datum=i, write=write, read=read)) > } > > This is my first try using R for statistics so I'm sure this isn't the best > solution. > The for-loop does it's job, but as I said is really slow. My data is for 21 > days and 1 line per second. > Is there a better way to select the date-ranges instead of a for-loop? The > line where I select all values for "day" seems to be the heaviest. Any > idea? > > Kind regards, > > Benjamin > > PS: I attached some sample data, in case you want to try for yourself. > > __ > 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] Associating the day of week to a daily xts object
Victor, The "weekdays" function will return the days of the week (as a character vector of names) that a given vector of dates (Date or POSIXct) fall on. These can then be converted into numbers using a look-up table/vector. See below for an example using the sample_matrix data included with the xts package. ## require(xts) data(sample_matrix) dates <- as.POSIXct(rownames(sample_matrix), format = "%Y-%m-%d") dayLookup <- 1:7 names(dayLookup) <- c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun") datesDays <- dayLookup[weekdays(dates, abbreviate = TRUE)] print(datesDays) ## >From here, you can just add the "datesDays" vector as an additional column to the original data, e.g. xts(cbind(sample_matrix, dow =datesDays), order.by = dates). HTH, Francisco On Mon, Mar 7, 2011 at 9:05 AM, Victor wrote: > I have the following xts objetct "temp" > > > str(temp) > An xts object from 2010-12-26 to 2011-03-05 containing: > Data: num [1:70, 1] 2.95 0.852 -0.139 1.347 2.485 ... > - attr(*, "dimnames")=List of 2 > ..$ : NULL > ..$ : chr "t_n" > Indexed by objects of class: [POSIXct,POSIXt] TZ: GMT > xts Attributes: > NULL > > > > temp > t_n > 2010-12-26 2.950 > 2010-12-27 0.8520833 > 2010-12-28 -0.1390625 > ... > > I would like to associate another column with the day of week in the form > of 1=Mon, 2=Tue, ..., 7=Sun > in order to obtain: > > >newtemp > > t_n dow > 2010-12-26 2.9507 > 2010-12-27 0.85208331 > 2010-12-28 -0.13906252 > .. > > How could make this in the shortest (and elegant?) way? > > Ciao from Rome > Vittorio > > __ > 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] Writing Rd files
Hi Nipesh, Try "\cr". See section 2.2 of the "Writing R Extensions" document for more information. HTH, - Francisco On Sun, Mar 6, 2011 at 2:03 PM, Nipesh Bajaj wrote: > Hi all, I have created a package and now into writing it's help files. > However I am having problem on, how to put a 'new line' in any > statement of the help file? For example please consider following: > > \title{ > This is a new function and this function will calculate the mean. > } > > However I want to write it is this way: > > \title{ > This is a new function & and, > this function will calculate the mean. > } > > I have tried using: "\n" or "\\" however could not achieve what I > want. Can somebody please help me on how to tell that, break here and > go to the next line? > > Thanks, > > __ > 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] Generic mixup?
Hi Nick, There is a difference between the newer S4 generic functions/methods and S3 generic functions. See for example: methods("print") showMethods("show") ?Methods HTH, -Francisco On Fri, Mar 4, 2011 at 10:14 AM, Nick Sabbe wrote: > Hello list. > > > > This is from an R session (admittedly, I'm still using R 2.11.1): > > > print > > function (x, ...) > > UseMethod("print") > > > > > showMethods("print") > > > > Function "print": > > > > > > Don't the two results contradict each other? Or do I have a terrible > misunderstanding of what comprises a generic function? > > > > Thx, > > > > Nick Sabbe > > -- > > ping: nick.sa...@ugent.be > > link: <http://biomath.ugent.be/> http://biomath.ugent.be > > wink: A1.056, Coupure Links 653, 9000 Gent > > ring: 09/264.59.36 > > > > -- Do Not Disapprove > > > > >[[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. > [[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] Species accumulation curves with differential sampling effort
Hello! I'm a PhD student working with coral reef fish diversity in Mexico. I want to do species accumulation curves but I have differential sampling effort for each "sample". Do you know or have developed an R script that consider the differential effort of each sample? PRIMER and other programs use each pack of species as a replica without the possibility of telling the program that the first group of species (sample) was collected during a 30 minutes dive and the second group (sample) was taking after 60 minutes dive and so on. I would really appreciate any kind of help! Best regards, Vanessa Francisco <">< [[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] Efficient nested loops
Dear R community, I am working with huge arrays, so I spend a lot of time computing. This is my code: for (x in 1:dim(variable)[1]){ for (y in 1:dim(variable)[2]){ for (z in 1:dim(variable)[3]){ result <- max(variable[x,y,z,]) } } } Is there a more efficient procedure to do this task? Thanks in advance! [[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 to use GMerrorsar () spdep
Hi Folks, I am trying to work with spatial logit model and I stopped my research in a Paper of Pinkse that make estimation of spatial models by GMM. I search in R communite and I find that GMerrorsar function in spdep package is used to it, but I have tryied so many times to make my model by this function and i did not figure it out. Some body knows how to work with this function? If so, can send me an example? I will be very glad for help. Thanks in advance, Francisco Gildemir Ferreira da Silva __ 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] How can I get the principal components after a varimax rotation using Varimax function?
Dear R users, I have some question about Varimax rotation of the loadings obtained from a PCA. Imagine X is a field where the rows are the observations and the columns are the variables. I obtain the loadings (L), the principal components (PC) and the percentage of the variance for each mode (PVar) from PCA. I would like to do a varimax rotation in order to find an easiest physical interpretation of the results. A way to do this task is to use the varimax function. I apply this function to the loadings as: L_varimax = L %*% varimax(L)$rotmat Now I would like to calculate the rotated principal components. It is possible to calculate it as pcs_varimax = scale(X) %*% L_varimax? If I do it, the standard desviation of the rotated principal components are not equal to 1. How can I solve it? How can I get the percentage of the variance for each rotated principal component? Thanks in advance [[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] Spatial Logit Model
Hello everybody, I am trying to do a spatial logit model with the spatially autoregressive error structure - SAE. Right now, I did an implementation and i would like to check the results. Somebody knows with there is a package in R that estimates parameters to spatial logit model with SAE structure? or Somebody else had made some implementations of this model and can help me? Thanks a lot in advances, Francisco Gildemir Ferreira da Silva __ 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] Levelplot
This would do levelplot(m, xlab="", ylab="") F zac...@lmb.uni-muenchen.de wrote: Dear mailing list, I am trying to find out, how do a levelplot without labels on the x- and y-axis. The labels=FALSE does not work...can anyone help me? m <- matrix(1:25, ncol=5) levelplot(m, labels=F) Regards, Benedikt __ 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] Spam of an Adjacency Matrix
Hi folks, I have a matrix of 3 columns and 17 lines that represents a graph or a adjacency matrix. I have also a vector of 30 elements with some of the nodes of the graph repeated. seems like: 1. matrix that represents a graph: 1 2 1 1 3 1 1 4 1 2 1 1 2 4 1 3 1 1 4 1 1 4 2 1 2. vector of nodes repeated: 1 1 1 1 1 2 2 3 3 1 1 4 4 4 4 Well, I have tried to do a code to take the repeated nodes and make a spam of the matrix above in a adjacency matrix representation and I get not sucess. The code seems like: mat_spam<-function(mat_rev,vet_adj){ i<-1 while (i<6) { j<-1 while (j<17) { m<-1 while (m<6) { if (mat_rev[i+1,m]==mat_rev[i, m+1]){ mat_rev[i+1,m+1]=0} else{if (mat_rev[i,m+1]==vet_adj[j,1] & mat_rev[i,m+2]==vet_adj[j,2]){ mat_rev[i+1,m+1]=1} else{mat_rev[i+1,m+1]=0}} m<-m+1} j<-j+1} i<-i+1} return(mat_rev)} Anyone has a sugestion how can I do this? The principal question is how to spam an adjacency matrix to a biggest one by repeating nodes that i am representing by a vector. Thanks in advances Francisco Gildemir Ferreira da Silva [[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] PDF text to work with maps.
Hi all folks, I would like to know if somebody has a PDF text with the first steps to use maps in R, like: insert maps, open maps, create adjacency matrix, make the moran index, etc. If somebody can send me a material of it I will be very glad and thankful Thanks a lot, Gildemir Silba __ 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] R-help Digest, Vol 86, Issue 26
Hi all folks, I would like to know if somebody has a PDF text with the first steps to use maps in R, like: insert maps, open maps, create adjacency matrix, make the moran index, etc. If somebody can send me a material of it I will be very glad and thankful Thanks a lot, Gildemir Silba __ 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] plotting a function
Take a look at the examples in ?curve. As suggested by Erik, in the future please read the posting guide so you can get a more accurate response. Regards, Francisco Dr. Francisco J. Zagmutt Vose Consulting 1643 Spruce St., Boulder Boulder, CO, 80302 USA www.voseconsulting.com Jin wrote: Hello Dear, I am trying to plot a function to see a minimum point (actually, using "optim"). For example, 1. y=f(x) 2. x has a range 3. plot(x,y) to see a point x minimizing y I tried "plot(x,y)", but it made an error. Am I doing a right way? Thanks, Jin __ 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] Problems if optimization
What's up fellows... I am a begginer in R and i am trying to find the parameters of one likelihood function, but when i otimize it, always appers a error or advertisement and the solve does not occur. The problem seems like that: "lMix<-function(pars,y){ beta1<-pars[1] beta2<-pars[2] beta3<-pars[3] beta4<-pars[4] beta5<-pars[5] alfa1<-pars[6] Fsp<-log(1/(1+exp(beta1*y[,10]+beta2*y[,3]+beta3*y[,3]+beta4*y[,5]+beta5*y[,6]+alfa1*y[,11]))) Frp<-log(1/(1+exp(beta1*y[,10]+beta2*y[,3]+beta3*y[,3]+beta4*y[,5]+beta5*y[,6]))) logl<- sum((y[,15]*Fsp)+(y[,19]*Frp)) return(-logl) } optim(c(1,1,1,1,1,1), llMix, y=Mix, method="CG") Erro em optim(c(1, 1, 1, 1, 1, 1), llMix, y = Mix, method = "CG") : Função não pode ser calculada nos parâmetros iniciais optimize(c(1,1,1,1,1,1),llMix,y=Mix, method="CG") Erro em min(interval) : invalid 'type' (closure) of argument optim(llMix,pars=c(0,0,0,0,1,1),y=Mix, method="CG") Erro em optim(llMix, pars = c(0, 0, 0, 0, 1, 1), y = Mix, method = "CG") : cannot coerce type 'closure' to vector of type 'double' " What is it happen? Any one can fix it for me or give me suggestions? I have also tried the maxLik package and it does not make the first interaction on the maximization process. Thanks a lot, Gildemir [[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] Error analysis for circular data
Dear R users, I would like to know if it is possible to calculate the Mean Error (ME), the Root Mean-squared error (RMSE) and absolute error (MAE) for two temporal series of directional data. Where Can I get documentation about it? Thanks in advance [[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] Doubt about cluster analysis
Dear R community, I'm a beginner with Cluster Analysis. I would like to know if there is a criterion to select the best set of clusters to do this analysis. Thanks in advance, [[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] MLE in R
Hello Liang, Besides looking at ?optim, you may also want to look into this nice working example www.mayin.org/ajayshah/KB/R/documents/mle/mle.html Regards, Francisco Francisco J. Zagmutt Vose Consulting 1643 Spruce St., Boulder Boulder, CO, 80302 USA www.voseconsulting.com Liang Wang wrote: Hi, dear R users I am a newbie in R and I need to use the method of meximum likelihood to fit a Weibull distribution to my survival data. I use "optim" as follows: optim(c(1.14,0.25),weibull.like,mydata=mydata,method="L-BFGS-B",hessian = TRUE) My question is: how do I setup the constraints that the two parametrs of Weibull to be pisotive? Many thanks! Any comments are greatly appreciated! Steven [[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: barchart() {lattice}
Hello Xin, Take a look at the examples under ?print.trellis Using your original example, you could use: require(lattice) p1=barchart(yield ~ variety | site, data = barley, groups = year, layout = c(1,6), ylab = "Barley Yield (bushels/acre)", scales = list(x = list(abbreviate = TRUE, minlength = 5))) p2=barchart(yield ~ variety | site, data = barley, groups = year, layout = c(1,6), stack = TRUE, ylab = "", scales = list(x = list(rot = 45))) #I removed the legend and ylab to make it look a bit better print(p1, split=c(1,1,2,1), more=TRUE) print(p2, split=c(2,1,2,1)) I hope this helps, Francisco Francisco J. Zagmutt Vose Consulting 1643 Spruce St., Boulder Boulder, CO, 80302 USA www.voseconsulting.com Xin Ge wrote: Hi All, I'm trying par(mfrow(c(1,2))) with barchart(), but its not working. Can I display two or more barcharts on a same page using some other function? I'm using following code --- where barchart() part is taken from help manual. library(lattice) par(mfrow=c(1,2)) barchart(yield ~ variety | site, data = barley, groups = year, layout = c(1,6), ylab = "Barley Yield (bushels/acre)", scales = list(x = list(abbreviate = TRUE, minlength = 5))) barchart(yield ~ variety | site, data = barley, groups = year, layout = c(1,6), stack = TRUE, auto.key = list(points = FALSE, rectangles = TRUE, space = "right"), ylab = "Barley Yield (bushels/acre)", scales = list(x = list(rot = 45))) par(mfrow=c(1,1)) Thanks, Xin [[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] Off topic - Compendium of distributions
This is not an R related posting but I thought it would be interesting for readers of this list. Apologies for any cross-posting Dear all Our company Vose Software has just made a very comprehensive “Compendium of Distributions” available for free online at www.vosesoftware.com/content/ebook.pdf. The document explains the thinking behind and the uses of 76 distributions and gives plots of the distributions with different parameter values. It also has a section to give you a more intuitive understanding of formulas for things like density, moments, etc. and gives lists of possible candidate distributions for different types of problems like waiting time, stock price movements, expert estimates, etc. We hope you will find it useful! Regards, Francisco Francisco J. Zagmutt Senior Risk Analysis Consultant Vose Consulting 1643 Spruce St., Boulder Boulder, CO, 80302 USA francisco(at)voseconsulting(com) www.voseconsulting.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] Doubt about CCA and PCA
Dear R community, I'm working with PCA and CCA methods, and I have a theoretical question. Why is it necesary to have more temporal values than variables when the CCA O PCA are going to be used? Could you advise to me some any paper about it? Thanks in advance, [[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] LondonR meeting - 21st of July
Thank you to everyone for showing such an interest in the next LondonR meeting. Below is the agenda for the meeting: LondonR meeting - 21st July 2009 Time:4pm - 7pm Venue: The Wall 45 Old Broad Street London ECN 1HU Tel: 020 7588 4845 Agenda: * 4pm - Richard Pugh Introduction * 4.15pm - Matthew Dowle Data.table: Higher speed time series queries * 4.45pm - Francisco Gochez (To be confirmed) An overview of RMetrics 2009 * 5.15pm - Brandon Whitcher Developing S4 Classes for Medical Imaging Data: Initial Experience(s) * 5.45pm - Richard Weeks Thoughts from useR 2009 - Rennes * 6.15pm - General discussion * Networking and feedback The evening is a free event. Please e-mail Sarah Lewis at sle...@mango-solutions.com if you would like to attend. We hope the evening will be as successful as the last one. We look forward to meeting all of the attendees. Best Regards, Francisco Gochez mangosolutions data analysis that delivers T: +44 (0)1249 767700 F: +44 (0)1249 767707 Unit 2 Greenways Business Park Bellinger Close Chippenham Wilts SN15 1BN UK [[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] Second LondonR Group meeting
Hello, Mango is pleased to announce that the next LondonR Group meeting will be held on 21st July. At this meeting we plan to have a feedback session from people who attended useR! in Rennes. To help us plan this, if you are attending at Rennes next month and would be interested in giving some short feedback on an aspect of your choice please let us know. DateTuesday 21st July 2009 Time4pm to 7pm Venue The Wall 45 Old Broad St London EC2N 1HU Tel 020 7588 4845 To register for this event please send an email to market...@mango-solutions.com I hope you can attend and look forward to meeting you. Kind regards, Francisco mangosolutions S & R Consulting and Training T: +44 (0)1249 767700 F: +44 (0)1249 767707 M: +44 (0)7966 062462 Unit 2 Greenways Business Park Bellinger Close Chippenham SN15 1BN UK [[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] please, a help or information
Estimado Fernando, I am afraid that we will not be able to help you until you provide a reproducible example. See the information at the bottom of my email about the posting guide. Also, you may be interested in joining the Spanish R-help list here https://stat.ethz.ch/mailman/listinfo/r-help-es Regards, Francisco Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com __ R-help@r-project.org mailing list PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. bernardo lagos alvarez wrote: Hi users R list I have been running some codes for some simulations. Such intensive use runif, nlm and simpleError("...") functions no more. But still remain the error that has been greatly discussed in the list of users of R, the problem !crashes RGui.¡ are sorry for the inconvenience." here are the details of the crash: AppName: rgui.exe AppVer: 2.81.47281.0 ModName: r.dll ModVer: 2.81.47281.0 Offset: 000c1939 Anybody know a version of R, that is not making such errors? I am waiting for your assistance as soon as possible. Thanks you by your atention. Here datails of windows and an installed R: --- Nombre de host:TOSHIBA-USER Versi¢n del sistema operativo: 5.1.2600 Service Pack 3 Compilaci¢n 2600 Configuraci¢n del sistema operativo: Estaci¢n de trabajo independiente Tipo de compilaci¢n del sistema operativo: Multiprocessor Free Propiedad de: Bernardo Organizaci¢n registrada: Id. del producto: 76460-OEM-0011903-00111 Fecha de instalaci¢n original: 25/09/2006, 08:06:26 p.m. Tiempo de actividad del sistema: 0 d¡as, 0 horas, 58 minutos, 26 segundos Fabricante del sistema:TOSHIBA Modelo el sistema: TECRA A7 Tipo de sistema: X86-based PC Procesador(es):1 Procesadores instalados. [01]: x86 Family 6 Model 14 Stepping 8 GenuineIntel ~1829 Mhz Versi¢n del BIOS: TOSINV - 604 Directorio de sistema: C:\WINDOWS\system32 Dispositivo de inicio: \Device\HarddiskVolume1 Configuraci¢n regional del sistema:0c0a Idioma:040A Zona horaria: N/D Cantidad total de memoria f¡sica: 1,022 MB Memoria f¡sica disponible: 396 MB Memoria virtual: tama¤o m ximo:2,048 MB Memoria virtual: disponible: 2,006 MB - release of R _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 8.1 year 2008 month 12 day22 svn rev47281 language R version.string R version 2.8.1 (2008-12-22) -- Bernardo. University of Concepción. Chile. __ 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] Surface for R outside of R
Hi Koj, I just completed a Windows application using the batch approach and it works very well. In our case, we used VBA for Excel to call different batch files that execute R code, but you can do the same from any platform. Here is a simple step-by-step example on how to make the batch files work (in Windows): 1. Write your R code and save it the directory where you want your results. Here is a very simple example using a file called "normal.r". This code takes 100 random samples from a standard normal and writes the results to a file called "data.txt" write(rnorm(100),"data.txt") q(save="no", runLast = F) 2. Make sure that the console finds the path to R without changing your path environment variable. The easiest way to do this is to download the batchfiles that Gabor Grothendieck has kindly provided in this site http://code.google.com/p/batchfiles/ and put one or all of the files in the same directory where your batch and R files will reside. The only file I have needed so far in different WinXP and Vista machines was "Rcmd.bat" but you may need to use others. 3. Open a text editor and create a file with a .bat or .cmd extension. Here is an example of a file called "TestBatch.CMD" that runs the code in "normal.r" @ECHO OFF title Random number generation ECHO Taking 100 random samples from a N(0,1) Rcmd BATCH --slave normal.r Log.txt ECHO Analysis done. See the file data.txt for results PING 1.1.1.1 -n 1 -w 1000 >NUL @ECHO OFF The main line here is "Rcmd BATCH --slave normal.r Log.txt" which tells the console to run your R code, and to write outputs printed in the console (i.e. errors) to the Log.txt file. The PING argument adds a short lag before the console closes, so the user can see what is echoed in the console. The other lines should be pretty self-explanatory. 4. Execute the "TestBatch.CMD" file (by hand or via your GUI) and watch the results. You should now see two new files ("data.txt" and "Log.txt") in the same directory as the batch file. I hope this helps. Let me know if you have any other questions. Regards, Francisco Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com Hans-Peter Suter wrote: want an analysis. The best case were a surface of e.g. 8 Buttons, each click leads to start a specific R file. My outputs are JPEG or CSV, so I don`t need the output inside of R. Could anyone can give me some recommendations, what could be a solution (e. g. Java)? Is such a solution possible? What about batch scripts which would be called by your GUI? (maybe you could even skip the gui...) (see ?Rscript (unix alike) or Windows FAQ, 2.12) -- Regards, Hans-Peter __ 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] Course announcement: R for Financial Data Analysis in New York (July 16-18th)
Dear userRs, Mango Solutions are pleased to announce that we will deliver a 3-day introductory R course focused on financial analysis in New York on the 16-18th of July. The course topics are as follows: * Introduction * The R Environment * Data Objects * Functions * Important R Functions * Traditional Graphics * Basic R Statistical and Mathematical Functions used in the Finance Industry * Time Series Analysis, Manipulation and Visualization * Performance and Risk Measurement As with all our courses, attendees are provided with comprehensive training manuals complete with detailed examples and laminated tip sheets for future reference. The cost of this course is $1800 for commercial attendees and $1000 for academic attendees. Should you wish to book a place on this course or have any questions please contact us at train...@mango-solutions.com. Further information is available on our web-site at http://www.mango-solutions.com/services/rtraining/r_finance.html Kind regards, Francisco mango solutions S & R Consulting and Training +44 (0)1249 767 700 __ 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] memory.limit
Hello Stephen, I can confirm that I get the same behavior in my Windows machine. Here is a summary: > memory.limit() [1] 2046 > memory.limit(2092) Error in trunc(.Internal(memory.size(size))) : Non-numeric argument to mathematical function > memory.limit() [1] 2092 As you described, the function reports an error but it indeed modified the memory allocation limit. This must be somehow related to the modification to memory.limit() described in the release notes for Windows R version 2.9.0 (http://cran.r-project.org/bin/windows/base/CHANGES.R-2.9.0) > version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 9.0 year 2009 month 04 day17 svn rev48333 language R version.string R version 2.9.0 (2009-04-17) All: is this a bug or are we missing something? Regards, Francisco __ Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA franci...@voseconsulting.com www.voseconsulting.com Derek Stephen Elmerick wrote: I ran the memory limit function in R 2.9.0 and received the 'error' below. The memory appears to update correctly, so there's probably no implication beyond cosmetic; however, thought I would make sure since the function as written did not generate the same error in my 2.8.0 version of R. Thanks memory.limit(size=4095) Error in trunc(.Internal(memory.size(size))) : Non-numeric argument to mathematical function memory.limit() [1] 4095 [[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@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] split a character variable into several character variable by a character
Hello Mao, If the popcode variable has a fixed number of characters (i.e each entry has 9 characters), you can use a simple call to substr: dat<-read.table("clipboard", header=T)#Read from your email varleft<-substr(dat$popcode,0,6) varright<-substr(dat$popcode,8,9) datnew<-data.frame(dat,varleft,varright) > datnew popcode codetot p3need varleft varright 1 BCPy01-01 BCPy01-01-1 100. BCPy01 01 2 BCPy01-01 BCPy01-01-2 100. BCPy01 01 3 BCPy01-01 BCPy01-01-3 100. BCPy01 01 4 BCPy01-02 BCPy01-02-1 92.5926 BCPy01 02 5 BCPy01-02 BCPy01-02-1 100. BCPy01 02 6 BCPy01-02 BCPy01-02-2 92.5926 BCPy01 02 7 BCPy01-02 BCPy01-02-2 100. BCPy01 02 8 BCPy01-02 BCPy01-02-3 92.5926 BCPy01 02 9 BCPy01-02 BCPy01-02-3 100. BCPy01 02 10 BCPy01-03 BCPy01-03-1 100. BCPy01 03 You can use a similar construction for codetot. I hope this helps, Francisco Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA franci...@voseconsulting.com www.voseconsulting.com Mao Jianfeng wrote: Dear, R-lister, I have a dataframe like the followed. And, I want to split a character variable ("popcode", or "codetot") into several new variables. For example, split "BCPy01-01" (in popcode) into "BCPy01" and "01". I need to know how to do that. I have tried strsplit() and substring() functions. But, I still can not perform the spliting. Any advice? Thanks in advance. df1: popcode codetot p3need BCPy01-01 BCPy01-01-1 100. BCPy01-01 BCPy01-01-2 100. BCPy01-01 BCPy01-01-3 100. BCPy01-02 BCPy01-02-1 92.5926 BCPy01-02 BCPy01-02-1 100. BCPy01-02 BCPy01-02-2 92.5926 BCPy01-02 BCPy01-02-2 100. BCPy01-02 BCPy01-02-3 92.5926 BCPy01-02 BCPy01-02-3 100. BCPy01-03 BCPy01-03-1 100. Regards, Mao Jian-feng [[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@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] Output an RWeka model via sink
Thanks. Just figured out what happened: out of quota. -Original Message- From: Coen van Hasselt Sent: Saturday, March 28, 2009 9:05 PM To: Francisco Javier Perez Caballero Cc: r-h...@stat.math.ethz.ch Subject: Re: [R] Output an RWeka model via sink I tried your code and it seems to work fine; the file contains the expected output. I used R 2.8.1 on WinXP. On Sun, Mar 29, 2009 at 07:48, Francisco Javier Perez Caballero wrote: > When I sink the output of an RWeka model to a text file, the output file > appears empty: > > library(RWeka) > model = LogitBoost(Species~.,data=iris) > print(model) > sink("output.txt") > print(model) > #file output.txt is created, but it is blank > sink() > > > Am I doing anything wrong? > [[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@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] Output an RWeka model via sink
When I sink the output of an RWeka model to a text file, the output file appears empty: library(RWeka) model = LogitBoost(Species~.,data=iris) print(model) sink("output.txt") print(model) #file output.txt is created, but it is blank sink() Am I doing anything wrong? [[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] bigglm() results different from glm()
That's a much cleaner solution! It would be nice if biglm takes the defaults from options(digits), but off course we can also just use print() as you pointed out. Thanks again for your replies and for making this library available to the community. Francisco -- Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com Thomas Lumley wrote: There's a digits= argument to the print method. a <- bigglm(ff,data=trees, chunksize=10, sandwich=TRUE) print(summary(a),digits=5) Large data regression model: bigglm(ff, data = trees, chunksize = 10, sandwich = TRUE) Sample size = 31 Coef (95% CI) SE p (Intercept) -6.63162 -8.08729 -5.17595 0.72783 0 log(Girth) 1.98265 1.87132 2.09398 0.05567 0 log(Height) 1.11712 0.73339 1.50085 0.19186 0 Sandwich (model-robust) standard errors I suppose I should make it take its default from options(digits)-3 or something. -thomas On Tue, 17 Mar 2009, Francisco J. Zagmutt wrote: Dear Thomas and John, Thanks for your prompt reply and for looking at your code to explain these differences. I see you replied very late at night, so I am sorry if my little problem kept you awake! As you pointed out, the model indeed converges (as reported in model$converged) once I specify a larger number of iterations. A very minor comment: it seems that the reporting of the estimates in summary.biglm() is not taking the parameters from options(digits). For example, using the same data and models as before: require(biglm) options(digits=8) dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log"), maxit=20) summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.3019509 0.0014147 1627.21 < 2.2e-16 *** ttment2 0.4052002 0.0018264 221.86 < 2.2e-16 *** summary(m1big) Coef (95% CI)SE p (Intercept) 2.302 2.299 2.305 0.001 0 ttment2 0.405 0.402 0.409 0.002 0 To get more digits I can extract the point estimates using coef(m1big), but after looking at str(m1big) the only way I could figure to extract the p-values was: summary(m1big)$mat[,"p"] (Intercept) ttment2 0 0 Is there a way I can get summary.biglm() to report more digits directly? Thanks again, Francisco -- Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com Thomas Lumley wrote: Yes, the slow convergence is easier to get with the log link. Overshooting the correct coefficient vector has more dramatic effects on the fitted values and weights, and in this example the starting value of (0,0) is a long way from the truth. The same sort of thing happens in the Cox model, where there are real data sets that will cause numeric overflow in a careless implementation. It might be worth trying to guess better starting values: saving an iteration or two would be useful with large data sets. -thomas On Tue, 17 Mar 2009, John Fox wrote: Dear Francisco, I was able to duplicate the problem that you reported, and in addition discovered that the problem seems to be peculiar to the poisson family. lm(y~ttment, data=dat) and biglm(y~ttment, data=dat) produce identical results, as do glm(y~ttment, data=dat) and bigglm(y~ttment, data=dat). Another example, with the binomial family: set.seed(12345) dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) dat$y0 <- ifelse(dat$y > 12, 1, 0) m1b <- glm(y0~ttment, data=dat, family=binomial) m1bigb <- bigglm(y0~ttment , data=dat, family=binomial()) coef(m1b) (Intercept) ttment2 -1.33508 2.34263 coef(m1bigb) (Intercept) ttment2 -1.33508 2.34263 deviance(m1b) [1] 109244 deviance(m1bigb) [1] 109244 I'm copying this message to Thomas Lumley, who's the author and maintainer of the biglm package. Regards, John -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Francisco J. Zagmutt Sent: March-16-09 10:26 PM To: r-h...@stat.math.ethz.ch Subject: [R] bigglm() results different from glm() Dear all, I am using the bigglm package to fit a few GLM's to a large dataset (3 million rows, 6 columns). While trying to fit a Poisson GLM I noticed that the coefficient estimates were very different from what I obtained when estimating the model on a smaller dataset using glm(), I wrote a very basic toy example to compare the results of bigglm() against a glm() call. Consider the following code: > require(biglm) > options(digits=6, scipen=3, contrasts = c("contr.treatment", "contr.poly")) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) >
Re: [R] bigglm() results different from glm()
Dear Thomas and John, Thanks for your prompt reply and for looking at your code to explain these differences. I see you replied very late at night, so I am sorry if my little problem kept you awake! As you pointed out, the model indeed converges (as reported in model$converged) once I specify a larger number of iterations. A very minor comment: it seems that the reporting of the estimates in summary.biglm() is not taking the parameters from options(digits). For example, using the same data and models as before: > require(biglm) > options(digits=8) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) > m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log"), maxit=20) > summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.3019509 0.0014147 1627.21 < 2.2e-16 *** ttment2 0.4052002 0.0018264 221.86 < 2.2e-16 *** > summary(m1big) Coef (95% CI)SE p (Intercept) 2.302 2.299 2.305 0.001 0 ttment2 0.405 0.402 0.409 0.002 0 To get more digits I can extract the point estimates using coef(m1big), but after looking at str(m1big) the only way I could figure to extract the p-values was: > summary(m1big)$mat[,"p"] (Intercept) ttment2 0 0 Is there a way I can get summary.biglm() to report more digits directly? Thanks again, Francisco -- Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com Thomas Lumley wrote: Yes, the slow convergence is easier to get with the log link. Overshooting the correct coefficient vector has more dramatic effects on the fitted values and weights, and in this example the starting value of (0,0) is a long way from the truth. The same sort of thing happens in the Cox model, where there are real data sets that will cause numeric overflow in a careless implementation. It might be worth trying to guess better starting values: saving an iteration or two would be useful with large data sets. -thomas On Tue, 17 Mar 2009, John Fox wrote: Dear Francisco, I was able to duplicate the problem that you reported, and in addition discovered that the problem seems to be peculiar to the poisson family. lm(y~ttment, data=dat) and biglm(y~ttment, data=dat) produce identical results, as do glm(y~ttment, data=dat) and bigglm(y~ttment, data=dat). Another example, with the binomial family: set.seed(12345) dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) dat$y0 <- ifelse(dat$y > 12, 1, 0) m1b <- glm(y0~ttment, data=dat, family=binomial) m1bigb <- bigglm(y0~ttment , data=dat, family=binomial()) coef(m1b) (Intercept) ttment2 -1.33508 2.34263 coef(m1bigb) (Intercept) ttment2 -1.33508 2.34263 deviance(m1b) [1] 109244 deviance(m1bigb) [1] 109244 I'm copying this message to Thomas Lumley, who's the author and maintainer of the biglm package. Regards, John -Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Francisco J. Zagmutt Sent: March-16-09 10:26 PM To: r-h...@stat.math.ethz.ch Subject: [R] bigglm() results different from glm() Dear all, I am using the bigglm package to fit a few GLM's to a large dataset (3 million rows, 6 columns). While trying to fit a Poisson GLM I noticed that the coefficient estimates were very different from what I obtained when estimating the model on a smaller dataset using glm(), I wrote a very basic toy example to compare the results of bigglm() against a glm() call. Consider the following code: > require(biglm) > options(digits=6, scipen=3, contrasts = c("contr.treatment", "contr.poly")) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) > m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log")) > summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.303050.001411629 <2e-16 *** ttment2 0.404290.00183 221 <2e-16 *** Null deviance: 151889 on 9 degrees of freedom Residual deviance: 101848 on 8 degrees of freedom AIC: 533152 > summary(m1big) Large data regression model: bigglm(y ~ ttment, data = dat, family = poisson(link = "log")) Sample size = 10 Coef (95% CI)SE p (Intercept) 2.651 2.650 2.653 0.001 0 ttment2 4.346 4.344 4.348 0.001 0 > m1big$deviance [1] 287158986 Notice that the coefficients and deviance are quite different in the model estimated using bigglm(). If I change the chunk to seq(1000,1,1000) the estimates remain the same. Can someone help me understand what is causing thes
Re: [R] bigglm() results different from glm()
Dear Thomas and John, Thanks for your prompt reply and for looking at your code to explain these differences. I see you replied very late at night, so I am sorry if my little problem kept you awake! As you pointed out, the model indeed converges (as reported in model$converged) once I specify a larger number of iterations. A very minor comment: it seems that the reporting of the estimates in summary.biglm() is not taking the parameters from options(digits). For example, using the same data and models as before: > require(biglm) > options(digits=8) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) > m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log"), maxit=20) > summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.3019509 0.0014147 1627.21 < 2.2e-16 *** ttment2 0.4052002 0.0018264 221.86 < 2.2e-16 *** > summary(m1big) Coef (95% CI)SE p (Intercept) 2.302 2.299 2.305 0.001 0 ttment2 0.405 0.402 0.409 0.002 0 To get more digits I can extract the point estimates using coef(m1big), but after looking at str(m1big) the only way I could figure to extract the p-values was: > summary(m1big)$mat[,"p"] (Intercept) ttment2 0 0 Is there a way I can get summary.biglm() to report more digits directly? Thanks again, Francisco -- Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.com Thomas Lumley wrote: Yes, the slow convergence is easier to get with the log link. Overshooting the correct coefficient vector has more dramatic effects on the fitted values and weights, and in this example the starting value of (0,0) is a long way from the truth. The same sort of thing happens in the Cox model, where there are real data sets that will cause numeric overflow in a careless implementation. It might be worth trying to guess better starting values: saving an iteration or two would be useful with large data sets. -thomas On Tue, 17 Mar 2009, John Fox wrote: Dear Francisco, I was able to duplicate the problem that you reported, and in addition discovered that the problem seems to be peculiar to the poisson family. lm(y~ttment, data=dat) and biglm(y~ttment, data=dat) produce identical results, as do glm(y~ttment, data=dat) and bigglm(y~ttment, data=dat). Another example, with the binomial family: set.seed(12345) dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) dat$y0 <- ifelse(dat$y > 12, 1, 0) m1b <- glm(y0~ttment, data=dat, family=binomial) m1bigb <- bigglm(y0~ttment , data=dat, family=binomial()) coef(m1b) (Intercept) ttment2 -1.33508 2.34263 coef(m1bigb) (Intercept) ttment2 -1.33508 2.34263 deviance(m1b) [1] 109244 deviance(m1bigb) [1] 109244 I'm copying this message to Thomas Lumley, who's the author and maintainer of the biglm package. Regards, John -Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Francisco J. Zagmutt Sent: March-16-09 10:26 PM To: r-h...@stat.math.ethz.ch Subject: [R] bigglm() results different from glm() Dear all, I am using the bigglm package to fit a few GLM's to a large dataset (3 million rows, 6 columns). While trying to fit a Poisson GLM I noticed that the coefficient estimates were very different from what I obtained when estimating the model on a smaller dataset using glm(), I wrote a very basic toy example to compare the results of bigglm() against a glm() call. Consider the following code: > require(biglm) > options(digits=6, scipen=3, contrasts = c("contr.treatment", "contr.poly")) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) > m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log")) > summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.303050.001411629 <2e-16 *** ttment2 0.404290.00183 221 <2e-16 *** Null deviance: 151889 on 9 degrees of freedom Residual deviance: 101848 on 8 degrees of freedom AIC: 533152 > summary(m1big) Large data regression model: bigglm(y ~ ttment, data = dat, family = poisson(link = "log")) Sample size = 10 Coef (95% CI)SE p (Intercept) 2.651 2.650 2.653 0.001 0 ttment2 4.346 4.344 4.348 0.001 0 > m1big$deviance [1] 287158986 Notice that the coefficients and deviance are quite different in the model estimated using bigglm(). If I change the chunk to seq(1000,1,1000) the estimates remain the same. Can someone help me understand what is causing thes
[R] bigglm() results different from glm()
Dear all, I am using the bigglm package to fit a few GLM's to a large dataset (3 million rows, 6 columns). While trying to fit a Poisson GLM I noticed that the coefficient estimates were very different from what I obtained when estimating the model on a smaller dataset using glm(), I wrote a very basic toy example to compare the results of bigglm() against a glm() call. Consider the following code: > require(biglm) > options(digits=6, scipen=3, contrasts = c("contr.treatment", "contr.poly")) > dat=data.frame(y =c(rpois(5, 10),rpois(5, 15)), ttment=gl(2,5)) > m1 <- glm(y~ttment, data=dat, family=poisson(link="log")) > m1big <- bigglm(y~ttment , data=dat, family=poisson(link="log")) > summary(m1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.303050.001411629 <2e-16 *** ttment2 0.404290.00183 221 <2e-16 *** Null deviance: 151889 on 9 degrees of freedom Residual deviance: 101848 on 8 degrees of freedom AIC: 533152 > summary(m1big) Large data regression model: bigglm(y ~ ttment, data = dat, family = poisson(link = "log")) Sample size = 10 Coef (95% CI)SE p (Intercept) 2.651 2.650 2.653 0.001 0 ttment2 4.346 4.344 4.348 0.001 0 > m1big$deviance [1] 287158986 Notice that the coefficients and deviance are quite different in the model estimated using bigglm(). If I change the chunk to seq(1000,1,1000) the estimates remain the same. Can someone help me understand what is causing these differences? Here is my version info: > version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 8.1 year 2008 month 12 day22 svn rev47281 language R version.string R version 2.8.1 (2008-12-22) Many thanks in advance for your help, Francisco -- Francisco J. Zagmutt Vose Consulting 2891 20th Street Boulder, CO, 80304 USA www.voseconsulting.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] RODBC crashes connecting to Teradata
Hi, I'm trying to connect to a Teradata database via RODBC on a Linux 64 machine (Red Hat Enterprise Linux 5). The ODBC driver is properly configured and queries sent via unixODBC's isql tool work properly. However, this is what happens when I try to connect via RODBC: > library(RODBC) > conn = odbcConnect("thedsn", uid="theuid", pwd="thepass") *** caught segfault *** address (nil), cause 'memory not mapped' Traceback: 1: .Call(C_RODBCGetInfo, attr(stat, "handle_ptr")) 2: odbcDriverConnect(st, ...) 3: odbcConnect("thedsn", uid = "theuid", pwd = "thepass") Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Any idea of what is going on here? [[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] Problem with dir.create in windows servers
Hello. I've being having problems to create directories in a windows server environment . It seems that the recursive argument is not working properly on the intranet, as it does in a local path. For example: > dir.create("server/directory1/directory2") , this works fine, and creates the directory2, but If we want to create the directory3 and another directory called directory4 inside of it. >dir.create("server/directory1/directory3/directory4", recursive=T) fails, and don't create any directory. The ShowWarning doesn't show any message. Any help please. Fran. __ 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] Problem with dir.create in windows servers
Hello. I've being having problems to create directories in a windows server environment . It seems that the recursive argument is not working properly on the intranet, as it does in a local path. For example: > dir.create("server/directory1/directory2") , this works fine, and creates the directory2, but If we want to create the directory3 and another directory called directory4 inside of it. >dir.create("server/directory1/directory3/directory4", recursive=T) fails, and don't create any directory. The ShowWarning doesn't show any message. Any help please. Fran. __ 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 eigenvectors
Hi everybody, I have some problems with the function eigen. I have a square matrix and I want to calculate the eigenvalues and eigenvectors. I apply the function eigen and I get it, however when I solve the same problem in Statistica software, I realise that some eigenvectors are the opposite. How can I get the same values? Thanks in advance [[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] reshape matrices
Hello everyone, I need reshape an array. For example, if we have next array: > a <- c(1,2,3,4,5,6,7,8,9,10,11,12) > dim(a) <- c(2,2,3) > a , , 1 [,1] [,2] [1,]13 [2,]24 , , 2 [,1] [,2] [1,]57 [2,]68 , , 3 [,1] [,2] [1,]9 11 [2,] 10 12 I need to get next matrices: 1 2 3 4 5 6 7 8 9 10 11 12 1 3 2 4 5 7 6 8 9 11 10 12 It exist any function that can be able to do it? Thanks and sorry for my english. [[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] reshape matrices
Hello everyone, I need reshape an array. For example, if we have next array: > a <- c(1,2,3,4,5,6,7,8,9,10,11,12) > dim(a) <- c(2,2,3) > a , , 1 [,1] [,2] [1,]13 [2,]24 , , 2 [,1] [,2] [1,]57 [2,]68 , , 3 [,1] [,2] [1,]9 11 [2,] 10 12 I need to get next matrices: 1 2 3 4 5 6 7 8 9 10 11 12 1 3 2 4 5 7 6 8 9 11 10 12 Does any function exist that can be able to do it ? Thanks in advance and sorry for my english. [[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] Handle missing values
Hi Ted I think I understand what you say. Maybe I should reformat my data and then divide my single annual file in twelve monthly files. But maybe then I miss the possibility of making an anual analysis but retain more valid data, I have to think about it. Hope any expert in clustering can help Thanks Ted Francisco Pastor escribió: Hi everyone I am new to R and have a question about missing values. I am trying to do a cluster analysis of monthly temperatures and my data are 14 columns with spatial coordinates (lat,lon) and 12 monthly values: /lat - lon - temp1 - //temp2 - temp3 - - //temp12/ If I omit missing values (my missing values are 99.00) with /mydata <- na.omit(mydata)/ every row with a missing value (i.e. eleven good temperature values and one month missing) is deleted. I would like to retain all valid values for the k-means analysis but excluding. I've been trying and searching about na.omit, na.action, na.exclude but can't find the right point. Any help would be appreciated. -- ------- Francisco Pastor Meteorology department Fundación CEAM [EMAIL PROTECTED] http://www.ceam.es/ceamet - http://www.ceam.es Parque Tecnologico, C/ Charles R. Darwin, 14 46980 PATERNA (Valencia), Spain Tlf. 96 131 82 27 - Fax. 96 131 81 90 --- Usuario Linux registrado: 363952 __ 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] Handle missing values
Hi everyone I am new to R and have a question about missing values. I am trying to do a cluster analysis of monthly temperatures and my data are 14 columns with spatial coordinates (lat,lon) and 12 monthly values: /lat - lon - temp1 - //temp2 - temp3 - - //temp12/ If I omit missing values (my missing values are 99.00) with /mydata <- na.omit(mydata)/ every row with a missing value (i.e. eleven good temperature values and one month missing) is deleted. I would like to retain all valid values for the k-means analysis but excluding. I've been trying and searching about na.omit, na.action, na.exclude but can't find the right point. Any help would be appreciated. -- ------- Francisco Pastor Meteorology department Fundación CEAM [EMAIL PROTECTED] http://www.ceam.es/ceamet - http://www.ceam.es Parque Tecnologico, C/ Charles R. Darwin, 14 46980 PATERNA (Valencia), Spain Tlf. 96 131 82 27 - Fax. 96 131 81 90 --- Usuario Linux registrado: 363952 __ 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] I need help with eofs
Hello, I'm a beginner in R. I'm learning to use this fantastic program, but I have some problems in how to use it. First of all, I have a txt file witch I am able to load to the program. I'm very interested in PCA, and I have a lot of packages, but I haven't got the results that I want. I would like to get the EOFS and to export it in a txt file. I would be pleased if you could help me. Thank you, and sorry for my english, Regards Francisco Javier [[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] I need help with eofs
Hello, I'm a beginner in R. I'm learning to use this fantastic program, but I have some problems in how to use it. First of all, I have a txt file witch I am able to load to the program. I'm very interested in PCA, and I have a lot of packages, but I haven't got the results that I want. I have got the PCs series, but I need obtain the spatial patterns (EOFS) and to export it in a txt file. I would be pleased if you could help me. Thank you, and sorry for my english, Regards Francisco Javier [[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] Plotting barplot and scatterplot on same device - x-axis problem
You could try using the arguments type, lwd and lend in the call to plot(). For example: plot(1:10, type="h", lend=2, lwd=10) points((1:10)*1.1, type="b") type="h" creates vertical lines under each observation, lwd controls the thickness of the lines, and lend=2 draws square ends on the lines. I hope this helps, Francisco Thomas Pedersen wrote: Hi R-users I'm a relative newbie and uses R mostly for graphical purpose. I have a layout problem when plotting a scatterplot and a barplot using par(new=TRUE). The baseline of the x-axis is not positioned equal for the two plotting functions (see picture) and I have been unable to find out how this is changed. http://www.nabble.com/file/p18025066/pic.jpeg I have added the script if this is of interrest: par(mar = c(7, 4, 4, 2) + 0.1) barplot(rbind(data$Asn,data$Glu,data$NH3), beside=T, axes=TRUE, xlim=c(0,57), ylim=c(0,10)) par(new=TRUE) plot(1:14,data$Acidification.time, axes=FALSE, type="b", xlab="", ylab="", xlim=c(0.3,14.7), ylim=c(6,8)) axis(1,pos=6, labels=FALSE, at=c(0.3,1:14,14.7)) text(1:14, par("usr")[3], srt = 90, adj = 1, labels = data$Month, xpd = TRUE) axis(4,pos=14.7) All help will be greatly appreciated __ 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] Helpfiles in HTML browser
Go to your Rprofile file (in the etc directory) and add the following line: options(htmlhelp=TRUE) I hope this helps Francisco Knut Krueger wrote: > I forgot how to switch between Windows helpfiles and Browser helpfiles. > f.e ?lm should open the browser. > Maybe anybody could give me a hint? > > Regards Knut > > __ > 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@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] Scientific Notation
I believe the argument to format is "scientific" i.e. axis(2, at=at, labels=format(at, scientific=FALSE)) Best regards, Francisco Duncan Murdoch wrote: > On 9/26/2007 11:24 AM, Jacques Wagnor wrote: >> Dear List: >> >> Below is how I specify an axis: >> >> axis(2, at=c(0.5, 0.0005)) >> >> R displays the numbers in scientific notation. What >> argument/parameter should I use to tell R to display the numbers as >> specified rather than in scientific notation? > > Something like > > axis(2, at=c(0.5, 0.0005), labels=c("0.5", "0.0005")) > > is a way to be 100% sure of what will be displayed, but in this > particular instance, > > at <- c(0.5, 0.0005) > axis(2, at=at, labels=format(at, sci=FALSE)) > > comes close, and there may be some other format spec that gets exactly > what you want. > > Duncan Murdoch > > __ > 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@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] Line Graph - Greater than 2 variables on plot
You can also use the facilities in the lattice package. Using Jim´s data names: require(lattice) xyplot(A+B+C~D, data=wag, type="l", auto.key=list(points = FALSE, lines = TRUE, space = "bottom"), ylab="value", main="Three variable plot") Regards, Francisco Jim Lemon wrote: > Wayne Aldo Gavioli wrote: >> Hello all, >> >> I was wondering if anyone knew how to construct a multiple line graph on R, >> where there are 2 (or more) sets of data points plotted against some x axis >> of >> data, and you can draw a line on the graph connecting each set of data >> points. >> >> For example: >> >> A B C D >> 0.65662.11851.23205 >> 0.647 2.08651.232510 >> 0.65322.10601.228715 >> 0.64872.12901.231320 >> 0.65942.12851.234125 >> 0.65772.10701.234330 >> 0.65792.13451.234035 >> 0.67342.17051.236240 >> 0.675 2.18451.237245 >> 0.65922.15501.234050 >> 0.66472.17101.230555 >> >> >> >> Would there be a way: >> a) To graph all the points of data in sets A, B and C as Y coordinates on one >> graph, using the points in set D as the X-axis/coordinates for all 3 sets >> (A, B >> and C)? >> b) To be able to draw 3 lines on the graph that connect each set of data (1 >> line >> connects all the A points, one line connects all the B points, one line >> connects >> all the C points) >> >> >> I couldn't find anything in the examples or the help section about multiple >> lines on the same graph, only one line. >> > Hi Wayne, > > Assume your data is in a data frame named "wag": > > plot(wag$D,wag$A,main="Three variable plot",xlab="D",ylab="Value", > ylim=range(wag[c("A","B","C")]),type="l",col=2) > lines(wag$D,wag$B,type="l",pch=2,col=3) > lines(wag$D,wag$C,type="l",pch=3,col=4) > legend(25,1.9,c("A","B","C"),lty=1,col=2:4) > > 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. > __ 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.