Re: [R] Runif Help: same variable, 3 different parameters
Cheers Dieter. You are a true bro. -- View this message in context: http://r.789695.n4.nabble.com/Runif-Help-same-variable-3-different-parameters-tp3073894p3074041.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Runif Help: same variable, 3 different parameters
pythonomics wrote: > > So I am working on an economic model and I need to change the parameters > of the runif statement as time goes on. > > X <-runif(1:50,0,5) > X <-runif(51:100,100,150) > X <-runif(100:T, 1,2) > > T=1000 c(runif(1:50,0,5), runif(51:100,100,150),runif(100:T, 1,2)) Dieter -- View this message in context: http://r.789695.n4.nabble.com/Runif-Help-same-variable-3-different-parameters-tp3073894p3074036.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Error in calcCurveGrob(x, x$debug) : End points must not be identical
Hi Bryan Hanson wrote: Hi All... I haven¹t found mention of this error anywhere. I'm trying to draw spline curves using grid graphics. Most of the time, I have no problems, but I have some data sets that give the error in the subject line. I'm not sure which end points are identical, but the end points passed to the function are definitely not identical. I get ... which(tst$x.st == tst$x.end & tst$y.st == tst$y.end) [1] 11 14 Also, the viewport you have set up will not show any of the curves because its "native" scale is 0 to 1. The following at least makes the curves visible ... grid.newpage() vp <- viewport(width = 0.9, height = 0.9, x = 0.5, y = 0.5, xscale=c(-10, 10), yscale=c(-15, 15)) pushViewport(vp) grid.rect(gp = gpar(lty = "dashed", col = "gray")) grid.points(0.5, 0.5, pch = 20, gp = gpar(cex = 0.5)) subset <- -c(11, 14) grid.curve(tst$x.st[subset], tst$y.st[subset], tst$x.end[subset], tst$y.end[subset], default.units = "native") Paul p.s. Intrigued to know what sort of image you are producing, if you are able to share. Any assistance appreciated! Bryan tst <- structure(list(x.st = c(-1, -2, -3, -1, -1.5, -3, -1.5, -1.5, -8, -1, -1.5, -1, -1.5, -2, -1.5, -2, -1, -1.5, -2), y.st = c(1.73205080756888, 3.46410161513776, 5.19615242270663, 1.73205080756888, 2.59807621135332, 5.19615242270663, 2.59807621135332, 2.59807621135332, 13.8564064605510, 1.73205080756888, 2.59807621135332, 1.73205080756888, 2.59807621135332, 3.46410161513776, 2.59807621135332, 3.46410161513776, 1.73205080756888, 2.59807621135332, 3.46410161513776), x.end = c(-6.5, -6.5, -6.5, -4, -4, -4, -1.5, -1, -1, -1.5, -1.5, -2, -2, -2, -3, -3, -8, -8, -8), y.end = c(-11.2583302491977, -11.2583302491977, -11.2583302491977, -6.92820323027551, -6.92820323027551, -6.92820323027551, -2.59807621135332, 1.73205080756888, 1.73205080756888, 2.59807621135332, 2.59807621135332, 3.46410161513776, 3.46410161513776, 3.46410161513776, 5.19615242270663, 5.19615242270663, 13.8564064605510, 13.8564064605510, 13.8564064605510 ), grp = c(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3), lwd = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 3, 1, 1, 1, 2, 5, 2)), .Names = c("x.st", "y.st", "x.end", "y.end", "grp", "lwd"), row.names = 34:52, class = "data.frame") grid.newpage() vp <- viewport(width = 0.9, height = 0.9, x = 0.5, y = 0.5) pushViewport(vp) grid.rect(gp = gpar(lty = "dashed", col = "gray")) grid.points(0.5, 0.5, pch = 20, gp = gpar(cex = 0.5)) grid.curve(tst$x.st, tst$y.st, tst$x.end, tst$y.end, default.units = "native") ### sessionInfo() R version 2.12.0 (2010-10-15) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] datasets tools grid grDevices graphics utils [7] stats methods base other attached packages: [1] gridExtra_0.7 GGally_0.2.2 xtable_1.5-6 [4] mvbutils_2.5.1 ggplot2_0.8.8 proto_0.3-8 [7] reshape_0.8.3 ChemoSpec_1.46 seriation_1.0-2 [10] colorspace_1.0-1 TSP_1.0-1 R.utils_1.5.3 [13] R.oo_1.7.4 R.methodsS3_1.2.1 rgl_0.92.794 [16] lattice_0.19-13mvoutlier_1.4 plyr_1.2.1 [19] RColorBrewer_1.0-2 chemometrics_1.0 som_0.3-5 [22] robustbase_0.5-0-1 rpart_3.1-46 pls_2.1-0 [25] pcaPP_1.8-3mvtnorm_0.9-92 nnet_7.3-1 [28] mclust_3.4.6 MASS_7.3-8 lars_0.9-7 [31] gclus_1.3 cluster_1.13.1 e1071_1.5-24 [34] class_7.3-2 __ 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. -- Dr Paul Murrell Department of Statistics The University of Auckland Private Bag 92019 Auckland New Zealand 64 9 3737599 x85392 p...@stat.auckland.ac.nz http://www.stat.auckland.ac.nz/~paul/ __ 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] Can't read data coded in Cyrillic
To whom it may concern, I have a database (an SPSS .sav file) with some vectors containing strings of words in Cyrillic. My console, however, doesn't seem to read Cyrillics. If I type ÐÑÐ¸Ð²ÐµÑ into the console it shows Ãðèâåò as I'm typing. Please note that courier new, my console's default font can encode Cyrillics elsewhere on my computer (e.g. MSWord and Chrome, from which this email is being sent). Please help me to understand why R can't do this. Thanks, Dan -- Daniel Winetsky MD, MS Candidate Stanford University School of Medicine [[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] How to this SAS transport file in R?
> -Original Message- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > project.org] On Behalf Of zhiji19 > Sent: Sunday, December 05, 2010 9:26 PM > To: r-help@r-project.org > Subject: [R] How to this SAS transport file in R? > > > Dear All, > > I try to read the SAS transport file in R, but it shows error. Please > help! > I am using R 2.11.1 > > library(foreign) > download.file("http://isites.harvard.edu/fs/docs/icb.topic35387.files/d > emo_c.xpt","C:/Desktop/demo_c.xpt") > sasxport <- read.xport("C:/Desktop/demo_c.xpt") > > Error in lookup.xport(file) : file not in SAS transfer format > > I think You need to use mode='wb' in your download statement download.file("http://isites.harvard.edu/fs/docs/icb.topic35387.files/demo_c.xpt";, "C:/Desktop/demo_c.xpt", mode='wb') Hope this is helpful, Dan Daniel J. Nordlund Washington State Department of Social and Health Services Planning, Performance, and Accountability Research and Data Analysis Division Olympia, WA 98504-5204 __ 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] How to this SAS transport file in R?
I've had good success with the read.xport function in the SASxport (not foreign) package. But couldn't you just use something like mydata<-read.xport('http/demo_c.xpt') The transport file looks like it's one of the NHANES demographic files. -- View this message in context: http://r.789695.n4.nabble.com/How-to-this-SAS-transport-file-in-R-tp3073969p3074009.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to this SAS transport file in R?
You didn't tell us your OS. That is a binary file, and from ?download.file: mode: character. The mode with which to write the file. Useful values are ‘"w"’, ‘"wb"’ (binary), ‘"a"’ (append) and ‘"ab"’. Only used for the ‘"internal"’ method. Only on Windows are the binary modes different, so it looks like you forgot mode="wb". On Sun, 5 Dec 2010, someone ashamed of her real name wrote: Dear All, I try to read the SAS transport file in R, but it shows error. Please help! I am using R 2.11.1 library(foreign) download.file("http://isites.harvard.edu/fs/docs/icb.topic35387.files/demo_c.xpt","C:/Desktop/demo_c.xpt";) sasxport <- read.xport("C:/Desktop/demo_c.xpt") Error in lookup.xport(file) : file not in SAS transfer format [...] PLEASE do read the posting guide http://www.R-project.org/posting-guide.html Yes, that does mean you. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595__ 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 to this SAS transport file in R?
Dear All, I try to read the SAS transport file in R, but it shows error. Please help! I am using R 2.11.1 library(foreign) download.file("http://isites.harvard.edu/fs/docs/icb.topic35387.files/demo_c.xpt","C:/Desktop/demo_c.xpt";) sasxport <- read.xport("C:/Desktop/demo_c.xpt") Error in lookup.xport(file) : file not in SAS transfer format -- View this message in context: http://r.789695.n4.nabble.com/How-to-this-SAS-transport-file-in-R-tp3073969p3073969.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] Runif Help: same variable, 3 different parameters
So I am working on an economic model and I need to change the parameters of the runif statement as time goes on. Ex. X <-runif(1:50,0,5) X <-runif(51:100,100,150) X <-runif(100:T, 1,2) Not sure how to go about entering this in to R properly. -- View this message in context: http://r.789695.n4.nabble.com/Runif-Help-same-variable-3-different-parameters-tp3073894p3073894.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] Runif Help: same variable, 3 different parameters
So I am working on an economic model and I need to change the parameters of the runif statement as time goes on. Ex. X <-runif(1:50,0,5) X <-runif(51:100,100,150) X <-runif(100:T, 1,2) Not sure how to go about entering this in to R properly. -- View this message in context: http://r.789695.n4.nabble.com/Runif-Help-same-variable-3-different-parameters-tp3073893p3073893.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to catch both warnings and errors?
On Dec 5, 2010, at 7:35 PM, Marius Hofert wrote: On 2010-12-06, at 01:07 , David Winsemius wrote: On Dec 5, 2010, at 3:13 PM, Marius Hofert wrote: Dear expeRts, I am struggling with warning/error handling. I would like to call a function which can produce either a) normal output b) a warning c) an error Since the function is called several (thousand) times in a loop, I would like to proceed "quietly" and collect the warnings and errors [to deal with them at a later point]. I do not see the function warnings() being used below: ?warnings It delivers the stored warnings with different default behavior for interactive and non-interactive sessions. I have seen constructs with tryCatch (which can deal with errors) and with withCallingHandlers (which can deal with warnings), but I cannot figure out how to catch *both* warnings and errors. Below is a minimal example of the function that is called several times in a large loop. The function should catch warnings and errors; the former work fine, but with the latter I do not know how to proceed. Made some changes in you code but don't know if it is what you were hoping for: f <- function(x){ ## warnings w.list <- NULL # init warning w.handler <- function(w){ # warning handler warn <- simpleWarning(w$message, w$call) # build warning # first change here w.list <<- c(w.list, paste(warnings(), collapse = " ")) # save warning invokeRestart("muffleWarning") } ## errors e.list <- NULL # init errors # not sure this is good idea e.handler <- function(e){ # error handler err <- c(e.list, e$message, e$call) # save error return( err) } ## execute command # wrapped cal in try() res <- withCallingHandlers(try(log(x)), warning = w.handler, error = e.handler) ## return result with warnings and errors list(result = res, warning = w.list, error = e.list) } Dear David, many thanks for your help. If I call your code with f(-1) and f("a"), I obtain: f(-1) $result [1] NaN $warning [1] "" $error NULL => The problem is that the warning is not given. f("a") Error in log(x) : Non-numeric argument to mathematical function $result [1] "Error in log(x) : Non-numeric argument to mathematical function \n" attr(,"class") [1] "try-error" $warning NULL $error NULL => The problem is that the error message is printed to the R console instead of suppressed (setting silent = TRUE didn't help either). Further, the $error component is empty (the error message should appear there -- if possible) Sorry. I was being misled by warnings that existed at my global level into thinking i had succeeded. This modification will suppress warning and error but will not populate the lists as we had hoped: f <- function(expr){ ## warnings w.list <- NULL # init warning w.handler <- function(w){ # warning handler warn <- c(w$message, w$call) # build warning # first change here muffleWarning <<- c(w.list, warn, collapse = " ") # save warning invokeRestart("muffleWarning")} ## errors e.list <- NULL # init errors # not sure this is good idea e.handler <- function(e){ # error handler e.list <<- c(e.list, e$message, e$call) # save error return( e.list) } ## execute command # wrapped cal in try() res <- withCallingHandlers(try(expr, silent=TRUE), warning = w.handler, error = e.handler) ## return result with warnings and errors list(result = res, warning = w.list, error = e.list) } > test <- f(log(-1)) > test $result [1] NaN $warning NULL $error NULL > test <- f(log("a")) > test $result [1] "Error in log(\"a\") : Non-numeric argument to mathematical function\n" attr(,"class") [1] "try-error" $warning NULL $error NULL Do you know a solution? Cheers, Marius David Winsemius, MD West Hartford, CT __ 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] How to get lasso fit coefficient(given penalty tuning parameter \lambda) using lars package
Oh, I just realize that I can actually just use this: beta = coef(lassofit, s=5, mode="lambda") Michael On Sun, Dec 5, 2010 at 8:10 PM, michael wrote: > Hi, all, > > I am using the lars package for lasso estimate. So I get a lasso > fit first: > > lassofit = lars(x,y,type ="lasso",normalize=T, intercept=T) > Then I want to get coefficient with respect to a certain value of \lambda > (the tuning parameter), I know lars has three mode options c("step", > "fraction", "norm"), but can I use the \lambda value instead of the norm > which is the L1 norm of coefficients? > > Thanks, > > Michael > [[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 to get lasso fit coefficient(given penalty tuning parameter \lambda) using lars package
Hi, all, I am using the lars package for lasso estimate. So I get a lasso fit first: lassofit = lars(x,y,type ="lasso",normalize=T, intercept=T) Then I want to get coefficient with respect to a certain value of \lambda (the tuning parameter), I know lars has three mode options c("step", "fraction", "norm"), but can I use the \lambda value instead of the norm which is the L1 norm of coefficients? Thanks, Michael [[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] Hmisc label function applied to data frame
Hello again, I have found that if I use sapply, I do not get a warning, i.e., lbl1 = sapply(d[,var1],label) works correctly and gives no warning. I'm sorry this did not occur to me earlier, my apologies! --Krishna On Thu, Dec 2, 2010 at 11:36 AM, Krishna Tateneni wrote: > Hello, > > I'm attempting to create a data frame with correlations between every pair > of variables in a data frame, so that I can then sort by the value of the > correlation coefficient and see which pairs of variables are most strongly > correlated. > > The sm2vec function in the corpcor library works very nicely as shown here: > > library(Hmisc) > library(corpcor) > > # Create example data > x1 = runif(50) > x2 = runif(50) > x3 = runif(50) > d = data.frame(x1=x1,x2=x2,x3=x3) > label(d$x1) = "Variable x1" > label(d$x2) = "Variable x2" > label(d$x3) = "Variable x3" > > # Get correlations > cormat = cor(d) > > # Get vector form of lower triangular elements > cors = sm2vec(cormat,diag=F) > inds = sm.index(cormat,diag=F) > > # Create a data frame > var1 = dimnames(cormat)[[1]][inds[,1]] > var2 = dimnames(cormat)[[2]][inds[,2]] > lbl1 = label(d[,var1]) > lbl2 = label(d[,var2]) > cor_df = data.frame(Var1=lbl1,Var2=lbl2,Cor=cors) > > The issue that I run into is when trying to get the labels in lbl1 and > lbl2. I get the warning: > > In mapply(FUN = label, x = x, default = default, MoreArgs = list(self = > TRUE), : > longer argument not a multiple of length of shorter > > My usage of label seems ambiguous since the data frame could also a label > attached to it, aside from labels attached to variables within the data > frame. However, the code above does work, with the warning. Aside from > using a loop to get the label of one variable at a time, is there any other > way of getting the labels for all variables in the data frame? > > Also, if there is a better way to achieve my goal of getting the > correlations between all variable pairs, I'd love to know. > > Thanks in advance for any responses! > > --Krishna > [[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] How to catch both warnings and errors?
On 2010-12-06, at 01:07 , David Winsemius wrote: > > On Dec 5, 2010, at 3:13 PM, Marius Hofert wrote: > >> Dear expeRts, >> >> I am struggling with warning/error handling. >> >> I would like to call a function which can produce either >> a) normal output >> b) a warning >> c) an error >> >> Since the function is called several (thousand) times in a loop, I would like >> to proceed "quietly" and collect the warnings and errors [to deal with them >> at a >> later point]. > > I do not see the function warnings() being used below: > > ?warnings > > It delivers the stored warnings with different default behavior for > interactive and non-interactive sessions. >> >> I have seen constructs with tryCatch (which can deal with errors) and with >> withCallingHandlers (which can deal with warnings), but I cannot figure out >> how >> to catch *both* warnings and errors. Below is a minimal example of the >> function >> that is called several times in a large loop. The function should catch >> warnings >> and errors; the former work fine, but with the latter I do not know how to >> proceed. >> > > > Made some changes in you code but don't know if it is what you were hoping > for: > > f <- function(x){ > ## warnings > w.list <- NULL # init warning > w.handler <- function(w){ # warning handler > warn <- simpleWarning(w$message, w$call) # build warning > # first change here > w.list <<- c(w.list, paste(warnings(), collapse = " ")) # save warning > invokeRestart("muffleWarning") > } > ## errors > e.list <- NULL # init errors # not sure this is good idea > e.handler <- function(e){ # error handler >err <- c(e.list, e$message, e$call) # save error >return( err) >} > ## execute command > # wrapped cal in try() > res <- withCallingHandlers(try(log(x)), warning = w.handler, error = > e.handler) > ## return result with warnings and errors > list(result = res, warning = w.list, error = e.list) > } > Dear David, many thanks for your help. If I call your code with f(-1) and f("a"), I obtain: > f(-1) $result [1] NaN $warning [1] "" $error NULL => The problem is that the warning is not given. > f("a") Error in log(x) : Non-numeric argument to mathematical function $result [1] "Error in log(x) : Non-numeric argument to mathematical function\n" attr(,"class") [1] "try-error" $warning NULL $error NULL => The problem is that the error message is printed to the R console instead of suppressed (setting silent = TRUE didn't help either). Further, the $error component is empty (the error message should appear there -- if possible) Do you know a solution? Cheers, Marius __ 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] How to catch both warnings and errors?
On Dec 5, 2010, at 3:13 PM, Marius Hofert wrote: Dear expeRts, I am struggling with warning/error handling. I would like to call a function which can produce either a) normal output b) a warning c) an error Since the function is called several (thousand) times in a loop, I would like to proceed "quietly" and collect the warnings and errors [to deal with them at a later point]. I do not see the function warnings() being used below: ?warnings It delivers the stored warnings with different default behavior for interactive and non-interactive sessions. I have seen constructs with tryCatch (which can deal with errors) and with withCallingHandlers (which can deal with warnings), but I cannot figure out how to catch *both* warnings and errors. Below is a minimal example of the function that is called several times in a large loop. The function should catch warnings and errors; the former work fine, but with the latter I do not know how to proceed. Made some changes in you code but don't know if it is what you were hoping for: f <- function(x){ ## warnings w.list <- NULL # init warning w.handler <- function(w){ # warning handler warn <- simpleWarning(w$message, w$call) # build warning # first change here w.list <<- c(w.list, paste(warnings(), collapse = " ")) # save warning invokeRestart("muffleWarning") } ## errors e.list <- NULL # init errors # not sure this is good idea e.handler <- function(e){ # error handler err <- c(e.list, e$message, e$call) # save error return( err) } ## execute command # wrapped cal in try() res <- withCallingHandlers(try(log(x)), warning = w.handler, error = e.handler) ## return result with warnings and errors list(result = res, warning = w.list, error = e.list) } The function should *always* return the list with the three components. How can I achieve this? Cheers, Marius Ps: How can I get the warning/error message in a nice string (as it is printed in a normal warning/error message)? The approach shown below is somehow ugly. Basically, I had trouble converting the call w$call to a string. ## based on http://tolstoy.newcastle.edu.au/R/help/04/06/0217.html ## and http://www.mail-archive.com/r-help@r-project.org/msg81380.html f <- function(x){ ## warnings w.list <- NULL # init warning w.handler <- function(w){ # warning handler warn <- simpleWarning(w$message, w$call) # build warning w.list <<- c(w.list, paste(warn, collapse = " ")) # save warning invokeRestart("muffleWarning") } ## errors e.list <- NULL # init errors e.handler <- function(e){ # error handler err <- simpleError(e$message, e$call) e.list <<- c(e.list, paste(err, collapse = " ")) # save error } ## execute command res <- withCallingHandlers(log(x), warning = w.handler, error = e.handler) ## return result with warnings and errors list(result = res, warning = w.list, error = e.list) } f(1) f(-1) f("a") __ 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 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Error in calcCurveGrob(x, x$debug) : End points must not be identical
Thanks Paul... I wasn't sure in what way the end points were identical and was testing differently, and failing to catch it. I'll have to check my computations now to see why I'm getting curves that meet this condition, as I don't think I should be. But I'm back on track now. The problem with the viewport being wrong was limited to the toy example I made. I'll send you a graphic directly so you can see what I'm working on. Thanks again for the correct "test" for identical endpoints. Should have been able to see that one myself! Bryan * Bryan Hanson Professor of Chemistry & Biochemistry DePauw University, Greencastle IN USA On 12/5/10 5:17 PM, "Paul Murrell" wrote: > Hi > > Bryan Hanson wrote: >> Hi All... I haven¹t found mention of this error anywhere. I'm trying to >> draw spline curves using grid graphics. Most of the time, I have no >> problems, but I have some data sets that give the error in the subject line. >> I'm not sure which end points are identical, but the end points passed to >> the function are definitely not identical. > > I get ... > > which(tst$x.st == tst$x.end & tst$y.st == tst$y.end) > [1] 11 14 > > Also, the viewport you have set up will not show any of the curves > because its "native" scale is 0 to 1. The following at least makes the > curves visible ... > > grid.newpage() > vp <- viewport(width = 0.9, height = 0.9, x = 0.5, y = 0.5, > xscale=c(-10, 10), yscale=c(-15, 15)) > pushViewport(vp) > grid.rect(gp = gpar(lty = "dashed", col = "gray")) > grid.points(0.5, 0.5, pch = 20, gp = gpar(cex = 0.5)) > subset <- -c(11, 14) > grid.curve(tst$x.st[subset], tst$y.st[subset], > tst$x.end[subset], tst$y.end[subset], > default.units = "native") > > Paul > > p.s. Intrigued to know what sort of image you are producing, if you are > able to share. > >> >> Any assistance appreciated! Bryan >> >> tst <- >> structure(list(x.st = c(-1, -2, -3, -1, -1.5, -3, -1.5, -1.5, >> -8, -1, -1.5, -1, -1.5, -2, -1.5, -2, -1, -1.5, -2), y.st = >> c(1.73205080756888, >> 3.46410161513776, 5.19615242270663, 1.73205080756888, 2.59807621135332, >> 5.19615242270663, 2.59807621135332, 2.59807621135332, 13.8564064605510, >> 1.73205080756888, 2.59807621135332, 1.73205080756888, 2.59807621135332, >> 3.46410161513776, 2.59807621135332, 3.46410161513776, 1.73205080756888, >> 2.59807621135332, 3.46410161513776), x.end = c(-6.5, -6.5, -6.5, >> -4, -4, -4, -1.5, -1, -1, -1.5, -1.5, -2, -2, -2, -3, -3, -8, >> -8, -8), y.end = c(-11.2583302491977, -11.2583302491977, -11.2583302491977, >> -6.92820323027551, -6.92820323027551, -6.92820323027551, -2.59807621135332, >> 1.73205080756888, 1.73205080756888, 2.59807621135332, 2.59807621135332, >> 3.46410161513776, 3.46410161513776, 3.46410161513776, 5.19615242270663, >> 5.19615242270663, 13.8564064605510, 13.8564064605510, 13.8564064605510 >> ), grp = c(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, >> 3, 3), lwd = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 3, 1, 1, 1, >> 2, 5, 2)), .Names = c("x.st", "y.st", "x.end", "y.end", "grp", >> "lwd"), row.names = 34:52, class = "data.frame") >> >> grid.newpage() >> vp <- viewport(width = 0.9, height = 0.9, x = 0.5, y = 0.5) >> pushViewport(vp) >> grid.rect(gp = gpar(lty = "dashed", col = "gray")) >> grid.points(0.5, 0.5, pch = 20, gp = gpar(cex = 0.5)) >> grid.curve(tst$x.st, tst$y.st, tst$x.end, tst$y.end, >> default.units = "native") >> >> ### >>> sessionInfo() >> R version 2.12.0 (2010-10-15) >> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) >> >> locale: >> [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 >> >> attached base packages: >> [1] datasets tools grid grDevices graphics utils >> [7] stats methods base >> >> other attached packages: >> [1] gridExtra_0.7 GGally_0.2.2 xtable_1.5-6 >> [4] mvbutils_2.5.1 ggplot2_0.8.8 proto_0.3-8 >> [7] reshape_0.8.3 ChemoSpec_1.46 seriation_1.0-2 >> [10] colorspace_1.0-1 TSP_1.0-1 R.utils_1.5.3 >> [13] R.oo_1.7.4 R.methodsS3_1.2.1 rgl_0.92.794 >> [16] lattice_0.19-13mvoutlier_1.4 plyr_1.2.1 >> [19] RColorBrewer_1.0-2 chemometrics_1.0 som_0.3-5 >> [22] robustbase_0.5-0-1 rpart_3.1-46 pls_2.1-0 >> [25] pcaPP_1.8-3mvtnorm_0.9-92 nnet_7.3-1 >> [28] mclust_3.4.6 MASS_7.3-8 lars_0.9-7 >> [31] gclus_1.3 cluster_1.13.1 e1071_1.5-24 >> [34] class_7.3-2 >> >> __ >> 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,
Re: [R] lm() and interactions in model formula for x passed as matrix
Thanks for the replies. I was just thinking that, for a two variable example, doing X<-cbind(x1,x2,x1*x2) lm(y~X) would work. So maybe that's what I'll do. This also allows me to pick and choose which interactions to include. Cheers Bill On Sun, Dec 5, 2010 at 8:19 PM, William Simpson wrote: > Suppose I have x variables x1, x2, x3 (however in general I don't know > how many x variables there are). I can do > X<-cbind(x1,x2,x3) > lm(y ~ X) > This fits the no-interaction model with b0, b1, b2, b3. > > How can I get lm() to fit the model that includes interactions when I > pass X to lm()? For my example, > lm(y~x1*x2*x3) > I am looking for something along the lines of > lm(y~X ...) > where ... is some extra stuff I need to fill in. > > Thanks for any help. > Bill > __ 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] The behaviour of read.csv().
On 05/12/2010 2:14 PM, Rolf Turner wrote: On 6/12/2010, at 3:00 AM, Duncan Murdoch wrote: I was going to suggest using DIF rather than CSV. It contains more internal information about the file (including the type of each entry), but has the disadvantage of being less readable, even though it is ascii. I don't think DIF is really the answer. My colleagues are familiar with the *.csv concept; they have never heard of ``DIF''. As I have said, we have had but few problems using *.csv. Better the devil you know ... Furthermore I have to deal with data provided by various sources ``external'' to the research project that I work for. I have to use the data that these sources provide, in the format in which they provide it. If they give me *.csv files I count myself lucky. Finally, there seems to be no ``write.DIF'' function, i.e. there is no way to produce *.DIF output, as far as I can tell. Hence it would not seem practical to use *.DIF as a data exchange standard. Sure, those are good points. However, in putting together a little demo, I found a couple of bugs in the R implementation of read.DIF, and it looks as though it ignores the internal type information. Sigh. As of r53778, the bugs I noticed should be fixed. read.DIF now respects the internal type information, so it will keep character strings like "001" as type character (unless you ask it to change the type). What does ``r53778'' mean? Revision 53778 from the version control system. When you start R-patched or R-devel it will print this in the startup message, e.g. R version 2.13.0 Under development (unstable) (2010-12-05 r53775) ^^ (from just before I saved the changes). Duncan Murdoch cheers, Rolf __ 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] lm() and interactions in model formula for x passed as matrix
Hi Bill, If you can put all (and only) your variables into a dataframe, (for example: X <- data.frame(y, x1, x2, x3) ) then another alternative to David's solution would be: lm(y ~ .^3, data = X) '.' will expand to every column except y, and then the ^3 will get you up to 3-way interactions. Cheers, Josh On Sun, Dec 5, 2010 at 12:19 PM, William Simpson wrote: > Suppose I have x variables x1, x2, x3 (however in general I don't know > how many x variables there are). I can do > X<-cbind(x1,x2,x3) > lm(y ~ X) > This fits the no-interaction model with b0, b1, b2, b3. > > How can I get lm() to fit the model that includes interactions when I > pass X to lm()? For my example, > lm(y~x1*x2*x3) > I am looking for something along the lines of > lm(y~X ...) > where ... is some extra stuff I need to fill in. > > Thanks for any help. > Bill > > __ > 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. > -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.com/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lm() and interactions in model formula for x passed as matrix
On Dec 5, 2010, at 3:19 PM, William Simpson wrote: Suppose I have x variables x1, x2, x3 (however in general I don't know how many x variables there are). I can do X<-cbind(x1,x2,x3) lm(y ~ X) This fits the no-interaction model with b0, b1, b2, b3. How can I get lm() to fit the model that includes interactions when I pass X to lm()? For my example, lm(y~x1*x2*x3) I am looking for something along the lines of lm(y~X ...) where ... is some extra stuff I need to fill in. The formula syntax in R allows you to specify interactions with the "^" operator but some testing makes me think you cannot use either y ~ .^3 or y ~ X^3 with matrix data arguments, here assuming you only want interaction up to third order. Assuming you know how to use do.call("cbind", varlist) perhaps: form = as.formula( paste("y ~ (", paste(colnames(X), collapse="+"), ")^3", sep="") ) lm(form) --- output-- Call: lm(formula = form) Coefficients: (Intercept) x1 x2 x3x1:x2 x1:x3 -0.383296-0.333429 0.003976 0.332982-0.001130 0.100698 x2:x3 x1:x2:x3 0.366745 0.122111 -- David Winsemius, MD West Hartford, CT __ 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] "less than or equal to" glyph
Hi, I suspect this has to do with your locale, as it seems to work on my system. Take a look at: Sys.getlocale() # see yours ?Sys.setlocale # documentation If you're are talking about something to add to a graph (e.g., in the title or labelling a point) try: sd_label <- substitute(LSD~(P <= 0.05) == n, list(n = 5)) plot(1:10, main = sd_label) HTH, Josh On Sun, Dec 5, 2010 at 11:25 AM, romzero wrote: > > If i insert \u2264 inside the text, like this > > lsd_label <- "LSD (P \u2264 0.05) = " > > the result is > > "LSD (P = 0.05) = " > > instead > > "LSD (P ≤ 0.05) = " > > how can i solve this problem? > > Thanks. > -- > View this message in context: > http://r.789695.n4.nabble.com/less-than-or-equal-to-glyph-tp3073557p3073557.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. > -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] lm() and interactions in model formula for x passed as matrix
Suppose I have x variables x1, x2, x3 (however in general I don't know how many x variables there are). I can do X<-cbind(x1,x2,x3) lm(y ~ X) This fits the no-interaction model with b0, b1, b2, b3. How can I get lm() to fit the model that includes interactions when I pass X to lm()? For my example, lm(y~x1*x2*x3) I am looking for something along the lines of lm(y~X ...) where ... is some extra stuff I need to fill in. Thanks for any help. Bill __ 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 to catch both warnings and errors?
Dear expeRts, I am struggling with warning/error handling. I would like to call a function which can produce either a) normal output b) a warning c) an error Since the function is called several (thousand) times in a loop, I would like to proceed "quietly" and collect the warnings and errors [to deal with them at a later point]. I have seen constructs with tryCatch (which can deal with errors) and with withCallingHandlers (which can deal with warnings), but I cannot figure out how to catch *both* warnings and errors. Below is a minimal example of the function that is called several times in a large loop. The function should catch warnings and errors; the former work fine, but with the latter I do not know how to proceed. The function should *always* return the list with the three components. How can I achieve this? Cheers, Marius Ps: How can I get the warning/error message in a nice string (as it is printed in a normal warning/error message)? The approach shown below is somehow ugly. Basically, I had trouble converting the call w$call to a string. ## based on http://tolstoy.newcastle.edu.au/R/help/04/06/0217.html ## and http://www.mail-archive.com/r-help@r-project.org/msg81380.html f <- function(x){ ## warnings w.list <- NULL # init warning w.handler <- function(w){ # warning handler warn <- simpleWarning(w$message, w$call) # build warning w.list <<- c(w.list, paste(warn, collapse = " ")) # save warning invokeRestart("muffleWarning") } ## errors e.list <- NULL # init errors e.handler <- function(e){ # error handler err <- simpleError(e$message, e$call) e.list <<- c(e.list, paste(err, collapse = " ")) # save error } ## execute command res <- withCallingHandlers(log(x), warning = w.handler, error = e.handler) ## return result with warnings and errors list(result = res, warning = w.list, error = e.list) } f(1) f(-1) f("a") __ 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] Problem storing lm() model in a list
Hi Tal, That is correct. In the code, I actually call the same exact data.frame that was used to create the object. I understand I can call predict() without a data.frame to get this result, but I would like to predict other datasets as well. Thanks, Harold On Dec 5, 2010, at 7:18 AM, Tal Galili wrote: > Hi Harold, > I think the error stems from the data.frame you are entering into the predict > function. > Your data.frame must have the EXACT column name as the one used for the > creation of the model object. > Is that the case in your code? > > > Best, > Tal > > Contact > Details:--- > Contact me: tal.gal...@gmail.com | 972-52-7275845 > Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | > www.r-statistics.com (English) > -- > > > > > On Sat, Dec 4, 2010 at 9:30 AM, Harold Pimentel > wrote: > Hi all, > > I recently wrote some code to store a number of polynomial models in a list > and return this list. The model is returned fine, but when I make a > subsequent call to predict(), I have an error. The code for polyModel > selection is listed at the end of the email. Here is an example of the error: > > > poly.fit <- polyModelSelection(x,y,10,F) > > for (d in 1:4) { > + lm.model <- poly.fit$models[[d]] > + curve(predict(lm.model,data.frame(x)),add=T,lty=d+1) > + } > Error: variable 'poly(x, d, raw = T)' was fitted with type "nmatrix.1" but > type "nmatrix.10" was supplied > > Does anyone have any ideas? > > > > Thanks, > > Harold > > > > > polyModelSelection <- function(x,y,maxd,add.dim=F) { >dim.mult <- 0 >if (add.dim) { >dim.mult = 1 >} >bestD <- 1 >bestError <- 1 >loss <- 1 >lm.models <- vector("list",maxd) >for (d in 1:maxd) { >lm.mod <- lm(y ~ 0 + poly(x,d,raw=T)) >lm.models[[d]] <- lm.mod >loss[d] <- modelError(lm.mod,data.frame(x),y)+d*dim.mult >if (d == 1){ >bestError <- loss[d] >} >if (loss[d] < bestError) { >bestError <- loss[d] >bestD <- d >} ># print(loss[d]) >} >list(loss=loss,models=lm.models,best=c(bestD)) > } > > modelError <- function(model,x,y) { >sum((y-predict(model,x))^2) > } > __ > 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] "less than or equal to" glyph
If i insert \u2264 inside the text, like this lsd_label <- "LSD (P \u2264 0.05) = " the result is "LSD (P = 0.05) = " instead "LSD (P ≤ 0.05) = " how can i solve this problem? Thanks. -- View this message in context: http://r.789695.n4.nabble.com/less-than-or-equal-to-glyph-tp3073557p3073557.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] HOW to use the survivalROC to get optimal cut-off values?
On Dec 5, 2010, at 11:33 AM, David Winsemius wrote: On Dec 5, 2010, at 11:14 AM, petre...@unina.it wrote: I have the same problem of a prevous request HOW to use the survivalROC (or another library in R) to get optimal cut-off values? I want to use the time-dependent survivalROC package.according to the,reference material,it only gives a set of ordered cut-off values .eg. Optimality specification requires some sort of valuation of incorrect decisions. If you are willing to defend a choice that a false positive has exactly the same loss as a false negative, which is generally not the case in medical decision-making, then the point on the ROC curve which is closest to the upper left-hand corner is "optimal". Snipped comments about small nubers of discrete values as that was a misunderstanding on my part. data(mayo) str(mayo) attach(mayo) ROC. 1 = survivalROC (Stime =time,status=censor,marker=mayoscore4,predict.time=365,lambda=0.05) str(ROC.1) plot(ROC.1$FP, ROC.1$TP, type="l", xlim=c(0,1), ylim=c(0,1), xlab=paste( "FP", "\n", "AUC = ",round(ROC.1$AUC,3)), ylab="TP",main="Mayoscore 4, Method = NNE \n Year = 1") abline(0,1) List of 6 $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ... *only 5 values * $ TP : num [1:313] 1 0.999 0.999 0.999 0.998 ... $ FP : num [1:313] 1 0.997 0.993 0.99 0.987 ... $ predict.time: num 365 $ Survival: num 0.93 $ AUC : num 0.888 so i dont know how to use the survivalROC to get optimal cut-off values?(only 5 values) I really do not understand what you are asking when you ask about "5 values". (I also see that another student GY QIAN in SHANGHAI,CHINA asked almost the same odd question on rhelp about 6 weeks ago making me wonder if there is some online class you are both taking for which this is homework?) There are 312 finite values for each of cut.values, FP, and TP. Is someone, your course instructor or collaborator, asking for a comparison at 5 selected cutpoints? I generally use Hmisc::describe for quick looks: > describe(ROC.1$cut.values[-1]) ROC.1$cut.values[-1] n missing uniqueMean .05 .10 .25 .50 . 75 .90 .95 312 0 312 6.538 5.320 5.459 5.842 6.312 6.949 7.966 8.675 lowest : 4.581 4.900 4.926 4.932 4.946, highest: 9.510 9.568 10.185 10.479 10.629 You could programmatically ask what values of the index (1:313) minimizes the sum of FN and (1-TP) = FP. And then you can use to index the TP and "FP" values. I guess it is common practice to call 1- specificity the "false positive rate" but I for one find that very confusing. even misleading, since it's not clear from the term what the denominator for such a "rate" really should be. (It's also not really a rate since no time is involved in the calculation.) At any rate, as it were: > with(ROC.1, which.min(1-TP+ FP)) [1] 259 > with(ROC.1, points(FP[259], TP[259], cex=3, col="red" )) -- David Winsemius, MD West Hartford, CT __ 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] Problem storing lm() model in a list
Your problem arises because predict() is using the value d=10 whenever it evaluates poly(x,d,raw=TRUE). The details are that this is the value that d in your function polyModelSelection had when the function finished. That value is stored in the environment environment(poly.fit$models[[i]]$terms) for i in 1:10 (which is the same environment for all 10 entries). It contains the following objects > objects(envir=environment(poly.fit$models[[1]]$terms)) [1] "add.dim" "bestD" "bestError" "d" [5] "dim.mult" "formula" "lm.mod""lm.models" [9] "loss" "maxd" "x" "y" and predict will look there for objects named in the formula but not found in its newdata= argument. Use get("d", envir=...) to see that d is 10 in all of them. I don't know the best way to avoid the problem. One way is to put a literal value for degree into your call to poly(), which also has the advantage that your printouts show what it is. In the following I used substitute() to put the explicit value of degree into your formula and used do.call when calling lm so the formula itself (not the name "formula") is put into the call component of the lm object. This seems kludgy to me, but I think it gets the job done. polyModelSelection <- function (x, y, maxd, add.dim = F) { dim.mult <- 0 if (add.dim) { dim.mult = 1 } bestD <- 1 bestError <- 1 loss <- 1 lm.models <- vector("list", maxd) for (d in 1:maxd) { # next assignment added formula <- substitute(y ~ 0 + poly(x, degree = d, raw = TRUE), list(d = d)) # change lm(...) to do.call("lm", list(...)) lm.mod <- do.call("lm", list(formula)) lm.models[[d]] <- lm.mod loss[d] <- modelError(lm.mod, data.frame(x), y) + d * dim.mult if (d == 1) { bestError <- loss[d] } if (loss[d] < bestError) { bestError <- loss[d] bestD <- d } } list(loss = loss, models = lm.models, best = c(bestD)) } Another way around the problem might be to replace the environment stored in those terms compononents with emptyenv(). Then none of the variables in the function polyModelSelection could be used. You could do that with for(d in seq_along(poly.fit$models)) { environment(poly.fit$models[[1]]$terms) <- emptyenv() } I like the first was better because I can easily see what d is. This sort of issue arises whenever we use functions that use nonstandard evaluation rules, like the modelling functions, subset(), library(), and help(). They are easy to call from the top level command line but harder to deal with when called with variables for arguments. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com > -Original Message- > From: r-help-boun...@r-project.org > [mailto:r-help-boun...@r-project.org] On Behalf Of Harold Pimentel > Sent: Friday, December 03, 2010 11:31 PM > To: r-help@r-project.org > Subject: [R] Problem storing lm() model in a list > > Hi all, > > I recently wrote some code to store a number of polynomial > models in a list and return this list. The model is returned > fine, but when I make a subsequent call to predict(), I have > an error. The code for polyModel selection is listed at the > end of the email. Here is an example of the error: > > > poly.fit <- polyModelSelection(x,y,10,F) > > for (d in 1:4) { > + lm.model <- poly.fit$models[[d]] > + curve(predict(lm.model,data.frame(x)),add=T,lty=d+1) > + } > Error: variable 'poly(x, d, raw = T)' was fitted with type > "nmatrix.1" but type "nmatrix.10" was supplied > > Does anyone have any ideas? > > > > Thanks, > > Harold > > > ## > ## > > polyModelSelection <- function(x,y,maxd,add.dim=F) { > dim.mult <- 0 > if (add.dim) { > dim.mult = 1 > } > bestD <- 1 > bestError <- 1 > loss <- 1 > lm.models <- vector("list",maxd) > for (d in 1:maxd) { > lm.mod <- lm(y ~ 0 + poly(x,d,raw=T)) > lm.models[[d]] <- lm.mod > loss[d] <- modelError(lm.mod,data.frame(x),y)+d*dim.mult > if (d == 1){ > bestError <- loss[d] > } > if (loss[d] < bestError) { > bestError <- loss[d] > bestD <- d > } > # print(loss[d]) > } > list(loss=loss,models=lm.models,best=c(bestD)) > } > > modelError <- function(model,x,y) { > sum((y-predict(model,x))^2) > } > __ > 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
Re: [R] The behaviour of read.csv().
On Dec 5, 2010, at 2:14 PM, Rolf Turner wrote: On 6/12/2010, at 3:00 AM, Duncan Murdoch wrote: As of r53778, the bugs I noticed should be fixed. read.DIF now respects the internal type information, so it will keep character strings like "001" as type character (unless you ask it to change the type). What does ``r53778'' mean? I assumed it was a version sequence number: http://cran.r-project.org/src/base-prerelease/ cheers, Rolf David Winsemius, MD West Hartford, CT __ 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] The behaviour of read.csv().
On 6/12/2010, at 3:00 AM, Duncan Murdoch wrote: >> I was going to suggest using DIF rather than CSV. It contains more >> internal information about the file (including the type of each entry), >> but has the disadvantage of being less readable, even though it is ascii. I don't think DIF is really the answer. My colleagues are familiar with the *.csv concept; they have never heard of ``DIF''. As I have said, we have had but few problems using *.csv. Better the devil you know ... Furthermore I have to deal with data provided by various sources ``external'' to the research project that I work for. I have to use the data that these sources provide, in the format in which they provide it. If they give me *.csv files I count myself lucky. Finally, there seems to be no ``write.DIF'' function, i.e. there is no way to produce *.DIF output, as far as I can tell. Hence it would not seem practical to use *.DIF as a data exchange standard. >> >> However, in putting together a little demo, I found a couple of bugs in >> the R implementation of read.DIF, and it looks as though it ignores the >> internal type information. Sigh. > > As of r53778, the bugs I noticed should be fixed. read.DIF now respects > the internal type information, so it will keep character strings like > "001" as type character (unless you ask it to change the type). What does ``r53778'' mean? cheers, Rolf __ 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] Checking for orthogonal contrasts
On Dec 3, 2010, at 17:17 , S Ellison wrote: > David, > > Thanks for the comments. > > I think, though, that I have found the answer to my own post. > >> Question: How would one check, in R, that [contrasts .. are >> 'orthogonal in the row-basis of the model matrix'] for a particular >> fitted linear model object? > > ?lm illustrates the use of crossprod() for probing the orthogonality of > a model matrix. If I understand correctly, the necessary condition is > essentially that all between-term off-diagonal elements of crossprod(m) > are zero if the contrasts are orthogonal, where 'term' refers to the > collection of columns related to a single term in the model formula. > > Example: > > y<-rnorm(27) > g <- gl(3, 9) > h <- gl(3,3,27) > > m1 <- model.matrix(y~g*h, contrasts = list(g="contr.sum", > h="contr.sum")) > crossprod(m1) > > #Compare with > m2 <- model.matrix(y~g*h, contrasts = list(g="contr.treatment", > h="contr.treatment")) > crossprod(m2) > #Note the nonzero off-diagonal elements between, say, g and h or > g, h and the various gi:hj elements > > > That presumably implies that one could test a linear model explicitly > for contrast orthogonality (and, implicitly, balanced design?) using > something like > > model.orthogonal.lm <- function(l) { > #l is a linear model > m <- model.matrix(l) > a <- attr(m, "assign") > a.outer <- outer(a, a, FUN="!=") > m.xprod <- crossprod(m) > all( m.xprod[a.outer] == 0 ) > } > > l1 <- lm(y~g*h, contrasts = list(g="contr.sum", h="contr.sum")) > > l2 <- lm(y~g*h, contrasts = list(g="contr.treatment", > h="contr.treatment")) > > model.orthogonal.lm(l1) > #TRUE > > model.orthogonal.lm(l2) > #FALSE > > Not sure how it would work on balanced incomplete block designs, > though. I'll have to try it. You'll find that the block and treatment terms are NOT orthogonal. That's where all the stuff about "efficiency factors" and "recovery of interblock information" comes from. > > Before I do, though, a) do I have the stats right? and b) this now > seems so obvious that someone must already have done it somewhere... ? a) basically, yes, I think you do b) yes, many, but there is an amazing amount of sloppily thought out "folklore" going around, including the common misconception that somehow sum-to-zero contrasts are inherently better than the other types. What does seem to be the case is just that they have computational advantages in completely balanced designs, because they then imply orthogonality of COLUMNS of the design matrix. That in turn means that you can construct the sum of squares for each model term based on its own columns only. In unbalanced designs, they just tend to give incorrect results... -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] HOW to use str() after the survivalROC (or another library in R) to get optimal cut-off values
Hi Mario, On Sun, Dec 5, 2010 at 10:11 AM, wrote: > > I have the same problem of a previous request > > HOW to use the survivalROC (or another library in R) to get optimal cut-off > values? > > I want to use the time-dependent survivalROC package.according to > the,reference > material,it only gives a partial set of ordered cut-off values .eg. I'm not exactly sure how to parse that sentence. Do you mean you are trying to follow the reference, but only get a partial set of cut-off values? Or do you mean that according to the reference, it only gives a partial set of ordered cut-off values? > > > data(mayo) > str(mayo) > attach(mayo) > ROC.1=survivalROC(Stime=time,status=censor,marker=mayoscore4,predict.time=365,lambda=0.05) > str(ROC.1) this prints the str()ucture of the ROC.1 object to the console so you can see how the object is stored in R > > plot(ROC.1$FP, ROC.1$TP, type="l", xlim=c(0,1), ylim=c(0,1), xlab=paste( > "FP", "\n", > "AUC = ",round(ROC.1$AUC,3)), ylab="TP",main="Mayoscore 4, Method = NNE \n > Year = 1") > abline(0,1) > > List of 6 > $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ... *only 5 values > > * $ TP : num [1:313] 1 0.999 0.999 0.999 0.998 ... > > $ FP : num [1:313] 1 0.997 0.993 0.99 0.987 ... > $ predict.time: num 365 > $ Survival : num 0.93 > > $ AUC : num 0.888 This is the output you should have gotten from str(ROC.1) > > In particular I'm unable to see and print all values after str(ROC.1) > and NOT only the 5 values listed above in the example in order to to get > optimal cut-off values Is this what you are looking for? ROC.1[["cut.values"]] or are you have some other problem? > > Thank you very much! > > > > Mario Petretta > Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche > Facoltà di Medicina e Chirurgia > Università di Napoli Federico II > 081 - 7462233 > > > > Mario Petretta > Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche > Facoltà di Medicina e Chirurgia > Università di Napoli Federico II > 081 - 7462233 > > __ > 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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] HOW to use str() after the survivalROC (or another library in R) to get optimal cut-off values
I have the same problem of a previous request HOW to use the survivalROC (or another library in R) to get optimal cut-off values? I want to use the time-dependent survivalROC package.according to the,reference material,it only gives a partial set of ordered cut-off values .eg. data(mayo) str(mayo) attach(mayo) ROC.1=survivalROC(Stime=time,status=censor,marker=mayoscore4,predict.time=365,lambda=0.05) str(ROC.1) plot(ROC.1$FP, ROC.1$TP, type="l", xlim=c(0,1), ylim=c(0,1), xlab=paste( "FP", "\n", "AUC = ",round(ROC.1$AUC,3)), ylab="TP",main="Mayoscore 4, Method = NNE \n Year = 1") abline(0,1) List of 6 $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ... *only 5 values * $ TP : num [1:313] 1 0.999 0.999 0.999 0.998 ... $ FP : num [1:313] 1 0.997 0.993 0.99 0.987 ... $ predict.time: num 365 $ Survival: num 0.93 $ AUC : num 0.888 In particular I'm unable to see and print all values after str(ROC.1) and NOT only the 5 values listed above in the example in order to to get optimal cut-off values Thank you very much! Mario Petretta Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche Facoltà di Medicina e Chirurgia Università di Napoli Federico II 081 - 7462233 Mario Petretta Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche Facoltà di Medicina e Chirurgia Università di Napoli Federico II 081 - 7462233 __ 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] Converting numbers into words
On Sun, Dec 5, 2010 at 11:50 AM, Thomas Levine wrote: > Example data > > desk=data.frame( > deskchoice=c('mid','mid','left','bookdrop','mid','bookdrop') > ) > > -- > > I like doing stuff like the line below, especially when I'm using Sweave. > > print(paste('Within the observation period,',nrow(desk), > 'patrons approached the circulation desk.')) > > > -- > > But what if I want to put it at the beginning of a sentence? > > print(sum(desk$deskchoice=='bookdrop'),'persons', > 'used the book drop. Everyone else interacted with a staff member.') > > Is there a pretty way to change the result of > sum(desk$deskchoice=='bookdrop') > from "2" to "Two"? > > -- > > And what if the number is one? > > print(sum(desk$deskchoice=='bookdrop'), > c('person','persons')[as.numeric(sum(desk$deskchoice=='bookdrop')!=1)+1], > 'used the book drop. Everyone else interacted with a staff member.') > > Is there a prettier way of choosing between "person" and "persons"? Using John Fox's numbers2words found here: http://tolstoy.newcastle.edu.au/R/help/05/04/2715.html and capwords found in the examples section of ?toupper try this: n <- sum(desk$deskchoice == "bookdrop") paste(capwords(numbers2words(n)), if (n == 1) "person" else "people", "used the book drop...") -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.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] Converting numbers into words
Example data desk=data.frame( deskchoice=c('mid','mid','left','bookdrop','mid','bookdrop') ) -- I like doing stuff like the line below, especially when I'm using Sweave. print(paste('Within the observation period,',nrow(desk), 'patrons approached the circulation desk.')) -- But what if I want to put it at the beginning of a sentence? print(sum(desk$deskchoice=='bookdrop'),'persons', 'used the book drop. Everyone else interacted with a staff member.') Is there a pretty way to change the result of sum(desk$deskchoice=='bookdrop') from "2" to "Two"? -- And what if the number is one? print(sum(desk$deskchoice=='bookdrop'), c('person','persons')[as.numeric(sum(desk$deskchoice=='bookdrop')!=1)+1], 'used the book drop. Everyone else interacted with a staff member.') Is there a prettier way of choosing between "person" and "persons"? -- Thanks Tom __ 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] HOW to use the survivalROC to get optimal cut-off values?
On Dec 5, 2010, at 11:14 AM, petre...@unina.it wrote: I have the same problem of a prevous request HOW to use the survivalROC (or another library in R) to get optimal cut-off values? I want to use the time-dependent survivalROC package.according to the,reference material,it only gives a set of ordered cut-off values .eg. Optimality specification requires some sort of valuation of incorrect decisions. If you are willing to defend a choice that a false positive has exactly the same loss as a false negative, which is generally not the case in medical decision-making, then the point on the ROC curve which is closest to the upper left-hand corner is "optimal". Having only 5 values is getting pretty close to violating the presumption of ROC analysis that the result be at least pseudo- continuous. I have see quite a few "ROC curves in this situation that do not have a clear winner ( now assuming equal cost for FP and FN which I already said was usually a faulty assumption) because the closest point on the curve was in the middle of the line segment between two of the points. I'm not sure that the typical practice of plotting ROC curves with slanting line segments is valid. There is no information between those discrete points. You should probably be using a table rather than a curve in this situation. data(mayo) str(mayo) attach(mayo) ROC. 1 = survivalROC (Stime =time,status=censor,marker=mayoscore4,predict.time=365,lambda=0.05) str(ROC.1) plot(ROC.1$FP, ROC.1$TP, type="l", xlim=c(0,1), ylim=c(0,1), xlab=paste( "FP", "\n", "AUC = ",round(ROC.1$AUC,3)), ylab="TP",main="Mayoscore 4, Method = NNE \n Year = 1") abline(0,1) List of 6 $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ... *only 5 values * $ TP : num [1:313] 1 0.999 0.999 0.999 0.998 ... $ FP : num [1:313] 1 0.997 0.993 0.99 0.987 ... $ predict.time: num 365 $ Survival: num 0.93 $ AUC : num 0.888 so i dont know how to use the survivalROC to get optimal cut-off values?(only 5 values) Thank you very much! Mario Petretta Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche Facoltà di Medicina e Chirurgia Università di Napoli Federico II 081 - 7462233 __ 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 __ 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 to use the survivalROC to get optimal cut-off values?
I have the same problem of a prevous request HOW to use the survivalROC (or another library in R) to get optimal cut-off values? I want to use the time-dependent survivalROC package.according to the,reference material,it only gives a set of ordered cut-off values .eg. data(mayo) str(mayo) attach(mayo) ROC.1=survivalROC(Stime=time,status=censor,marker=mayoscore4,predict.time=365,lambda=0.05) str(ROC.1) plot(ROC.1$FP, ROC.1$TP, type="l", xlim=c(0,1), ylim=c(0,1), xlab=paste( "FP", "\n", "AUC = ",round(ROC.1$AUC,3)), ylab="TP",main="Mayoscore 4, Method = NNE \n Year = 1") abline(0,1) List of 6 $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ... *only 5 values * $ TP : num [1:313] 1 0.999 0.999 0.999 0.998 ... $ FP : num [1:313] 1 0.997 0.993 0.99 0.987 ... $ predict.time: num 365 $ Survival: num 0.93 $ AUC : num 0.888 so i dont know how to use the survivalROC to get optimal cut-off values?(only 5 values) Thank you very much! Mario Petretta Dipartimento di Medicina Clinica Scienze Cardiovascolari e Immunologiche Facoltà di Medicina e Chirurgia Università di Napoli Federico II 081 - 7462233 __ 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] merging two vectors
On Dec 5, 2010, at 9:57 AM, ram basnet wrote: Dear R users, It may be very simple but it is being difficult for me. I have two vectors with some common string. And, i want to combine into a vector in such a way that it includes string from both vectors and make a unique. For example: x <- paste(rep("A",5),1:5,sep = ".") x [1] "A.1" "A.2" "A.3" "A.4" "A.5" y <- paste(rep("A"),3:7, sep = ".") y [1] "A.3" "A.4" "A.5" "A.6" "A.7" Now, I want to combine these two vectors in the following way: "A.1" "A.2" "A.3" "A.4" "A.5" "A.6" "A.7" I tried with merge(), but not able to get as I want. ?union -- David Winsemius, MD West Hartford, CT __ 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] Problem storing lm() model in a list
Hi Harold, I think the error stems from the data.frame you are entering into the predict function. Your data.frame must have the EXACT column name as the one used for the creation of the model object. Is that the case in your code? Best, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Sat, Dec 4, 2010 at 9:30 AM, Harold Pimentel wrote: > Hi all, > > I recently wrote some code to store a number of polynomial models in a list > and return this list. The model is returned fine, but when I make a > subsequent call to predict(), I have an error. The code for polyModel > selection is listed at the end of the email. Here is an example of the > error: > > > poly.fit <- polyModelSelection(x,y,10,F) > > for (d in 1:4) { > + lm.model <- poly.fit$models[[d]] > + curve(predict(lm.model,data.frame(x)),add=T,lty=d+1) > + } > Error: variable 'poly(x, d, raw = T)' was fitted with type "nmatrix.1" but > type "nmatrix.10" was supplied > > Does anyone have any ideas? > > > > Thanks, > > Harold > > > > > > polyModelSelection <- function(x,y,maxd,add.dim=F) { >dim.mult <- 0 >if (add.dim) { >dim.mult = 1 >} >bestD <- 1 >bestError <- 1 >loss <- 1 >lm.models <- vector("list",maxd) >for (d in 1:maxd) { >lm.mod <- lm(y ~ 0 + poly(x,d,raw=T)) >lm.models[[d]] <- lm.mod >loss[d] <- modelError(lm.mod,data.frame(x),y)+d*dim.mult >if (d == 1){ >bestError <- loss[d] >} >if (loss[d] < bestError) { >bestError <- loss[d] >bestD <- d >} ># print(loss[d]) >} >list(loss=loss,models=lm.models,best=c(bestD)) > } > > modelError <- function(model,x,y) { >sum((y-predict(model,x))^2) > } > __ > 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] merging two vectors
unique(c(x,y)) On Sun, Dec 5, 2010 at 8:27 PM, ram basnet wrote: > Dear R users, > > It may be very simple but it is being difficult for me. > I have two vectors with some common string. And, i want to combine into a > vector in such a way that it includes string from both vectors and make a > unique. > > For example: > > x <- paste(rep("A",5),1:5,sep = ".") > x > [1] "A.1" "A.2" "A.3" "A.4" "A.5" > > y <- paste(rep("A"),3:7, sep = ".") > y > [1] "A.3" "A.4" "A.5" "A.6" "A.7" > > Now, I want to combine these two vectors in the following way: > > "A.1" "A.2" "A.3" "A.4" "A.5" "A.6" "A.7" > > I tried with merge(), but not able to get as I want. > > Is there any way to do this in R ? > > If it is, it will be great. > Thanks in advance. > > Regards, > Ram Kumar Basnet. > Wageningen, Netherland. > > > > [[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] what is this averaging function called ?, has R a built in function for it ?
Hi there, I didn't fully follow your code, but allow me to ask: *What is your input - and what do you want to output to be?* Is you input a vector (array) of points? Is your output some sort of averaging of the points according to order? If you searching for a function to "smooth" a sequance of points? Have a look at this: http://en.wikipedia.org/wiki/Lowess In R you can use ?loess Cheers, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Sun, Dec 5, 2010 at 12:11 AM, madr wrote: > > I know little of statistics and have created this function out of > intuition. > But since this algorithm is so basic I wonder what is the proper name of > this function and is it build in R. > > here is some code in PHP to illustrate what the function is doing, it uses > some function I created but the meaning is obvious: > > #get csv file and interchange rows with columns to get two arrays > $csv = aic(getcsv(file_get_contents("out.csv"))); > #now those arrays contained in one bigger array are sorted > array_multisort($csv[0],SORT_NUMERIC,$csv[1],SORT_NUMERIC); > > #second array is created and values that will be put on x or 0 axis are > made > unique with every y or 1 > # value is going into array under x/0 it will be used after to make mean > arithmetic, geometric or harmonic > foreach ($csv[0] as $k=>$x) { >$sum[$x][] = $csv[1][$k]; > } > > #the x values are put on other array for later use > $x = array_keys($sum); > $rang = $sum = array_values($sum); > > #and here is the key feature, to smooth the line the function looks for (in > this case) 500 values above and beond given value > # if they exist of course, the search stops when search goes outside the > array > # the search also stop when number of gathered values goes beyond 500 or > next value that would be added will be making > # this value more than 500, you can imagine that there could be a large > spike in data and this would be affecting points near > # if this precaution haven't been conceived > foreach ($rang as $k=>&$v) { >if (!($k % 100)) echo $k.' '; >$up = $down = array(); >$walk = 0; >while (true) { >++$walk; >if (isset($sum[$k-$walk]) and > count($v)+count($up)+count($sum[$k-$walk])<500) >$up = array_merge($up,$sum[$k-$walk]); >else break; >} >$walk = 0; >while (true) { >++$walk; >if (isset($sum[$k+$walk]) and > count($v)+count($down)+count($sum[$k+$walk])<500) >$down = array_merge($down,$sum[$k+$walk]); >else break; >} > >$rang[$k] = array_merge($up,$rang[$k],$down); ># after gathering data for given point it makes a mean, in this case > arithmetic >$rang[$k] = array_sum($rang[$k])/count($rang[$k]); > } > # now the array with x values can be added and fipped array is ready to go > to a file > $csv = aic(array($x,$rang)); > > # in php this is awfully slow but I like it because it is sensitive for the > densiti of the data and to not goes away in strange > # directions when data density becomes very low > -- > View this message in context: > http://r.789695.n4.nabble.com/what-is-this-averaging-function-called-has-R-a-built-in-function-for-it-tp3072826p3072826.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] merging two vectors
On 05/12/2010 9:57 AM, ram basnet wrote: Dear R users, It may be very simple but it is being difficult for me. I have two vectors with some common string. And, i want to combine into a vector in such a way that it includes string from both vectors and make a unique. For example: x<- paste(rep("A",5),1:5,sep = ".") x [1] "A.1" "A.2" "A.3" "A.4" "A.5" y<- paste(rep("A"),3:7, sep = ".") y [1] "A.3" "A.4" "A.5" "A.6" "A.7" Now, I want to combine these two vectors in the following way: "A.1" "A.2" "A.3" "A.4" "A.5" "A.6" "A.7" I tried with merge(), but not able to get as I want. Is there any way to do this in R ? unique(c(x,y)) 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] merging two vectors
Dear R users, It may be very simple but it is being difficult for me. I have two vectors with some common string. And, i want to combine into a vector in such a way that it includes string from both vectors and make a unique. For example: x <- paste(rep("A",5),1:5,sep = ".") x [1] "A.1" "A.2" "A.3" "A.4" "A.5" y <- paste(rep("A"),3:7, sep = ".") y [1] "A.3" "A.4" "A.5" "A.6" "A.7" Now, I want to combine these two vectors in the following way: "A.1" "A.2" "A.3" "A.4" "A.5" "A.6" "A.7" I tried with merge(), but not able to get as I want. Is there any way to do this in R ? If it is, it will be great. Thanks in advance. Regards, Ram Kumar Basnet. Wageningen, Netherland. [[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] foreach question
Hello Group, I am experimenting with parallel processing ... here is part of my code. require(snow) require(doSNOW) require(foreach) cl.tmp = makeCluster(rep("localhost",4), type="SOCK") registerDoSNOW(cl.tmp) foreach(i=1:NROW(sDat),.packages="gdata",.verbose=TRUE) %dopar% { do something 1 do something 2 do something 3 } I ran this individually without a parallel setup and it works perfectly Now when I do the parallel invocation the error being thrown is that the functions which I created are not being passed to the parallel instances. do something 1,2,3 above use my functions which I created in the code itself. Is there someway that I can have then pass to the parallel routines? Is this the right approach? Thanks and have a nice weekend. S __ 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] The behaviour of read.csv().
On 03/12/2010 7:08 AM, Duncan Murdoch wrote: On 02/12/2010 9:59 PM, Rolf Turner wrote: On 3/12/2010, at 3:48 PM, David Scott wrote: On 03/12/10 14:33, Duncan Murdoch wrote: I think the fill=TRUE option arrived about 10 years ago, in R 1.2.0. The comment in the NEWS file suggests it was in response to some strange csv file coming out of Excel. The real problem with the CSV format is that there really isn't a well defined standard for it. The first RFC about it was published in 2005, and it doesn't claim to be authoritative. Excel is kind of a standard, but it does some very weird things. (For example: enter the string 01 into a field. To keep the leading 0, you need to type it as '01. Save the file, read it back: goodbye 0. At least that's what a website I was just on says about Excel, and what OpenOffice does.) I've been burned so many times by storing data in .csv files, that I just avoid them whenever I can. Absolutely agree with this Duncan. Playing around with .csv files is like playing with some sort of unstable explosive. I also avoid them as much as possible. Where I work, everybody but me uses (yeuuccchhh!!!) Excel or SPSS. If we are to share data sets, *.csv files seem to be the most efficacious, if not the only, way to go. I was going to suggest using DIF rather than CSV. It contains more internal information about the file (including the type of each entry), but has the disadvantage of being less readable, even though it is ascii. However, in putting together a little demo, I found a couple of bugs in the R implementation of read.DIF, and it looks as though it ignores the internal type information. Sigh. As of r53778, the bugs I noticed should be fixed. read.DIF now respects the internal type information, so it will keep character strings like "001" as type character (unless you ask it to change the type). Duncan Murdoch Duncan Murdoch So far, we've had very few problems. The one that started off this thread is the only one I can think of that related to the *.csv format. At least *.csv files have the virtue of being ASCII files, whence if things go wrong it is at least possible to dig into them with a text editor and figure out just what the problem is. cheers, Rolf __ 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 time varying covariate-unfold function
Hello All, I am trying to use the unfold function in RcmdrPlugin.survival library, which converts the survival data with time varying covariates to the counting process notation. The problem is somehow, the event indicator created is not correct. Below is the data, I am trying to convert: CASE TRT FAILTIME FAILCENS SEX AGE IGG0 IGG28 IGG42 IGG84 IGG364 26003 A 11.203321 43.953508080 320 NA IGG0 to IGG364 are my time varying covariates and FAILCENS = 2 indicates observed times. So, when I tried to use, igg2.long<-unfold(igg.2[3,],time="FAILTIME",event="FAILCENS",cov=9:13,cov.na mes="iggtitre"), I obtained the following output: start stop FAILCENS.time CASE TRT FAILTIME FAILCENS SURVTIME SURVCENS SEX 3.1 01 0 26003 A 11.20331 17.54411 1 3.2 12 0 26003 A 11.20331 17.54411 1 3.3 23 0 26003 A 11.20331 17.54411 1 3.4 34 0 26003 A 11.20331 17.54411 1 AGE IGGrel iggtitre 3.1 43.9535 NA0 3.2 43.9535 NA 80 3.3 43.9535 NA 80 3.4 43.9535 NA 320 > As it can be seen, this data was observed, but my censoring indicator (FAILCENS.time) for the last start stop line should have been one, which is not. Could somebody help me with this problem? I did try using the Rossi data also, which is an example in the unfold function. Even there, the arrest.time are all zero, even when the event happened. Thanks a lot for your time. Mathangi [[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] Automatic reprinting in base graphics when resizing device - how?
Dear R-users, I want to write a class that reprints automatically when the device is resized. Using grid graphics this can be achieved e.g. by setting the drawDetails method for a grob which is called automatically when the device has changed. Grid graphics are too slow for my application and I need to do the same in base graphics. How can I automatically update (reprint) the device when it changes in size? Thanks, Mark Mark Heckmann Blog: www.markheckmann.de R-Blog: http://ryouready.wordpress.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help with diff(sqrt()) function in terms of time series
Addi Wei wrote: > > year10 is the time series data set below > 11.64 > 11.50 > 11.49 > > ...I tried my best to produce actual code. > Please post str(year10). Dieter -- View this message in context: http://r.789695.n4.nabble.com/Help-with-diff-sqrt-function-in-terms-of-time-series-tp3072877p3073097.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] grep for strings
Hi Santosh, you may also try require(stringr) x <- "abcdefghijklcd" str_locate_all(x,"cd") On Sun, Dec 5, 2010 at 1:21 PM, Petr Savicky wrote: > On Sun, Dec 05, 2010 at 08:04:08AM +0530, Santosh Srinivas wrote: > > I am trying to find the function where I can search for a pattern in a > > text string (I thought I could use grep for this but no :(). > > > > > x > > [1] "abcdefghijkl" > > > > I want to find the positions (i.e. equivalent of nchar) for "cd" and > > in case there are multiple hits .. then the results as a array > [[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.