Re: [R] Negative exponential fit
On 30.11.2011 07:25, Indrajit Sengupta wrote: Agreed, nobody has deemed it default, but is there any other such package that you can think of for this purpose. Not that I know. Suggesting it is perfect. I'd just not call any package a default if it does not ship with R. Best, Uwe Regards, From: Uwe Liggeslig...@statistik.tu-dortmund.de To: Indrajit Senguptaindra_cali...@yahoo.com Cc: r-help@r-project.orgr-help@r-project.org; rch4r...@geneseo.edu Sent: Tuesday, November 29, 2011 7:33 PM Subject: Re: [R] Negative exponential fit On 29.11.2011 07:06, Indrajit Sengupta wrote: What have you tried so far - can you explain? fitdistrplus package is the default package for fitting distributions. It is a contributed packages, and perhaps it is a good one (I do not know), but calling it the *default* ... who defined that? Best, Uwe Ligges Regards, Indrajit From: rch4r...@geneseo.edu To: r-help@r-project.org Sent: Tuesday, November 29, 2011 8:39 AM Subject: [R] Negative exponential fit We need help We are doing a project for a statistical class in and we are looking at world record times in different running events over time. We are trying to fit the data with a negative exponential but we just cant seem to get a function that works properly. we have on our x-axis the date and on the y-axis the time(in seconds). So as you can imagine, the times have decreased and appear to be approaching a limit. Any ideas for a nls function that would work for us would be greatly appreciated. Rob -- View this message in context: http://r.789695.n4.nabble.com/Negative-exponential-fit-tp4117889p4117889.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] Non parametric, repeated-measures, factorial ANOVA
On 11/29/2011 10:31 PM, RobH2011 wrote: Hi I have data from an experiment that used a repeated-measures factorial 2x2 design (i.e. each participant contributed data to both levels of both factors). I need a non-parametric version of the repeated-measures factorial ANOVA to analyse the data. SPSS only has non-parametric tests for one-way ANOVAs but I have been told that the test I need can be implemented using the R software. Unfortunantly I haven't been able to work out how to do this test in R. I've got the Wilcox book 'Robust Estimation and Hypothesis testing' but to be honest I don't have enough of a background in stats to understand a great deal of it, and I wasn't able to work out how to do the test I want. Can anyone let me know how I can do such a test in R? I presume this is a test that people frequently need to use and that there will be some code already written that can implement the test. Hi Rob, The nparLD package, based on the work of Professor Brunner, might do what you want. 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] Parameters setting in functions optimization
Le 11/30/2011 2:09 AM, Florent D. a écrit : Thanks for your answer ! I also think your last write-up for LogLiketot (using a single argument par) is the correct approach if you want to feed it to optim(). I'm not dedicated to optim() fonction. I just want to optimise my two parameters and the loglikelihood result, and if there's a better fonction for that, I wish I could use it. So now you have a problem with log(LikeGi(l, i, par[1], par[2])) for some values of par[1] and par[2]. Where is LikeGi coming from? a package or is it your own function? My own function, otherwise it would be simplier to discuss about my problems. You could add some print statements (if you are familiar with browser() it is even better) so you may see what values of par are causing trouble. I'm not familiar, but I'll search about browser(). If the function with par is correct, any idea of what I've made with this : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) ? On Tue, Nov 29, 2011 at 1:15 PM, Diane Bailleul diane.baill...@u-psud.fr wrote: Good afternoon everybody, I'm quite new in functions optimization on R and, whereas I've read lot's of function descriptions, I'm not sure of the correct settings for function like optimx and nlminb. I'd like to minimize my parameters and the loglikelihood result of the function. My parameters are a mean distance of dispersion and a proportion of individuals not assigned, coming from very far away. The function LikeGi reads external tables and it's working as I want (I've got a similar model on Mathematica). My final function is LogLiketot : LogLiketot- function(dist,ms) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } dist is the mean dispersal distance (0, lots of meters) and ms the proportion of individuals (0-1). Of course, I want them to be as low as possible. I'd tried to enter the initials parameters as indicated in the tutorials : optim(c(40,0.5), fn=LogLiketot) Error in 1 - ms : 'ms' is missing But ms is 0.5 ... So I've tried this form : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) with different values for the two parameters : par fvalues method fns grs itns conv KKT1 KKT2 xtimes 219.27583, 25.37964 2249.698BFGS 12 8 NULL0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL0 TRUE TRUE 66.3 The first line is not possible but as I've not constrained the optimization ... but the second line would be a very good result ! Then, searching for another similar cases, I've tried to change my function form: LogLiketot- function(par) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,par[1],par[2]))) } } return(-sum(res)) } where dist=par[1] and ms=par[2] And I've got : optimx(c(40,0.5), fn=LogLiketot) par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083BFGS 29 10 NULL0 TRUE NA 92.03 1 39.7372199, 0.9778101 1064.083 Nelder-Mead 53 NA NULL0 TRUE NA 70.83 And I've got now a warning message : In log(LikeGi(l, i, par[1], par[2])) : NaNs produced (which are very bad results in that case) Anyone with previous experiences in optimization of several parameters could indicate me the right way to enter the initial parameters in this kind of functions ? Thanks a lot for helping me ! Diane -- Diane Bailleul Doctorante Université Paris-Sud 11 - Faculté des Sciences d'Orsay Unité Ecologie, Systématique et Evolution Département Biodiversité, Systématique et Evolution UMR 8079 - UPS CNRS AgroParisTech Porte 320, premier étage, Bâtiment 360 91405 ORSAY CEDEX FRANCE (0033) 01.69.15.56.64 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Diane Bailleul Doctorante Université Paris-Sud 11 - Faculté des Sciences d'Orsay Unité Ecologie, Systématique et Evolution Département Biodiversité, Systématique et Evolution UMR 8079 - UPS CNRS AgroParisTech Porte 320, premier étage, Bâtiment 360 91405 ORSAY CEDEX FRANCE (0033) 01.69.15.56.64 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Storing the linear model object
Good morning, Normally, if you save your R environment, some objects, like your lm, must be saved in with. Otherwise, have you tried save() or dput() functions ? Le 11/30/2011 6:10 AM, arunkumar a écrit : Hi Please let me know if we can store the linear model object in the data base and retrive the object and output from them Data- read.csv(C:/FE and RE.csv) Formula=Y~X2+X3+X4 lmobject = lm(formula=Formula,data=Data) can i store the lm object in the database and and is it possible to retrive it and get the summary information -- View this message in context: http://r.789695.n4.nabble.com/Storing-the-linear-model-object-tp4121944p4121944.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. -- Diane Bailleul Doctorante Université Paris-Sud 11 - Faculté des Sciences d'Orsay Unité Ecologie, Systématique et Evolution Département Biodiversité, Systématique et Evolution UMR 8079 - UPS CNRS AgroParisTech Porte 320, premier étage, Bâtiment 360 91405 ORSAY CEDEX FRANCE (0033) 01.69.15.56.64 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Replacing a value in a dataframe
hi I have data like this in a dataframe Var Value Cheque X140FALSE X220FALSE X328TRUE I want to replace it FLASE with 0 and TRUE with 1. is there any method by which i can do without using LOOP -- View this message in context: http://r.789695.n4.nabble.com/Replacing-a-value-in-a-dataframe-tp4122188p4122188.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] Hiding message
On 30/11/11 20:22, Jeff Newmiller wrote: Don't do that. Use $ notation to refer to elements of lists/data frames explicitly. Or use with(). It is truly-ruly handy and effective. cheers, Rolf Turner --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.us Basics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Aldomichael.v.claw...@gmail.com wrote: It is from attatching and detaching data frames and moving them between functions, there is no library calls in the loop Michael Weylandt wrote What does your code look like that you see such messages repeatedly? I only see those message upon loading a package: maybe move the library() call outside the loop. Also, the suppressPackageStartupMessages() function can hide them, but you should first figure out for yourself why they are there. Michael On Tue, Nov 29, 2011 at 9:39 PM, Aldolt;michael.v.clawson@gt; wrote: So I am receiving messages in a loop The following object(s) are masked from... I tried triedoptions(warn=-1) and options(show.error.message:=FALSE) but it still tells me The following object(s) are masked from... repeatedly, and my problem is that it masks the print statement from my function Any ideas? Thanks -- View this message in context: http://r.789695.n4.nabble.com/Hiding-message-tp4121625p4121625.html Sent from the R help mailing list archive at Nabble.com. __ R-help@ mailing list https://stat.ethz.ch/mailman/listinfo/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@ mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://r.789695.n4.nabble.com/Hiding-message-tp4121625p4121905.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Replacing a value in a dataframe
On 30.11.2011 09:16, arunkumar wrote: hi I have data like this in a dataframe Var Value Cheque X140FALSE X220FALSE X328TRUE I want to replace it FLASE with 0 and TRUE with 1. is there any method by which i can do without using LOOP dataframe$Cheque - as.integer(Cheque) Uwe Ligges -- View this message in context: http://r.789695.n4.nabble.com/Replacing-a-value-in-a-dataframe-tp4122188p4122188.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Invalid number of components, ncomp
Error in mvr(Kd_nM ~ qsar, ncomp = 6, data = my, validation = CV, method = kernelpls) : Invalid number of components, ncomp How I can fix this? [[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] Generalized singular value decomposition
Hello, I would like to perform a generalized singular value decomposition with R. The only possibility I found is GSVD that is based on LAPACK/BLAS. Are there other possibilities too? If not, has anybody used LAPACK/BLAS under Windows XP? How can I install them? Following [1] did not help. I hope this is the right place for my question. Thank you very much! Oana Tomescu [1] http://sites.google.com/site/jivsoft/Home/r-blas-interface __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Parameters setting in functions optimization
With optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) where LogLiketot- function(dist,ms) { res - NULL for(i in 1:nrow(pop5)) { for(l in 1:nrow(freqvar)) { res - c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } I think it will do something like this for the first call to LogLiketot: LogLiketot(c(30,50), ms=c(0.4,0.5)) which is obviously not the usage you had in mind. Also, I see you say the results for the bad usage above: par fvalues method fns grs itns conv KKT1 KKT2 xtimes 219.27583, 25.37964 2249.698BFGS 12 8 NULL0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL0 TRUE TRUE 66.3 look very good but you do not comment about the results for the correct usage of optimx: par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083BFGS 29 10 NULL0 TRUE NA 92.03 1 39.7372199, 0.9778101 1064.083 Nelder-Mead 53 NA NULL0 TRUE NA 70.83 Do you realize optimx is trying to _minimize_ your function? See that the fvalues from the correct usage are much better (smaller) than you first (bad) usage. On Wed, Nov 30, 2011 at 4:16 AM, Diane Bailleul diane.baill...@u-psud.fr wrote: Le 11/30/2011 2:09 AM, Florent D. a écrit : Thanks for your answer ! I also think your last write-up for LogLiketot (using a single argument par) is the correct approach if you want to feed it to optim(). I'm not dedicated to optim() fonction. I just want to optimise my two parameters and the loglikelihood result, and if there's a better fonction for that, I wish I could use it. So now you have a problem with log(LikeGi(l, i, par[1], par[2])) for some values of par[1] and par[2]. Where is LikeGi coming from? a package or is it your own function? My own function, otherwise it would be simplier to discuss about my problems. You could add some print statements (if you are familiar with browser() it is even better) so you may see what values of par are causing trouble. I'm not familiar, but I'll search about browser(). If the function with par is correct, any idea of what I've made with this : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) ? On Tue, Nov 29, 2011 at 1:15 PM, Diane Bailleul diane.baill...@u-psud.fr wrote: Good afternoon everybody, I'm quite new in functions optimization on R and, whereas I've read lot's of function descriptions, I'm not sure of the correct settings for function like optimx and nlminb. I'd like to minimize my parameters and the loglikelihood result of the function. My parameters are a mean distance of dispersion and a proportion of individuals not assigned, coming from very far away. The function LikeGi reads external tables and it's working as I want (I've got a similar model on Mathematica). My final function is LogLiketot : LogLiketot- function(dist,ms) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } dist is the mean dispersal distance (0, lots of meters) and ms the proportion of individuals (0-1). Of course, I want them to be as low as possible. I'd tried to enter the initials parameters as indicated in the tutorials : optim(c(40,0.5), fn=LogLiketot) Error in 1 - ms : 'ms' is missing But ms is 0.5 ... So I've tried this form : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) with different values for the two parameters : par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 19.27583, 25.37964 2249.698 BFGS 12 8 NULL 0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL 0 TRUE TRUE 66.3 The first line is not possible but as I've not constrained the optimization ... but the second line would be a very good result ! Then, searching for another similar cases, I've tried to change my function form: LogLiketot- function(par) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,par[1],par[2]))) } } return(-sum(res)) } where dist=par[1] and ms=par[2] And I've got : optimx(c(40,0.5), fn=LogLiketot) par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083 BFGS 29 10 NULL 0 TRUE NA 92.03 1 39.7372199, 0.9778101 1064.083 Nelder-Mead 53 NA NULL 0 TRUE NA 70.83 And I've got now a warning message : In log(LikeGi(l, i, par[1], par[2])) : NaNs produced (which are very bad results in that case) Anyone with previous experiences in optimization of several parameters could indicate me the right way to enter the initial parameters in this kind of functions ? Thanks a lot for helping me ! Diane -- Diane Bailleul Doctorante Université Paris-Sud 11 - Faculté des
Re: [R] Replacing a value in a dataframe
er... Uwe, shouldn't that be, e.g. dataframe$Cheque - as.integer(dataframe$Cheque) ## or building on Rolf's suggestion dataframe - within(dataframe, Cheque - as.integer(Cheque)) While I am at it, is there any practical difference in efficiency between these two approaches? -- Bert 2011/11/30 Uwe Ligges lig...@statistik.tu-dortmund.de: On 30.11.2011 09:16, arunkumar wrote: hi I have data like this in a dataframe Var Value Cheque X1 40 FALSE X2 20 FALSE X3 28 TRUE I want to replace it FLASE with 0 and TRUE with 1. is there any method by which i can do without using LOOP dataframe$Cheque - as.integer(Cheque) Uwe Ligges -- View this message in context: http://r.789695.n4.nabble.com/Replacing-a-value-in-a-dataframe-tp4122188p4122188.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Negative exponential fit
Am 29.11.2011 07:06, schrieb Indrajit Sengupta: What have you tried so far - can you explain? fitdistrplus package is the default package for fitting distributions. Regards, Indrajit From: rch4 r...@geneseo.edu To: r-help@r-project.org Sent: Tuesday, November 29, 2011 8:39 AM Subject: [R] Negative exponential fit We need help We are doing a project for a statistical class in and we are looking at world record times in different running events over time. We are trying to fit the data with a negative exponential but we just cant seem to get a function that works properly. we have on our x-axis the date and on the y-axis the time(in seconds). So as you can imagine, the times have decreased and appear to be approaching a limit. Any ideas for a nls function that would work for us would be greatly appreciated. Rob -- View this message in context: http://r.789695.n4.nabble.com/Negative-exponential-fit-tp4117889p4117889.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]] I do not want to shed out any doubt as to the merits of pkg fitdistrplus, in particular for censored data, but seconding Uwe's reply later in this thread, you may also want to check out pkg distrMod on CRAN --- M. Kohl, P. Ruckdeschel (2010): R Package distrMod: S4 Classes and Methods for Probability Models. Journal of Statistical Software, 35(10), 1-27. URL http://www.jstatsoft.org/v35/i10/. which offers quite some additional flexibility for model fitting---including new models (built on distributions which have no name but instead are image distributions under arithmetic transformations of existing ones, see example M2 in the cited ref) and (nonlinear) transformations of the parameter (see Example p.15, cited ref). Best regards, Peter __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] SPSS - R
At 16:41 27/11/2011, Kristi Shoemaker wrote: Hi John, Your assumptions are correct and those examples were very helpful, thanks! I think I'm almost there, but I'm screwing up something with the within-subjects factor (the example has two, I only have one). See below. Why am I not seeing my within-subjects factors in the ANOVA report? I am by no means an expert in this but my reading of the documentation suggests that the idata and idesign parameters are for multivariate models and that is not what you specified. Dat2 Subj Age Sex Time HippV 1 s01 32.9 1 0w 6.50098 2 s01 32.9 1 6w 6.91793 3 s02 35.1 0 0w 7.32480 4 s02 35.1 0 6w 7.56012 5 s03 34.4 0 0w 6.51385 6 s03 34.4 0 6w 6.56875 9 s05 39.9 1 0w 6.92855 10 s05 39.9 1 6w 6.94926 11 s06 29.5 1 0w 6.99383 12 s06 29.5 1 6w 7.10568 13 s07 45.9 1 0w 6.94380 14 s07 45.9 1 6w 7.08190 15 s08 20.3 1 0w 7.76881 16 s08 20.3 1 6w 7.72725 17 s09 26.9 0 0w 5.37566 18 s09 26.9 0 6w 5.74887 21 s11 22.0 0 0w 7.12992 22 s11 22.0 0 6w 7.16237 23 s12 31.0 1 0w 6.70629 24 s12 31.0 1 6w 6.80872 25 s13 50.1 1 0w 7.22649 26 s13 50.1 1 6w 7.58900 27 s14 22.2 0 0w 5.97577 28 s14 22.2 0 6w 5.80801 29 s16 20.4 1 0w 7.99554 30 s16 20.4 1 6w 8.09260 31 s20 24.0 0 0w 7.01014 32 s20 24.0 0 6w 6.87821 33 s21 32.4 0 0w 5.90883 34 s21 32.4 0 6w 5.95392 37 s24 22.2 0 0w 6.15474 38 s24 22.2 0 6w 6.00906 41 s27 22.2 1 0w 7.88765 42 s27 22.2 1 6w 7.76038 49 s35 24.3 0 0w 6.05998 50 s35 24.3 0 6w 6.07399 51 s36 23.5 0 0w 7.83182 52 s36 23.5 0 6w 7.60268 53 s38 59.7 1 0w 7.39672 54 s38 59.7 1 6w 6.98291 55 s39 40.5 0 0w 7.31330 56 s39 40.5 0 6w 7.50559 57 s40 24.2 1 0w 8.54958 58 s40 24.2 1 6w 8.65016 59 s41 23.6 1 0w 7.76049 60 s41 23.6 1 6w 7.58946 61 s42 53.3 0 0w 7.03388 62 s42 53.3 0 6w 7.48384 63 s43 34.4 0 0w 6.86967 64 s43 34.4 0 6w 6.81076 65 s44 44.8 0 0w 7.33779 66 s44 44.8 0 6w 7.86175 67 s45 40.2 1 0w 6.55963 68 s45 40.2 1 6w 6.52577 71 s47 26.5 0 0w 7.09418 72 s47 26.5 0 6w 6.92850 73 s48 22.7 0 0w 6.77078 74 s48 22.7 0 6w 6.67289 75 s50 36.1 1 0w 7.47208 76 s50 36.1 1 6w 7.55876 library(car) Loading required package: MASS Loading required package: nnet Loading required package: survival Loading required package: splines timepoints - factor(c(0w, 6w), + levels=c(0w,6w)) idata=data.frame(timepoints) mod.ok - lm(HippV ~ Age*Sex,data=Dat2) (av.ok - Anova(mod.ok, idata=idata, idesign=~timepoints)) Anova Table (Type II tests) Response: HippV Sum Sq Df F valuePr(F) Age0.0234 1 0.071 0.7909202 Sex5.3082 1 16.128 0.0001781 *** Age:Sex4.7772 1 14.514 0.0003478 *** Residuals 18.4314 56 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 From: John Fox j...@mcmaster.ca Cc: r-help@r-project.org Sent: Saturday, November 26, 2011 7:19 PM Subject: RE: [R] SPSS - R Dear Kristi, I assume that this is a repeated-measures ANOVA with one within-subjects factor (Time) and two between-subjects factors (Age and Sex, which are crossed). If Age is numeric, and not a factor, then the type-III tests that you requested don't test sensible hypotheses. In any event, if my guess is right about the design, then you can use the Anova() function in the car package for an equivalent analysis. See the repeated-measures example in ?Anova (for the O'Brien and Kaiser data). You've already had an answer to the more general question. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Kristi Shoemaker Sent: November-26-11 11:08 AM To: r-help@r-project.org Subject: [R] SPSS - R I'm an SPSS user trying to make the transition to R. Can someone help me translate the following SPSS code into R?: GLM Total_tp1 Total_tp2 WITH Age Sex /WSFACTOR=Time 2 Repeated /METHOD=SSTYPE(3) /CRITERIA=ALPHA(.05) /WSDESIGN= Time /DESIGN= Age Sex Age*Sex. Also. can anyone recommend any resources to help SPSS users learn to things in R? Thanks, -kristi [[alternative HTML version deleted]] [[alternative HTML version deleted]] Michael Dewey i...@aghmed.fsnet.co.uk http://www.aghmed.fsnet.co.uk/home.html __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Replacing a value in a dataframe
On 30.11.2011 12:26, Bert Gunter wrote: er... Uwe, shouldn't that be, e.g. dataframe$Cheque- as.integer(dataframe$Cheque) Sure, thanks. ## or building on Rolf's suggestion dataframe- within(dataframe, Cheque- as.integer(Cheque)) While I am at it, is there any practical difference in efficiency between these two approaches? Well, just profile it. The latter has some overhead, of course: d - data.frame(a=c(TRUE, FALSE)) system.time(for(i in 1:1e4) {d - data.frame(a=c(TRUE, FALSE)); d$a - as.integer(d$a)}) system.time(for(i in 1:1e4) {d - data.frame(a=c(TRUE, FALSE)); d - within(d, a - as.integer(a))}) Uwe Ligges -- Bert 2011/11/30 Uwe Liggeslig...@statistik.tu-dortmund.de: On 30.11.2011 09:16, arunkumar wrote: hi I have data like this in a dataframe Var Value Cheque X140FALSE X220FALSE X328TRUE I want to replace it FLASE with 0 and TRUE with 1. is there any method by which i can do without using LOOP dataframe$Cheque- as.integer(Cheque) Uwe Ligges -- View this message in context: http://r.789695.n4.nabble.com/Replacing-a-value-in-a-dataframe-tp4122188p4122188.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Storing the linear model object
If you mean storing it within one session, it's just like any other object m - lm(...) and it can be put in lists, modified, whatever just like any other object. Michael On Nov 30, 2011, at 4:23 AM, Diane Bailleul diane.baill...@u-psud.fr wrote: Good morning, Normally, if you save your R environment, some objects, like your lm, must be saved in with. Otherwise, have you tried save() or dput() functions ? Le 11/30/2011 6:10 AM, arunkumar a écrit : Hi Please let me know if we can store the linear model object in the data base and retrive the object and output from them Data- read.csv(C:/FE and RE.csv) Formula=Y~X2+X3+X4 lmobject = lm(formula=Formula,data=Data) can i store the lm object in the database and and is it possible to retrive it and get the summary information -- View this message in context: http://r.789695.n4.nabble.com/Storing-the-linear-model-object-tp4121944p4121944.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. -- Diane Bailleul Doctorante Université Paris-Sud 11 - Faculté des Sciences d'Orsay Unité Ecologie, Systématique et Evolution Département Biodiversité, Systématique et Evolution UMR 8079 - UPS CNRS AgroParisTech Porte 320, premier étage, Bâtiment 360 91405 ORSAY CEDEX FRANCE (0033) 01.69.15.56.64 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] All combinations
Dear all, I would like something simple to do in R that I do not know how I should search for it. Let's say that I have a list of a-c(1,2,3,4,5) b-(6,7,8) and I want to get back all their possible cominations like 1,6 1,7 1,8 2,6 2,7 2,8 3,6 3,7 3,8 and so on. How I can do that? B.R Alex [[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] All combinations
expand.grid() This one is admittedly rather hard to find... Michael On Nov 30, 2011, at 7:15 AM, Alaios ala...@yahoo.com wrote: Dear all, I would like something simple to do in R that I do not know how I should search for it. Let's say that I have a list of a-c(1,2,3,4,5) b-(6,7,8) and I want to get back all their possible cominations like 1,6 1,7 1,8 2,6 2,7 2,8 3,6 3,7 3,8 and so on. How I can do that? B.R Alex [[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] Generalized singular value decomposition
Hi, Googling for R-help gsvd leads to a number of interesting entries on the mailing list. It seems that GSVD is present in the lapack version included with R and can be called (see the mailing list entries). https://stat.ethz.ch/pipermail/r-help/2004-August/056713.html http://r.789695.n4.nabble.com/problem-with-generalized-singular-value-decomposition-using-LAPACK-td834930.html http://r.789695.n4.nabble.com/Generalized-SVD-td798791.html regards, Paul On 11/30/2011 10:51 AM, Oana Tomescu wrote: Hello, I would like to perform a generalized singular value decomposition with R. The only possibility I found is GSVD that is based on LAPACK/BLAS. Are there other possibilities too? If not, has anybody used LAPACK/BLAS under Windows XP? How can I install them? Following [1] did not help. I hope this is the right place for my question. Thank you very much! Oana Tomescu [1] http://sites.google.com/site/jivsoft/Home/r-blas-interface __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Paul Hiemstra, Ph.D. Global Climate Division Royal Netherlands Meteorological Institute (KNMI) Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39 P.O. Box 201 | 3730 AE | De Bilt tel: +31 30 2206 494 http://intamap.geo.uu.nl/~paul http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Notation
Hi Silvano, Your function is not part of the standard R installation. After some googling I found out, but next time: PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. The betareg package provides manual pages and a vignette (tutorial) which should provide the information you are looking for. http://cran.r-project.org/web/packages/betareg/index.html cheers, Paul On 11/29/2011 04:50 PM, Silvano wrote: Hi, what's mean / in command: betareg(inf~Grupo/Sexo, data=dados) it's a effect nested? -- Silvano Cesar da Costa Departamento de Estatística Universidade Estadual de Londrina Fone: 3371-4346 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Paul Hiemstra, Ph.D. Global Climate Division Royal Netherlands Meteorological Institute (KNMI) Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39 P.O. Box 201 | 3730 AE | De Bilt tel: +31 30 2206 494 http://intamap.geo.uu.nl/~paul http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in log
I'd suggest you do some leg-work and figure out why you are getting values 1. If your algorithm is motivated by some approximation then a min() or pmin() *might* be the right fix, but if there are no approximations you may need to start debugging properly to see why you are getting an out of bounds value. Since there's no random number generation involved, I'd hesitate to just throw out the result without knowing its source. Also keep in mind the limitations of floating point arithmetic if you expect alpha*d^beta to be small. Michael On Nov 29, 2011, at 6:58 PM, Sarah Goslee sarah.gos...@gmail.com wrote: On Tue, Nov 29, 2011 at 6:55 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: yes, log of negative number is undefined and R also do the same and produces NaNs. Here I want to reject the value of exp(-alpha*d^(-beta)) when greater than 1, and want to run the loop otherwise. Thanks Then just add another if() statement checking for that condition. On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com wrote: Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. Wait... you're complaining that you can't take the natural log of a negative number in R? You can't do that anywhere. What do you expect to happen? The log of a negative number IS NaN. Sarah On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: I have following code: loglikelihood - function(alpha,beta= 0.1){ loglh-0 d-0 p-0 k-NULL data-read.table(epidemic.txt,header = TRUE) attach(data, warn.conflicts = F) k -which(inftime==1) d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta) p-1 - exp(-alpha*d) for(i in 1:100){ if(i!=k){ if(inftime[i]==0){ loglh-loglh +log(1-p[i]) } if(inftime[i]==2){ loglh-loglh + log(p[i]) } } } return(loglh) } Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: No, that,s not a problem Michael, I have following code: loglikelihood - function(alpha,beta= 0.1){ loglh-0 d-0 p-0 k-NULL data-read.table(epidemic.txt,header = TRUE) attach(data, warn.conflicts = F) k -which(inftime==1) d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta) p-1 - exp(-alpha*d) for(i in 1:100){ if(i!=k){ if(inftime[i]==0){ loglh-loglh +log(1-p[i]) } if(inftime[i]==2){ loglh-loglh + log(p[i]) } } } return(loglh) } Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. -- Sarah Goslee http://www.functionaldiversity.org __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in log
The loglikelihood() looks ok and gives some value. But I am using this function for the simulated annealing, generating the random samples from uniform proposal density., The codes are as follows epiannea - function(T0 = 1, N = 500,beta = 0.1,x0 = 0.1, rho = 0.90, eps = 0.1, loglikelihood, data){ moving - 1 count - 0 Temp - T0 alpha - x0 while(moving 0){ mmoving - 0 for(i in 1:N){ r - alpha + runif(1,-eps,eps) if(r 0){ a - min(1,exp((loglikelihood(alpha,beta)-loglikelihood(r,beta))/Temp)) } if(runif(1) a){ moving - moving +1 alpha - r } } Temp - Temp*rho count - count + 1 } #plot(a,loglikelihood(alpha,beta), type = l) fmin - loglikelihood(alpha, beta) return(c(fmin, alpha, count)) } epiannea(loglikelihood=loglikelihood) Some times it shows the warnings that logp[i] produces NaNs, that means p[i] is negative, but it should not be that because p[i] is the probabiliy and can't be negative. Some times the code runs but never stop. On Wed, Nov 30, 2011 at 7:34 AM, R. Michael Weylandt michael.weyla...@gmail.com michael.weyla...@gmail.com wrote: I'd suggest you do some leg-work and figure out why you are getting values 1. If your algorithm is motivated by some approximation then a min() or pmin() *might* be the right fix, but if there are no approximations you may need to start debugging properly to see why you are getting an out of bounds value. Since there's no random number generation involved, I'd hesitate to just throw out the result without knowing its source. Also keep in mind the limitations of floating point arithmetic if you expect alpha*d^beta to be small. Michael On Nov 29, 2011, at 6:58 PM, Sarah Goslee sarah.gos...@gmail.com wrote: On Tue, Nov 29, 2011 at 6:55 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: yes, log of negative number is undefined and R also do the same and produces NaNs. Here I want to reject the value of exp(-alpha*d^(-beta)) when greater than 1, and want to run the loop otherwise. Thanks Then just add another if() statement checking for that condition. On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com wrote: Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. Wait... you're complaining that you can't take the natural log of a negative number in R? You can't do that anywhere. What do you expect to happen? The log of a negative number IS NaN. Sarah On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: I have following code: loglikelihood - function(alpha,beta= 0.1){ loglh-0 d-0 p-0 k-NULL data-read.table(epidemic.txt,header = TRUE) attach(data, warn.conflicts = F) k -which(inftime==1) d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta) p-1 - exp(-alpha*d) for(i in 1:100){ if(i!=k){ if(inftime[i]==0){ loglh-loglh +log(1-p[i]) } if(inftime[i]==2){ loglh-loglh + log(p[i]) } } } return(loglh) } Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: No, that,s not a problem Michael, I have following code: loglikelihood - function(alpha,beta= 0.1){ loglh-0 d-0 p-0 k-NULL data-read.table(epidemic.txt,header = TRUE) attach(data, warn.conflicts = F) k -which(inftime==1) d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta) p-1 - exp(-alpha*d) for(i in 1:100){ if(i!=k){ if(inftime[i]==0){ loglh-loglh +log(1-p[i]) } if(inftime[i]==2){ loglh-loglh + log(p[i]) } } } return(loglh) } Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. -- Sarah Goslee http://www.functionaldiversity.org __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.htmlhttp://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
Re: [R] Generalized singular value decomposition
Thank you for your help! The GSVD (SVDgen) in PTAk does not perform a generalized singular value decomposition for 2 matrices, only for 1. I should have mentioned this - sorry. Are there maybe other packages? I have also found the last 2 links, but I look for a way to use LAPACK with Windows - this issue is unfortunately not addressed in these posts. Can anyone help me with LAPACK and Windows XP? I have noticed that there are already 2 dlls in the R/bin folder (Rblas.dll and Rlapack.dll). But I still can not load the libraries. I guess that I still have to install/do something ... I would appreciate any help! Regards, Oana Am 30.11.2011 13:18, schrieb Paul Hiemstra: Hi, Googling for R-help gsvd leads to a number of interesting entries on the mailing list. It seems that GSVD is present in the lapack version included with R and can be called (see the mailing list entries). https://stat.ethz.ch/pipermail/r-help/2004-August/056713.html http://r.789695.n4.nabble.com/problem-with-generalized-singular-value-decomposition-using-LAPACK-td834930.html http://r.789695.n4.nabble.com/Generalized-SVD-td798791.html regards, Paul On 11/30/2011 10:51 AM, Oana Tomescu wrote: Hello, I would like to perform a generalized singular value decomposition with R. The only possibility I found is GSVD that is based on LAPACK/BLAS. Are there other possibilities too? If not, has anybody used LAPACK/BLAS under Windows XP? How can I install them? Following [1] did not help. I hope this is the right place for my question. Thank you very much! Oana Tomescu [1] http://sites.google.com/site/jivsoft/Home/r-blas-interface __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Invito a collegarsi su LinkedIn
Vorrei aggiungerti alla mia rete professionale su LinkedIn. -giovanni giovanni parrinello Senior Statistical Consultant at University of Brescia Italy Confirm that you know giovanni parrinello: https://www.linkedin.com/e/-92ffr2-gvmdq8pr-1/isd/5085379353/o5O4tnc8/?hs=falsetok=0BYzKxkO2MnB01 -- You are receiving Invitation to Connect emails. Click to unsubscribe: http://www.linkedin.com/e/-92ffr2-gvmdq8pr-1/vxMuFwS2PpRnFjm0VsifityWSt9K_bXjxepMJ1/goo/R-help%40stat%2Emath%2Eethz%2Ech/20061/I1773723218_1/?hs=falsetok=0xPSvS6c2MnB01 (c) 2011 LinkedIn Corporation. 2029 Stierlin Ct, Mountain View, CA 94043, USA. [[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] glm step() fkt
Hi volks, i have a question about the step() fkt. Is there a possibility to save the last model generated from this method. I have a loop and so i generate 100 different models with the step fkt and i want to know which model is the most common. CODE: ... missStep - numeric(100) for (j in 1:100) { trainindex - sample(c(1:462),300) train - data[trainindex,] test - data[-trainindex,] mod.step - step(mod, direction=both,trace=T) } ... So i want to store each iteration the best model. Thx in advanced people, greetings -- View this message in context: http://r.789695.n4.nabble.com/glm-step-fkt-tp4122588p4122588.html Sent from the R help mailing list archive at Nabble.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.
[R] how to call a function for each row
Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp4122906p4122906.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] Non parametric, repeated-measures, factorial ANOVA
Hi Thanks for your reply. The package appears to be for the analysis of longitudinal designs, and requires a time factor to be input on the relevant . My experiment is just a standard repeated-measures design, and I just need a direct, non-parametricc equivalent of the Factorial repeated-measures ANOVA. It occurs to me that the terminology I am using might be a problem - when I say repeated measures I mean that it is a within subjects design, so that each subject does each condition, counterbalanced so that order/time effects are not an issue. Thanks Rob -- View this message in context: http://r.789695.n4.nabble.com/Non-parametric-repeated-measures-factorial-ANOVA-tp4118816p4122501.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] Replace columns in a data.frame randomly splitted
Dear community, I'm working with the data.frame attached ( http://r.789695.n4.nabble.com/file/n4122926/df1.xls df1.xls ), let's call it df1. I typed: df1- read.xls(C:/... dir .../df1.xls,colNames= TRUE, rowNames= TRUE) Then I splited randomly df1 using splitdf function (http://gettinggeneticsdone.blogspot.com/2011/03/splitting- dataset-revisited-keeping.html) So now, I have df1 divided in another 2 dataframes: trainset and testset. (They both have the same row names and column names as df1) I'd like to change df1$v1 and df1$v2 values in the trainset and testset by the ones in df2 ( http://r.789695.n4.nabble.com/file/n4122926/df2.xls df2.xls ) How can I explain R to identify in trainset/testset$v1 and $v2 the rownames and replace the values by the corresponding ones in df2? Thanks in advance, u...@host.com -- View this message in context: http://r.789695.n4.nabble.com/Replace-columns-in-a-data-frame-randomly-splitted-tp4122926p4122926.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] Problems with Raster and clim.pact packages with large netcdf files (2.7G) in x64 bit R
Emmanuel Jjunju ejjunju at gmail.com writes: I normally use the raster or clim.pact pckages to read netcdf (.nc) files. This has always worked out for me until this weekend every time i try to read a .nc file i get the following error ... Could you please: 1. Consider moving this thread to R-sig-geo, where it may get more response? 2. Try to reproduce this just using the functions in the ncdf package, as this is the package that both clim.pact and raster use for reading NetCDF files. 3. If you have access to a non-Windows platform, see if the fault is the same, which would point to something in ncdf/64-bit, or if not, would point to an issue in the Windows 64-bit NetCDF external dependency. 4. Since this is also file-dependent, see if any other similarly large NetCDF files with a similar structure. There is some traffic on NetCDF lists suggesting possible issues with posix.c Hope this helps, Roger I have confirmed that the problem is with R x64 Bit regardless of whether using Rstudio or R GUI or any other GUI. I repeat: No problems with 32 bit which is strange given that 32 bit restricts my memory to 3583 !! and if i load more data in the workspace i will run out of memeory. My machine has 16GB of RAM __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 call a function for each row
Read ?apply This is easiest: df - matrix(c(1,2,3,3,4,6,5,6,1), 3) apply(df, 1, function(x) 6*x[1]+7*x[2]+8*x[3]) But it's much more efficient to do it with matrix multiplication. In keeping with the best of tradition, this is left as an exercise to the reader. Michael On Wed, Nov 30, 2011 at 8:10 AM, arunkumar akpbond...@gmail.com wrote: Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp4122906p4122906.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 call a function for each row
will this do it: x - read.table(text = X1 X2 X3 + 1 3 5 + 2 4 6 + 3 6 1, header = TRUE) x X1 X2 X3 1 1 3 5 2 2 4 6 3 3 6 1 apply(x, 1, function(a) 6 * a[1] + 7 * a[2] + 8 * a[3]) [1] 67 88 68 On Wed, Nov 30, 2011 at 8:10 AM, arunkumar akpbond...@gmail.com wrote: Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp4122906p4122906.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. -- Jim Holtman Data Munger Guru What is the problem that you are trying to solve? Tell me what you want to do, not how you want to do 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] how to call a function for each row
Homework. ? I would say that this indiicates that you need to open the R tutorials and start reading. -- Bert On Wed, Nov 30, 2011 at 5:10 AM, arunkumar akpbond...@gmail.com wrote: Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp4122906p4122906.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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] glm step() fkt
Put them in a list: ModelList - vector(list, 100) ModelList[[i]] - mod.step - step(mod, direction=both,trace=T) Then come back and use sapply() to do whatever you want to the set of models to compare/count/etc them Michael On Wed, Nov 30, 2011 at 6:12 AM, Schrabauke bern...@yahoo.de wrote: Hi volks, i have a question about the step() fkt. Is there a possibility to save the last model generated from this method. I have a loop and so i generate 100 different models with the step fkt and i want to know which model is the most common. CODE: ... missStep - numeric(100) for (j in 1:100) { trainindex - sample(c(1:462),300) train - data[trainindex,] test - data[-trainindex,] mod.step - step(mod, direction=both,trace=T) } ... So i want to store each iteration the best model. Thx in advanced people, greetings -- View this message in context: http://r.789695.n4.nabble.com/glm-step-fkt-tp4122588p4122588.html Sent from the R help mailing list archive at Nabble.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. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Negative exponential fit
Thanks Peter, will take a look at this package. Regards, Indrajit From: Peter Ruckdeschel peter.ruckdesc...@web.de To: r-h...@stat.math.ethz.ch Sent: Wednesday, November 30, 2011 5:03 PM Subject: Re: [R] Negative exponential fit I do not want to shed out any doubt as to the merits of pkg fitdistrplus, in particular for censored data, but seconding Uwe's reply later in this thread, you may also want to check out pkg distrMod on CRAN --- M. Kohl, P. Ruckdeschel (2010): R Package distrMod: S4 Classes and Methods for Probability Models. Journal of Statistical Software, 35(10), 1-27. URL http://www.jstatsoft.org/v35/i10/. which offers quite some additional flexibility for model fitting---including new models (built on distributions which have no name but instead are image distributions under arithmetic transformations of existing ones, see example M2 in the cited ref) and (nonlinear) transformations of the parameter (see Example p.15, cited ref). Best regards, Peter __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[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 call a function for each row
Isn't this even easier? X1 - c(1:3) X2 - c(3, 4, 6) X3 - c(5, 6, 1) Y - 6*X1 + 7*X2 + 8*X3 Y [1] 67 88 68 Or if you really need a function: MakeY - function(x, y, z) 6*x + 7*y + 8*z MakeY(X1, X2, X3) [1] 67 88 68 -- David L Carlson Associate Professor of Anthropology Texas AM University College Station, TX 77843-4352 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of R. Michael Weylandt Sent: Wednesday, November 30, 2011 8:17 AM To: arunkumar Cc: r-help@r-project.org Subject: Re: [R] how to call a function for each row Read ?apply This is easiest: df - matrix(c(1,2,3,3,4,6,5,6,1), 3) apply(df, 1, function(x) 6*x[1]+7*x[2]+8*x[3]) But it's much more efficient to do it with matrix multiplication. In keeping with the best of tradition, this is left as an exercise to the reader. Michael On Wed, Nov 30, 2011 at 8:10 AM, arunkumar akpbond...@gmail.com wrote: Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp41 22906p4122906.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] SPSS - R
Dear Michael, The original poster and I ended up pursuing this question off-list and, I think, resolving the issue, discovering in the process that SPSS apparently reports type-II tests that are really type-III tests. (Yes, one fits a multivariate linear model prior to calling Anova() with the idata and idesign arguments.) Best, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Michael Dewey Sent: November-30-11 6:49 AM To: Kristi Shoemaker; r-help@r-project.org Subject: Re: [R] SPSS - R At 16:41 27/11/2011, Kristi Shoemaker wrote: Hi John, Your assumptions are correct and those examples were very helpful, thanks! I think I'm almost there, but I'm screwing up something with the within-subjects factor (the example has two, I only have one). See below. Why am I not seeing my within-subjects factors in the ANOVA report? I am by no means an expert in this but my reading of the documentation suggests that the idata and idesign parameters are for multivariate models and that is not what you specified. Dat2 Subj Age Sex Time HippV 1 s01 32.9 1 0w 6.50098 2 s01 32.9 1 6w 6.91793 3 s02 35.1 0 0w 7.32480 4 s02 35.1 0 6w 7.56012 5 s03 34.4 0 0w 6.51385 6 s03 34.4 0 6w 6.56875 9 s05 39.9 1 0w 6.92855 10 s05 39.9 1 6w 6.94926 11 s06 29.5 1 0w 6.99383 12 s06 29.5 1 6w 7.10568 13 s07 45.9 1 0w 6.94380 14 s07 45.9 1 6w 7.08190 15 s08 20.3 1 0w 7.76881 16 s08 20.3 1 6w 7.72725 17 s09 26.9 0 0w 5.37566 18 s09 26.9 0 6w 5.74887 21 s11 22.0 0 0w 7.12992 22 s11 22.0 0 6w 7.16237 23 s12 31.0 1 0w 6.70629 24 s12 31.0 1 6w 6.80872 25 s13 50.1 1 0w 7.22649 26 s13 50.1 1 6w 7.58900 27 s14 22.2 0 0w 5.97577 28 s14 22.2 0 6w 5.80801 29 s16 20.4 1 0w 7.99554 30 s16 20.4 1 6w 8.09260 31 s20 24.0 0 0w 7.01014 32 s20 24.0 0 6w 6.87821 33 s21 32.4 0 0w 5.90883 34 s21 32.4 0 6w 5.95392 37 s24 22.2 0 0w 6.15474 38 s24 22.2 0 6w 6.00906 41 s27 22.2 1 0w 7.88765 42 s27 22.2 1 6w 7.76038 49 s35 24.3 0 0w 6.05998 50 s35 24.3 0 6w 6.07399 51 s36 23.5 0 0w 7.83182 52 s36 23.5 0 6w 7.60268 53 s38 59.7 1 0w 7.39672 54 s38 59.7 1 6w 6.98291 55 s39 40.5 0 0w 7.31330 56 s39 40.5 0 6w 7.50559 57 s40 24.2 1 0w 8.54958 58 s40 24.2 1 6w 8.65016 59 s41 23.6 1 0w 7.76049 60 s41 23.6 1 6w 7.58946 61 s42 53.3 0 0w 7.03388 62 s42 53.3 0 6w 7.48384 63 s43 34.4 0 0w 6.86967 64 s43 34.4 0 6w 6.81076 65 s44 44.8 0 0w 7.33779 66 s44 44.8 0 6w 7.86175 67 s45 40.2 1 0w 6.55963 68 s45 40.2 1 6w 6.52577 71 s47 26.5 0 0w 7.09418 72 s47 26.5 0 6w 6.92850 73 s48 22.7 0 0w 6.77078 74 s48 22.7 0 6w 6.67289 75 s50 36.1 1 0w 7.47208 76 s50 36.1 1 6w 7.55876 library(car) Loading required package: MASS Loading required package: nnet Loading required package: survival Loading required package: splines timepoints - factor(c(0w, 6w), + levels=c(0w,6w)) idata=data.frame(timepoints) mod.ok - lm(HippV ~ Age*Sex,data=Dat2) (av.ok - Anova(mod.ok, idata=idata, idesign=~timepoints)) Anova Table (Type II tests) Response: HippV Sum Sq Df F valuePr(F) Age0.0234 1 0.071 0.7909202 Sex5.3082 1 16.128 0.0001781 *** Age:Sex4.7772 1 14.514 0.0003478 *** Residuals 18.4314 56 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 From: John Fox j...@mcmaster.ca Cc: r-help@r-project.org Sent: Saturday, November 26, 2011 7:19 PM Subject: RE: [R] SPSS - R Dear Kristi, I assume that this is a repeated-measures ANOVA with one within-subjects factor (Time) and two between-subjects factors (Age and Sex, which are crossed). If Age is numeric, and not a factor, then the type-III tests that you requested don't test sensible hypotheses. In any event, if my guess is right about the design, then you can use the Anova() function in the car package for an equivalent analysis. See the repeated-measures example in ?Anova (for the O'Brien and Kaiser data). You've already had an answer to the more general question. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On
Re: [R] Generalized singular value decomposition
tom wrote Thank you for your help! The GSVD (SVDgen) in PTAk does not perform a generalized singular value decomposition for 2 matrices, only for 1. I should have mentioned this - sorry. Are there maybe other packages? I have also found the last 2 links, but I look for a way to use LAPACK with Windows - this issue is unfortunately not addressed in these posts. Can anyone help me with LAPACK and Windows XP? I have noticed that there are already 2 dlls in the R/bin folder (Rblas.dll and Rlapack.dll). But I still can not load the libraries. I guess that I still have to install/do something ... I don't think you have to install something. You have not told us how you are trying to load the dll's. You can use the GSVD function from http://r.789695.n4.nabble.com/problem-with-generalized-singular-value-decomposition-using-LAPACK-td834930.html but you will need the exact path to the Lapack dll in the dyn.load call. Something like dyn.load(C:\\R\\.\\libRlapack.dll). You do not need to dyn.load the Blas dll. I have tried this GSVD function with Mac OS X with success. Berend -- View this message in context: http://r.789695.n4.nabble.com/Generalized-singular-value-decomposition-tp4122542p4123335.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] Parameters setting in functions optimization
optimx does allow you to use bounds. The default is using only methods from optim(), but even though I had a large hand in those methods, and they work quite well, there are other tools available within optimx that should be more appropriate for your problem. For example, the current version of optimx should work quite well with lower and upper bounds specified, and you can see which methods work well by putting control=list(all.methods=TRUE). all.methods would be overkill for general use of course. I am hoping to have a new version of optimx up on R-forge within a week -- there are so many options to check -- that traps the NaNs etc. if it can, and also allows parameter and function scaling as well as several other new features. This is all experimental at the moment, but initial results are promising. This is less about better methods than about trapping errors and bad scaling etc. However, if you are able to share your script and data, I'll be happy to use it as a test and report back to you if you can communicate it to me off-list. Best, JN On 11/30/2011 06:00 AM, r-help-requ...@r-project.org wrote: Message: 68 Date: Tue, 29 Nov 2011 19:15:43 +0100 From: Diane Bailleul diane.baill...@u-psud.fr To: r-help@r-project.org Subject: [R] Parameters setting in functions optimization Message-ID: 4ed5214f.7030...@u-psud.fr Content-Type: text/plain; charset=ISO-8859-1; format=flowed Good afternoon everybody, I'm quite new in functions optimization on R and, whereas I've read lot's of function descriptions, I'm not sure of the correct settings for function like optimx and nlminb. I'd like to minimize my parameters and the loglikelihood result of the function. My parameters are a mean distance of dispersion and a proportion of individuals not assigned, coming from very far away. The function LikeGi reads external tables and it's working as I want (I've got a similar model on Mathematica). My final function is LogLiketot : LogLiketot- function(dist,ms) { res - NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res - c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } dist is the mean dispersal distance (0, lots of meters) and ms the proportion of individuals (0-1). Of course, I want them to be as low as possible. I'd tried to enter the initials parameters as indicated in the tutorials : optim(c(40,0.5), fn=LogLiketot) Error in 1 - ms : 'ms' is missing But ms is 0.5 ... So I've tried this form : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) with different values for the two parameters : par fvalues method fns grs itns conv KKT1 KKT2 xtimes 219.27583, 25.37964 2249.698BFGS 12 8 NULL0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL0 TRUE TRUE 66.3 The first line is not possible but as I've not constrained the optimization ... but the second line would be a very good result ! Then, searching for another similar cases, I've tried to change my function form: LogLiketot- function(par) { res - NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res - c(res, pop5[i,l]*log(LikeGi(l,i,par[1],par[2]))) } } return(-sum(res)) } where dist=par[1] and ms=par[2] And I've got : optimx(c(40,0.5), fn=LogLiketot) par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083BFGS 29 10 NULL0 TRUE NA 92.03 1 39.7372199, 0.9778101 1064.083 Nelder-Mead 53 NA NULL0 TRUE NA 70.83 And I've got now a warning message : In log(LikeGi(l, i, par[1], par[2])) : NaNs produced (which are very bad results in that case) Anyone with previous experiences in optimization of several parameters could indicate me the right way to enter the initial parameters in this kind of functions ? Thanks a lot for helping me ! Diane -- Diane Bailleul Doctorante Universit? Paris-Sud 11 - Facult? des Sciences d'Orsay Unit? Ecologie, Syst?matique et Evolution D?partement Biodiversit?, Syst?matique et Evolution UMR 8079 - UPS CNRS AgroParisTech Porte 320, premier ?tage, B?timent 360 91405 ORSAY CEDEX FRANCE (0033) 01.69.15.56.64 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] All combinations
On Nov 30, 2011, at 7:18 AM, R. Michael Weylandt wrote: expand.grid() This one is admittedly rather hard to find... Well, it is linked from the `combn` help page. And it is the likely to be first or second in a search with ??combinations since it is in pkg:base and at least on my interface the displayed hits are sorted by alpha(pakgname), so I would disagree with it being hard to find. Other ideas After replacing the missing `c` function: outer(a,b,FUN=paste, sep=,) [,1] [,2] [,3] [1,] 1,6 1,7 1,8 [2,] 2,6 2,7 2,8 [3,] 3,6 3,7 3,8 [4,] 4,6 4,7 4,8 [5,] 5,6 5,7 5,8 Perhaps not what the OP asked for, but then exactly what did the OP ask for, anyway? Perhaps this? (Or not.) as.data.frame(sapply(a, function(x){ sapply(b, function(y) c(x,y),simplify=FALSE)}) ) V1 V2 V3 V4 V5 1 1, 6 2, 6 3, 6 4, 6 5, 6 2 1, 7 2, 7 3, 7 4, 7 5, 7 3 1, 8 2, 8 3, 8 4, 8 5, 8 Interesting how print() handles data.frame columns of lists, don't you think? And then, of course, building it from scratch: matrix(c( rep(a, length(b)), rep(b, each=length(a))), ncol=2) [,1] [,2] [1,]16 [2,]26 [3,]36 [4,]46 [5,]56 [6,]17 [7,]27 [8,]37 [9,]47 [10,]57 [11,]18 [12,]28 [13,]38 [14,]48 [15,]58 -- David. Michael On Nov 30, 2011, at 7:15 AM, Alaios ala...@yahoo.com wrote: Dear all, I would like something simple to do in R that I do not know how I should search for it. Let's say that I have a list of a-c(1,2,3,4,5) b-(6,7,8) and I want to get back all their possible cominations like 1,6 1,7 1,8 2,6 2,7 2,8 3,6 3,7 3,8 and so on. How I can do that? B.R Alex David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Parameters setting in functions optimization
Dear Florent, I know that I'm asking to optim to minimize my values, and that the results with a lower fvalue are best supported than those with a higher fvalue. My comment was just from a data point of view. I'd like the lower ms (second parameter) as possible, as well as the fvalue. So a ms of 0.97 (i.e. that 97% of my individuals are coming from outside the experiment plot) is very disapointing. Dear John, I can send you my data by email. It's very kind of you to offer to use my data as test for your new and coming optimx ! thank you. Le 11/30/2011 12:08 PM, Florent D. a écrit : With optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) where LogLiketot- function(dist,ms) { res- NULL for(i in 1:nrow(pop5)) { for(l in 1:nrow(freqvar)) { res- c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } I think it will do something like this for the first call to LogLiketot: LogLiketot(c(30,50), ms=c(0.4,0.5)) which is obviously not the usage you had in mind. Also, I see you say the results for the bad usage above: par fvalues method fns grs itns conv KKT1 KKT2 xtimes 219.27583, 25.37964 2249.698BFGS 12 8 NULL0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL0 TRUE TRUE 66.3 look very good but you do not comment about the results for the correct usage of optimx: par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083BFGS 29 10 NULL0 TRUE NA 92.03 1 39.7372199, 0.9778101 1064.083 Nelder-Mead 53 NA NULL0 TRUE NA 70.83 Do you realize optimx is trying to _minimize_ your function? See that the fvalues from the correct usage are much better (smaller) than you first (bad) usage. On Wed, Nov 30, 2011 at 4:16 AM, Diane Bailleul diane.baill...@u-psud.fr wrote: Le 11/30/2011 2:09 AM, Florent D. a écrit : Thanks for your answer ! I also think your last write-up for LogLiketot (using a single argument par) is the correct approach if you want to feed it to optim(). I'm not dedicated to optim() fonction. I just want to optimise my two parameters and the loglikelihood result, and if there's a better fonction for that, I wish I could use it. So now you have a problem with log(LikeGi(l, i, par[1], par[2])) for some values of par[1] and par[2]. Where is LikeGi coming from? a package or is it your own function? My own function, otherwise it would be simplier to discuss about my problems. You could add some print statements (if you are familiar with browser() it is even better) so you may see what values of par are causing trouble. I'm not familiar, but I'll search about browser(). If the function with par is correct, any idea of what I've made with this : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) ? On Tue, Nov 29, 2011 at 1:15 PM, Diane Bailleul diane.baill...@u-psud.frwrote: Good afternoon everybody, I'm quite new in functions optimization on R and, whereas I've read lot's of function descriptions, I'm not sure of the correct settings for function like optimx and nlminb. I'd like to minimize my parameters and the loglikelihood result of the function. My parameters are a mean distance of dispersion and a proportion of individuals not assigned, coming from very far away. The function LikeGi reads external tables and it's working as I want (I've got a similar model on Mathematica). My final function is LogLiketot : LogLiketot- function(dist,ms) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,dist,ms))) } } return(-sum(res)) } dist is the mean dispersal distance (0, lots of meters) and ms the proportion of individuals (0-1). Of course, I want them to be as low as possible. I'd tried to enter the initials parameters as indicated in the tutorials : optim(c(40,0.5), fn=LogLiketot) Error in 1 - ms : 'ms' is missing But ms is 0.5 ... So I've tried this form : optimx(c(30,50),ms=c(0.4,0.5), fn=LogLiketot) with different values for the two parameters : par fvalues method fns grs itns conv KKT1 KKT2 xtimes 219.27583, 25.37964 2249.698BFGS 12 8 NULL0 TRUE TRUE 57.5 1 29.6787861, 0.1580298 2248.972 Nelder-Mead 51 NA NULL0 TRUE TRUE 66.3 The first line is not possible but as I've not constrained the optimization ... but the second line would be a very good result ! Then, searching for another similar cases, I've tried to change my function form: LogLiketot- function(par) { res- NULL for(i in 1:nrow(pop5)){ for(l in 1:nrow(freqvar)){ res- c(res, pop5[i,l]*log(LikeGi(l,i,par[1],par[2]))) } } return(-sum(res)) } where dist=par[1] and ms=par[2] And I've got : optimx(c(40,0.5), fn=LogLiketot) par fvalues method fns grs itns conv KKT1 KKT2 xtimes 2 39.9969607, 0.9777634 1064.083BFGS
Re: [R] Negative exponential fit
On Nov 29, 2011, at 23:19 , Ben Bolker wrote: rch4 rch4 at geneseo.edu writes: We need help We are doing a project for a statistical class in and we are looking at world record times in different running events over time. We are trying to fit the data with a negative exponential but we just cant seem to get a function that works properly. we have on our x-axis the date and on the y-axis the time(in seconds). So as you can imagine, the times have decreased and appear to be approaching a limit. Any ideas for a nls function that would work for us would be greatly appreciated. Rob I disagree with the other solutions posted here: think you're looking not for a distribution, but for the change over time. Not that this is anywhere near my areas of expertise, but wouldn't you want to be even more careful than that? I mean, surely the record time is nondecreasing, and one would expect that the time between records to carry information about the issue (e.g., in a stable situation, it should increase as a lower limit is being approached)? You could start with fit1 - lm(log(time)~I(date-date[1])) where the intercept will be the *log* of the intercept (value on the first date) and the slope will be the exponential coefficient. If you need to be more careful about your statistical assumptions (e.g. if the variance appears to be homogeneous on the original scale but not on the log scale) then something like fit2 - nls(exp(logint)*exp(-r*(date-date[1])), start=...) should work. You need to set the starting values appropriately -- the values from the linear fit above should be pretty good. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Negative exponential fit
On 11-11-30 11:32 AM, peter dalgaard wrote: On Nov 29, 2011, at 23:19 , Ben Bolker wrote: rch4 rch4 at geneseo.edu writes: We need help We are doing a project for a statistical class in and we are looking at world record times in different running events over time. We are trying to fit the data with a negative exponential but we just cant seem to get a function that works properly. we have on our x-axis the date and on the y-axis the time(in seconds). So as you can imagine, the times have decreased and appear to be approaching a limit. Any ideas for a nls function that would work for us would be greatly appreciated. Rob I disagree with the other solutions posted here: think you're looking not for a distribution, but for the change over time. Not that this is anywhere near my areas of expertise, but wouldn't you want to be even more careful than that? I mean, surely the record time is nondecreasing, and one would expect that the time between records to carry information about the issue (e.g., in a stable situation, it should increase as a lower limit is being approached)? All that seems reasonable. In addition, the lower limit is not zero (which my answer assumed). However, the OP can't get a negative exponential fit to work in the first place, so they should probably be starting with something simple ... To deal with the lower limit is not zero problem they can just add a parameter: fit2 - nls(minval+exp(logint)*exp(-r*(date-date[1])), start=...) You could start with fit1 - lm(log(time)~I(date-date[1])) where the intercept will be the *log* of the intercept (value on the first date) and the slope will be the exponential coefficient. If you need to be more careful about your statistical assumptions (e.g. if the variance appears to be homogeneous on the original scale but not on the log scale) then something like fit2 - nls(exp(logint)*exp(-r*(date-date[1])), start=...) should work. You need to set the starting values appropriately -- the values from the linear fit above should be pretty good. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Replace columns in a data.frame randomly splitted
Two ways immediately come to mind: 1) Change the values before splitting 2) Use the same seed for the two splits so the rows match exactly and then just do the changes directly Any more automated solution will depend on whether your data has rownames or not. Michael PS - As a general rule, most folks don't like to download unknown xls files -- it's much easier to create a plain text representation of R data using the dput() command. On Wed, Nov 30, 2011 at 8:15 AM, agent dunham crossp...@hotmail.com wrote: Dear community, I'm working with the data.frame attached ( http://r.789695.n4.nabble.com/file/n4122926/df1.xls df1.xls ), let's call it df1. I typed: df1- read.xls(C:/... dir .../df1.xls,colNames= TRUE, rowNames= TRUE) Then I splited randomly df1 using splitdf function (http://gettinggeneticsdone.blogspot.com/2011/03/splitting- dataset-revisited-keeping.html) So now, I have df1 divided in another 2 dataframes: trainset and testset. (They both have the same row names and column names as df1) I'd like to change df1$v1 and df1$v2 values in the trainset and testset by the ones in df2 ( http://r.789695.n4.nabble.com/file/n4122926/df2.xls df2.xls ) How can I explain R to identify in trainset/testset$v1 and $v2 the rownames and replace the values by the corresponding ones in df2? Thanks in advance, u...@host.com -- View this message in context: http://r.789695.n4.nabble.com/Replace-columns-in-a-data-frame-randomly-splitted-tp4122926p4122926.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Non-finite finite-difference value error in eha's aftreg
On Wed, Nov 16, 2011 at 3:06 PM, Milan Bouchet-Valat nalimi...@club.frwrote: Hi list! I'm getting an error message when trying to fit an accelerated failure time parametric model using the aftreg() function from package eha: Error in optim(beta, Fmin, method = BFGS, control = list(trace = as.integer(printlevel)), : non-finite finite-difference value [2] This only happens when adding four specific covariates at the same time in the model (see below). I understand that kind of problem can come from a too high correlations between my covariates, but is there anything I can do to avoid it? Yes, remove one (or more) covariate(s). The design matrix is almost certainly singular. Does something need to be improved in aftreg.fit? Yes, error messages. I'll look at it. Göran My data set is constituted of 34,505 observations (years) of 2,717 individuals, which seems reasonable to me to fit a complex model like that (covariates are all factors with less than 10 levels). I can send it by private mail if somebody wants to help debugging this. The details of the model and errors follow, but feel free to ask for more testing. I'm using R 2.13.1 (x86_64-redhat-linux-gnu), eha 2.0-5 and survival 2.36-9. Thanks for your help! m -aftreg(Surv(start, end, event) ~ homo1 + sexego + dipref1 ++ t.since.school.q, +data=ms, dist=loglogistic, id=ident) Error in optim(beta, Fmin, method = BFGS, control = list(trace = as.integer(printlevel)), : non-finite finite-difference value [2] Calls: aftreg - aftreg.fit - aftp0 - optim traceback() 4: optim(beta, Fmin, method = BFGS, control = list(trace = as.integer(printlevel)), hessian = TRUE) 3: aftp0(printlevel, ns, nn, id, strata, Y, X, offset, dis, means) 2: aftreg.fit(X, Y, dist, strats, offset, init, shape, id, control, center) 1: aftreg(Surv(start, end, event) ~ homo1 + sexego + dipref1 + t.since.school.q, data = ms, dist = loglogistic, id = ident) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Göran Broström [[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] install multtest and preprocessCore packages in Bioconductor library
Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] install multtest and preprocessCore packages in Bioconductor library
As the posting guide strongly suggests, I think the first step is to update to R 2.14: there have been changes in package design and I'm not sure back-compatibility will run that far. Doing so should let you download fresh binaries straight from BioC. I'm not sure why you would get compilation errors if you are just installing a binary -- are you sure you're doing it right? What are your commands? I think follow-up should be directed to the dedicated BioConductor list, but I assume they will tell you the same thing: 1) Update R and try to install prepackaged libraries; 2) install local copies of the binaries; 3) only compile from source if needed (you probably don't have a fortran compiler which seems to be part of the problem) Michael On Wed, Nov 30, 2011 at 12:57 PM, UyenThao Nguyen ungu...@tethysbio.com wrote: Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] install multtest and preprocessCore packages in Bioconductor library
My apologies: the posting guide doesn't actually suggest updating to the most recent stable version -- though perhaps it would be a good addition to the before posting section -- but I still think that's the necessary fix to your problem. Michael On Wed, Nov 30, 2011 at 1:08 PM, R. Michael Weylandt michael.weyla...@gmail.com wrote: As the posting guide strongly suggests, I think the first step is to update to R 2.14: there have been changes in package design and I'm not sure back-compatibility will run that far. Doing so should let you download fresh binaries straight from BioC. I'm not sure why you would get compilation errors if you are just installing a binary -- are you sure you're doing it right? What are your commands? I think follow-up should be directed to the dedicated BioConductor list, but I assume they will tell you the same thing: 1) Update R and try to install prepackaged libraries; 2) install local copies of the binaries; 3) only compile from source if needed (you probably don't have a fortran compiler which seems to be part of the problem) Michael On Wed, Nov 30, 2011 at 12:57 PM, UyenThao Nguyen ungu...@tethysbio.com wrote: Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] install multtest and preprocessCore packages in Bioconductor library
On Nov 30, 2011, at 12:57 PM, UyenThao Nguyen wrote: Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) You should read the directions about how to start with BioConductor. Generally, you need to install BiocLite or some such. I forget the exact specifics but they are all in the accompanying documentation. And if you need to install from sources, see comment below. 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran This suggests that in addition to not reading the Posting Guide (where you are asked not to append unrelated questions to unrelated threads AND you are asked to post Mac problems to the Mac list) that you also failed to read the Mac sections of the Installation and Administration Guide where you would have learned that you need to install Xcode. Please study that document more thouroughly if you plan to install from sources ... although in this case that should not be necessary since the BioC repository has binaries. collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/ library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/ Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. There is no L drive. Macs don't refer to their drives or partitions by letters. The dash indicates that L argument to a system command was not found. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls help
chuck.01 wrote datum - structure(list(Y = c(415.5, 3847.8325, 1942.83325, 1215.223, 950.142857325, 2399.585, 804.75, 579.5, 841.70825, 494.053571425 ), X = c(1.081818182, 0.492727273, 0.756363636, 0.896363636, 1.518181818, 0.49917, 1.354545455, 1.61, 1.706363636, 1.063636364 )), .Names = c(Y, X), row.names = c(NA, -10L), class = data.frame) with(datum, plot(Y~X)) As you can see there is a non-linear association between X and Y, and I would like to fit an appropriate model. I was thinking an exponential decay model might work well. I tried the following (a and k starting values are based off of a lm() fit), but get an error. fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1867)) Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates Try plot(datum$X,datum$Y) For more complex cases, plot the initial function you are trying to fit, but in this case it is easy to see that k is more in the order of 2. So try with k=2. Dieter -- View this message in context: http://r.789695.n4.nabble.com/nls-help-tp4123876p4124164.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] problem during installing bayesQR package for R 2.14 version
1. Start R (as administrator i.e. right click on the R icon on your desktop and select run as administrator). 2 Select Packages form menu bar 3 Select Install package(s) from local zip files.. and negotiate to your zip file. When you want to use the package you must enter library(bayesQR) I would strongly recomment that you read section 4 of the windows FAQ. This provides a lot of solutions to various problems that occur when one tries to install packages. This explains, in particular, how to install packages on systems where you do not have full admin rights. Best Regads John On 30 November 2011 05:31, kalam narendar reddy narendarcse...@gmail.com wrote: @John hello sir, i have R 2.14 version.and i have downloaded bayesQR package zip file as u said from following link http:// http://cran.r-project.org/web/packages/bayesQR/index.ht ml my OS is Windows7.i have downloaded Windows binary: bayesQR_1.3.zip file from above link.I am new to R. So please tell me what is the next step i have to do inorder to install the bayesQR package.pls reply me as quickly as possible.i am not getting the what the R manual for package installation for windows is a thanks in advance with regards *Kalam Narendar Reddy,* Masters In technology, University Of HYderabad, -- John C Frain Economics Department Trinity College Dublin Dublin 2 Ireland www.tcd.ie/Economics/staff/frainj/home.html mailto:fra...@tcd.ie mailto:fra...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in log
You still haven't provided anything reproducible since we can't get to your data file (also you used the-you-really-shouldn't-use-this-function function attach) but here's what I'd guess: Your say the problem occurs when exp(-alpha*d^(-beta)) 1 Of course, this happens when alpha*d^(-beta) 0. It seems like alpha , beta 0, so you need to ask yourself why d 0 : I can't immediately see why this would happen in your code. The easiest way to get at this it just insert a print(d) statement or, even better, a stopifnot(d 0) with options(error = recover). This will open an interactive debugger and you can do a post-mortem to see exactly what happened. This is really all I can see without minimal reproducible code, which you aren't providing. You may also want to work on vectorizing your code to make it faster for long simulations (and arguably easier to read, though it's a matter of taste) E.g., for(i in 1:100){ if(i!=k){ if(inftime[i]==0){ loglh-loglh +log(1-p[i]) } if(inftime[i]==2){ loglh-loglh + log(p[i]) } } } could become something like (very , very untested): loglh - loglh + sum((ifelse(inftime == 0, log(1-p), 0) + ifelse(inftime == 2, log(p), 0))[-k]) which will be much faster. This takes all the elements and vectorizes the two if statements, throws out the k case that you don't want, then sums them and adds to loglh. It should be much faster than your loop since it's all vectorized. Finally, within the loglikelihood() function, you don't need to pre-define d,p,k: it suffices to initialize them at the values you calculate. Hope something in here helps, but my bet is that your problem is coming from the unprovided data and/or something with attach. Michael On Wed, Nov 30, 2011 at 7:43 AM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: The loglikelihood() looks ok and gives some value. But I am using this function for the simulated annealing, generating the random samples from uniform proposal density., The codes are as follows epiannea - function(T0 = 1, N = 500,beta = 0.1,x0 = 0.1, rho = 0.90, eps = 0.1, loglikelihood, data){ moving - 1 count - 0 Temp - T0 alpha - x0 while(moving 0){ mmoving - 0 for(i in 1:N){ r - alpha + runif(1,-eps,eps) if(r 0){ a - min(1,exp((loglikelihood(alpha,beta)-loglikelihood(r,beta))/Temp)) } if(runif(1) a){ moving - moving +1 alpha - r } } Temp - Temp*rho count - count + 1 } #plot(a,loglikelihood(alpha,beta), type = l) fmin - loglikelihood(alpha, beta) return(c(fmin, alpha, count)) } epiannea(loglikelihood=loglikelihood) Some times it shows the warnings that logp[i] produces NaNs, that means p[i] is negative, but it should not be that because p[i] is the probabiliy and can't be negative. Some times the code runs but never stop. On Wed, Nov 30, 2011 at 7:34 AM, R. Michael Weylandt michael.weyla...@gmail.com michael.weyla...@gmail.com wrote: I'd suggest you do some leg-work and figure out why you are getting values 1. If your algorithm is motivated by some approximation then a min() or pmin() *might* be the right fix, but if there are no approximations you may need to start debugging properly to see why you are getting an out of bounds value. Since there's no random number generation involved, I'd hesitate to just throw out the result without knowing its source. Also keep in mind the limitations of floating point arithmetic if you expect alpha*d^beta to be small. Michael On Nov 29, 2011, at 6:58 PM, Sarah Goslee sarah.gos...@gmail.com wrote: On Tue, Nov 29, 2011 at 6:55 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: yes, log of negative number is undefined and R also do the same and produces NaNs. Here I want to reject the value of exp(-alpha*d^(-beta)) when greater than 1, and want to run the loop otherwise. Thanks Then just add another if() statement checking for that condition. On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com wrote: Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After generating the out put we can omit it, but the problem is different. Wait... you're complaining that you can't take the natural log of a negative number in R? You can't do that anywhere. What do you expect to happen? The log of a negative number IS NaN. Sarah On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: I have following code: loglikelihood - function(alpha,beta= 0.1){ loglh-0 d-0 p-0 k-NULL data-read.table(epidemic.txt,header = TRUE) attach(data, warn.conflicts = F) k -which(inftime==1) d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta) p-1 - exp(-alpha*d) for(i in 1:100){ if(i!=k){ if(inftime[i]==0){
[R] SAS to R: I would like to replicate a statistical analysis performed in SAS in R.
Hello everybody, A statistician performed an analysis in SAS for me which I would like to replicate in R. I have however problems in figuring out the R code to do that. As I understood it was a covariance regression model. In the analysis, baseline was used as covariate and autoregressive (1) as covariance structure. The model included baseline, session, group and interaction between baseline and group. The degree of freedom was corrected by Kenward-Roger method. Study design: - 2 Groups: group A: n=10 subjects, group B: n=10 - 5 motor test: baseline; learning sessions immediately, 6h, 24h, and 30d after memory task - dependent variable: response time (RT) Research question: Do the two groups learn differently across time? I would appreciate a lot any help on the R commands or the choice of statistical approach. Thanks! Marianne, Columbia, SC [[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] R help
Just throwing this out there: there's no probability distribution in the equation you gave, so there's no context for an MLE: what is the likelihood for m? Knowing a little bit about the NS model, it seems like it makes more sense to use nls() to fit all four parameters to the data at once. Something like this: nls(y ~ a + b*(1-exp(-m/tau))/(m/tau)+c*((1-exp(-m/tau))/(m/tau)-exp(-m/tau)), data = data.frame(y = y, tau = tau), start = list(a = 1, b = 1, c = 1, m = 1)) Note that as presented, this doesn't work, but it's hopefully a good start (I'm not very good with nls). However, the smartest thing to do is to do your homework before jumping into a project: library(YieldCurve) Nelson.Siegel(y, tau, c(48, 60, 84)) Michael On Tue, Nov 29, 2011 at 1:13 PM, Heh Ness imos...@yahoo.com wrote: I have a model like this (Nelson and Siegel 1987) img src=http://r.789695.n4.nabble.com/file/n4120161/31f188c684764cd431619dbb59fed5ae.png; border=0/ where tau and y are the maturities and yields, respectively, given to me in my data file.. y-c(4.863,5.662,6.41,6.864,7.153,7.352,7.409,7.474,7.503,7.644,7.676,7.701,7.674,7.668,7.665,7.741,7.743,7.742) tau-c(1,3,6,9,12,15,18,21,24,30,36,48,60,72,84,96,108,120) I firstly need to find the MLE of m which maximises the likelihood function and then I can easily find the b1, b2 and b3 constants using this m value via least squares estimation.. But does anybody know how I can go abouts finding the MLE of m and if you could help with providing r code for it, I would appreciate that a lot. I have been pulling my hair out for the past week now trying to do it :) [[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] nls help
Hello, I have data like the following: datum - structure(list(Y = c(415.5, 3847.8325, 1942.83325, 1215.223, 950.142857325, 2399.585, 804.75, 579.5, 841.70825, 494.053571425 ), X = c(1.081818182, 0.492727273, 0.756363636, 0.896363636, 1.518181818, 0.49917, 1.354545455, 1.61, 1.706363636, 1.063636364 )), .Names = c(Y, X), row.names = c(NA, -10L), class = data.frame) with(datum, plot(Y~X)) As you can see there is a non-linear association between X and Y, and I would like to fit an appropriate model. I was thinking an exponential decay model might work well. I tried the following (a and k starting values are based off of a lm() fit), but get an error. fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1867)) Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates I have never attempted to fit a non-linear model before, and thus the model may be inappropriately specified, or it is also possible that I have no idea what I am doing. Would someone please offer some advice. Thanks. Chuck -- View this message in context: http://r.789695.n4.nabble.com/nls-help-tp4123876p4123876.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help about fitdistr funtion
Hi, I have a variable X classified in a lot of groups and I need to run the [fitdistr] funtion for each group. I tried with the [by] or the [tapply] funtions because my data is organize in two columns (variable and the groups), but neither of these command work. If somebody have a tip to help me up I really appreciate it. thanks, [[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] Package field missing from repositories
Hi all; After an overnight Ubuntu upgrade to 10.10, I had to reinstall some R packages. But the fields package appears to be missing from my usual repository, and a few of the other repositories I've tested: install.packages('fields',lib='/usr/local/lib/R/site-library',repos='http://cran.stat.sfu.ca') Warning message: In getDependencies(pkgs, dependencies, available, lib) : package ‘fields’ is not available The package 'fields' *does* appear in the list of available packages at my local repository: http://cran.stat.sfu.ca/, and I've installed the dependencies ('spam'). Any suggestions? Thanks! -- Joseph Shea Post-Doctoral Fellow Departments of Geography University of British Columbia and University of Northern British Columbia Phone: 250-960-5840 Email: joseph.s...@geog.ubc.ca __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] install multtest and preprocessCore packages in Bioconductor library
Hi Nguyen, Subject: [R] install multtest and preprocessCore packages in Bioconductor library Date: Wed, 30 Nov 2011 09:57:36 -0800 From: UyenThao Nguyen ungu...@tethysbio.com To: r-help r-help@r-project.org CC: uth.ngu...@ucdavis.edu uth.ngu...@ucdavis.edu Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen You don't specify what commands you are using to install these packages. Bioconductor recommends installing packages with biocLite(), like so: source(http://bioconductor.org/biocLite.R;) biocLite(c(multtest, preprocessCore)) Please send the output of those commands, as well as the output of sessionInfo() so we can help solve your problem. Thanks Dan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Nomogram with stratified cph in rms package, how to get failure probability
Hello, I am using Dr. Harrell's rms package to make a nomogram. I was able to make a beautiful one. However, I want to change 5-year survival probability to 5-year failure probability. I couldn’t get hazard rate from Hazard(f1) because I used cph for the model. Here is my code: library(rms) f1 - cph(Surv(retime,dfs) ~ age+her2+t_stage+n_stage+er+grade+cytcyt+Cyt_PCDK2 , data=data11, surv=T, x=T, y=T, time.inc=5) surv- Survival(f1) haz- Hazard(f1) Here is the Error in UseMethod(Hazard) : no applicable method for 'Hazard' applied to an object of class c('cph', 'Design', 'coxph') surv10 - function(lp) surv(10,lp) surv5 - function(lp) surv(5,lp) quant - Quantile(f1) at.surv - c(0.1, 0.3, 0.5, 0.7, 0.9) at.med1-c(2,3,4, 5,6,7,8, 10,15,20,25, 30) par(cex=0.8) nom- nomogram(f1, conf.int=F, fun=list(surv5, surv10), funlabel=c('5-Year Survival Probability', '10-Year Survival Probability' ), lp=F, fun.at=c(at.surv, at.surv),label.every=1, force.label=FALSE, cex.axis=0.8, verbose=TRUE, cex.var=0.8) How can I show failure probability instead of survival probability? I would very much appreciate any assistance in this matter. Thank you Very much. -- View this message in context: http://r.789695.n4.nabble.com/Nomogram-with-stratified-cph-in-rms-package-how-to-get-failure-probability-tp4123249p4123249.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] How can I pick a matrix from a function? (Out Product of Gradient)
Hi all, I would like to use optim() to estimate the equation by the log-likelihood function and gradient function which I had written. I try to use OPG(Out Product of Gradient) to calculate the Hessian matrix since sometime Hessian matrix is difficult to calculate. Thus I want to pick the Gradient matrix from the gradient function. Moreover, could R show the process of calculation on gradient function I written by using optim( ) or other commands? Thanks, Yanghao X - cbind(rep(1,n),sex,age,yrmarry,children,rating) dy - (mydata$y0)*1 # * # Probit model: log-likelihood # * fprobit - *function*(beta,y,X) { n - length(y) k - ncol(X) b - beta[1:k] kk - 2*y-1 z - X %*% b L - log(pnorm(kk*z)) L - sum(L) return(-L) } # * # Probit model: analytic gradient # * gprobit - *function*(b,y,X) { kappa=2*y-1 z - kappa*(X %*% b) imr - kappa*dnorm(z)/pnorm(z) G-matrix(ncol=ncol(X),nrow=nrow(X)) *for*(i *in* 1:nrow(imr)){ G[i,]-imr[i,]*X[i,] } g - apply(G,2,sum) return(-g) } #** # How can I pick the G matrix from this function? #** For initial value# xx - cbind(sex,age,yrmarry,children,rating) reg1 - lm(dy~xx) (b0 - reg1$coef) ## (mle - optim(b0,fprobit,gr=gprobit, method=BFGS,hessian=T,y=dy,X=X,control=list(trace=T))) [[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] formula for calculating the survival probability for nomogram
Hi, I used Dr. Harrell's rms package to make a nomogram. Below is my code for nomogram and calculate total points and probability *in original data set* used for building the nomogram. *My question is how I get the formula for calculating the survival probability for this nomogram. Then I can use this formula to do validation by using other data set. * f1 - cph(Surv(retime,dfs) ~ age+her2+t_stage+n_stage+er+cytcyt+Cyt_PCDK2 , data=data11, surv=T, x=T, y=T, time.inc=5) surv- Survival(f1) surv10 - function(lp) surv(10,lp) surv5 - function(lp) surv(5,lp) quant - Quantile(f1) at.surv - c(0.1, 0.3, 0.5, 0.7, 0.9) nom- nomogram(f1, conf.int=F, fun=list(surv5, surv10), funlabel=c('5-Year Survival Probability', '10-Year Survival Probability' ), lp=F, fun.at=c(at.surv, at.surv)) # total points ## age.p - lm(nom[[1]]$points ~ nom[[1]]$age) tot.points - (age*age.p$coef[2]+age.p$coef[1])+ (her2=='Positive')*nom[[2]]$points[2] + (t_stage=='T2/T3')*nom[[3]]$points[1]+(n_stage=='N1/N2')*nom[[4]]$points[2]+(er=='Negative')*nom[[5]]$points[1] + (cytcyt=='Positive')*nom[[6]]$points[2]+(Cyt_PCDK2=='Positive')*nom[[7]]$points[2] summary(tot.points) data22 - data.frame(retime=retime, dfs=dfs, tot.points=tot.points) fit - cph(Surv(retime, dfs) ~ tot.points,data=data22 ,surv=T, x=T, y=T, time.inc=5) survff - survfit.cph(fit, newdata=data22) times -survff$time probb - survff$surv # 5 year ## indd - times==times[times=5][length(times[times=5])] probb2 - probb[indd,] ## 5-year surv prob for each patient I would very much appreciate any assistance in this matter. Thank you Very much. Min -- View this message in context: http://r.789695.n4.nabble.com/formula-for-calculating-the-survival-probability-for-nomogram-tp4124387p4124387.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] Problem installing the package tkrplot
Hello people. My problem is that when I try to install the package tkrplot, I got the next problem: install.packages(tkrplot) Installing package(s) into ‘/home/marcos/R/i686-pc-linux-gnu-library/2.13’ (as ‘lib’ is unspecified) probando la URL 'http://www.vps.fmvz.usp.br/CRAN/src/contrib/tkrplot_0.0-23.tar.gz' Content type 'application/x-gzip' length 39037 bytes (38 Kb) URL abierta == downloaded 38 Kb * installing *source* package ‘tkrplot’ ... configure: creating ./config.status config.status: creating src/Makevars ** libs gcc -I/usr/share/R/include -I/usr/include/tcl8.5 -I/usr/include/tcl8.5 -fpic -std=gnu99 -O3 -pipe -g -c tcltkimg.c -o tcltkimg.o tcltkimg.c:2:16: fatal error: tk.h: No existe el fichero o el directorio compilation terminated. make: *** [tcltkimg.o] Error 1 ERROR: compilation failed for package ‘tkrplot’ * removing ‘/home/marcos/R/i686-pc-linux-gnu-library/2.13/tkrplot’ Warning in install.packages : installation of package 'tkrplot' had non-zero exit status The downloaded packages are in ‘/tmp/RtmpHRfZ1k/downloaded_packages’ In this moment I have R version 2.13.2 (2011-09-30). I stayed thankfull for the any help. --- Marcos Amaris Gonzalez __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Bootstrapping for residuals
Hi, Ive been trying to write a program for bootstrapping residuals in R but without much success. A lecturer wants to predict the performance of students in an end-of-year physics exam, y. The lecturer has the students results from a mid-term physics exam, x, and a mid-term biology exam, z. He proposes the following linear model for the end-of-year exam result yi = α + βxi + γzi + qi, where q is the error. Y is a matrix of the data and we have y=first column of the data and X=second 2 columns(the x z data) Now I need to write a program for obtaining bootstrap estimates, i have: x=scan(data) Y=matrix(x,ncol=3,byrow=T) y=Y[,1] X=Y[,2:3] ls=lsfit(X,y) beta=ls$coef yest=beta[1]+beta[2]*X[,1]+beta[3]*X[,2] res=y-yest boot=function(X,res,beta,b) { n=24 output=matrix(0,ncol=2,nrow=b) for(i in 1:b) { error=sample(res,n,replace=T) ystar=beta[1]+beta[2]*X[,1]+beta[3]*X[,2]+error ls=lsfit(X,ystar) output[i,]=ls$coef } output } I think the first 8 lines are right but my function might be wrong? Any help? -- View this message in context: http://r.789695.n4.nabble.com/Bootstrapping-for-residuals-tp4123657p4123657.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] nls help
It's a scaling problem: If you do this: datum - datum[order(datum$X),] with(datum, plot(Y~X)) with(datum, lines(X, 3400*exp(-1867*X))) you'll see that your initial guess is just so far gone that the nls() optimizer can't handle it. If you try a more reasonable initial guess it works fine: fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1.867)) Michael On Wed, Nov 30, 2011 at 12:14 PM, chuck.01 charliethebrow...@gmail.com wrote: Hello, I have data like the following: datum - structure(list(Y = c(415.5, 3847.8325, 1942.83325, 1215.223, 950.142857325, 2399.585, 804.75, 579.5, 841.70825, 494.053571425 ), X = c(1.081818182, 0.492727273, 0.756363636, 0.896363636, 1.518181818, 0.49917, 1.354545455, 1.61, 1.706363636, 1.063636364 )), .Names = c(Y, X), row.names = c(NA, -10L), class = data.frame) with(datum, plot(Y~X)) As you can see there is a non-linear association between X and Y, and I would like to fit an appropriate model. I was thinking an exponential decay model might work well. I tried the following (a and k starting values are based off of a lm() fit), but get an error. fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1867)) Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates I have never attempted to fit a non-linear model before, and thus the model may be inappropriately specified, or it is also possible that I have no idea what I am doing. Would someone please offer some advice. Thanks. Chuck -- View this message in context: http://r.789695.n4.nabble.com/nls-help-tp4123876p4123876.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How can I pick a matrix from a function? (Out Product of Gradient)
Hi all, I would like to use optim() to estimate the equation by the log-likelihood function and gradient function which I had written. I try to use OPG(Out Product of Gradient) to calculate the Hessian matrix since sometime Hessian matrix is difficult to calculate. Thus I want to pick the Gradient matrix from the gradient function. Moreover, could R show the process of calculation on gradient function I written by using optim( ) or other commands? Thanks, Yanghao X - cbind(rep(1,n),sex,age,yrmarry,children,rating) dy - (mydata$y0)*1 # * # Probit model: log-likelihood # * fprobit - function(beta,y,X) { n - length(y) k - ncol(X) b - beta[1:k] kk - 2*y-1 z - X %*% b L - log(pnorm(kk*z)) L - sum(L) return(-L) } # * # Probit model: analytic gradient # * gprobit - function(b,y,X) { kappa=2*y-1 z - kappa*(X %*% b) imr - kappa*dnorm(z)/pnorm(z) G-matrix(ncol=ncol(X),nrow=nrow(X)) for(i in 1:nrow(imr)){ G[i,]-imr[i,]*X[i,] } g - apply(G,2,sum) return(-g) } For initial value# xx - cbind(sex,age,yrmarry,children,rating) reg1 - lm(dy~xx) (b0 - reg1$coef) ## (mle - optim(b0,fprobit,gr=gprobit, method=BFGS,hessian=T,y=dy,X=X,control=list(trace=T))) -- View this message in context: http://r.789695.n4.nabble.com/How-can-I-pick-a-matrix-from-a-function-Out-Product-of-Gradient-tp4123965p4123965.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] help about fitdistr funtion
I like to use split() to split your data into groups, then run lapply() to use fitdistr on each element of the list. E.g, df - data.frame(X = rnorm(500), ID = sample(letters[1:5], 500, TRUE)) temp - split(df$X, df$ID) lapply(temp, fitdistr, normal) Though it's just as easy with tapply(): tapply(df$X, df$ID, fitdistr, normal) Michael On Wed, Nov 30, 2011 at 10:00 AM, Carlos Javier Rincon Rodriguez cjrinc...@unal.edu.co wrote: Hi, I have a variable X classified in a lot of groups and I need to run the [fitdistr] funtion for each group. I tried with the [by] or the [tapply] funtions because my data is organize in two columns (variable and the groups), but neither of these command work. If somebody have a tip to help me up I really appreciate it. thanks, [[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] install multtest and preprocessCore packages in Bioconductor library
Thank you All for your prompt helps. Dear Dan, Before I used source(http://bioconductor.org/biocLite.R;) biocLite(Biobase) # GET BIOCONDUCTOR PACKAGE FROM SOURCE then, I went to Install packages in Rgui pointing to BioConductor to load these packages in. That was why I got the error messages. Now, I added biocLite(c(multtest,preprocessCore)) and saw that the packages were saved in /var/folders/9m/9mgg3bdhHcmLfSBpO7I6CTI/-Tmp-//Rtmpylgrz5/downloaded_packages, which I am not sure where. # biocLite(c(multtest,preprocessCore)) Using R version 2.10.1, biocinstall version 2.5.11. Installing Bioconductor version 2.5 packages: [1] multtest preprocessCore Please wait... trying URL 'http://cran.cnr.Berkeley.edu/bin/macosx/leopard/contrib/2.10/multtest_2.4.0.tgz' Content type 'application/x-gzip' length 1694933 bytes (1.6 Mb) opened URL == downloaded 1.6 Mb trying URL 'http://www.bioconductor.org/packages/2.5/bioc/bin/macosx/leopard/contrib/2.10/preprocessCore_1.8.0.tgz' Content type 'application/x-gzip' length 235382 bytes (229 Kb) opened URL == downloaded 229 Kb The downloaded packages are in /var/folders/9m/9mgg3bdhHcmLfSBpO7I6CTI/-Tmp-//Rtmpylgrz5/downloaded_packages Then, I did library(multtest) ?multtest No documentation for 'multtest' in specified packages and libraries: you could try '??multtest' I don't know how to load these functions after they were saved? I am not a computer person, so I don't know what to do when things don't work. Thanks again, Nguyen On Nov 30, 2011, at 10:09 AM, Dan Tenenbaum wrote: Hi Nguyen, Subject: [R] install multtest and preprocessCore packages in Bioconductor library Date: Wed, 30 Nov 2011 09:57:36 -0800 From: UyenThao Nguyen ungu...@tethysbio.commailto:ungu...@tethysbio.com To: r-help r-help@r-project.orgmailto:r-help@r-project.org CC: uth.ngu...@ucdavis.edumailto:uth.ngu...@ucdavis.edu uth.ngu...@ucdavis.edumailto:uth.ngu...@ucdavis.edu Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package preprocessCore * removing /Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore * restoring previous /Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen You don't specify what commands you are using to install these packages. Bioconductor recommends installing packages with biocLite(), like so: source(http://bioconductor.org/biocLite.R;) biocLite(c(multtest, preprocessCore)) Please send the output of those commands, as well as the output of sessionInfo() so we can help solve your problem. Thanks Dan __ R-help@r-project.orgmailto:R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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.orgmailto:R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] install multtest and preprocessCore packages in Bioconductor library
Hi Nguyen, On Wed, Nov 30, 2011 at 12:11 PM, UyenThao Nguyen ungu...@tethysbio.com wrote: Thank you All for your prompt helps. Dear Dan, Before I used source(http://bioconductor.org/biocLite.R;) biocLite(Biobase) # GET BIOCONDUCTOR PACKAGE FROM SOURCE then, I went to Install packages in Rgui pointing to BioConductor to load these packages in. That was why I got the error messages. Now, I added biocLite(c(multtest,preprocessCore)) and saw that the packages were saved in /var/folders/9m/9mgg3bdhHcmLfSBpO7I6CTI/-Tmp-//Rtmpylgrz5/downloaded_packages, which I am not sure where. # biocLite(c(multtest,preprocessCore)) Using R version 2.10.1, biocinstall version 2.5.11. Installing Bioconductor version 2.5 packages: [1] multtest preprocessCore Please wait... trying URL 'http://cran.cnr.Berkeley.edu/bin/macosx/leopard/contrib/2.10/multtest_2.4.0.tgz' Content type 'application/x-gzip' length 1694933 bytes (1.6 Mb) opened URL == downloaded 1.6 Mb trying URL 'http://www.bioconductor.org/packages/2.5/bioc/bin/macosx/leopard/contrib/2.10/preprocessCore_1.8.0.tgz' Content type 'application/x-gzip' length 235382 bytes (229 Kb) opened URL == downloaded 229 Kb The downloaded packages are in /var/folders/9m/9mgg3bdhHcmLfSBpO7I6CTI/-Tmp-//Rtmpylgrz5/downloaded_packages Then, I did library(multtest) The fact that this produced no error means that multtest is successfully installed. ?multtest No documentation for 'multtest' in specified packages and libraries: you could try '??multtest' To get help, do this: help(package=multtest) I don't know how to load these functions after they were saved? I am not a computer person, so I don't know what to do when things don't work. You have already loaded the package; library(packagename) is the way to do it. This and the help() bit above should get you started. If you have further questions about these specific packages, please post to the bioconductor list, after reading our posting guide. http://bioconductor.org/help/mailing-list/posting-guide/ Thanks, Dan Thanks again, Nguyen On Nov 30, 2011, at 10:09 AM, Dan Tenenbaum wrote: Hi Nguyen, Subject: [R] install multtest and preprocessCore packages in Bioconductor library Date: Wed, 30 Nov 2011 09:57:36 -0800 From: UyenThao Nguyen ungu...@tethysbio.com To: r-help r-help@r-project.org CC: uth.ngu...@ucdavis.edu uth.ngu...@ucdavis.edu Hi All, I've tried to install these multtest and preprocessCore packages in Mac, but kept getting error messages. I tried to load the packages using 2 ways: 1. Installed from BioConductor (sources) 2. And installed from BioConductor (binaries) Both ways, I got these error messages: For preprocessCore: ld: warning: directory '/usr/local/lib' following -L not found ld: library not found for -lgfortran collect2: ld returned 1 exit status make: *** [preprocessCore.so] Error 1 ERROR: compilation failed for package ‘preprocessCore’ * removing ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ * restoring previous ‘/Library/Frameworks/R.framework/Versions/2.10/Resources/library/preprocessCore’ For multtest: ld: warning: directory '/usr/local/lib' following -L not found ** arch - x86_64 Looks like the error got to do with L drive, which I don't know how to fix it. Please help. Any suggestions will be greatly appreciated. Thank you very much, Nguyen You don't specify what commands you are using to install these packages. Bioconductor recommends installing packages with biocLite(), like so: source(http://bioconductor.org/biocLite.R;) biocLite(c(multtest, preprocessCore)) Please send the output of those commands, as well as the output of sessionInfo() so we can help solve your problem. Thanks Dan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] SAS to R: I would like to replicate a statistical analysis performed in SAS in R.
Marianne Stephan mariannestephan at hotmail.com writes: A statistician performed an analysis in SAS for me which I would like to replicate in R. I have however problems in figuring out the R code to do that. As I understood it was a covariance regression model. In the analysis, baseline was used as covariate and autoregressive (1) as covariance structure. The model included baseline, session, group and interaction between baseline and group. The degree of freedom was corrected by Kenward-Roger method. Study design: - 2 Groups: group A: n=10 subjects, group B: n=10 - 5 motor test: baseline; learning sessions immediately, 6h, 24h, and 30d after memory task - dependent variable: response time (RT) Research question: Do the two groups learn differently across time? I would appreciate a lot any help on the R commands or the choice of statistical approach. Thanks! Marianne, Columbia, SC [[alternative HTML version deleted]] You would be better off submitting this to the r-sig-mixed-models mailing list. If possible, posting your data somewhere, or a reasonable facsimile of it, would improve your chances of getting an answer. Kenward-Roger has not been available in R until recently (see a long thread referenced at http://glmm.wikidot.com/faq on degrees of freedom, but has recently been implemented in (an experimental version of?) the 'doBy' package. You might also get more answers if you show what you have tried so far. Ben Bolker __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Bootstrapping for residuals
On Nov 30, 2011, at 11:31 AM, bubbles1990 wrote: Hi, Ive been trying to write a program for bootstrapping residuals in R but without much success. A lecturer wants to predict the performance of students in an end-of- year physics exam, y. The lecturer has the students results from a mid-term physics exam, x, and a mid-term biology exam, z. He proposes the following linear model for the end-of-year exam result yi = α + βxi + γzi + qi, where q is the error. Y is a matrix of the data and we have y=first column of the data and X=second 2 columns(the x z data) Now I need to write a program for obtaining bootstrap estimates, i have: x=scan(data) Y=matrix(x,ncol=3,byrow=T) y=Y[,1] X=Y[,2:3] ls=lsfit(X,y) beta=ls$coef yest=beta[1]+beta[2]*X[,1]+beta[3]*X[,2] res=y-yest boot=function(X,res,beta,b) { n=24 output=matrix(0,ncol=2,nrow=b) for(i in 1:b) { error=sample(res,n,replace=T) ystar=beta[1]+beta[2]*X[,1]+beta[3]*X[,2]+error ls=lsfit(X,ystar) output[i,]=ls$coef } output } I think the first 8 lines are right but my function might be wrong? Any help? You should ask your instructor or teaching assistant for help. -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Package field missing from repositories
On Nov 30, 2011, at 1:15 PM, joseph.s...@geog.ubc.ca wrote: Hi all; After an overnight Ubuntu upgrade to 10.10, I had to reinstall some R packages. But the fields package appears to be missing from my usual repository, and a few of the other repositories I've tested: install.packages('fields',lib='/usr/local/lib/R/site- library',repos='http://cran.stat.sfu.ca') Warning message: In getDependencies(pkgs, dependencies, available, lib) : package ‘fields’ is not available The package 'fields' *does* appear in the list of available packages at my local repository: http://cran.stat.sfu.ca/, and I've installed the dependencies ('spam'). Any suggestions? Maybe you should post a bit more about your setup, as the Posting Guide suggests. -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls help
Ah, I see, thank you both. Michael Weylandt wrote It's a scaling problem: If you do this: datum - datum[order(datum$X),] with(datum, plot(Y~X)) with(datum, lines(X, 3400*exp(-1867*X))) you'll see that your initial guess is just so far gone that the nls() optimizer can't handle it. If you try a more reasonable initial guess it works fine: fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1.867)) Michael On Wed, Nov 30, 2011 at 12:14 PM, chuck.01 lt;CharlieTheBrown77@gt; wrote: Hello, I have data like the following: datum - structure(list(Y = c(415.5, 3847.8325, 1942.83325, 1215.223, 950.142857325, 2399.585, 804.75, 579.5, 841.70825, 494.053571425 ), X = c(1.081818182, 0.492727273, 0.756363636, 0.896363636, 1.518181818, 0.49917, 1.354545455, 1.61, 1.706363636, 1.063636364 )), .Names = c(Y, X), row.names = c(NA, -10L), class = data.frame) with(datum, plot(Y~X)) As you can see there is a non-linear association between X and Y, and I would like to fit an appropriate model. I was thinking an exponential decay model might work well. I tried the following (a and k starting values are based off of a lm() fit), but get an error. fit - nls(Y ~ a*exp(-k * X), datum, start=c(a=3400, k=1867)) Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates I have never attempted to fit a non-linear model before, and thus the model may be inappropriately specified, or it is also possible that I have no idea what I am doing. Would someone please offer some advice. Thanks. Chuck -- View this message in context: http://r.789695.n4.nabble.com/nls-help-tp4123876p4123876.html Sent from the R help mailing list archive at Nabble.com. __ R-help@ mailing list https://stat.ethz.ch/mailman/listinfo/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@ mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://r.789695.n4.nabble.com/nls-help-tp4123876p4124772.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] Error adding Bigmemory package
I am trying to add the bigmemory packages but I get the following error message: In file included from bigmemory.cpp:14:0: ../inst/include/bigmemory/isna.hpp: In function 'bool neginf(double)': ../inst/include/bigmemory/isna.hpp:22:57: error: 'isinf' was not declared in this scope make: *** [bigmemory.o] Error 1 ERROR: compilation failed for package 'bigmemory' * removing '/usr/local/lib/R/library/bigmemory' * restoring previous '/usr/local/lib/R/library/bigmemory' The downloaded packages are in '/tmp/RtmpNrkt2a/downloaded_packages' Updating HTML index of packages in '.Library' Warning message: In install.packages(bigmemory) : installation of package 'bigmemory' had non-zero exit status install.packages(gsl) Error in install.packages(gsl) : object 'gsl' not found install.packages(gnulib) Error in install.packages(gnulib) : object 'gnulib' not found I am using the GNU compiler version 4.5.1 from blastwave on a Solaris 10 box. I can't seem to find how to compile bigmembory successfully. I can complile and add other packages. roger -- View this message in context: http://r.789695.n4.nabble.com/Error-adding-Bigmemory-package-tp4124863p4124863.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] Calculating marginal effects in R
The relevant code is contained in the following working paper: http://ideas.repec.org/p/ucn/wpaper/201122.html -- View this message in context: http://r.789695.n4.nabble.com/Calculating-marginal-effects-in-R-tp887488p4124996.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] Bootstrapping for residuals
I study a part-time long distance learning course so only really have access to online sources for help, i'm sure ill crack it eventually. I understand how to do it with 2 variables, its just having 3 that is confusing me -- View this message in context: http://r.789695.n4.nabble.com/Bootstrapping-for-residuals-tp4123657p4125325.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] forecasting linear regression from lagged variable
I'm currently working with some time series data with the xts package, and would like to generate a forecast 12 periods into the future. There are limited observations, so I am unable to use an ARIMA model for the forecast. Here's the regression setup, after converting everything from zoo objects to vectors. hire.total.lag1 - lag(hire.total, lag=-1, na.pad=TRUE) lm.model - lm(hire.total ~ case.total + hire.total.lag1) hire.total is a monthly historical time series from Jan 2010 to Oct 2011. hire.total.lag1 is the same time series, lagged 1 period backwards. case.total is a monthly historical time series from Jan 2010 to Oct 2011, plus forecasts forward another 12 periods. I'd like to use this model to forecast hire.total for the next period, and use each successive prediction of hire.total as the lag1 observation for the next prediction. I have enough observed values for case.total to forecast out as far as I need. I might be able to construct this using a loop, but I have a feeling it will be messy and slow. Any suggestions? -- View this message in context: http://r.789695.n4.nabble.com/forecasting-linear-regression-from-lagged-variable-tp4125091p4125091.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] forecasting linear regression from lagged variable
hi: that model is what the political scientists called an LDV ( lagged dependent variable model ), the marketing people call a koyck distributed model and the economists call it an ADL(1,0) or ADL(0,1). ( i forget ). You have to be careful about the error term because, if it's not white noise, then the coefficient on hire.total.lag1 is biased. ( there's a whole literature on this. google for regression with lagged dependent variable and a lot of things should come up ). In that case, you need to build the likelihood and use optim or optimx because lm is not the correct approach. as far as prediction, most people assume an impulse at x ( you're case.total variable ) and zero afterwards but, since you have them and don't need to assume an implulse, then I would just do the prediction by brute force as you were planning to do. let me know off line if you want papers about below from many different perspectives ( econ, marketing, poli sci ) and also, send it to R-Sig-Finance because it's more relevant there so more people will respond probably. but apologize for cross-posting. mark On Wed, Nov 30, 2011 at 4:40 PM, AaronB aa...@communityattributes.comwrote: I'm currently working with some time series data with the xts package, and would like to generate a forecast 12 periods into the future. There are limited observations, so I am unable to use an ARIMA model for the forecast. Here's the regression setup, after converting everything from zoo objects to vectors. hire.total.lag1 - lag(hire.total, lag=-1, na.pad=TRUE) lm.model - lm(hire.total ~ case.total + hire.total.lag1) hire.total is a monthly historical time series from Jan 2010 to Oct 2011. hire.total.lag1 is the same time series, lagged 1 period backwards. case.total is a monthly historical time series from Jan 2010 to Oct 2011, plus forecasts forward another 12 periods. I'd like to use this model to forecast hire.total for the next period, and use each successive prediction of hire.total as the lag1 observation for the next prediction. I have enough observed values for case.total to forecast out as far as I need. I might be able to construct this using a loop, but I have a feeling it will be messy and slow. Any suggestions? -- View this message in context: http://r.789695.n4.nabble.com/forecasting-linear-regression-from-lagged-variable-tp4125091p4125091.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem installing the package tkrplot
On Nov 30, 2011, at 15:31 , Marcos Amaris Gonzalez wrote: Hello people. My problem is that when I try to install the package tkrplot, I got the next problem: You seem to be missing the tk.h file. This is, most likely, in a tcl/tk development package that you haven't installed. This is part of your Linux distribution, not R. (You're not telling us which Linux you are running, and it can't be pinpointed any further without that information, but look for something like tk-devel and tcl-devel.) install.packages(tkrplot) Installing package(s) into ‘/home/marcos/R/i686-pc-linux-gnu-library/2.13’ (as ‘lib’ is unspecified) probando la URL 'http://www.vps.fmvz.usp.br/CRAN/src/contrib/tkrplot_0.0-23.tar.gz' Content type 'application/x-gzip' length 39037 bytes (38 Kb) URL abierta == downloaded 38 Kb * installing *source* package ‘tkrplot’ ... configure: creating ./config.status config.status: creating src/Makevars ** libs gcc -I/usr/share/R/include -I/usr/include/tcl8.5 -I/usr/include/tcl8.5 -fpic -std=gnu99 -O3 -pipe -g -c tcltkimg.c -o tcltkimg.o tcltkimg.c:2:16: fatal error: tk.h: No existe el fichero o el directorio compilation terminated. make: *** [tcltkimg.o] Error 1 ERROR: compilation failed for package ‘tkrplot’ * removing ‘/home/marcos/R/i686-pc-linux-gnu-library/2.13/tkrplot’ Warning in install.packages : installation of package 'tkrplot' had non-zero exit status The downloaded packages are in ‘/tmp/RtmpHRfZ1k/downloaded_packages’ In this moment I have R version 2.13.2 (2011-09-30). I stayed thankfull for the any help. --- Marcos Amaris Gonzalez __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Bootstrapping for residuals
On Nov 30, 2011, at 11:31 AM, bubbles1990 wrote: Hi, Ive been trying to write a program for bootstrapping residuals in R but without much success. A lecturer wants to predict the performance of students in an end-of- year physics exam, y. The lecturer has the students results from a mid-term physics exam, x, and a mid-term biology exam, z. He proposes the following linear model for the end-of-year exam result yi = α + βxi + γzi + qi, where q is the error. Y is a matrix of the data and we have y=first column of the data and X=second 2 columns(the x z data) Now I need to write a program for obtaining bootstrap estimates, You replied but included no context, a common and not appreciated (in both senses of that verb) failing in posts from Nabble. My question would be: bootstrap estimate ... of what? -- David. i have: x=scan(data) Y=matrix(x,ncol=3,byrow=T) y=Y[,1] X=Y[,2:3] ls=lsfit(X,y) beta=ls$coef yest=beta[1]+beta[2]*X[,1]+beta[3]*X[,2] res=y-yest boot=function(X,res,beta,b) { n=24 output=matrix(0,ncol=2,nrow=b) for(i in 1:b) { error=sample(res,n,replace=T) ystar=beta[1]+beta[2]*X[,1]+beta[3]*X[,2]+error ls=lsfit(X,ystar) output[i,]=ls$coef } output } I think the first 8 lines are right but my function might be wrong? Any help? -- View this message in context: http://r.789695.n4.nabble.com/Bootstrapping-for-residuals-tp4123657p4123657.html Sent from the R help mailing list archive at Nabble.com. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] forecasting linear regression from lagged variable
On Wed, Nov 30, 2011 at 4:40 PM, AaronB aa...@communityattributes.com wrote: I'm currently working with some time series data with the xts package, and would like to generate a forecast 12 periods into the future. There are limited observations, so I am unable to use an ARIMA model for the forecast. Here's the regression setup, after converting everything from zoo objects to vectors. hire.total.lag1 - lag(hire.total, lag=-1, na.pad=TRUE) lm.model - lm(hire.total ~ case.total + hire.total.lag1) hire.total is a monthly historical time series from Jan 2010 to Oct 2011. hire.total.lag1 is the same time series, lagged 1 period backwards. case.total is a monthly historical time series from Jan 2010 to Oct 2011, plus forecasts forward another 12 periods. I'd like to use this model to forecast hire.total for the next period, and use each successive prediction of hire.total as the lag1 observation for the next prediction. I have enough observed values for case.total to forecast out as far as I need. I might be able to construct this using a loop, but I have a feeling it will be messy and slow. Any suggestions? Its not clear that you can't use an arima model. You would use the n.ahead= argument of predict.Arima. The dyn and dynlm packages handle lagged lm's of zoo objects but you will have to do a predict and then append the prediction to the data and repeat in a loop as you describe. If you want to be careful about error terms then its more complex. -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Invalid number of components, ncomp
I'm not familiar with the package and I don't currently have it loaded, but looking at the examples in the documentation here http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/pls.pcr/html/mvr.html it seems that you need ncomp = 1:6. Saying you want to do multivariate regression on only a single variable is probably the source of the error. Does that help? Michael On Wed, Nov 30, 2011 at 5:10 AM, Vytautas Rakeviius vytautas1...@yahoo.com wrote: Error in mvr(Kd_nM ~ qsar, ncomp = 6, data = my, validation = CV, method = kernelpls) : Invalid number of components, ncomp How I can fix this? [[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 to call a function for each row
if you want to store the result in a column of your data.frame: within(df, Y - 6*X1+7*X2+8*X3) On Wed, Nov 30, 2011 at 9:59 AM, David L Carlson dcarl...@tamu.edu wrote: Isn't this even easier? X1 - c(1:3) X2 - c(3, 4, 6) X3 - c(5, 6, 1) Y - 6*X1 + 7*X2 + 8*X3 Y [1] 67 88 68 Or if you really need a function: MakeY - function(x, y, z) 6*x + 7*y + 8*z MakeY(X1, X2, X3) [1] 67 88 68 -- David L Carlson Associate Professor of Anthropology Texas AM University College Station, TX 77843-4352 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of R. Michael Weylandt Sent: Wednesday, November 30, 2011 8:17 AM To: arunkumar Cc: r-help@r-project.org Subject: Re: [R] how to call a function for each row Read ?apply This is easiest: df - matrix(c(1,2,3,3,4,6,5,6,1), 3) apply(df, 1, function(x) 6*x[1]+7*x[2]+8*x[3]) But it's much more efficient to do it with matrix multiplication. In keeping with the best of tradition, this is left as an exercise to the reader. Michael On Wed, Nov 30, 2011 at 8:10 AM, arunkumar akpbond...@gmail.com wrote: Hi I have a data-frame which look like this X1 X2 X3 1 3 5 2 4 6 3 6 1 I want to apply a formula Y=6*X1 + 7*X2 + 8*X3 for every row Thanks in Advance -- View this message in context: http://r.789695.n4.nabble.com/how-to-call-a-function-for-each-row-tp41 22906p4122906.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] FIML with missing data in sem package
Is there a way to use full information maximum likelihood (FIML) to estimate missing data in the sem package? For example, suppose I have a dataset with complete information on X1-X3, but missing data (MAR) on X4. Is there a way to use FIML in this case? I know lavaan and openmx allow you to do it, but I couldn't find anything in the documentation for the sem package. Thanks! -- Dustin Fife Graduate Student Quantitative Psychology University of Oklahoma [[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] Random Forests in R
I understand the original implementation of Random Forest was done in Fortran code. In the source files of the R implementation there is a note C wrapper for random forests: get input from R and drive the Fortran routines.. I'm far from an expert on this...does that mean that the implementation in R is through calls to C functions only (not Fortran)? So, would knowing C be enough to understand this code, or Fortran is also necessary? Thanks for your help Axel. [[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] R logo in eps formt
G'day all, Sorry if this message has been posted before, but searching for R is always difficult... I was hoping for a copy of the logo in eps format? Can I do this from R, or is one available for download? cheers Ben __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] R endnote entry
I know citation() gives the R citation to be used in publications. Has anyone put this into endnote nicely? I'm not very experienced with endnote, and the way I have it at the momeny the 'R Development Core Team' becomes R. D. C. T. etc. Cheers. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R endnote entry
I have R entered as a computer program, and in the Programmer Name field I write it as: R Development Core Team, Notice the comma after Team. That seems to be the key to getting Endnote to treat the whole thing as a surname that doesn't get abbreviated. On Wed, Nov 30, 2011 at 10:56 PM, Matt Cooper mattcst...@gmail.com wrote: I know citation() gives the R citation to be used in publications. Has anyone put this into endnote nicely? I'm not very experienced with endnote, and the way I have it at the momeny the 'R Development Core Team' becomes R. D. C. T. etc. Cheers. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ___ Luke Miller Postdoctoral Researcher Marine Science Center Northeastern University Nahant, MA (781) 581-7370 x318 [[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] Random Forests in R
On Wed, Nov 30, 2011 at 7:48 PM, Axel Urbiz axel.ur...@gmail.com wrote: I understand the original implementation of Random Forest was done in Fortran code. In the source files of the R implementation there is a note C wrapper for random forests: get input from R and drive the Fortran routines.. I'm far from an expert on this...does that mean that the implementation in R is through calls to C functions only (not Fortran)? So, would knowing C be enough to understand this code, or Fortran is also necessary? I haven't seen the C and Fortran code for Random Forest but I understand the note to say that R code calls some C functions that pre-process (possibly re-format etc) the data, then call the actual Random Forest method that's written in Fortran, then possibly post-process the output and return it to R. It would imply that to understand the actual Random Forest code, you will have to read the Fortran source code. Best, Peter __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] RTisean executable problem in STLperArea
Hi, I’m trying to use the STLperArea function in the ndvits package. I can run the sample data (“SLPSAs_full”) without any problem. However, when I come to run the function on my own data, I get the following message. Waiting to confirm page change... Error in .checkPath(path) : no TISEAN executables found in that directory. Please set a proper TISEAN executables path using 'setTISEANpath' I’ve tried setting the TISEANpath to the directory where the RTisean library is on my computer. Could someone please suggest what the problem might be? Many thanks, Louise -- View this message in context: http://r.789695.n4.nabble.com/RTisean-executable-problem-in-STLperArea-tp4125382p4125382.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] R code
Hi everybody, I am unable to resolve this error using the following for loop. Would appreciate help. The same loop works with for(i in 1:92) strangely. I checked the .raw input file and everything is kosher, even Line 547 mentioned in the error message. I wonder if there is any problem with the paste function. Thanks very much in advance. ** for(i in 1:93) { inputdta- read.table(file=gsub( ,,paste(JSR_network_[,i,].RAW), fixed=TRUE), header=TRUE) colnames(inputdta) - c(jsrid,firmid,jsrstart,injsr) print(nrow(inputdta)) } *** Message: Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : line 547 did not have 4 elements -- View this message in context: http://r.789695.n4.nabble.com/R-code-tp4125904p4125904.html Sent from the R help mailing list archive at Nabble.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] Average of Two Matrices
What if you have over 50 matrices and you don't want to write them all out one-by-one? I know there's something really quite simple, but I haven't found it yet!... -- View this message in context: http://r.789695.n4.nabble.com/Average-of-Two-Matrices-tp860672p4126615.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] error in ZIGP package and optim
Dear all: I am trying to use the response.zigp and est.zigp from the ZIGP package. Firstly, I generated a simple data via response.zigp and then I try to fit them by est.zigp function, this is the code: - library(ZIGP) #try-ZIGP-1 n-100 X-matrix(c(rep(1,n),runif(n*2,0,1)),n,3) W-matrix(c(rep(1,n),rbinom(n,1,0.5)),n,2) Z-matrix(c(rep(1,n),rbinom(n,1,0.5)),n,2) beta-c(1,-2,1) alpha-c(1,1.5) gamma-c(-1,2) y-response.zigp(X,W,Z,beta,alpha,gamma) x1-as.vector(X[,1]) x2-as.vector(X[,2]) x3-as.vector(X[,3]) w1-as.vector(W[,1]) w2-as.vector(W[,2]) z1-as.vector(Z[,1]) z2-as.vector(Z[,2]) fm.X-~x1+x2+x3 fm.W-~1 fm.Z-~z1+z2 est.zigp(y,fm.X,fm.W,fm.Z,Offset=rep(1,length(y)),init=FALSE) - Then I get the *error:* Error in optim(par = start.delta, fn = loglikelihood.zigp, gr = gradient, : initial value in 'vmmin' is not finite Any help you can give me would be greatly appreciated. -- View this message in context: http://r.789695.n4.nabble.com/error-in-ZIGP-package-and-optim-tp4127053p4127053.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] ordinal data simulation
Dear all, I am doing an ordinal data simulation. I have a question regarding the cut off values between categories. Suppose I have three categories, if I do regression, there should be two cut off values. I find some simulation code for the ordinal data. However, they usually only generate a uniform random number and compare the probability with the random number. Could any one be so kind to tell me the rational behind it. I don't know why we need a random number rather than just two fixed cut off values. Thanks very much! Myckel __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Average of Two Matrices
Well, you're rather stuck writing them all out unless you have them in some other data structure. ## three dimensional array (50 4 x 4 matrices) x - array(rnorm(16 * 50), dim = c(4, 4, 50)) apply(x, 1:2, mean) [,1][,2] [,3][,4] [1,] -0.09460574 0.01572077 -0.166194625 -0.28176540 [2,] -0.26649044 -0.23546332 -0.003559647 -0.09857819 [3,] 0.12657883 0.10979631 0.368414189 -0.13653997 [4,] -0.02827496 0.07634844 -0.088057890 0.06208034 ## list of length 50 each element is 4 x 4 matrix y - rep(list(matrix(rnorm(16), 4, 4)), 50) length(y) [1] 50 Reduce(`+`, y)/length(y) [,1][,2][,3] [,4] [1,] -1.8334660 0.67812999 0.20159314 1.6163501 [2,] -0.6716821 -0.23942474 -0.50482638 -0.5765309 [3,] 0.7991496 0.38674047 -0.02128386 -1.1868209 [4,] -1.1654429 0.06773386 0.13538268 0.3847375 See ?array ?apply ?Reduce for documentation Hope this helps, Josh On Wed, Nov 30, 2011 at 7:22 PM, T.D. Rudolph tylerdrudo...@gmail.com wrote: What if you have over 50 matrices and you don't want to write them all out one-by-one? I know there's something really quite simple, but I haven't found it yet!... -- View this message in context: http://r.789695.n4.nabble.com/Average-of-Two-Matrices-tp860672p4126615.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Ph.D. Student, Health Psychology Programmer Analyst II, ATS Statistical Consulting Group University of California, Los Angeles https://joshuawiley.com/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Generalized singular value decomposition
Am 30.11.2011 16:13, schrieb Berend Hasselman: tom wrote I don't think you have to install something. You have not told us how you are trying to load the dll's. You can use the GSVD function from http://r.789695.n4.nabble.com/problem-with-generalized-singular-value-decomposition-using-LAPACK-td834930.html but you will need the exact path to the Lapack dll in the dyn.load call. Something like dyn.load(C:\\R\\.\\libRlapack.dll). This was just what I was looking for, thank you very much! When I read the post before I did not notice the line with the dyn.load. You do not need to dyn.load the Blas dll. I have tried this GSVD function with Mac OS X with success. I tried the GSVD with the folowing matrices: A = matrix(rnorm(100,100,20),10,10) B = matrix(rnorm(200,10,3),20,10) res = GSVD(A,B) The resulted U and V have all entries equal to zero. This seems a little wired since A=U*E1*Q' and B=V*E2*Q'. If U and V are zero then A and B should also be zero. Are my example matrices a special case for the decomposition? I just started working with it, so I am a completely newbie ... Thank you again very much for your help! (I was already thinking about doing the decomposition with Matlab and loading the result into R.) Regards, Oana Berend -- View this message in context: http://r.789695.n4.nabble.com/Generalized-singular-value-decomposition-tp4122542p4123335.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] RTisean executable problem in STLperArea
It is one of the rather confusing things about the RTisean package. It does not say what your subject says: it says 'TISEAN executables'. The DESCRIPTION *should* have a SystemRequirements line (and this has been reported), It does say It requires that you have the Tisean-3.0.1 algorithms available on your computer. (It means the Tisean-3.0.1 binary programs ) Have you? If so, that is what you should be seting setTISEANpath to point to. On Wed, 30 Nov 2011, lglew wrote: Hi, I’m trying to use the STLperArea function in the ndvits package. I can run the sample data (“SLPSAs_full”) without any problem. However, when I come to run the function on my own data, I get the following message. Waiting to confirm page change... Error in .checkPath(path) : no TISEAN executables found in that directory. Please set a proper TISEAN executables path using 'setTISEANpath' I’ve tried setting the TISEANpath to the directory where the RTisean library is on my computer. Could someone please suggest what the problem might be? Many thanks, Louise -- View this message in context: http://r.789695.n4.nabble.com/RTisean-executable-problem-in-STLperArea-tp4125382p4125382.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. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595__ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.