Re: [R] Functional data analysis - problems with smoothing
The error message is pretty clear: regardless of what *you* think, R says that 'isi' is not numeric. Are you sure that 'isi' is not a *factor* object? I'm willing to bet that it is. Use str() to check your data. -Peter Ehlers Benjamin Cheah wrote: Hi all, I'm having major issues with smoothing my functional data I'm referring to Jim Ramsay's examples in his books. The following error message keeps appearing, despite all my data being numeric can anyone kindly offer any suggestions? isi - vector of argument values - i.e. the independent variable of the curves rlz - data array TMSfdPar - functional data parameter. TMSfdPar = fdPar(TMSbasis, 4, 0.01) TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar) Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric. I don't understand why the error message keeps popping up. I've tried playing around with Jim Ramsay's datasets and I think my data is organised in a similar manner to his, so can't understand what's going on oh, the frustration! Thanks in advance, Ben (the amateur R data analyst and statistician.) __ [[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] coxph() and survfit()
Dear All, I have a question regarding the output of survfit() when I supply a Cox model. Lets say for example: library(survival) fit - coxph(Surv(time, status == 2) ~ factor(spiders), data = pbc) fit # HR for spiders is significant newdata - data.frame(spiders = factor(0:1)) sf - survfit(fit, newdata = newdata) sum.sf - summary(sfit, times = c(2000, 2500, 3000)) # survival estimates for the yes/no spiders # and the 3 follow up times sum.sf$surv # corresponding lower limits of the # 95% CI sum.sf$low # corresponding upper limits of the # 95% CI sum.sf$up we observe that the 95% CIs overlap!! How is this possible since the HR for spiders is significant. Regards, Mura __ 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] cut2 once, bin twice...
Hello, I'm using the Hmisc cut2 function to bin a set of data. It produces bins that I like with results like this: [96,270]:171 [69, 96): 54 [49, 69): 40 [35, 49): 28 [28, 35): 14 [24, 28): 8 (Other) : 48 I would like to take a second set of data, and assign it to bins based on factors defined by my call to cut 2. Does anyone know how I can do this? Thank you, -S -- View this message in context: http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26020736.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] cut2 once, bin twice...
sdanzige wrote: I'm using the Hmisc cut2 function to bin a set of data. It produces bins that I like with results like this: [96,270]:171 [69, 96): 54 [49, 69): 40 [35, 49): 28 [28, 35): 14 [24, 28): 8 (Other) : 48 I would like to take a second set of data, and assign it to bins based on factors defined by my call to cut 2. It used to be quite tricky, but on popular request Brian Ripley has added an example how to extract the intervals using regular expression on the bottom of the examples for cut (note:cut in base, not cut2 in Hmisc). If someone knows of an easier way, please correct me. How about adding this information as attribute to the standard cut? Dieter -- View this message in context: http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26022244.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] twoord.plot y lab size
On 10/22/2009 10:26 PM, Jacob Kasper wrote: I am using twoord.plot and my Y axis units on the left y overlap. I tried using cex.axis in my par command but that only adjusted the x units, not the y. cex.axis in twoord.plot did not help. How do I adjust Y units in the twoord.plot? example: twoord.plot (lx=myear, ly=z, rx=myear, ry=sta,xlim=c(1985,2010), xlab=Year,ylab=Individuals, rylab=# Stations) Hi Jacob, While there is no way to do this in the current version, the upcoming version 2.7-2 will have an axislab.cex argument that allows control of the axis and tick labels. Thanks for the idea. Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Bayesian regression stepwise function?
Charles C. Berry wrote: On Thu, 22 Oct 2009, Ben Bolker wrote: Allan.Y wrote: Hi everyone, I am wondering if there exists a stepwise regression function for the Bayesian regression model. I tried googling, but I couldn't find anything. I know step function exists for regular stepwise regression, but nothing for Bayes. Why? That seems so ... un-Bayesian ... If 'fools rush in where angels fear to tread', then Bayesians 'jump' in where frequentists fear to 'step'... Seriously, there are Bayesian regression approaches that priorize the model size (sometimes only implicitly by assigning a prior for the inclusion of each candidate regressors). Then they 'jump' between models of different sizes. On CRAN, Package qtlbim (which is specialized to a particular genetics problem) implements one such, I think. Package bqtl does not implement the jumping approach, but does explore a model space with differing numbers of regressors for the same (qtl) problem. Perhaps the closest to a general purpose 'stepwise flavored' Bayesian regression is implemented in Package BMA, which IIRC borrows step() for some of its work. But CRAN now has more packages than my cortex has neurons, so there are probably more packages that do something like this. Try RSiteSearch(jump regression, restric='functions') and start reading. library(sos) jr - ???'jump regression' summary(jr) # Downloaded 22 links in 14 packages. # packages with 3 help pages matching 'jump regression': # monomvn, CalciOMatic, polydect; # qtl has 1 ... bs - ???'Bayesian step' bsw - ???'Bayesian stepwise' bs. - bs|bsw summary(bs.) Total number of matches: 119 Downloaded 115 links in 63 packages. Packages with at least 2 matches using search pattern 'Bayesian+step | Bayesian+stepwise': Package Count MaxScore TotalScoreDate 1 DPpackage 82 11 2009-03-17 06:38:42 2 BMA 61 6 2009-05-01 12:21:01 3 ensembleBMA 51 5 2009-07-01 07:16:43 4 MMIX 43 8 2009-08-16 05:47:37 5mclust 42 8 2009-07-14 05:39:10 6adlift 41 4 2009-02-01 06:31:16 7 tgp 34 8 2009-07-29 11:14:36 8 geoRglm 34 6 2009-10-19 19:21:57 9 G1DBN 32 5 2008-01-24 16:27:29 10 bayesm 31 3 2008-06-14 15:31:09 11 bqtl 31 3 2009-01-15 07:02:43 12 qtlbim 2 12 13 2009-07-29 10:50:13 13 BAYSTAR 2 11 12 2009-10-19 19:17:18 14 gap 24 5 2009-01-15 07:03:58 15 MSBVAR 24 5 2009-07-29 10:49:07 16bayesSurv 23 5 2008-12-15 10:06:16 17 bnlearn 22 3 2009-09-01 08:25:43 18 spBayes 22 3 2009-03-28 07:33:41 19 dlm 21 2 2009-07-01 07:16:28 20GOFSN 21 2 2009-07-01 07:47:10 21 MCMChybridGP 21 2 2009-10-01 05:27:43 22 mgcv 21 2 2009-08-24 16:03:12 23nlreg 21 2 2007-11-02 13:14:56 24 nlt 21 2 2009-03-02 07:30:34 25 pARtial 21 2 2006-11-01 02:59:23 26plink 21 2 2009-07-01 07:19:44 27 qtl 21 2 2009-09-17 15:46:05 28 timsac 21 2 2009-03-28 07:27:28 # This took more longer to describe it than run. # Hope this helps. # Spencer HTH, Chuck -- View this message in context: http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p26015081.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. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:cbe...@tajo.ucsd.eduUC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ 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. -- Spencer Graves, PE, PhD President and Chief Operating Officer Structure Inspection and Monitoring, Inc. 751 Emerson Ct. San José, CA 95126 ph: 408-655-4567 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help
Re: [R] Functional data analysis - problems with smoothing
Hi, Ben: Which chapter in which of Ramsay's books? For his most recent book (Functional Data Analysis with R and Matlab, with Giles Hooker and your's truly), there is one script file for each chapter with names like fdarm-ch01.R, ... fdarm-ch11.R. These script files reproduce nearly all the examples in the book including all but one of the 76 figures (according to what I wrote for back cover). The preface tells you how to find those script files: system.file('scripts', package='fda') [1] c:/Users/sgraves/R/R-2.9.1/library/fda/scripts Obviously, the answer is the path on my installation; the answer you get should tell you where to look on the computer (and start-up sequence) you are using. This same scripts directory contains similar partial scripts for Functional Data Analysis and Applied Functional Data Analysis. If you'd care to send me updates to any of those partial scripts, I can see that they are included in the package to make those examples easier for others in the future. Hope this helps. Spencer Graves Peter Ehlers wrote: The error message is pretty clear: regardless of what *you* think, R says that 'isi' is not numeric. Are you sure that 'isi' is not a *factor* object? I'm willing to bet that it is. Use str() to check your data. -Peter Ehlers Benjamin Cheah wrote: Hi all, I'm having major issues with smoothing my functional data I'm referring to Jim Ramsay's examples in his books. The following error message keeps appearing, despite all my data being numeric can anyone kindly offer any suggestions? isi - vector of argument values - i.e. the independent variable of the curves rlz - data array TMSfdPar - functional data parameter. TMSfdPar = fdPar(TMSbasis, 4, 0.01) TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar) Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric. I don't understand why the error message keeps popping up. I've tried playing around with Jim Ramsay's datasets and I think my data is organised in a similar manner to his, so can't understand what's going on oh, the frustration! Thanks in advance, Ben (the amateur R data analyst and statistician.) __ [[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. -- Spencer Graves, PE, PhD President and Chief Operating Officer Structure Inspection and Monitoring, Inc. 751 Emerson Ct. San José, CA 95126 ph: 408-655-4567 __ 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] PDF too large, PNG bad quality
On 10/23/2009 06:07 AM, Lasse Kliemann wrote: I wish to save a scatter plot comprising approx. 2 million points in order to include it in a LaTeX document. Using 'pdf(...)' produces a file of size about 20 MB, which is useless. Using 'cairo_pdf(...)' produces a smaller file, around 3 MB. This is still too large. Not only that the document will be too large, but also PDF viewers choke on this. Moreover, Cairo has problems with text: by default text looks ugly, like scaled bitmaps. After hours of trying different settings, I discovered that choosing a different font family can help, e.g.: 'par(family=Mono)'. This gives good-looking text. Yet, the problem with the file size remains. There exists the hint to produdc EPS instead and then convert to PDF using 'epstopdf'. The resulting PDF files are slightly smaller, but still too large, and PDF viewers still don't like it. So I gave PNG a try. PNG files are much smaller and PDF viewers have no trouble with them. However, fonts look ugly. The same trick that worked for Cairo PDF has no effect for PNG. When I view the PNGs with a dedicated viewer like 'qiv', even the fonts look good. But not when included in LaTeX; I simply use '\includegraphics{...}' and run the document through 'pdflatex'. I tried both, creating PNG with 'png(...)' and converting from PDF to PNG using 'convert' from ImageMagick. So my questions are: - Is there a way to produce sufficiently lean PDFs directly in R, even when the plot comprises several million points? - How to produce a PNG that still looks nice when included in a LaTeX PDF document? Any hints will be greatly appreciated. Hi Lasse, I may be right off the track, but I can't make much sense of 2 million points in a scatterplot. If you are interested in the density of points within the plot, you could compute this using something like the bkde2 function in the KernSmooth package and then plot that using something like image. 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] data file with columns of unequal length
I am running an expt that presents a point process input x and measures a point process output y. The times of each event are recorded. The lengths of the data records of x and y are necessarily different, and can be different by a factor of 10. I would like to save these data after each experiment as a file with two columns, one for x and one for y. However, R dataframes require columns of equal length. One solution is to fill the empty places in y with NAs so it has the same length as x. I view that as unsatisfactory (there are in reality no missing values). Another possibility is to store x and y in separate files. I also view that as unsatisfactory (it is too easy to lose track of the y file corresponding to a given x file). Can anyone suggest a way to deal with this situation? Thanks very much 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] data file with columns of unequal length
On 10/23/2009 07:58 PM, William Simpson wrote: I am running an expt that presents a point process input x and measures a point process output y. The times of each event are recorded. The lengths of the data records of x and y are necessarily different, and can be different by a factor of 10. I would like to save these data after each experiment as a file with two columns, one for x and one for y. However, R dataframes require columns of equal length. One solution is to fill the empty places in y with NAs so it has the same length as x. I view that as unsatisfactory (there are in reality no missing values). Another possibility is to store x and y in separate files. I also view that as unsatisfactory (it is too easy to lose track of the y file corresponding to a given x file). Can anyone suggest a way to deal with this situation? Hi Bill, xy-list(x=1:10,y=1:100) Note that this cheerfully ignores how you are going to figure out which x goes with which y(s). 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] standard errors in bbmle
Hello, Mle2 is a little unforthcoming in the matter of standard errors? Is there a way to ask the program to supply standard errors along with estimates in cases when it doesn't print them 'voluntarily'? regards, s [[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] data file with columns of unequal length
Jim Lemon wrote: On 10/23/2009 07:58 PM, William Simpson wrote: I am running an expt that presents a point process input x and measures a point process output y. The times of each event are recorded. The lengths of the data records of x and y are necessarily different, and can be different by a factor of 10. I would like to save these data after each experiment as a file with two columns, one for x and one for y. However, R dataframes require columns of equal length. One solution is to fill the empty places in y with NAs so it has the same length as x. I view that as unsatisfactory (there are in reality no missing values). Another possibility is to store x and y in separate files. I also view that as unsatisfactory (it is too easy to lose track of the y file corresponding to a given x file). Can anyone suggest a way to deal with this situation? Hi Bill, xy-list(x=1:10,y=1:100) Note that this cheerfully ignores how you are going to figure out which x goes with which y(s). As I understand it, they don't come in pairs anyway. For the same reason a data frame is just the wrong kind of data structure. If you don't want separate data files, you can use one file with two columns where the second column is (say) 1 for the x and 2 for the y. (There are other options, like concatenating the x and the y with some sort of separator inbetween, but it easily gets painful to read them back in.) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - (p.dalga...@biostat.ku.dk) FAX: (+45) 35327907 __ 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] standard errors in bbmle
alexander russell wrote: Hello, Mle2 is a little unforthcoming in the matter of standard errors? Is there a way to ask the program to supply standard errors along with estimates in cases when it doesn't print them 'voluntarily'? regards, s What did you do? Which cases? Did you look at the examples? (as in, example(mle2) ) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - (p.dalga...@biostat.ku.dk) FAX: (+45) 35327907 __ 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] data file with columns of unequal length
The way you do it is to compute the cross-intensity function (you can google this; a key name is David Brillinger). The general problem is that of system identification for point processes. Bill On Fri, Oct 23, 2009 at 10:31 AM, Jim Lemon j...@bitwrit.com.au wrote: On 10/23/2009 07:58 PM, William Simpson wrote: I am running an expt that presents a point process input x and measures a point process output y. The times of each event are recorded. The lengths of the data records of x and y are necessarily different, and can be different by a factor of 10. I would like to save these data after each experiment as a file with two columns, one for x and one for y. However, R dataframes require columns of equal length. One solution is to fill the empty places in y with NAs so it has the same length as x. I view that as unsatisfactory (there are in reality no missing values). Another possibility is to store x and y in separate files. I also view that as unsatisfactory (it is too easy to lose track of the y file corresponding to a given x file). Can anyone suggest a way to deal with this situation? Hi Bill, xy-list(x=1:10,y=1:100) Note that this cheerfully ignores how you are going to figure out which x goes with which y(s). Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] data file with columns of unequal length
As I understand it, they don't come in pairs anyway. Correct. For the same reason a data frame is just the wrong kind of data structure. If you don't want separate data files, you can use one file with two columns where the second column is (say) 1 for the x and 2 for the y. Could you explain this further? I don't really get it. My xs and ys have different lengths. How would I read them in? Ideally I would just read in the x column and y columns separately x-read.file(file.dat, column1) y-read.file(file.dat, column2) But I know of no way to do that... [Oh, by the way, I never mentioned this, the x and ys are event times and are in ascending order (by time of occurrence).] (There are other options, like concatenating the x and the y with some sort of separator inbetween, but it easily gets painful to read them back in.) Thanks 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] data file with columns of unequal length
Thanks Jim. BTW the times in x and y are in ascending order (time of occurrence). If I do it this way, how do I actually read the data in and store in the file? Toy code, please. Bill Hi Bill, xy-list(x=1:10,y=1:100) Note that this cheerfully ignores how you are going to figure out which x goes with which y(s). Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] data file with columns of unequal length
William Simpson wrote: As I understand it, they don't come in pairs anyway. Correct. For the same reason a data frame is just the wrong kind of data structure. If you don't want separate data files, you can use one file with two columns where the second column is (say) 1 for the x and 2 for the y. Could you explain this further? I don't really get it. My xs and ys have different lengths. How would I read them in? Same format as when you have data from two different groups. Don't really know how to make it clearer than that. Like the sleep dataset (OK, so they are actually paired...). Ideally I would just read in the x column and y columns separately x-read.file(file.dat, column1) y-read.file(file.dat, column2) But I know of no way to do that... [Oh, by the way, I never mentioned this, the x and ys are event times and are in ascending order (by time of occurrence).] (There are other options, like concatenating the x and the y with some sort of separator inbetween, but it easily gets painful to read them back in.) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - (p.dalga...@biostat.ku.dk) FAX: (+45) 35327907 __ 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] data file with columns of unequal length
OK thanks, I look at sleep and get it Bill On Fri, Oct 23, 2009 at 12:21 PM, Peter Dalgaard p.dalga...@biostat.ku.dk wrote: William Simpson wrote: As I understand it, they don't come in pairs anyway. Correct. For the same reason a data frame is just the wrong kind of data structure. If you don't want separate data files, you can use one file with two columns where the second column is (say) 1 for the x and 2 for the y. Could you explain this further? I don't really get it. My xs and ys have different lengths. How would I read them in? Same format as when you have data from two different groups. Don't really know how to make it clearer than that. Like the sleep dataset (OK, so they are actually paired...). Ideally I would just read in the x column and y columns separately x-read.file(file.dat, column1) y-read.file(file.dat, column2) But I know of no way to do that... [Oh, by the way, I never mentioned this, the x and ys are event times and are in ascending order (by time of occurrence).] (There are other options, like concatenating the x and the y with some sort of separator inbetween, but it easily gets painful to read them back in.) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - (p.dalga...@biostat.ku.dk) FAX: (+45) 35327907 __ 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] data file with columns of unequal length
On 10/23/2009 10:07 PM, William Simpson wrote: Thanks Jim. BTW the times in x and y are in ascending order (time of occurrence). If I do it this way, how do I actually read the data in and store in the file? Toy code, please. Hi Bill, This seems a bit like some heartbeat data that I had to deal with some years ago. There were two output streams, one the time of each R wave and the other an asynchronous stimulus. In that case, I had to work out code to interdigitate the signals and do some processing of the resulting data. It sounds like you have two files, each recording times of input and output respectively one value to a line. If so, my first guess is: x-as.vector(read.table(input_times.dat)) y-as.vector(read.table(output_times.dat)) xy-list(x,y) # to store this list in a file write.csv(xy,file=xy_23_10_2009.csv) Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] data file with columns of unequal length
On 10/23/2009 10:07 PM, William Simpson wrote: Thanks Jim. BTW the times in x and y are in ascending order (time of occurrence). If I do it this way, how do I actually read the data in and store in the file? Toy code, please. Hi Bill, This seems a bit like some heartbeat data that I had to deal with some years ago. There were two output streams, one the time of each R wave and the other an asynchronous stimulus. In that case, I had to work out code to interdigitate the signals and do some processing of the resulting data. It sounds like you have two files, each recording times of input and output respectively one value to a line. If so, my first guess is: x-as.vector(read.table(input_times.dat)) y-as.vector(read.table(output_times.dat)) xy-list(x,y) # to store this list in a file # oops, bit too quick on the keyboard dput(xy,file=xy_23_10_2009.R) Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] arima crashes too
Barry Rowlingson wrote: However, arima crashes for this: arima(c(1.71, 1.78, 1.95, 1.59, 2.13), order=c(1,0,0)) I'm not getting what I'd call 'crashes' with your arma or arima examples- I get an error message and a warning: arma(c(2.01, 2.22, 2.09, 2.17, 2.42), order=c(1,0)) Error in AA %*% t(X) : requires numeric/complex matrix/vector arguments In addition: Warning message: In ar.ols(x, order.max = k, aic = FALSE, demean = FALSE, intercept = include.intercept) : model order: 2singularities in the computation of the projection matrixresults are only valid up to model order1 You've not told us what you get, Because I get a message in Portuguese. But the meaning is precisely that one: Error in AA %*% t(X). and the phrase 'crash' normally means some kind of memory error that *terminates* a running R session. Ah, maybe then crash is not the correct word. It stops running whatever it is running. Maybe abort is the best word? Are you really crashing R such that it terminates? In which case, what version number/platform etc? I mean that, if I run a loop, it doesn't finish. Or, more catastrophically, if I am running a loop and saving data to an open file, it terminates the loop and does not close the file. Reproducible example: test.arima - function() { lets.crash.arima - c(71, 78, 95, 59) # , 113 for (x in 90:120) { reg - arima(c(lets.crash.arima, x), order = c(1,0,0)) cat(ok for x =, x, \n) } cat(close file and prepare a nice summary\n) return(arima passed the test) } test.arima() As you can see, the loop aborts, the function never returns, with potentially nasty effects (namely: I have to finish R with q() to close the files and examine them). Alberto 'Darth Albmont' Monteiro __ 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] opposite estimates from zeroinfl() and hurdle()
Dear all, A question related to the following has been asked on R-help before, but I could not find any answer to it. Input will be much appreciated. I got an unexpected sign of the slope parameter associated with a covariate (diam) using zeroinfl(). It led me to compare the estimates given by zeroinfl() and hurdle(): The (significant) negative estimate here is surprising, given the biology of the species: summary(zeroinfl(bnl ~ 1| diam, dist = poisson, data = valdaekar, EM = TRUE)) Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(|z|) (Intercept) 3.74604 0.02635 142.2 2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(|z|) (Intercept) 21.7510 7.6525 2.842 0.00448 ** diam -1.1437 0.3941 -2.902 0.00371 ** Number of iterations in BFGS optimization: 1 Log-likelihood: -582.8 on 3 Df The hurdle model gives the same estimates, but with opposite (and expected) signs of the parameters: summary(hurdle(bnl ~ 1| diam, dist = poisson, data = valdaekar)) Count model coefficients (truncated poisson with log link): Estimate Std. Error z value Pr(|z|) (Intercept) 3.74604 0.02635 142.2 2e-16 *** Zero hurdle model coefficients (binomial with logit link): Estimate Std. Error z value Pr(|z|) (Intercept) -21.7510 7.6525 -2.842 0.00448 ** diam 1.1437 0.3941 2.902 0.00371 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 8 Log-likelihood: -582.8 on 3 Df Why is this so? thanks, Tord Windows NT, R 2.8.1, pcsl 1.03 -- Tord Snäll Department of Ecology / Swedish Species Information Centre Swedish University of Agricultural Sciences (SLU) P.O. 7044, SE-750 07 Uppsala, Sweden Office/Mobile/Fax +46-18-672612/+46-76-7662612/+46-18-673537 www.ekol.slu.se/staff_tordsnall www.artdata.slu.se/personal/fototsn.asp __ 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] a problem about integrate function in R .thank you .
I don't seem to get a problem with this. Have you tried a Monte Carlo approach to verify that you are getting incorrect answers? For me, I get when the upper is 1 that integrate(e2, lower = 0, upper = 1) -0.2820948 with absolute error 5e-05 sum(e2(runif(1)))/1 [1] -0.2825667 which seems to be consistent, while for upper = 0.75 I get 0.75*sum(e2(runif(1, min=0, max=0.75)))/1 [1] -0.2333506 integrate(e2,lower=0,upper = 0.75) -0.2341178 with absolute error 7.8e-05 On Oct 23, 2:52 pm, fuzuo xie xiefu...@gmail.com wrote: e2 - function(x) { out - 0*x for(i in 1:length(x)) out[i] -integrate(function(y) qnorm(y),lower=0,upper=x[i])$value out } integrate(e2,lower=0, upper=a)$value above is my code , when a is small , say a0.45 the result is right . however , when a0.5 the result is incorrect . why ? thank you . [[alternative HTML version deleted]] __ r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guidehttp://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] arima crashes too
On Fri, Oct 23, 2009 at 12:32 PM, Alberto Monteiro albm...@centroin.com.br wrote: I mean that, if I run a loop, it doesn't finish. Or, more catastrophically, if I am running a loop and saving data to an open file, it terminates the loop and does not close the file. Reproducible example: test.arima - function() { lets.crash.arima - c(71, 78, 95, 59) # , 113 for (x in 90:120) { reg - arima(c(lets.crash.arima, x), order = c(1,0,0)) cat(ok for x =, x, \n) } cat(close file and prepare a nice summary\n) return(arima passed the test) } test.arima() As you can see, the loop aborts, the function never returns, with potentially nasty effects (namely: I have to finish R with q() to close the files and examine them). If you're doing anything in a loop that has the potential to fail because of singularities or other conditions when your model can't be fitted, you need to stick what you are doing in a 'try' clause. This lets you trap errors and do something with them. Plenty of examples in help(try) or this from me: for(i in 1:10){ print(solve(matrix(c(3,3,3,i),2,2))) } This stops the loop at i=3. Now stick it in a try() clause: for(i in 1:10){ print(try(solve(matrix(c(3,3,3,i),2,2 } and it gives a warning and carries on. If you want your code to do something with the failure cases then the help for try() tells you what to look for. I'm not sure why your arima produces an error, but I'm assuming the numbers are such that the model can't be fitted. I don't really know what arima is doing. Barry __ 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] reset par() within plot layout
Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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] reset par() within plot layout
Maybe something like this? #Desired result is a layout of 2 plots: one red and one black par(mfrow=c(2,1)) par(col=red) plot(1:100) par(col=black) plot(1:10) On Fri, Oct 23, 2009 at 7:26 AM, Janke ten Holt j.c.ten.h...@rug.nl wrote: Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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. -- Stephen Sefick Let's not spend our time and resources thinking about things that are so little or so large that all they really do for us is puff us up and make us feel like gods. We are mammals, and have not exhausted the annoying little problems of being mammals. -K. Mullis __ 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] Functional data analysis - problems with smoothing
Hi Peter and Spencer (and everyone else out there), Thank you for your prompt reply to my post re. FDA - FYI - i'm very new to R and FDA. Currently half way through a biostats masters degree in Sydney, Australia! Absolutely loving it and I hope to include R computing on my future CV! Just to clarify things: When I entered what you suggested, this is what I obtained with my own dataset str(isi) num [1:14] 1 1.5 2 2.5 3 3.5 4 5 7 10 ... Similarly, for Jim Ramsay's dataset (chapter five of Use R book): str(age) num [1:31] 1 1.25 1.5 1.75 2 3 4 5 6 7 ... Doesn't this mean that my data object, 'isi' is numeric? Also, I was looking through Jim Ramsay's datasets - I'm fairly sure that my vector, 'isi' was similarly organised to his vector, 'age' in chapter 5 of his new Use R FDA book, so I'm fairly certain that the data is numeric but it obviously isn't Hope this email clarified things - apologies if it hasn't - still very new to the terminology and overall feel of R. The R user network is terrific! Thank you for caring! Ben From: Peter Ehlers ehl...@ucalgary.ca Cc: R-help@r-project.org Sent: Fri, 23 October, 2009 5:07:15 PM Subject: Re: [R] Functional data analysis - problems with smoothing The error message is pretty clear: regardless of what *you* think, R says that 'isi' is not numeric. Are you sure that 'isi' is not a *factor* object? I'm willing to bet that it is. Use str() to check your data. -Peter Ehlers Benjamin Cheah wrote: Hi all, I'm having major issues with smoothing my functional data I'm referring to Jim Ramsay's examples in his books. The following error message keeps appearing, despite all my data being numeric can anyone kindly offer any suggestions? isi - vector of argument values - i.e. the independent variable of the curves rlz - data array TMSfdPar - functional data parameter. TMSfdPar = fdPar(TMSbasis, 4, 0.01) TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar) Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric. I don't understand why the error message keeps popping up. I've tried playing around with Jim Ramsay's datasets and I think my data is organised in a similar manner to his, so can't understand what's going on oh, the frustration! Thanks in advance, Ben (the amateur R data analyst and statistician.) __ [[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.
Re: [R] reset par() within plot layout
OK, maybe I oversimplified the example, because you are of course correct... I am working on a function, which I feed data to plot in layouts of specified dimensions. I want to be able to set any par() variables for each plot in the layout. Because the next plot does not know which par() variables were changed in the previous plot, I want to reset all par() variables after each plot, rather than simply change back one specific par() variable. So I really need to reset par() without 'messing with' my layout. stephen sefick wrote: Maybe something like this? #Desired result is a layout of 2 plots: one red and one black par(mfrow=c(2,1)) par(col=red) plot(1:100) par(col=black) plot(1:10) On Fri, Oct 23, 2009 at 7:26 AM, Janke ten Holt j.c.ten.h...@rug.nl wrote: Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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] Wavelets
Hi all, I am trying to do wavelets and I got an error message saying The length of data is not a power of 2 Is there a way of handing that? or should the data length be exactly the power of 2? I am using R version 2.9.2 (2009-08-24) The is library(wavethresh). wds - wd(ds$v,filter.number=1) 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] Wavelets
I believe that you can extend it to the power of 2 by padding it with zeros. I don't remember if it is at the begining or the end of the time series. hth stephen On Fri, Oct 23, 2009 at 7:56 AM, Ashta sewa...@gmail.com wrote: Hi all, I am trying to do wavelets and I got an error message saying The length of data is not a power of 2 Is there a way of handing that? or should the data length be exactly the power of 2? I am using R version 2.9.2 (2009-08-24) The is library(wavethresh). wds - wd(ds$v,filter.number=1) 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. -- Stephen Sefick Let's not spend our time and resources thinking about things that are so little or so large that all they really do for us is puff us up and make us feel like gods. We are mammals, and have not exhausted the annoying little problems of being mammals. -K. Mullis __ 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] Wavelets
hi all, I don't use the package wavethresh, but I suppose that should be possible to fill of ZEROS (zero padding) the time series, such that its length has a 2^n elements. Another way is to cut (remove) some elements, so that, the new length has 2^n elements. You can get more info. @ A practical guide to Wavelet Analysis C. torrence and G. P.Compo Bulleting of the American Meteorology. Society. 1998 I hope this help you, -- Josue Polanco On Fri, Oct 23, 2009 at 2:56 PM, Ashta sewa...@gmail.com wrote: Hi all, I am trying to do wavelets and I got an error message saying The length of data is not a power of 2 Is there a way of handing that? or should the data length be exactly the power of 2? I am using R version 2.9.2 (2009-08-24) The is library(wavethresh). wds - wd(ds$v,filter.number=1) 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. -- Josué Mosés Polanco Martínez Correo-e alternativo jom...@linuxmail.org It is a wasted day unless you have learned something new and made someone smile -Mark Weingartz. __ 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] twoord.plot y lab size
For now, you might be able to use par(las=1). You may have to make enough margin room with par(mar=...) and you'll probably want to set ylab= and add the y-label with mtext(). -Peter Ehlers Jim Lemon wrote: On 10/22/2009 10:26 PM, Jacob Kasper wrote: I am using twoord.plot and my Y axis units on the left y overlap. I tried using cex.axis in my par command but that only adjusted the x units, not the y. cex.axis in twoord.plot did not help. How do I adjust Y units in the twoord.plot? example: twoord.plot (lx=myear, ly=z, rx=myear, ry=sta,xlim=c(1985,2010), xlab=Year,ylab=Individuals, rylab=# Stations) Hi Jacob, While there is no way to do this in the current version, the upcoming version 2.7-2 will have an axislab.cex argument that allows control of the axis and tick labels. Thanks for the idea. 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.
Re: [R] How to find moving averages within each subgroup of a data frame
Dear all, If I have the following data frame: set.seed(21) df1 - data.frame(col1=c(rep('a',5), rep('b',5), rep('c',5)), col4=rnorm(1:15)) col1 col4 1 a 0.793013171 2 a 0.522251264 3 a 1.74641 4 a -1.271336123 5 a 2.197389533 6 b 0.433130777 7 b -1.570199630 8 b -0.934905667 9 b 0.063493345 10b -0.002393336 11c -2.276781240 12c 0.757412225 13c -0.548405554 14c 0.172549478 15c 0.562853068 How can i make a 2 point moving average within each group? i.e. col1col4SMA a 0.793013171 NA a 0.522251264 0.657632218 a 1.74641 1.134236753 a -1.2713361230.237443059 a 2.197389533 0.463026705 b 0.433130777 NA b -1.57019963 -0.568534427 b -0.934905667-1.252552649 b 0.063493345 -0.435706161 b -0.0023933360.030550005 c -2.27678124 NA c 0.757412225 -0.759684508 c -0.5484055540.104503336 c 0.172549478 -0.187928038 c 0.562853068 0.367701273 From what i've read, i was thinking it would be something along the lines of: aggregate(df1$col4, by=list(df1$col1), function(x) {filter(x, rep(1/2,2), sides=1 )} ) Error in aggregate.data.frame(as.data.frame(x), ...) : 'FUN' must always return a scalar But this tells me (i think) that aggregate should only return a single value per group. So what i need, i guess, is something that takes all the values in a given group, and returns a vector of the same length. Not sure which function to use for that. Thanks in advance, C xx P.S. on a side note, is it possible to extract the values which aggregate passes to the function(x) in my example above? __ 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] data frame is killing me! help
Steve Lianoglou-6 wrote: Hi, On Oct 22, 2009, at 2:35 PM, bbslover wrote: Usage data(gasoline) Format A data frame with 60 observations on the following 2 variables. octane a numeric vector. The octane number. NIR a matrix with 401 columns. The NIR spectrum and I see the gasoline data to see below NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR.1694 nm NIR.1696 nm NIR.1698 nm NIR.1700 nm 1 1.242645 1.250789 1.246626 1.250985 1.264189 1.244678 1.245913 1.221135 2 1.189116 1.223242 1.253306 1.282889 1.215065 1.225211 1.227985 1.198851 3 1.198287 1.237383 1.260979 1.276677 1.218871 1.223132 1.230321 1.208742 4 1.201066 1.233299 1.262966 1.272709 1.211068 1.215044 1.232655 1.206696 5 1.259616 1.273713 1.296524 1.299507 1.226448 1.230718 1.232864 1.202926 6 1.24109 1.262138 1.288401 1.291118 1.229769 1.227615 1.22763 1.207576 7 1.245143 1.265648 1.274731 1.292441 1.218317 1.218147 1.73 1.200446 8 1.222581 1.245782 1.26002 1.290305 1.221264 1.220265 1.227947 1.188174 9 1.234969 1.251559 1.272416 1.287405 1.211995 1.213263 1.215883 1.196102 look at this NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR. 1694 nm NIR.1696 nm NIR.1698 nm NIR.1700 nm how can I add letters NIR to my variable, because my 600 independents never have NIR as the prefix. however, it is needed to model the plsr. for example aa=plsr(y~NIR, data=data ,), the prefix NIR is necessary, how can I do with it? I'm not really sue that I'm getting you, but if your problem is that the column names of your data.frame don't match the variable names you'd like to use in your formula, just change the colnames of your data.frame to match your formula. BTW - I have no idea where to get this gasoline data set, so I'm just imagining: eg. colnames(gasoline) - c('put', 'the', 'variable', 'names', 'that', 'you', 'want', 'here') -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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. thanks for you. but the numbers of indenpendence are so many, it is not easy to identify them one by one, is there some better way? -- View this message in context: http://www.nabble.com/data-frame-is-killing-me%21-help-tp26015079p26024985.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] opposite estimates from zeroinfl() and hurdle()
Dear all, A question related to the following has been asked on R-help before, but I could not find any answer to it. Input will be much appreciated. I got an unexpected sign of the slope parameter associated with a covariate (diam) using zeroinfl(). It led me to compare the estimates given by zeroinfl() and hurdle(): The (significant) negative estimate here is surprising, given the biology of the species: summary(zeroinfl(bnl ~ 1| diam, dist = poisson, data = valdaekar, EM = TRUE)) Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(|z|) (Intercept) 3.746040.02635 142.2 2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(|z|) (Intercept) 21.7510 7.6525 2.842 0.00448 ** diam -1.1437 0.3941 -2.902 0.00371 ** Number of iterations in BFGS optimization: 1 Log-likelihood: -582.8 on 3 Df The hurdle model gives the same estimates, but with opposite (and expected) signs of the parameters: summary(hurdle(bnl ~ 1| diam, dist = poisson, data = valdaekar)) Count model coefficients (truncated poisson with log link): Estimate Std. Error z value Pr(|z|) (Intercept) 3.746040.02635 142.2 2e-16 *** Zero hurdle model coefficients (binomial with logit link): Estimate Std. Error z value Pr(|z|) (Intercept) -21.7510 7.6525 -2.842 0.00448 ** diam 1.1437 0.3941 2.902 0.00371 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 8 Log-likelihood: -582.8 on 3 Df Why is this so? thanks, Tord Windows NT, R 2.8.1, pcsl 1.03 __ 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] Wavelets
hi all, On Fri, Oct 23, 2009 at 3:05 PM, stephen sefick ssef...@gmail.com wrote: I believe that you can extend it to the power of 2 by padding it with zeros. I don't remember if it is at the begining or the end of the time series. hth at the end of the time series Cheers, -- Josue Polanco stephen On Fri, Oct 23, 2009 at 7:56 AM, Ashta sewa...@gmail.com wrote: Hi all, I am trying to do wavelets and I got an error message saying The length of data is not a power of 2 Is there a way of handing that? or should the data length be exactly the power of 2? I am using R version 2.9.2 (2009-08-24) The is library(wavethresh). wds - wd(ds$v,filter.number=1) 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. -- Stephen Sefick Let's not spend our time and resources thinking about things that are so little or so large that all they really do for us is puff us up and make us feel like gods. We are mammals, and have not exhausted the annoying little problems of being mammals. -K. Mullis __ 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. -- Josué Mosés Polanco Martínez Correo-e alternativo jom...@linuxmail.org It is a wasted day unless you have learned something new and made someone smile -Mark Weingartz. __ 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] access elements of a named list using a factor
Hi I have a factor 'f' and a named list 'jj'. I want names(jj) to match up with levels(f). How do I use levels(f) to access elements of jj? f - factor(c(pigs,pigs,slugs)) f [1] pigs pigs slugs Levels: pigs slugs jj - list(pigs=1:10,slugs=1:3) My attempts to produce jj$pigs all give errors: jj$levels(f)[1] Error: attempt to apply non-function do.call($,jj,levels(f)[1]) Error in if (quote) { : argument is not interpretable as logical $(jj,levels(f)[1]) Error in jj$levels(f)[1] : invalid subscript type 'language' -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ 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] Memory Problems with CSV and Survey Objects
I'm working with a 350MB CSV file on a server that has 3GB of RAM, yet I'm hitting a memory error when I try to store the data frame into a survey design object, the R object that stores data for complex sample survey data. When I launch R, I execute the following line from Windows: C:\Program Files\R\R-2.9.1\bin\Rgui.exe --max-mem-size=2047M Anything higher, and I get an error message saying the maximum has been set to 2047M. Here are the commands: library(survey) #this step takes more than five minutes data08-read.csv(data08.csv,header=TRUE,nrows=210437) object.size(data08) #329877112 bytes #Looking at Windows Task Manager, Mem Usage for Rgui.exe is already 659,632K brr.dsgn -svrepdesign( data = data08 , repweights = data08[, grep( ^repwgt , colnames( data08)) ], type = BRR , combined.weights = TRUE , weights = data08$mainwgt ) #Error: cannot allocate vector of size 254.5 Mb #The survey design object does not get created. #This also causes Windows Task Manager, Mem Usage to spike to 1,748,136K #And here are some memory diagnostics memory.limit() [1] 2047 memory.size() [1] 1449.06 gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 131148 3.6 593642 15.9 15680924 418.8 Vcells 45479988 347.0 173526492 1324.0 220358611 1681.3 A description of the survey package can be found here: http://faculty.washington.edu/tlumley/survey/ I tried creating a work-around by using the database-backed survey objects (DB SO), included in the survey package to conserve memory on larger datasets like this one. Unfortunately, I don't think the survey package supports database connections for replicate weight designs yet, since I've only been able to get a database connection working after creating a svydesign object and not a svrepdesign object - and also because neither the DB SO website nor the svrepdesign help page make any mention of those parameters. The DB SOs are described in detail here: http://faculty.washington.edu/tlumley/survey/svy-dbi.html Any advice would be truly appreciated. Thanks, Anthony Damico [[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] access elements of a named list using a factor
do you mean: f - factor(c(pigs, pigs, slugs)) jj - list(pigs = 1:10, slugs = 1:3) jj[levels(f)[1]] jj[[levels(f)[1]]] Best, Dimitris Robin Hankin wrote: Hi I have a factor 'f' and a named list 'jj'. I want names(jj) to match up with levels(f). How do I use levels(f) to access elements of jj? f - factor(c(pigs,pigs,slugs)) f [1] pigs pigs slugs Levels: pigs slugs jj - list(pigs=1:10,slugs=1:3) My attempts to produce jj$pigs all give errors: jj$levels(f)[1] Error: attempt to apply non-function do.call($,jj,levels(f)[1]) Error in if (quote) { : argument is not interpretable as logical $(jj,levels(f)[1]) Error in jj$levels(f)[1] : invalid subscript type 'language' -- Dimitris Rizopoulos Assistant Professor Department of Biostatistics Erasmus University Medical Center Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands Tel: +31/(0)10/7043478 Fax: +31/(0)10/7043014 __ 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] Inserting rows
Hi all, I have the data set df with three varaibles, x1 x2 x3 1 2 5 2 4 1 5 6 0 1 1 2 I want to insert more rows ( eg, 3 rows with value filled with zeros) 1 2 5 2 4 1 5 6 6 1 1 2 0 0 0 0 0 0 0 0 0 Can any body help me out? 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] access elements of a named list using a factor
Hello Dimitris thanks for this. It works! I guess I was fixated on the dollar sign. I must confess that I don't really understand any of the error messages below. Can anyone help me interpret them? rksh Dimitris Rizopoulos wrote: do you mean: f - factor(c(pigs, pigs, slugs)) jj - list(pigs = 1:10, slugs = 1:3) jj[levels(f)[1]] jj[[levels(f)[1]]] Best, Dimitris Robin Hankin wrote: Hi I have a factor 'f' and a named list 'jj'. I want names(jj) to match up with levels(f). How do I use levels(f) to access elements of jj? f - factor(c(pigs,pigs,slugs)) f [1] pigs pigs slugs Levels: pigs slugs jj - list(pigs=1:10,slugs=1:3) My attempts to produce jj$pigs all give errors: jj$levels(f)[1] Error: attempt to apply non-function do.call($,jj,levels(f)[1]) Error in if (quote) { : argument is not interpretable as logical $(jj,levels(f)[1]) Error in jj$levels(f)[1] : invalid subscript type 'language' -- Robin K. S. Hankin Uncertainty Analyst University of Cambridge 19 Silver Street Cambridge CB3 9EP 01223-764877 __ 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] Inserting rows
Try this: rbind(df, matrix(0, 3, 3, dimnames = list(NULL, names(df On Fri, Oct 23, 2009 at 11:44 AM, Ashta sewa...@gmail.com wrote: Hi all, I have the data set df with three varaibles, x1 x2 x3 1 2 5 2 4 1 5 6 0 1 1 2 I want to insert more rows ( eg, 3 rows with value filled with zeros) 1 2 5 2 4 1 5 6 6 1 1 2 0 0 0 0 0 0 0 0 0 Can any body help me out? 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. -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40 S 49° 16' 22 O __ 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] data frame is killing me! help
bbslover wrote: Steve Lianoglou-6 wrote: Hi, On Oct 22, 2009, at 2:35 PM, bbslover wrote: Usage data(gasoline) Format A data frame with 60 observations on the following 2 variables. octane a numeric vector. The octane number. NIR a matrix with 401 columns. The NIR spectrum and I see the gasoline data to see below NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR.1694 nm NIR.1696 nm NIR.1698 nm NIR.1700 nm 1 1.242645 1.250789 1.246626 1.250985 1.264189 1.244678 1.245913 1.221135 2 1.189116 1.223242 1.253306 1.282889 1.215065 1.225211 1.227985 1.198851 3 1.198287 1.237383 1.260979 1.276677 1.218871 1.223132 1.230321 1.208742 4 1.201066 1.233299 1.262966 1.272709 1.211068 1.215044 1.232655 1.206696 5 1.259616 1.273713 1.296524 1.299507 1.226448 1.230718 1.232864 1.202926 6 1.24109 1.262138 1.288401 1.291118 1.229769 1.227615 1.22763 1.207576 7 1.245143 1.265648 1.274731 1.292441 1.218317 1.218147 1.73 1.200446 8 1.222581 1.245782 1.26002 1.290305 1.221264 1.220265 1.227947 1.188174 9 1.234969 1.251559 1.272416 1.287405 1.211995 1.213263 1.215883 1.196102 look at this NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR. 1694 nm NIR.1696 nm NIR.1698 nm NIR.1700 nm how can I add letters NIR to my variable, because my 600 independents never have NIR as the prefix. however, it is needed to model the plsr. for example aa=plsr(y~NIR, data=data ,), the prefix NIR is necessary, how can I do with it? I'm not really sue that I'm getting you, but if your problem is that the column names of your data.frame don't match the variable names you'd like to use in your formula, just change the colnames of your data.frame to match your formula. BTW - I have no idea where to get this gasoline data set, so I'm just imagining: eg. colnames(gasoline) - c('put', 'the', 'variable', 'names', 'that', 'you', 'want', 'here') -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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. thanks for you. but the numbers of indenpendence are so many, it is not easy to identify them one by one, is there some better way? You don't need to identify anything. What you need to do is read the help page for the function you want to use, so you (at the very least) know how to use the function. library(pls) data(gasoline) fit - plsr(octane~NIR, data=gasoline, validation = CV) summary(fit) Data: X dimension: 60 401 Y dimension: 60 1 Fit method: kernelpls Number of components considered: 53 VALIDATION: RMSEP Cross-validated using 10 random segments. (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps CV 1.5431.372 0.3827 0.2522 0.2347 0.2455 0.2281 adjCV1.5431.367 0.3740 0.2497 0.2360 0.2407 0.2243 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps CV 0.2311 0.2352 0.24550.25340.27370.28140.2832 adjCV 0.2257 0.2303 0.23950.24730.26460.27050.2726 14 comps 15 comps 16 comps 17 comps 18 comps 19 comps 20 comps CV 0.29130.29320.29850.31370.32890.33230.3391 adjCV0.28080.28210.28630.30080.31410.31720.3228 21 comps 22 comps 23 comps 24 comps 25 comps 26 comps 27 comps CV 0.34760.33840.33160.32130.31550.31180.3062 adjCV0.33070.32170.31540.30570.30020.29640.2908 28 comps 29 comps 30 comps 31 comps 32 comps 33 comps 34 comps CV 0.30330.30340.30740.30830.30940.30870.3105 adjCV0.28810.28810.29170.29260.29360.29290.2946 35 comps 36 comps 37 comps 38 comps 39 comps 40 comps 41 comps CV 0.31080.31060.31050.31040.31040.31050.3105 adjCV0.29490.29470.29460.29450.29450.29450.2946 42 comps 43 comps 44 comps 45 comps 46 comps 47 comps 48 comps CV 0.31050.31050.31050.31050.31050.31050.3105 adjCV0.29460.29460.29460.29460.29460.29460.2946 49 comps 50 comps 51 comps 52 comps 53 comps CV 0.31050.31050.31050.31050.3105 adjCV0.29460.29460.29460.29460.2946 TRAINING: % variance explained 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps X 70.9778.5686.15 95.496.1296.9797.32 98.1 octane31.9094.6697.71 98.098.6898.9399.06 99.1 9
Re: [R] reset par() within plot layout
Janke, Janke ten Holt schrieb: Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) look at par.ini: it's a list with all the argument-value pairs for par(). You might be able to solve your problem by removing the appropriate elements from par.ini before calling par(par.ini). Do the following to look which ones need to be kept for the layout: par() layout(matrix(1:2, nr=2)) par() Tom par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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] reset par() within plot layout
Tom Gottfried wrote: Janke, Janke ten Holt schrieb: Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) look at par.ini: it's a list with all the argument-value pairs for par(). You might be able to solve your problem by removing the appropriate elements from par.ini before calling par(par.ini). Do the following to look which ones need to be kept for the layout: par() layout(matrix(1:2, nr=2)) par() Tom par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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. Yes, that had occured to me too. So I tried: layout(matrix(1:2, nr=2)) par(no.readonly=TRUE) plot(1:10) par(no.readonly=TRUE) This has differences in fig mfg usr xaxp yaxp But even keeping these back does not solve my problem. So I figured there must be something else going on that I am unaware of... btw, your exact suggestion, par() layout(matrix(1:2, nr=2)) par() does not result in any differences. __ 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] Functional data analysis - problems with smoothing
Ben, I lose my bet! Good thing my opinion is never worth more than $.02. Looking at the code for smooth.spline() I see that, although 'argvals' is the first argument to the function, the error message is generated by the second argument. [** Spencer, is this intentional? **] So, Ben, try str(rlz) or class(rlz) and see whether it conforms to one of the structures described in ?smooth.basis. -Peter Ehlers Benjamin Cheah wrote: Hi Peter and Spencer (and everyone else out there), Thank you for your prompt reply to my post re. FDA - FYI - i'm very new to R and FDA. Currently half way through a biostats masters degree in Sydney, Australia! Absolutely loving it and I hope to include R computing on my future CV! Just to clarify things: When I entered what you suggested, this is what I obtained with my own dataset str(isi) num [1:14] 1 1.5 2 2.5 3 3.5 4 5 7 10 ... Similarly, for Jim Ramsay's dataset (chapter five of Use R book): str(age) num [1:31] 1 1.25 1.5 1.75 2 3 4 5 6 7 ... Doesn't this mean that my data object, 'isi' is numeric? Also, I was looking through Jim Ramsay's datasets - I'm fairly sure that my vector, 'isi' was similarly organised to his vector, 'age' in chapter 5 of his new Use R FDA book, so I'm fairly certain that the data is numeric but it obviously isn't Hope this email clarified things - apologies if it hasn't - still very new to the terminology and overall feel of R. The R user network is terrific! Thank you for caring! Ben From: Peter Ehlers ehl...@ucalgary.ca Cc: R-help@r-project.org Sent: Fri, 23 October, 2009 5:07:15 PM Subject: Re: [R] Functional data analysis - problems with smoothing The error message is pretty clear: regardless of what *you* think, R says that 'isi' is not numeric. Are you sure that 'isi' is not a *factor* object? I'm willing to bet that it is. Use str() to check your data. -Peter Ehlers Benjamin Cheah wrote: Hi all, I'm having major issues with smoothing my functional data I'm referring to Jim Ramsay's examples in his books. The following error message keeps appearing, despite all my data being numeric can anyone kindly offer any suggestions? isi - vector of argument values - i.e. the independent variable of the curves rlz - data array TMSfdPar - functional data parameter. TMSfdPar = fdPar(TMSbasis, 4, 0.01) TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar) Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric. I don't understand why the error message keeps popping up. I've tried playing around with Jim Ramsay's datasets and I think my data is organised in a similar manner to his, so can't understand what's going on oh, the frustration! Thanks in advance, Ben (the amateur R data analyst and statistician.) __ [[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-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] Inserting rows
?rbind df1 - data.frame(matrix(rep(0,9),nrow=3)) names(df1) - c(x1,x2,x3) rbind(df,df1) --- On Fri, 10/23/09, Ashta sewa...@gmail.com wrote: From: Ashta sewa...@gmail.com Subject: [R] Inserting rows To: R help r-help@r-project.org Received: Friday, October 23, 2009, 9:44 AM Hi all, I have the data set df with three varaibles, x1 x2 x3 1 2 5 2 4 1 5 6 0 1 1 2 I want to insert more rows ( eg, 3 rows with value filled with zeros) 1 2 5 2 4 1 5 6 6 1 1 2 0 0 0 0 0 0 0 0 0 Can any body help me out? 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. __ Looking for the perfect gift? Give the gift of Flickr! http://www.flickr.com/gift/ __ 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] Change positions of columns in data frame
Hi all, Probably a simple question, but I just can't find a simple answear in the older threads or anywhere else. I've added some new vectors as columns in a data frame using cbind(). As they're all put as the last columns inte the data frame, I would like to move them to specific positions. How do you do to change the position of a column in a data frame? I know I can use fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame with the given columns in the specified order, but there must be an easier way..? All the best, Joel _ Nya Windows 7 - Hitta en dator som passar dig! Mer information. http://windows.microsoft.com/shop [[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] Bayesian regression stepwise function?
If 'fools rush in where angels fear to tread', then Bayesians 'jump' in where frequentists fear to 'step'... Very nice, Chuck! Definitely one for my list of fortunes. Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Charles C. Berry Sent: Friday, October 23, 2009 1:18 AM To: Ben Bolker Cc: r-help@r-project.org Subject: Re: [R] Bayesian regression stepwise function? On Thu, 22 Oct 2009, Ben Bolker wrote: Allan.Y wrote: Hi everyone, I am wondering if there exists a stepwise regression function for the Bayesian regression model. I tried googling, but I couldn't find anything. I know step function exists for regular stepwise regression, but nothing for Bayes. Why? That seems so ... un-Bayesian ... If 'fools rush in where angels fear to tread', then Bayesians 'jump' in where frequentists fear to 'step'... Seriously, there are Bayesian regression approaches that priorize the model size (sometimes only implicitly by assigning a prior for the inclusion of each candidate regressors). Then they 'jump' between models of different sizes. On CRAN, Package qtlbim (which is specialized to a particular genetics problem) implements one such, I think. Package bqtl does not implement the jumping approach, but does explore a model space with differing numbers of regressors for the same (qtl) problem. Perhaps the closest to a general purpose 'stepwise flavored' Bayesian regression is implemented in Package BMA, which IIRC borrows step() for some of its work. But CRAN now has more packages than my cortex has neurons, so there are probably more packages that do something like this. Try RSiteSearch(jump regression, restric='functions') and start reading. HTH, Chuck -- View this message in context: http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p2601 5081.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. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ 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] PDF too large, PNG bad quality
Hello Lasse, Why not try this? (1) Create 20MB PDF from R (2) use convert command in linux, examples below: convert -resize 50% 20mbfile.pdf smallerfile.pdf convert -resize 75% 20mbfile.pdf image.png Ghostscript can help you as well for conversion! Using vector formats (pdf,ps,eps) are good for this purpose, as opposed to scalar (png, jpg, etc.). Best, Parthiban. 2009/10/23 Jim Lemon j...@bitwrit.com.au On 10/23/2009 06:07 AM, Lasse Kliemann wrote: I wish to save a scatter plot comprising approx. 2 million points in order to include it in a LaTeX document. Using 'pdf(...)' produces a file of size about 20 MB, which is useless. Using 'cairo_pdf(...)' produces a smaller file, around 3 MB. This is still too large. Not only that the document will be too large, but also PDF viewers choke on this. Moreover, Cairo has problems with text: by default text looks ugly, like scaled bitmaps. After hours of trying different settings, I discovered that choosing a different font family can help, e.g.: 'par(family=Mono)'. This gives good-looking text. Yet, the problem with the file size remains. There exists the hint to produdc EPS instead and then convert to PDF using 'epstopdf'. The resulting PDF files are slightly smaller, but still too large, and PDF viewers still don't like it. So I gave PNG a try. PNG files are much smaller and PDF viewers have no trouble with them. However, fonts look ugly. The same trick that worked for Cairo PDF has no effect for PNG. When I view the PNGs with a dedicated viewer like 'qiv', even the fonts look good. But not when included in LaTeX; I simply use '\includegraphics{...}' and run the document through 'pdflatex'. I tried both, creating PNG with 'png(...)' and converting from PDF to PNG using 'convert' from ImageMagick. So my questions are: - Is there a way to produce sufficiently lean PDFs directly in R, even when the plot comprises several million points? - How to produce a PNG that still looks nice when included in a LaTeX PDF document? Any hints will be greatly appreciated. Hi Lasse, I may be right off the track, but I can't make much sense of 2 million points in a scatterplot. If you are interested in the density of points within the plot, you could compute this using something like the bkde2 function in the KernSmooth package and then plot that using something like image. 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. [[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] Change positions of columns in data frame
dataframe xx x1 x2 x3 1 2 5 2 4 1 5 6 0 1 1 2 data.frame(xx$x2,xx$x1,xx$x3) # or Awkward but works --- On Fri, 10/23/09, Joel Fürstenberg-Hägg joel_furstenberg_h...@hotmail.com wrote: From: Joel Fürstenberg-Hägg joel_furstenberg_h...@hotmail.com Subject: [R] Change positions of columns in data frame To: r-help@r-project.org Received: Friday, October 23, 2009, 10:19 AM Hi all, Probably a simple question, but I just can't find a simple answear in the older threads or anywhere else. I've added some new vectors as columns in a data frame using cbind(). As they're all put as the last columns inte the data frame, I would like to move them to specific positions. How do you do to change the position of a column in a data frame? I know I can use fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame with the given columns in the specified order, but there must be an easier way..? All the best, Joel _ Nya Windows 7 - Hitta en dator som passar dig! Mer information. http://windows.microsoft.com/shop [[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. __ Be smarter than spam. See how smart SpamGuard is at giving ju __ 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] Nested logit in GEV family
Hi All, I know someone has posted similar messages before, but there is no reply. So I wonder whether there is a way to run a nested logit model in R. It is in the GEV family; however, the commands fitting GEV don't seem to work for nested logit. Or maybe I have to write down the log-liklihood function and let R maximize it? Thank you! Best wishes. Min -- Min Chen Graduate Student Department of Agricultural, Food, and Resource Economics 208 Cook Hall Michigan State University chenm...@msu.edu / chenmin0...@gmail.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] How to make R packages?
These two are for windows. I have a linux machine. Are there any other reference for linux? On Fri, Oct 23, 2009 at 12:08 AM, Jorge Ivan Velez jorgeivanve...@gmail.com wrote: Hi Peng, Also, take a look at the following links: http://robjhyndman.com/research/Rpackages_notes.pdf http://www.biostat.wisc.edu/~kbroman/Rintro/Rwinpack.html HTH, Jorge On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu wrote: I found the following document on making R packages. But it is old. I'm wondering if there is more current ones and hopefully more complete ones. http://biosun1.harvard.edu/courses/individual/bio271/lectures/L6/Rpkg.pdf __ 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] Documenting related data sets
Say I want to include the data from Minard's march on Moscow graphic in a package. They consist of three data frames, Minard.troops, Minard.cities, Minard.temp. I've determined that I can document them all in one .Rd file, Minard.Rd, that starts: \name{Minard} \Rdversion{1.1} \alias{Minard} \alias{Minard.cities} \alias{Minard.troops} \alias{Minard.temp} \docType{data} \title{ Data from Minard's famous graphic map of Napoleon's march on Moscow } ... \usage{ data(Minard.troops) data(Minard.cities) data(Minard.temp) } ... But, AFAICS, it seems I have to save each of these in a separate .RData / .rda file for them to be found by R CMD check. Is this correct, or is there some way to include them all in a Minard.RData, via save(Minard.troops, Minard.cities, Minard.temp, file=Minard.RData) -Michael -- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html Toronto, ONT M3J 1P3 CANADA __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Change positions of columns in data frame
DF[ ,names(DF)[c(your_preferred_order)]] -Peter Ehlers Joel Fürstenberg-Hägg wrote: Hi all, Probably a simple question, but I just can't find a simple answear in the older threads or anywhere else. I've added some new vectors as columns in a data frame using cbind(). As they're all put as the last columns inte the data frame, I would like to move them to specific positions. How do you do to change the position of a column in a data frame? I know I can use fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame with the given columns in the specified order, but there must be an easier way..? All the best, Joel _ Nya Windows 7 - Hitta en dator som passar dig! Mer information. http://windows.microsoft.com/shop [[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] standard errors in bbmle
Peter Dalgaard wrote: alexander russell wrote: Hello, Mle2 is a little unforthcoming in the matter of standard errors? Is there a way to ask the program to supply standard errors along with estimates in cases when it doesn't print them 'voluntarily'? regards, s What did you do? Which cases? Did you look at the examples? (as in, example(mle2) ) As the package author: +1. Perhaps you're looking for summary()? Ben Bolker -- View this message in context: http://www.nabble.com/standard-errors-in-bbmle-tp26023856p26028127.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] removing random effect from nlme or using varPower() in nls
Hi Kingsford, That's exactly what I was looking for! Thanks. Mike On Thu, 22 Oct 2009, Kingsford Jones wrote: help(gnls, pack=nlme) hth, Kingsford Jones On Thu, Oct 22, 2009 at 4:36 PM, Michael A. Gilchrist mi...@utk.edu wrote: Hello, I've been fitting a random effects model using nlme to some data, but I am discovering that the variation in my random effect is very small. As a result, I would like to replace it as a fixed effect (i.e. essentially fit the same model but with no random effect). As I understand it I could do this using nls(), but I'm using a number of options such as weights = varPower() which I am at a loss on how to implement in that framework. Is there a way to use nlme but with out a random effect? (a bit absurd, I realize, but I have the syntax working...) Alternatively, is there a way to use weights = varPower() with nls? Any help would be appreciated. Sincerely, Mike __ 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] cut2 once, bin twice...
On Fri, Oct 23, 2009 at 3:58 AM, Dieter Menne dieter.me...@menne-biomed.de wrote: sdanzige wrote: I'm using the Hmisc cut2 function to bin a set of data. It produces bins that I like with results like this: [96,270]:171 [69, 96): 54 [49, 69): 40 [35, 49): 28 [28, 35): 14 [24, 28): 8 (Other) : 48 I would like to take a second set of data, and assign it to bins based on factors defined by my call to cut 2. It used to be quite tricky, but on popular request Brian Ripley has added an example how to extract the intervals using regular expression on the bottom of the examples for cut (note:cut in base, not cut2 in Hmisc). If someone knows of an easier way, please correct me. How about adding this information as attribute to the standard cut? The strapply function in gsubfn can do it with a simpler regular expression since it extracts based on content rather than delimiters, which is what you want here: # create sample data library(gsubfn) set.seed(1) dat - seq(4, 7, by = 0.05) x - sample(dat, 30) . # use cut groups - cut(x, breaks = 10) # extract interval boundaries using strapply strapply(levels(groups), [[:digit:].]+, as.numeric, simplify = TRUE) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 4.0 4.3 4.6 4.9 5.2 5.5 5.8 6.1 6.4 6.7 [2,] 4.3 4.6 4.9 5.2 5.5 5.8 6.1 6.4 6.7 7.0 The above is from demo(gsubfn-cut) For more see the gsubfn home page at http://gsubfn.googlecode.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 find moving averages within each subgroup of a data frame
Try this: df1$SMA - ave(df1$col4, df1$col1, FUN = function(x) c(NA, (head(x, -1) + tail(x, -1))/2)) It would also be possible to convert it from long form to wide form using reshape (or read.zoo in the devel version of zoo), convert that to a zoo series and the use rollapply in the zoo package. On Thu, Oct 22, 2009 at 12:28 PM, clair.crossup...@googlemail.com clair.crossup...@googlemail.com wrote: Dear all, If I have the following data frame: set.seed(21) df1 - data.frame(col1=c(rep('a',5), rep('b',5), rep('c',5)), col4=rnorm(1:15)) col1 col4 1 a 0.793013171 2 a 0.522251264 3 a 1.74641 4 a -1.271336123 5 a 2.197389533 6 b 0.433130777 7 b -1.570199630 8 b -0.934905667 9 b 0.063493345 10 b -0.002393336 11 c -2.276781240 12 c 0.757412225 13 c -0.548405554 14 c 0.172549478 15 c 0.562853068 How can i make a 2 point moving average within each group? i.e. col1 col4 SMA a 0.793013171 NA a 0.522251264 0.657632218 a 1.74641 1.134236753 a -1.271336123 0.237443059 a 2.197389533 0.463026705 b 0.433130777 NA b -1.57019963 -0.568534427 b -0.934905667 -1.252552649 b 0.063493345 -0.435706161 b -0.002393336 0.030550005 c -2.27678124 NA c 0.757412225 -0.759684508 c -0.548405554 0.104503336 c 0.172549478 -0.187928038 c 0.562853068 0.367701273 From what i've read, i was thinking it would be something along the lines of: aggregate(df1$col4, by=list(df1$col1), function(x) {filter(x, rep(1/2,2), sides=1 )} ) Error in aggregate.data.frame(as.data.frame(x), ...) : 'FUN' must always return a scalar But this tells me (i think) that aggregate should only return a single value per group. So what i need, i guess, is something that takes all the values in a given group, and returns a vector of the same length. Not sure which function to use for that. Thanks in advance, C xx P.S. on a side note, is it possible to extract the values which aggregate passes to the function(x) in my example above? __ 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] reset par() within plot layout
Janke ten Holt wrote: Tom Gottfried wrote: Janke, Janke ten Holt schrieb: Dear list, I would like to produce a matrix of plots, where par() is reset after each plot (see below [simplified] example). When I use layout() to do so, I seem to also reset the layout. I have not been able to figure out how to prevent this from happening. Any help is greatly appreciated! Janke Example code: #Desired result is a layout of 2 plots: one red and one black layout(matrix(1:2, nr=2)) par.ini - par(no.readonly=TRUE) look at par.ini: it's a list with all the argument-value pairs for par(). You might be able to solve your problem by removing the appropriate elements from par.ini before calling par(par.ini). Do the following to look which ones need to be kept for the layout: par() layout(matrix(1:2, nr=2)) par() Tom par(col=red) plot(1:100) par(par.ini) plot(1:10) -- Janke ten Holt Dept. of Psychology/Sociology University of Groningen, the Netherlands __ 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. Yes, that had occured to me too. So I tried: layout(matrix(1:2, nr=2)) par(no.readonly=TRUE) plot(1:10) par(no.readonly=TRUE) This has differences in fig mfg usr xaxp yaxp But even keeping these back does not solve my problem. So I figured there must be something else going on that I am unaware of... I think that what you want is layout(matrix(1:2, nr=2)) opar - par(col=red) plot(1:100) par(opar) plot(1:10) -Peter Ehlers btw, your exact suggestion, par() layout(matrix(1:2, nr=2)) par() does not result in any differences. __ 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] How to make R packages?
Hi Pen, On Oct 23, 2009, at 10:31 AM, Peng Yu wrote: These two are for windows. I have a linux machine. Are there any other reference for linux? The steps are pretty much the same. For instance, in robjhyndman.com link, you can simply start from step 2 and go on. Also, the official Writing R Extensions covers this topic pretty thoroughly: http://cran.r-project.org/doc/manuals/R-exts.html#Creating-R-packages -steve On Fri, Oct 23, 2009 at 12:08 AM, Jorge Ivan Velez jorgeivanve...@gmail.com wrote: Hi Peng, Also, take a look at the following links: http://robjhyndman.com/research/Rpackages_notes.pdf http://www.biostat.wisc.edu/~kbroman/Rintro/Rwinpack.html HTH, Jorge On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu wrote: I found the following document on making R packages. But it is old. I'm wondering if there is more current ones and hopefully more complete ones. http://biosun1.harvard.edu/courses/individual/bio271/lectures/L6/Rpkg.pdf __ 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. -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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] reference for mahalanobis {stats}
The help on mahalanobis {stats} does not include any reference. I'm interested in understand why Mahalanobis is defined in its current way and how to use it. Could somebody point me a good book on this? I have looked through a few books, but they all give very light explanation on it. __ 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] cut2 once, bin twice...
Dieter Menne wrote: It used to be quite tricky, but on popular request Brian Ripley has added an example how to extract the intervals using regular expression on the bottom of the examples for cut (note:cut in base, not cut2 in Hmisc). Thank you, but the regular expression example doesn't seem to work correctly. labs-levels(df$p_bin) labs [1] 012345 [7] 6789 10 11 [13] 12 13 14 15 16 17 [19] 18 19 20 [21, 24) [24, 28) [28, 35) [25] [35, 49) [49, 69) [69, 96) [96,270] cbind(lower = as.numeric( sub(\\((.+),.*, \\1, labs) ), upper = as.numeric( sub([^,]*,([^]]*)\\], \\1, labs) )) Warning in cbind(lower = as.numeric(sub(\\((.+),.*, \\1, labs)), upper = as.numeric(sub([^,]*,([^]]*)\\], : NAs introduced by coercion Warning in cbind(lower = as.numeric(sub(\\((.+),.*, \\1, labs)), upper = as.numeric(sub([^,]*,([^]]*)\\], : NAs introduced by coercion lower upper [1,] 0 0 [2,] 1 1 [3,] 2 2 [4,] 3 3 [5,] 4 4 [6,] 5 5 [7,] 6 6 [8,] 7 7 [9,] 8 8 [10,] 9 9 [11,]1010 [12,]1111 [13,]1212 [14,]1313 [15,]1414 [16,]1515 [17,]1616 [18,]1717 [19,]1818 [20,]1919 [21,]2020 [22,]NANA [23,]NANA [24,]NANA [25,]NANA [26,]NANA [27,]NANA [28,]NA 270 -- Any ideas? Thank you, -S -- View this message in context: http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26027643.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] Change positions of columns in data frame
Joel - Suppose the columns are named x1, x2, x3, x4, and x5. You can use subscripting: x[c('x2','x4','x1','x3','x5')] or, to save typing the quotes subset(x,select=c(x2,x4,x1,x3,x5)) - Phil Spector Statistical Computing Facility Department of Statistics UC Berkeley spec...@stat.berkeley.edu On Fri, 23 Oct 2009, Joel Fürstenberg-Hägg wrote: Hi all, Probably a simple question, but I just can't find a simple answear in the older threads or anywhere else. I've added some new vectors as columns in a data frame using cbind(). As they're all put as the last columns inte the data frame, I would like to move them to specific positions. How do you do to change the position of a column in a data frame? I know I can use fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame with the given columns in the specified order, but there must be an easier way..? All the best, Joel _ Nya Windows 7 - Hitta en dator som passar dig! Mer information. http://windows.microsoft.com/shop [[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] cut2 once, bin twice...
sdanzige wrote: Thank you, but the regular expression example doesn't seem to work correctly. I wrote a regular expression that does seem to work, so I'll post it here for anyone else that needs it. labs-levels(df$p_bin) cbind(lower=as.numeric(sub([[(],,sub(,.*,,labs))), upper=as.numeric(sub([])],,sub([[(].*, *,,labs))) ) I fear my inelegance will peg me as a Windows programmer, but so be it... -S -- View this message in context: http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26028296.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] reference for mahalanobis {stats}
Hi, On Oct 23, 2009, at 11:36 AM, Peng Yu wrote: The help on mahalanobis {stats} does not include any reference. I'm interested in understand why Mahalanobis is defined in its current way and how to use it. Could somebody point me a good book on this? I have looked through a few books, but they all give very light explanation on it. The wikipedia page looks pretty good ... the Intuitive Explanation section should help: http://en.wikipedia.org/wiki/Mahalanobis_distance Googling also provides many more explanations: http://matlabdatamining.blogspot.com/2006/11/mahalanobis-distance.html HTH, -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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] Determining which file is newer
Colleagues, I wish to execute a task only if a particular file is newer than a second file. I can access the file modification dates using file.info()$mtime. This yields a time object (? POSIX). sizeisdir modemtime ctime Testfile4421FALSE 755 2009-10-23 08:59:09 2009-10-23 08:59:09 What is the simplest means to compare these two objects? R2.9.1 in OSX, Windows, Linux Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-866-PLessThan (1-866-753-7784) www.PLessThan.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] opposite estimates from zeroinfl() and hurdle()
Tord Snäll-4 wrote: Dear all, A question related to the following has been asked on R-help before, but I could not find any answer to it. Input will be much appreciated. I got an unexpected sign of the slope parameter associated with a covariate (diam) using zeroinfl(). It led me to compare the estimates given by zeroinfl() and hurdle(): [snip] The right thing to do in this case is to poke through the code of hurdle() and zeroinfl(), but a simple (?) demonstration shows that hurdle() and zeroinfl() are indeed reporting opposite values : hurdle reports -log(p/(1-p)) = -qlogis(p), where p is the probability of a zero count: z = rpois(500,lambda=3) z = (z[z0])[1:90] z = c(z,rep(0,10)) hurdle(z~1) ## -qlogis(0.1) ## zero coefficient always == -qlogis(0.1) zeroinfl reports log(p/(1-p)), where p is the zero-inflation: z = rpois(90,lambda=3) z = c(z,rep(0,10)) zeroinfl(z~1) ## qlogis(0.1) tmpf = function() { z = rpois(90,lambda=3) z = c(z,rep(0,10)) coef(zeroinfl(z~1))[2] } rr = replicate(1000,tmpf()) hist(rr,breaks=1000) summary(rr) qlogis(0.1) Perhaps it would be worth sending an e-mail to the package maintainers to request a note to this effect in the documentation, particularly if this a FAQ ... -- View this message in context: http://www.nabble.com/opposite-estimates-from-zeroinfl%28%29-and-hurdle%28%29-tp26024735p26029131.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] arima crashes too
Barry Rowlingson wrote: If you're doing anything in a loop that has the potential to fail because of singularities or other conditions when your model can't be fitted, you need to stick what you are doing in a 'try' clause. This lets you trap errors and do something with them. Plenty of examples in help(try) or this from me: for(i in 1:10){ print(solve(matrix(c(3,3,3,i),2,2))) } This stops the loop at i=3. Now stick it in a try() clause: for(i in 1:10){ print(try(solve(matrix(c(3,3,3,i),2,2 } and it gives a warning and carries on. If you want your code to do something with the failure cases then the help for try() tells you what to look for. I will try try. # irreristible pun But... I'm not sure why your arima produces an error, but I'm assuming the numbers are such that the model can't be fitted. I don't really know what arima is doing. That's a problem. In this case, I'm fitting an AR(1) time series using arima and arma. The documentation does not mention that we should be careful about what numbers are passed to it; on contrary, it's perfectly possible (in theory) that an AR(1) time series was created using an evil phi in x[i+1] - phi * x[i] + sigma * n01[i] and the regression should return the evil phi. As in this reproducible example: phi - 1.5 sigma - 0.5 n01 - c(-1, 0.3, -0.3, 0.2, -0.7) x - rep(0, 6) for (i in 1:5) x[i+1] - phi * x[i] + sigma * n01[i] arma(x, order=c(1,0)) arima(x, order=c(1,0,0)) It seems like arma does not force phi to be in the -1:1 range, because it correctly returns evil phi's, but arima kind of bayesianises the output, forcing the time series to be stationary. Probably a longer example may scrash/s stop arima but not arma. Non-reproducible example (almost surely scrashes/s stops arima but plays safely in arma): phi - 1.05 sigma - 0.5 n - 100 n01 - rnorm(n) x - rep(0, n) for (i in 1:n) x[i+1] - phi * x[i] + sigma * n01[i] arma(x, order=c(1,0)) arima(x, order=c(1,0,0)) Increasing phi (or n) will scrash/s stop both arma and arima: phi - 1.1 sigma - 0.5 n - 100 n01 - rnorm(n) x - rep(0, n) for (i in 1:n) x[i+1] - phi * x[i] + sigma * n01[i] arma(x, order=c(1,0)) arima(x, order=c(1,0,0)) Alberto 'Darth Albmont' Monteiro __ 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] Determining which file is newer
Have a look at your company name and your phone number. Nudge nudge, wink wink, say no more... /H On Fri, Oct 23, 2009 at 9:02 AM, Dennis Fisher fis...@plessthan.com wrote: Colleagues, I wish to execute a task only if a particular file is newer than a second file. I can access the file modification dates using file.info()$mtime. This yields a time object (? POSIX). size isdir mode mtime ctime Testfile 4421 FALSE 755 2009-10-23 08:59:09 2009-10-23 08:59:09 What is the simplest means to compare these two objects? R2.9.1 in OSX, Windows, Linux Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-866-PLessThan (1-866-753-7784) www.PLessThan.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@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] Determining which file is newer
On Fri, Oct 23, 2009 at 5:02 PM, Dennis Fisher fis...@plessthan.com wrote: Colleagues, I wish to execute a task only if a particular file is newer than a second file. I can access the file modification dates using file.info()$mtime. This yields a time object (? POSIX). size isdir mode mtime ctime Testfile 4421 FALSE 755 2009-10-23 08:59:09 2009-10-23 08:59:09 What is the simplest means to compare these two objects? Split the string representation of the object into date and time parts, then split those parts into day, month , year and hour minute second, taking care of internationalisation and localisation, then recombine all those parts into a number of seconds since Jan 1 1970, taking care of leap seconds and the slowing down of the earth due to the tides of the moon. Or just use the usual operators: file.info(/etc/passwd)$mtime = file.info(/etc/motd)$mtime [1] TRUE Barry __ 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] Change positions of columns in data frame
You can change the order of nearly everything just by applying an appropriate index and assigning to the original object. Something like data - data[,order_of_columns] Tom Joel Fürstenberg-Hägg schrieb: Hi all, Probably a simple question, but I just can't find a simple answear in the older threads or anywhere else. I've added some new vectors as columns in a data frame using cbind(). As they're all put as the last columns inte the data frame, I would like to move them to specific positions. How do you do to change the position of a column in a data frame? I know I can use fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame with the given columns in the specified order, but there must be an easier way..? All the best, Joel _ Nya Windows 7 - Hitta en dator som passar dig! Mer information. http://windows.microsoft.com/shop [[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. -- Technische Universität München Department für Pflanzenwissenschaften Lehrstuhl für Grünlandlehre Am Hochanger 1 85350 Freising / Germany Phone: ++49 (0)8161 715324 Fax: ++49 (0)8161 713243 email: tom.gottfr...@wzw.tum.de http://www.wzw.tum.de/gruenland __ 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 make R packages?
On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu pengyu...@gmail.com wrote: I found the following document on making R packages. But it is old. I'm wondering if there is more current ones and hopefully more complete ones. Friedrich Leisch made an excellent tutorial for making R packages, with a tutorial that works in Linux. You can find it here; http://epub.ub.uni-muenchen.de/6175/ with some additional commentary here: http://blog.revolution-computing.com/2009/08/creating-r-packages-a-tutorial-draft.html # David Smith -- David M Smith da...@revolution-computing.com VP of Community, REvolution Computing http://blog.revolution-computing.com Tel: +1 (206) 577-4778 x3203 (Palo Alto, CA, USA) Download REvolution R free: http://revolution-computing.com/downloads/revolution-r.php __ 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] Capturing lagged values?
Consider a data structure where header rows supply ID values (school and grade in the example below) that are not captured in subsequent rows df - read.table(textConnection( school grade studentid score 101 5 NA NA NA NA 123 21 NA NA 124 25 201 6 NA NA NA NA 231 33), header=TRUE) closeAllConnections() How do I capture the missing values for school and grade? Art * Art Burke Associate, Evaluation Program Education Northwest 101 SW Main St, Ste 500 Portland, OR 97204 Phone: 503.275.9592 art.bu...@educationnorthwest.org mailto:art.bu...@educationnorthwest.org http://educationnorthwest.org http://educationnorthwest.org/ We have recently changed our name to Education Northwest from Northwest Regional Educational Laboratory. Please note the new e-mail and Web addresses in the signature above. You may continue to find us on the Web at http://www.nwrel.org http://www.nwrel.org/ for the immediate future as well. [[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] coxph() and survfit()
On Thu, 22 Oct 2009, Mura Tamakou wrote: Dear All, I have a question regarding the output of survfit() when I supply a Cox model. Lets say for example: snipped: code that doesn't quite run we observe that the 95% CIs overlap!! How is this possible since the HR for spiders is significant. It's perfectly natural. To a good approximation, the p-value for the comparison will be significant when the point estimate for each group is outside the confidence interval for the other group. Suppose you had two point estimates with standard error equal to 1.0. The standard error of the difference would be 1.414, so the p-value would be less than 0.05 if the two point estimates differ by more than 1.96*1.414. The confidence intervals will overlap if the point estimates differ by less than 1.96*2. -thomas Thomas Lumley Assoc. Professor, Biostatistics tlum...@u.washington.eduUniversity of Washington, Seattle __ 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] ggplot2: stat_bin ..count.. with geom_text when NA is present
One for the ggplot2 gurus... I have a function which makes a plot just fine if the response vector (res in the example; fac1 is a factor) has no NA in it. It plots the data, then makes a little annotation at the bottom with the data counts using: p - p + geom_text(aes(x = fac1, y = min(res) - 0.1 * diff(range(res)), label = paste(n = , ..count.. , sep = )), color = black, size = 4.0, stat = bin) If there are NA in the res vector, I get warnings from stat_summary and geom_point about removing rows; these arise from an earlier part of the function and the points and error bars all plot. However, the count annotation does not appear on the plot when there are NA in res. Looking at the ggplot2 web site, there is a drop parameter for stat_bin. I inserted drop = TRUE several places in the snippet above and the function did not complain but still did not plot the counts. I looked at the function bin{ggplot2} which apparently does the work. There are some programming tricks there I'm not really familiar with, but generally it looks like it na.rm or na.omit's in several places, while the drop = TRUE is carried out as the last step. So, any suggestions about why the counts don't appear on my plot? I suppose I can always clean the data first, but it would be much more practical to do that in the background during the preparation of the plot. Thanks as always, Bryan * Bryan Hanson Acting Chair Professor of Chemistry Biochemistry DePauw University, Greencastle IN USA __ 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 Problems with CSV and Survey Objects
Yes, a 350Mb data frame is a bit big for 32-bit R to handle conveniently. As you note, the survey package doesn't yet do database-backed replicate-weight designs. You can get the same effect yourself without too much work. First, put the data into a database, such as SQLite. If you have the data frame read in then dbWriteTable will do it. Now, drop most of the variables, keeping the sampling weights, replicate weights, and a couple of other variables. Create a svrepdesign() with the reduced data set. When you want to do an analysis, use dbGetQuery() to load the variables you need for the analysis, and put them in the $variables component of the svrepdesign. That's exactly what the database-backed functions do for svydesign objects. [If you only ever want to use a small subset of the variables, it's even easier: drop all the extraneous variables and create a svrepdesign with the variables you want] -thomas On Fri, 23 Oct 2009, Anthony Damico wrote: I'm working with a 350MB CSV file on a server that has 3GB of RAM, yet I'm hitting a memory error when I try to store the data frame into a survey design object, the R object that stores data for complex sample survey data. When I launch R, I execute the following line from Windows: C:\Program Files\R\R-2.9.1\bin\Rgui.exe --max-mem-size=2047M Anything higher, and I get an error message saying the maximum has been set to 2047M. Here are the commands: library(survey) #this step takes more than five minutes data08-read.csv(data08.csv,header=TRUE,nrows=210437) object.size(data08) #329877112 bytes #Looking at Windows Task Manager, Mem Usage for Rgui.exe is already 659,632K brr.dsgn -svrepdesign( data = data08 , repweights = data08[, grep( ^repwgt , colnames( data08)) ], type = BRR , combined.weights = TRUE , weights = data08$mainwgt ) #Error: cannot allocate vector of size 254.5 Mb #The survey design object does not get created. #This also causes Windows Task Manager, Mem Usage to spike to 1,748,136K #And here are some memory diagnostics memory.limit() [1] 2047 memory.size() [1] 1449.06 gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 131148 3.6 593642 15.9 15680924 418.8 Vcells 45479988 347.0 173526492 1324.0 220358611 1681.3 A description of the survey package can be found here: http://faculty.washington.edu/tlumley/survey/ I tried creating a work-around by using the database-backed survey objects (DB SO), included in the survey package to conserve memory on larger datasets like this one. Unfortunately, I don't think the survey package supports database connections for replicate weight designs yet, since I've only been able to get a database connection working after creating a svydesign object and not a svrepdesign object - and also because neither the DB SO website nor the svrepdesign help page make any mention of those parameters. The DB SOs are described in detail here: http://faculty.washington.edu/tlumley/survey/svy-dbi.html Any advice would be truly appreciated. Thanks, Anthony Damico [[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. Thomas Lumley Assoc. Professor, Biostatistics tlum...@u.washington.eduUniversity of Washington, Seattle __ 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 do you know which functions are being debugged?
This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? Thanks, Andrew [[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] Bonferroni with unequal sample sizes
Hello- I have run an ANOVA on 4 treatments with unequal sample sizes (n=9,7,10 and 10). I want to determine where my sig. differences are between treatments using a Bonferroni test, and have run the code: pairwise.t.test(Wk16, Treatment, p.adf=bonf) I receive an error message stating that my arguments are of unequal length: Error in tapply(x, g, mean, na.rm = TRUE) : arguments must have same length Is there a way to run this test even with unequal sample sizes? Thanks. Stephanie Reiner Andrews Hall 410 Shipping: Rt. 1208 Greate Road Gloucester Point, VA 23062 work: 804-684-7869 cell: 443-286-4795 __ 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] Bonferroni with unequal sample sizes
On 10/23/2009 1:30 PM, slrei...@vims.edu wrote: Hello- I have run an ANOVA on 4 treatments with unequal sample sizes (n=9,7,10 and 10). I want to determine where my sig. differences are between treatments using a Bonferroni test, and have run the code: pairwise.t.test(Wk16, Treatment, p.adf=bonf) I receive an error message stating that my arguments are of unequal length: Error in tapply(x, g, mean, na.rm = TRUE) : arguments must have same length Is there a way to run this test even with unequal sample sizes? Unequal sample sizes are not causing the reported error. The problem is that 'Wk16' and 'Treatment' do not have the same number of observations. If the group sizes are 9, 7, 10, and 10, then 'Wk16' and 'Treatment' should each contain 36 observations. Thanks. Stephanie Reiner Andrews Hall 410 Shipping: Rt. 1208 Greate Road Gloucester Point, VA 23062 work: 804-684-7869 cell: 443-286-4795 __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 do you know which functions are being debugged?
list all objects first; use a loop (explicitly or not) to check whether (1) your objects are functions (2) functions are debugged f = function(x) x g = 1 x = ls() debug(f) sapply(x[sapply(x, function(i) is.function(get(i)))], isdebugged) f TRUE undebug(f) sapply(x[sapply(x, function(i) is.function(get(i)))], isdebugged) f FALSE Regards, Yihui -- Yihui Xie xieyi...@gmail.com Phone: 515-294-6609 Web: http://yihui.name Department of Statistics, Iowa State University 3211 Snedecor Hall, Ames, IA On Fri, Oct 23, 2009 at 12:28 PM, Andrew Yee y...@post.harvard.edu wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? Thanks, Andrew [[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] how do you know which functions are being debugged?
On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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.
Re: [R] how do you know which functions are being debugged?
Oops... I forgot to mention that 'envir' (or 'pos') should be specified in ls()/get() in my last reply if you are looking for debugged functions in environments other than .GlobalEnv. Regards, Yihui -- Yihui Xie xieyi...@gmail.com Phone: 515-294-6609 Web: http://yihui.name Department of Statistics, Iowa State University 3211 Snedecor Hall, Ames, IA On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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] how do you know which functions are being debugged?
I have often wanted to get such a list as well without having to write code myself. If a function which automatically iterated over all possible objects and listed out the ones that being debugged were available in the core it would be useful. Even if its slow this would be used at debug time and not a production time so it would not be too bad. Even a heuristic that checked likely locations, e.g. just the global environment, but not all locations might be useful. On Fri, Oct 23, 2009 at 2:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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] Bayesian regression stepwise function?
The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's not stepwise. Ntzoufras shows how to use WinBUGS to conduct Bayesian model selection, but again it's not stepwise Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of Statistical Software 7(7), 1--19. Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ. Raftery, A. E. (1995), 'Bayesian model selection in social research', Sociological Methodology 25, 111-163. Allan.Y all...@cmu.edu Sent by: r-help-boun...@r-project.org 10/22/2009 01:09 PM To r-help@r-project.org cc Subject [R] Bayesian regression stepwise function? Hi everyone, I am wondering if there exists a stepwise regression function for the Bayesian regression model. I tried googling, but I couldn't find anything. I know step function exists for regular stepwise regression, but nothing for Bayes. Thanks -- View this message in context: http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p26013725.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Prediction Error Calculation
Hello List, I am fitting a logistic regression model for some presence/absence type data. I have numerous covariates I am fitting to explain variation, and I am using AIC to rank models. However, I would like to report how well my best model (s) do at prediction. I have looked over the archives and the web and have come up with something that gives me what I think is the mean prediction error, BUT I am not sure of that. I am sort of unfamiliar with these types of statistics. Here is my code: metrics.global-glm(Type~MPI+IJI+ED+PRD+class2+class3+class5, family=binomial, data=metrics)## ##Type is my binary response 0 or 1 muhat-metrics.global$fitted.values ##assigns the fitted values a name muhat global.diag-glm.diag(metrics.global) ##creates a the diagnostic values cv.err-mean((metrics.global$y-muhat)^2/(1-global.diag$h)^2) ###calculates cv.err cv.err My main problem is I am unsure how to interpret what cv.err means for my model. I know that h is a leverage statistic for each observation. I would appreciate some interpretation clarification. Thank you. -- View this message in context: http://www.nabble.com/Prediction-Error-Calculation-tp26031236p26031236.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 do you know which functions are being debugged?
On 10/23/2009 2:27 PM, Gabor Grothendieck wrote: I have often wanted to get such a list as well without having to write code myself. If a function which automatically iterated over all possible objects and listed out the ones that being debugged were available in the core it would be useful. Even if its slow this would be used at debug time and not a production time so it would not be too bad. Even a heuristic that checked likely locations, e.g. just the global environment, but not all locations might be useful. If you look at the code in setBreakpoint, it does a search for functions in lots of places (but not everywhere). You could start from that. Generally speaking, it's quite hard to think of a way to express what the answer would look like to a perfectly general function like you're asking for. Suppose a function lives in the parent environment of the environment of an expression that was returned in the result of a call to lm that is a local variable in the function that called the function where you are asking the question: how would you return that as an answer??? In a language with pointers, you'd just return a pointer to it, but R doesn't have those. Duncan Murdoch On Fri, Oct 23, 2009 at 2:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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.
[R] help using R's linprog for LP
Hi, I am using R in one of my courses. I am trying to use R's linprog package to solve to formulate 2-class classification problem as Linear programming problem. For my formulation, I need to set to cvec to all 0s. I know the points are linearly separable so an optimal solution x does exist, which satisfies all the constraints. But given the constraints and setting cvec to all 0s it simply gives me an x of all 0s. How can I fix this? Medha __ 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 do you know which functions are being debugged?
On Fri, Oct 23, 2009 at 2:34 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 2:27 PM, Gabor Grothendieck wrote: I have often wanted to get such a list as well without having to write code myself. If a function which automatically iterated over all possible objects and listed out the ones that being debugged were available in the core it would be useful. Even if its slow this would be used at debug time and not a production time so it would not be too bad. Even a heuristic that checked likely locations, e.g. just the global environment, but not all locations might be useful. If you look at the code in setBreakpoint, it does a search for functions in lots of places (but not everywhere). You could start from that. Generally speaking, it's quite hard to think of a way to express what the answer would look like to a perfectly general function like you're asking for. Suppose a function lives in the parent environment of the environment of an expression that was returned in the result of a call to lm that is a local variable in the function that called the function where you are asking the question: how would you return that as an answer??? In a language with pointers, you'd just return a pointer to it, but R doesn't have those. I agree its problematic but it would still be useful since after you have debugged a number of functions its easy to forget which ones are being debugged. __ 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 do you know which functions are being debugged?
Thanks everyone for the responses. As a suggestion, wouldn't it be useful, that when you run isdebugged() without any parameters, it would just report back the functions that are being flagged for debugging? Andrew On Fri, Oct 23, 2009 at 2:12 PM, Yihui Xie xieyi...@gmail.com wrote: Oops... I forgot to mention that 'envir' (or 'pos') should be specified in ls()/get() in my last reply if you are looking for debugged functions in environments other than .GlobalEnv. Regards, Yihui -- Yihui Xie xieyi...@gmail.com Phone: 515-294-6609 Web: http://yihui.name Department of Statistics, Iowa State University 3211 Snedecor Hall, Ames, IA On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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. [[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 do you know which functions are being debugged?
On 10/23/2009 2:49 PM, Andrew Yee wrote: Thanks everyone for the responses. As a suggestion, wouldn't it be useful, that when you run isdebugged() without any parameters, it would just report back the functions that are being flagged for debugging? Could you give an example of what the answer would look like? Not all functions have names, not all names are unique, ... Duncan Murdoch Andrew On Fri, Oct 23, 2009 at 2:12 PM, Yihui Xie xieyi...@gmail.com wrote: Oops... I forgot to mention that 'envir' (or 'pos') should be specified in ls()/get() in my last reply if you are looking for debugged functions in environments other than .GlobalEnv. Regards, Yihui -- Yihui Xie xieyi...@gmail.com Phone: 515-294-6609 Web: http://yihui.name Department of Statistics, Iowa State University 3211 Snedecor Hall, Ames, IA On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 1:28 PM, Andrew Yee wrote: This is kind of a dumb question: I know you can use isdebugged() to find out if a specific function is flagged for debugging, but is there a way to list all the functions that are flagged for debugging? No, R doesn't keep any master list, it sets a flag in each one. So you could iterate over all visible objects to find the ones that have the flag set (and this is quite slow, there are a lot of places to look), but there would always be the chance that one was hiding somewhere you didn't look. 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.
[R] Data format for KSVM
Hi, I have a process using svm from the e1071 library. it works. I want to try using the KSVM library instead. The same data used wiht e1071 gives me an error with KSVM. My data is a data.frame. sample code: svm_formula - formula(y ~ a + B + C) svm_model - ksvm(formula, data=train_data, type=C-svc, kernel=rbfdot, C=1) I get the following error: object is not a matrix So I tried this: svm_model - ksvm(formula, data=as.matrix(train_data), type=C-svc, kernel=rbfdot, C=1, scaled=FALSE) Now I get this error: Error in model.fram.definition(data = list(v1 = c(1.1234, -2.3232: Object is not a matrix My data was previously scaled with the scale() function so that the mean is centered at 0. and the range is {-1,1} Can anyone provide some suggestions as to why I'm getting an error? Thanks! -N __ 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] Bayesian regression stepwise function?
The stepAIC() function in MASS can be used, with k = log(n), to implement your suggestion of quasi-Bayesian stepwise selection using the BIC criterion. Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of jlu...@ria.buffalo.edu Sent: Friday, October 23, 2009 2:31 PM To: Allan.Y Cc: r-help@r-project.org; r-help-boun...@r-project.org Subject: Re: [R] Bayesian regression stepwise function? The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's not stepwise. Ntzoufras shows how to use WinBUGS to conduct Bayesian model selection, but again it's not stepwise Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of Statistical Software 7(7), 1--19. Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ. Raftery, A. E. (1995), 'Bayesian model selection in social research', Sociological Methodology 25, 111-163. Allan.Y all...@cmu.edu Sent by: r-help-boun...@r-project.org 10/22/2009 01:09 PM To r-help@r-project.org cc Subject [R] Bayesian regression stepwise function? Hi everyone, I am wondering if there exists a stepwise regression function for the Bayesian regression model. I tried googling, but I couldn't find anything. I know step function exists for regular stepwise regression, but nothing for Bayes. Thanks -- View this message in context: http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p2601 3725.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-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 do you know which functions are being debugged?
Hi, On Oct 23, 2009, at 2:39 PM, Gabor Grothendieck wrote: On Fri, Oct 23, 2009 at 2:34 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 2:27 PM, Gabor Grothendieck wrote: I have often wanted to get such a list as well without having to write code myself. If a function which automatically iterated over all possible objects and listed out the ones that being debugged were available in the core it would be useful. Even if its slow this would be used at debug time and not a production time so it would not be too bad. Even a heuristic that checked likely locations, e.g. just the global environment, but not all locations might be useful. If you look at the code in setBreakpoint, it does a search for functions in lots of places (but not everywhere). You could start from that. Generally speaking, it's quite hard to think of a way to express what the answer would look like to a perfectly general function like you're asking for. Suppose a function lives in the parent environment of the environment of an expression that was returned in the result of a call to lm that is a local variable in the function that called the function where you are asking the question: how would you return that as an answer??? In a language with pointers, you'd just return a pointer to it, but R doesn't have those. I agree its problematic but it would still be useful since after you have debugged a number of functions its easy to forget which ones are being debugged. If this is something you just want to use on your own functions that are set to `debug`, why not just make write a debug wraps the base::debug and does some keeps track of which functions you're debugging in some global environment variable (of your creation). You should probably use something beside substitute here, but: R .DLIST - character() R debug - function(fun, text=, condition=NULL) { .DLIST - c(.DLIST, substitute(fun)) base::debug(fun, text=text, condition=condition) } R myfun - function(something) { a - length(something) b - sum(something) b } R debug(myfun) R debug(myfun) R .DLIST [[1]] myfun R myfun(1:10) debugging in: myfun(1:10) debug: { a - length(something) b - sum(something) b } Browse[2] n debug: a - length(something) Browse[2] n debug: b - sum(something) Browse[2] n debug: b Browse[2] n exiting from: myfun(1:10) [1] 55 You can define your own undebug function accordingly. Would that do what you want? -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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] reference for mahalanobis {stats}
On Fri, Oct 23, 2009 at 10:59 AM, Steve Lianoglou mailinglist.honey...@gmail.com wrote: Hi, On Oct 23, 2009, at 11:36 AM, Peng Yu wrote: The help on mahalanobis {stats} does not include any reference. I'm interested in understand why Mahalanobis is defined in its current way and how to use it. Could somebody point me a good book on this? I have looked through a few books, but they all give very light explanation on it. The wikipedia page looks pretty good ... the Intuitive Explanation section should help: http://en.wikipedia.org/wiki/Mahalanobis_distance Googling also provides many more explanations: http://matlabdatamining.blogspot.com/2006/11/mahalanobis-distance.html I have read the above links. But I still feel that a book might explain it in a more complete fashion. __ 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 do you know which functions are being debugged?
On Fri, Oct 23, 2009 at 2:52 PM, Duncan Murdoch murd...@stats.uwo.ca wrote: On 10/23/2009 2:49 PM, Andrew Yee wrote: Thanks everyone for the responses. As a suggestion, wouldn't it be useful, that when you run isdebugged() without any parameters, it would just report back the functions that are being flagged for debugging? Could you give an example of what the answer would look like? Not all functions have names, not all names are unique, ... One could use the formal argument list of the function, i.e. the output of args, when there is no notion of name since that is always available. __ 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 using R's linprog for LP
Hi, I found the reason. By default it puts a condition for x = 0. Is there a way to get rid of this condition? Medha On Fri, Oct 23, 2009 at 2:34 PM, Medha Atre medha.a...@gmail.com wrote: Hi, I am using R in one of my courses. I am trying to use R's linprog package to solve to formulate 2-class classification problem as Linear programming problem. For my formulation, I need to set to cvec to all 0s. I know the points are linearly separable so an optimal solution x does exist, which satisfies all the constraints. But given the constraints and setting cvec to all 0s it simply gives me an x of all 0s. How can I fix this? Medha __ 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 apply the Wilcoxon test to a hole table at once?
Hi, I have a data set: Dataset X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 1 user1 m 22 19 28 24 12 18 9 7 4 5 4 7 5 7 9 2 user2 f 25 19 23 18 18 15 6 8 6 6 7 10 7 7 7 3 user3 f 28 21 24 18 15 12 10 6 7 9 5 10 5 9 5 4 user4 f 26 19 26 21 12 18 6 6 5 1 3 8 6 5 6 5 user5 m 21 22 26 18 9 6 4 6 1 7 2 4 4 6 4 6 user6 m 24 8 25 12 18 12 7 8 4 1 4 6 7 5 6 ... 71 user71 m 18 4 10 6 3 6 9 5 10 8 4 5 6 5 5 I can apply the Wilcoxon test on each column one by one, but how to do this on the hole table at once? wilcox.test(X3 ~ X2, alternative=two.sided, data=Dataset) Wilcoxon rank sum test with continuity correction data: X3 by X2 W = 439, p-value = 0.1291 alternative hypothesis: true location shift is not equal to 0 I researched on this, but I can't find a solution. I would really appreciate any help. P.S. Excuse my lack of terminology :). Iurie Malai Moldova Pedagogical State University -- View this message in context: http://www.nabble.com/How-to-apply-the-Wilcoxon-test-to-a-hole-table-at-once--tp26030572p26030572.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] Reporting discriminat analysis
Dear all, I want to report the results of discriminant analysis in a way that would enable the readers to use my results (functions) for classification of their own cases. For this it should suffice to provide discriminant functions along with the cut-off values for distinguishing the groups. However, I cannot find these cut-off values. I use the function lda of the MASS package, although I have also looked at some other procedures for discriminant analysis, and equally failed. Thanks, Stats Wolf __ 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] Bayesian regression stepwise function?
Ravi Varadhan wrote: The stepAIC() function in MASS can be used, with k = log(n), to implement your suggestion of quasi-Bayesian stepwise selection using the BIC criterion. Ravi. Although many statisticians use BIC otherwise, it was only designed to compare two pre-specified models. Frank --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of jlu...@ria.buffalo.edu Sent: Friday, October 23, 2009 2:31 PM To: Allan.Y Cc: r-help@r-project.org; r-help-boun...@r-project.org Subject: Re: [R] Bayesian regression stepwise function? The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's not stepwise. Ntzoufras shows how to use WinBUGS to conduct Bayesian model selection, but again it's not stepwise Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of Statistical Software 7(7), 1--19. Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ. Raftery, A. E. (1995), 'Bayesian model selection in social research', Sociological Methodology 25, 111-163. Allan.Y all...@cmu.edu Sent by: r-help-boun...@r-project.org 10/22/2009 01:09 PM To r-help@r-project.org cc Subject [R] Bayesian regression stepwise function? Hi everyone, I am wondering if there exists a stepwise regression function for the Bayesian regression model. I tried googling, but I couldn't find anything. I know step function exists for regular stepwise regression, but nothing for Bayes. Thanks -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ 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.