Re: [R] Odp: Error in matrix, not ordered vectors or numerical value, and SIAR.
Thank-you Petr. I have consulted as many manuals and help pages, and search engines, and trying various things in SIAR, but continuous errors prevail.. I have tried changing all matrices to numeric, but the first column , of the two latter matrices, returns as NA. Any other suggestions? I had put in one of these messages more detail about data, see below. pond - c(A3, A3, A3, B3, B3, B3,C2, A7, A7) d13C - c(-22.07, -20.39, -21.11...) d15N - c(16.06, 17.13, 16.7, ...) FConsumerData - data.frame(pond, d13C, d15N) FTEFData - data.frame(tefsources, tefMeand13C, tefSDd13C, tefMeand15N, tefSDd15N) sources - c(Copepoda, Algal Paste, Amphipoda, Meand13C - c(-26.57, -20.41, -21.65,... SDd13C - c(1.78, 4.01, 1.46,... Meand15N - c(15.55, 9.58, 14.07... SDd15N - c(1.77, 1.34, 1.66,... FSourceData - data.frame(sources, Meand13C, SDd13C, Meand15N, SDd15N) tefsources - c(Copepoda, Algal Paste, Amphipoda,... tefMeand13C - c(0.4, 0.4, 0.4,... tefSDd13C - c(1.3, 1.3, 1.3, tefMeand15N - c(3.4, 3.4, 3.4, tefSDd15N - c(1.0, 1.0, 1.0, Hope that helps explain? ?? -- View this message in context: http://r.789695.n4.nabble.com/Error-in-matrix-not-ordered-vectors-or-numerical-value-and-SIAR-tp4024578p4030650.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] Errors in SIAR
On Nov 11, 2011 at 1:06am Alex DJ wrote: Alex, I haven't followed this thread closely (now split), but this most recent request for help from you is very difficult to make sense of. I think we need access to your data set. Further, there appear to be problems with the function you are using in package siar[1], based on trying to run demos in the package. (Though this may not be the issue.) Two things you can try are the following: First, the function you are using (siarmcmcdirichletv4) calls for matrices, not data frames (which you are using). Usually this throws an error, so maybe this is not the problem. Secondly, if you have a factor then you turn it into an ordered factor by doing ordered(myFactor) where myFactor is the name of the factor you want to convert to an ordered factor. You can call the levels of your factor what you like. Here, reading up on these functions, and trying the examples, would help you to help yourself. Have you even looked at ?factor. The fact that you are new to R really means that you should be reading everything about R that you can lay your hands on. This will help you to ask sensible questions about things that still puzzle you. And if you want to convert a factor into a numeric then you have to do something like as.numeric(as.character(myFactor)) This is all documented in R. Have you read the manualAn Introduction to R, which comes with R? There is also plenty of good documentation on how to use R on the web and indeed on R`s homepage. Regards, Mark. [1] R is case-sensitive: the package is called siar, not SIAR. Please respect the author's designation. - Mark Difford (Ph.D.) Research Associate Botany Department Nelson Mandela Metropolitan University Port Elizabeth, South Africa -- View this message in context: http://r.789695.n4.nabble.com/Errors-in-SIAR-tp4029804p4030667.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] Have write.csv write csv without observation number
Hello all, Have a little annoying issue with wreite.csv. It always writes out observation number. Is there a way to avoid this. Thanks, DL -- View this message in context: http://r.789695.n4.nabble.com/Have-write-csv-write-csv-without-observation-number-tp4030672p4030672.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] Generating the Ctrl-M character
The \r\n worked as the the diff can see no difference between my program's file and the file generated from the windows computer. For David, I use Emacs whitespace mode to see the spaces and the \n's. The \r does not show up in the whitespace mode. Maybe there is some way of turning it on. but THANK you guys. Ashim : ) On Fri, Nov 11, 2011 at 12:11 PM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote: No need to get defensive... I couldn't tell whether you remembered the details. This is starting to go OT, though. As to handing \n as expected, in the context of generating text files for their native platforms they do act as expected. However, in the context of people writing binary characters in a file in an attempt to make it suitable for some other platform, munging can happen, and people can be surprised by the results if they don't get the whole picture. Specifically, AFAIK \n is always stored as a single character (LF) in memory on UTF8 or ASCII systems. On Windows, when that character is written to a text file a \r is written just before it, creating the expected end-of-line marker for a text file on that OS. On the old Macintosh operating system a lone \r was the EOL marker, so \n was automatically translated to \r (but not on OSX, which uses UNIX conventions). From the perspective of the person working with native text files, this makes things just work. From the perspective of someone trying to obtain a particular but pattern in the output file this amounts to munging. The C programming standard library includes the concept of distinguishing between text files and binary files just so this native adaptation can occur or not, depending on the desired behavior. --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. 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. David Winsemius dwinsem...@comcast.net wrote: On Nov 10, 2011, at 11:55 PM, Jeff Newmiller wrote: Wow, Deadpan David. How about using the escape sequence \r? So shoot me. It wasn't documented in a manner that I recognized in the places I looked. ?character ?Syntax ?Constants And I did look at ?Quotes where \r is listed but did not know that it was == cntrl-M (if in fact it is.) I assumed (probably incorrectly) that \r was cntrl-R. There was proabably a time in the past when I could have told you a lot of the decimal and maybe even hexadecimal equivalents for cr, lf, beep, but those days are behind me. Keep in mind that Ctrl-M is used as the end-of-line character on some operating systems, so accomplishing this may not be portable, and you didn't specify your operating system. On the three main platforms (*nix, Windows/DOS, and Mac), \r should work, but \n may get munged. What? I very much doubt that any of those systems will not handle \n as expected. \r on the other hand I'm not so sure of. -- David. --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. 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. David Winsemius dwinsem...@comcast.net wrote: On Nov 10, 2011, at 9:35 PM, Ashim Kapoor wrote: Dear R-helpers, I want to append a Ctrl-M character to a string and then save it to a text file. mystring-This is a test. # How do I add a Ctrl-M to it in the end ?? cat(mystring,file=testfile) cntrl_m - intToUtf8(13) cat(cntrl_m,file=testfile) The resulting file seems to have a blank line in my editor. -- 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. David Winsemius, MD West Hartford, CT [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting
Re: [R] Have write.csv write csv without observation number
If you find yourself asking is there a way..?, put that question mark in front of the function you are wondering about and push return ?write.csv In this case, you will see that the write function has a parameter to say whether to write out row names ornot, which defaults to TRUE... row.names = TRUE On Nov 11, 2011, at 7:46 AM, dlarmour wrote: Hello all, Have a little annoying issue with wreite.csv. It always writes out observation number. Is there a way to avoid this. Thanks, DL -- View this message in context: http://r.789695.n4.nabble.com/Have-write-csv-write-csv-without-observation-number-tp4030672p4030672.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] Error with Source()
On 11-11-10 9:25 AM, ftonini wrote: Hi everybody, I started to receive a weird message in R that I have never seen before...also I haven't found anything on google or on this forum about it. Whenever I use the command source(...) to point to one of my scripts, I get the following message: Error in source(myfunctions.R) : myfunctions.R:884:9: unexpected symbol 883: 884: cond ^ I am using the same commands as I did in the past and it was working...I started to receive this error (not sure if it has to do with it or not) after trying to create a batch file to run one of my .R scripts with double-click. That batch file worked...but as soon as I use the source() command it does not work any more. Any help is appreciated! The message indicates that on line 884 of that script, the parser is finding a syntax error. A common way to generate that error is to fail to close parentheses, e.g. ( a cond gives a similar one. Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Gibbs sampler
On 11-11-10 1:28 PM, Gyanendra Pokharel wrote: I have the following code, gibbs-function(m,theta = 0.25, lambda =0.55, n =1){ alpha- 1.5 beta- 1.5 gamma- 1.5 x- array(0,c(m+1, 3)) x[1,1]- theta x[1,2]- lambda x[1,3]- n for(t in 2:(m+1)){ x[t,1]- rbinom(1, x[t-1,3], x[t-1,1]) x[t,2]-rbeta(1, x[t-1,1] + alpha, x[t-1,3] - x[t-1,1] + beta) Are you certain that x[t-1,3] - x[t-1,1] + beta will always be positive? If it is negative, rbeta will return NaN. Duncan Murdoch x[t,3]- rpois(1,(1 - x[t-1,1])*gamma) } x } gibbs(100) it returns only 1 or 2 values of theta, some times NaN, this may be if any theta is greater than 1, which is used as the probability for the next rbinom(), so returns NaN. Can some one suggest to solve this problem? [[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] beanplot without log scale?
Hi: Try the log = argument to beanplot(). Here's an example: library('beanplot') v - exp(runif(1000, 0, 10)) beanplot(v) beanplot(v, log = ') HTH, Dennis On Thu, Nov 10, 2011 at 1:27 PM, Twila Moon twi...@uw.edu wrote: Is it possible to force beanplot not to use a log scale? I want to be able to create multiple beanplots of different data on the same specified y-axis. When I try to force a ylim, I get the following... beanplot(reg5vel,ylim=(c(0,12000))) log=y selected Warning message: In plot.window(xlim = c(0.5, 7.5), ylim = c(0, 12000), log = y) : nonfinite axis limits [GScale(-inf,4.07918,2, .); log=1] Thanks, Twila [[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] performance of adaptIntegrate vs. integrate
baptiste auguie baptiste.auguie at googlemail.com writes: Dear list, [cross-posting from Stack Overflow where this question has remained unanswered for two weeks] I'd like to perform a numerical integration in one dimension, I = int_a^b f(x) dx where the integrand f: x in IR - f(x) in IR^p is vector-valued. integrate() only allows scalar integrands, thus I would need to call it many (p=200 typically) times, which sounds suboptimal. The cubature package seems well suited, as illustrated below, library(cubature) Nmax - 1e3 tolerance - 1e-4 integrand - function(x, a=0.01) c(exp(-x^2/a^2), cos(x)) adaptIntegrate(integrand, -1, 1, tolerance, 2, max=Nmax)$integral [1] 0.01772454 1.68294197 However, adaptIntegrate appears to perform quite poorly when compared to integrate. Consider the following example (one-dimensional integrand), library(cubature) integrand - function(x, a=0.01) exp(-x^2/a^2)*cos(x) Nmax - 1e3 tolerance - 1e-4 # using cubature's adaptIntegrate time1 - system.time(replicate(1e3, { a - adaptIntegrate(integrand, -1, 1, tolerance, 1, max=Nmax) }) ) # using integrate time2 - system.time(replicate(1e3, { b - integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax) }) ) time1 user system elapsed 2.398 0.004 2.403 time2 user system elapsed 0.204 0.004 0.208 a$integral [1] 0.0177241 b$value [1] 0.0177241 a$functionEvaluations [1] 345 b$subdivisions [1] 10 Somehow, adaptIntegrate was using many more function evaluations for a similar precision. Both methods apparently use Gauss-Kronrod quadrature, though ?integrate adds a Wynn's Epsilon algorithm. Could that explain the large timing difference? Cubature is astonishingly slow here though it was robust and accurate in most cases I used it. You may write to the maintainer for more information. Even a pure R implementation of adaptive Gauss-Kronrod as in pracma::quadgk is faster than cubature: library(pracma) time3 - system.time(replicate(1e3, { d - quadgk(integrand, -1, 1) }# 0.0177240954011303 ) ) # user system elapsed # 0.899 0.0020.897 If your functions are somehow similar and you are willing to dispense with the adaptive part, then compute Gaussian quadrature nodes xi and weights wi for the fixed interval and compute sum(wi * fj(xi)) for each function fj. This avoids recomputing nodes and weights for every call or function. Package 'gaussquad' provides such nodes and weights. Or copy the relevant calculations from quadgk (the usual (7,15)-rule). Regards, Hans Werner I'm open to suggestions of alternative ways of dealing with vector-valued integrands. Thanks. baptiste __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Building a statically-linked extension?
On Thu, 10 Nov 2011, Tyler Pirtle wrote: Hi there, I'm writing an R extension that has a C component that relies on two third party libraries that I'm bundling with the extension. I'd like to statically link my resulting extension so that I can rely on the bundled versions of the libraries I'm distributing, so there's two questions - 1, does R allow statically linked C extensions to be used at runtime? Yes 2, are there any standard ways of having R build my extension statically? Yes. If you want to know more, follow the posting guide: this is not an R-help topic, the missing 'at a minimum' information is vital Thanks, Tyler [[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. -- 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.
Re: [R] match first consecutive list of capitalized words in string
Thank you very much for your help. I'm using the second suggestion in my program and it works very well. Jonas -Original Message- From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com] Sent: Thu 11/10/2011 5:58 AM To: Richter-Dumke, Jonas Cc: r-help@r-project.org Subject: Re: [R] match first consecutive list of capitalized words in string On Tue, Nov 8, 2011 at 7:48 AM, Richter-Dumke, Jonas rich...@demogr.mpg.de wrote: Dear R-Helpers, this is my first post ever to a mailing list, so please feel free to point out any missunderstandings on my side regarding the conventions of this mailing list. My problem: Assuming the following character vector is given: names - c(filia Maria, vidua Joh Dirck Kleve (oo 02.02.1732), Bernardus Engelb Franciscus Linde j.u.Doktor referendarius sereniss Judex et gograven Rheinensis) Is there a regular expression matching the first consecutive list of capitalized words in a single characterstring (Maria, Joh Dirck Kleve, Bernardus Engelb Franciscus Linde)? This expression would very reliably seperate the person names from the additional information in my historic church register transcription. Try this. It matches a word boundary followed by zero or more of the parenthesized expression. That expression is an upper case letter followed by zero or more lower case letters followed by one or more spaces. Finally we match the last word which consists of an upper case letter followed by zero or more lower case letters and a word boundary. Note that it assumes R 2.14.0 or later: re - \\b([[:upper:]][[:lower:]]* +)*[[:upper:]][[:lower:]]*\\b regmatches(names, regexpr(re, names)) [1] Maria Joh Dirck Kleve [3] Bernardus Engelb Franciscus Linde -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com -- This mail has been sent through the MPI for Demographic ...{{dropped:10}} __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Win upgrade pb (virus)
I just upgraded my Win7 32bits installation to 2.14.0 after deinstalling 2.12.x First thing I moved the win-library from 2.12 to 2.14 and executed a update.packages(ask='graphics',checkBuilt=TRUE) (Swiss mirror). This aborts with the console message: Error in if (any(diff)) { : missing value where TRUE/FALSE needed And the antivirus (AVG) pops up complaining that colorspace.dll contains the virus Win32/Heur I cannot ignore the thread because update.packages has already crashed when I can reach the antivirus dialogbox to push the ignore button. To continue I had to remove the colorspace and ggplot2 packages. Is it a known problem? What else can I do (except removing the above packages or changing antivirus)? Thanks! mario -- Ing. Mario Valle Data Analysis and Visualization Group| http://www.cscs.ch/~mvalle Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60 v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Win upgrade pb (virus)
On Fri, 11 Nov 2011, Mario Valle wrote: I just upgraded my Win7 32bits installation to 2.14.0 after deinstalling 2.12.x First thing I moved the win-library from 2.12 to 2.14 and executed a update.packages(ask='graphics',checkBuilt=TRUE) (Swiss mirror). This aborts with the console message: Error in if (any(diff)) { : missing value where TRUE/FALSE needed And the antivirus (AVG) pops up complaining that colorspace.dll contains the virus Win32/Heur I cannot ignore the thread because update.packages has already crashed when I can reach the antivirus dialogbox to push the ignore button. To continue I had to remove the colorspace and ggplot2 packages. Is it a known problem? What else can I do (except removing the above packages or changing antivirus)? Yes (a known problem in AVG), and the standard advice is to change antivirus (or set exceptions, which I gather AVG does not allow). And please report the false positive to AVG: the more they hear about the more attention they will give it. The binary distribution has been checked by e.g. Sophos (which both winbuilder and ox.ac.uk use). Thanks! mario -- Ing. Mario Valle Data Analysis and Visualization Group| http://www.cscs.ch/~mvalle Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60 v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82 -- 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.
Re: [R] optim seems to be finding a local minimum
Ben Bolker bbolker at gmail.com writes: Simulated annealing and other stochastic global optimization methods are also possible solutions, although they may or may not work better than the many-starting-points solution -- it depends on the problem, and pretty much everything has to be tuned. Tabu search http://en.wikipedia.org/wiki/Tabu_search is another possibility, although I don't know much about it ... It is known that the Excel Solver has much improved during recent years. Still there are slightly better points such as myfunc(c(0.889764228112319, 94701144.5712312)) # 334.18844 restricting the domain to [0, 1] x [0, 10^9] for an evolutionary approach, for instance DEoptim::DEoptim(). Finding a global optimum in 2 dimensions is not so difficult. Here the scale of the second variable could pose a problem as small local minima might be overlooked easily. Hans Werner __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 package for segmentation with both continuous and categorical input variables XXXX
Th CRAN taskview on clustering has many clustering and mixture packages which may be useful for your precise situation: http://cran.r-project.org/web/views/Cluster.html hth, Ingmar On Thu, Nov 10, 2011 at 7:41 PM, Dan Abner dan.abne...@gmail.com wrote: Hello everyone, Can anyone suggest a decently documented (with good examples in the documentation) R package/function that performs segmentation (cluster, mixture modeling) of a population using both continuous and categorical input variables? Thank you, Dan [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Odp: Error in matrix, not ordered vectors or numerical value, and SIAR.
Re: [R] Odp: Error in matrix, not ordered vectors or numerical value, and SIAR. Thank-you Petr. I have consulted as many manuals and help pages, and search engines, and trying various things in SIAR, but continuous errors prevail.. I have tried changing all matrices to numeric, but the first column , of the two latter matrices, returns as NA. Any other suggestions? I had put in one of these messages more detail about data, see below. pond - c(A3, A3, A3, B3, B3, B3,C2, A7, A7) d13C - c(-22.07, -20.39, -21.11...) d15N - c(16.06, 17.13, 16.7, ...) FConsumerData - data.frame(pond, d13C, d15N) FTEFData - data.frame(tefsources, tefMeand13C, tefSDd13C, tefMeand15N, tefSDd15N) sources - c(Copepoda, Algal Paste, Amphipoda, Meand13C - c(-26.57, -20.41, -21.65,... SDd13C - c(1.78, 4.01, 1.46,... Meand15N - c(15.55, 9.58, 14.07... SDd15N - c(1.77, 1.34, 1.66,... FSourceData - data.frame(sources, Meand13C, SDd13C, Meand15N, SDd15N) tefsources - c(Copepoda, Algal Paste, Amphipoda,... tefMeand13C - c(0.4, 0.4, 0.4,... tefSDd13C - c(1.3, 1.3, 1.3, tefMeand15N - c(3.4, 3.4, 3.4, tefSDd15N - c(1.0, 1.0, 1.0, Hope that helps explain? Not much. FConsumerData is data frame with nonnumeric first column, the same is with FTEFData and FSourceData. If you change it to matrix you will get character matrix. Trying to change character values to numeric obviously leads to NA. How the poor R shall know what number it shall make from A3, B3, C2 or A7? I do not use siar package but in the help page it is clearly stated that input shall be ***numeric*** matrices. Not data frames, not character values not factors. If you do not know distinction among those objects you shall look into R intro document (at least first 6 chapters) which should came along with your installation. BTW, in help page there is out - siarmcmcdirichletv4(geese1demo,sourcesdemo,correctionsdemo,concdepdemo) It does not produce any error does it? If it does not produce any error it is sign that something is wrong with your data. Look at structure of your objects you are trying to evaluate and compare it with those demo objects. str(object) is the first option when something does not seem to be OK. Regards Petr ?? -- View this message in context: http://r.789695.n4.nabble.com/Error-in- matrix-not-ordered-vectors-or-numerical-value-and-SIAR-tp4024578p4030650.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] Building package problem
Hello Many thanks for the replies. I am note sure whether you've got what you meant (Prof. Ripley) but here is the output of Sys.getlocale() Sys.getlocale()[1] LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 Is that what you meant? Many thanks Ed On Tue, Nov 8, 2011 at 1:03 PM, Prof Brian Ripley rip...@stats.ox.ac.ukwrote: R CMD check is *not* 'building a package'. Nor is making a Windows binary package. 'Building a package' is creating a source tarball from a source directory. On Tue, 8 Nov 2011, Joshua Wiley wrote: Hi Ed, If the only error is in examples then this should work: R CMD check --no-examples foopkg should not have anything to do with vignettes (although those may also not run, who knows). As far as building a binary, look at: R CMD INSTALL --help which leads you to R CMD INSTALL --build foopkg And as for the hydroGOF issue, my guess is that the problem is your locale or timezone. But despite the posting guide, you failed to tell us. AFAIK CRAN only checks packages in English locales. (One thing we know is that in Columbia there was no midnight on one of the dates in that file. So hydroGOF really ought to be specifying a timezone when reading character datetimes.) HTH, Josh On Tue, Nov 8, 2011 at 4:35 AM, Eduardo M. A. M. Mendes emammen...@gmail.com wrote: Dear R-users I am trying to recompile a CRAN package on Windows 32. Rtools for 2.14 (that is the version I am running) and miktex were sucessfully installed on my machine. Problems: a) hydroGOF is a CRAN package, but R CMD check does not work on it. C:\Users\eduardo\Documents\R_**tests2R CMD check hydroGOF * using log directory 'C:/Users/eduardo/Documents/R_** tests2/hydroGOF.Rcheck' * using R version 2.14.0 (2011-10-31) * using platform: i386-pc-mingw32 (32-bit) * using session charset: ISO8859-1 * checking for file 'hydroGOF/DESCRIPTION' ... OK * checking extension type ... Package * this is package 'hydroGOF' version '0.3-2' * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for executable files ... OK * checking whether package 'hydroGOF' can be installed ... OK * checking installed package size ... OK * checking package directory ... OK * checking for portable file names ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... OK * checking whether the package can be loaded with stated dependencies ... OK * checking whether the package can be unloaded cleanly ... OK * checking whether the namespace can be loaded with stated dependencies ... OK * checking whether the namespace can be unloaded cleanly ... OK * checking for unstated dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of 'data' directory ... OK * checking data for non-ASCII characters ... OK * checking data for ASCII and uncompressed saves ... OK * checking examples ... ERROR Running examples in 'hydroGOF-Ex.R' failed The error most likely occurred in: ### Name: plot2 ### Title: Plotting 2 Time Series ### Aliases: plot2 ### Keywords: dplot ### ** Examples sim - 2:11 obs - 1:10 ## Not run: ##D plot2(sim, obs) ## End(Not run) ## # Loading daily streamflows of the Ega River (Spain), from 1961 to 1970 require(zoo) Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from 'package:base': as.Date, as.Date.numeric data(EgaEnEstellaQts) obs - EgaEnEstellaQts # Generating a simulated daily time series, initially equal to the observed se ries sim - obs # Randomly changing the first 2000 elements of 'sim', by using a normal distri bution # with mean 10 and standard deviation equal to 1 (default of 'rnorm'). sim[1:2000] - obs[1:2000] + rnorm(2000, mean=10) # Plotting 'sim' and 'obs' in 2 separate panels plot2(x=obs, y=sim) # Plotting 'sim' and 'obs' in the same window plot2(x=obs, y=sim, plot.type=single) Error in as.POSIXlt.character(x, tz, ...) : character string is not in a
Re: [R] A question on Programming
Christofer Bogaso bogaso.christofer at gmail.com writes: [...] Here my question is, is there any speed reduction if I put them within a function (I think there may be some speed reduction at least within for-loop, because that loop needs to call that function many times), relative to if I used that group of codes as-it-is in many places? [...] You did not asked for it, you may know it, or may not know it: if you use apply functions and other vector oriented functions, this can bring you a huge speedup, compared to a for-loop. Ciao, Oliver __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Building package problem
The posting guide asks for the output of sessionInfo() And what does Sys.timezone() say (it isn't always helpful). On Fri, 11 Nov 2011, Eduardo Mendes wrote: Hello Many thanks for the replies. I am note sure whether you've got what you meant (Prof. Ripley) but here is the output of Sys.getlocale() Sys.getlocale() [1] LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC _MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1 252 Is that what you meant? Many thanks Ed On Tue, Nov 8, 2011 at 1:03 PM, Prof Brian Ripley rip...@stats.ox.ac.uk wrote: R CMD check is *not* 'building a package'. Nor is making a Windows binary package. 'Building a package' is creating a source tarball from a source directory. On Tue, 8 Nov 2011, Joshua Wiley wrote: Hi Ed, If the only error is in examples then this should work: R CMD check --no-examples foopkg should not have anything to do with vignettes (although those may also not run, who knows). As far as building a binary, look at: R CMD INSTALL --help which leads you to R CMD INSTALL --build foopkg And as for the hydroGOF issue, my guess is that the problem is your locale or timezone. But despite the posting guide, you failed to tell us. AFAIK CRAN only checks packages in English locales. (One thing we know is that in Columbia there was no midnight on one of the dates in that file. So hydroGOF really ought to be specifying a timezone when reading character datetimes.) HTH, Josh On Tue, Nov 8, 2011 at 4:35 AM, Eduardo M. A. M. Mendes emammen...@gmail.com wrote: Dear R-users I am trying to recompile a CRAN package on Windows 32. Rtools for 2.14 (that is the version I am running) and miktex were sucessfully installed on my machine. Problems: a) hydroGOF is a CRAN package, but R CMD check does not work on it. C:\Users\eduardo\Documents\R_tests2R CMD check hydroGOF * using log directory 'C:/Users/eduardo/Documents/R_tests2/hydroGOF.Rcheck' * using R version 2.14.0 (2011-10-31) * using platform: i386-pc-mingw32 (32-bit) * using session charset: ISO8859-1 * checking for file 'hydroGOF/DESCRIPTION' ... OK * checking extension type ... Package * this is package 'hydroGOF' version '0.3-2' * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for executable files ... OK * checking whether package 'hydroGOF' can be installed ... OK * checking installed package size ... OK * checking package directory ... OK * checking for portable file names ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... OK * checking whether the package can be loaded with stated dependencies ... OK * checking whether the package can be unloaded cleanly ... OK * checking whether the namespace can be loaded with stated dependencies ... OK * checking whether the namespace can be unloaded cleanly ... OK * checking for unstated dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of 'data' directory ... OK * checking data for non-ASCII characters ... OK * checking data for ASCII and uncompressed saves ... OK * checking examples ... ERROR Running examples in 'hydroGOF-Ex.R' failed The error most likely occurred in:
Re: [R] optim seems to be finding a local minimum
Hans W Borchers hwborchers at googlemail.com writes: Ben Bolker bbolker at gmail.com writes: Simulated annealing and other stochastic global optimization methods are also possible solutions, although they may or may not work better than the many-starting-points solution -- it depends on the problem, and pretty much everything has to be tuned. Tabu search http://en.wikipedia.org/wiki/Tabu_search is another possibility, although I don't know much about it ... It is known that the Excel Solver has much improved during recent years. Still there are slightly better points such as myfunc(c(0.889764228112319, 94701144.5712312)) # 334.18844 restricting the domain to [0, 1] x [0, 10^9] for an evolutionary approach, for instance DEoptim::DEoptim(). Finding a global optimum in 2 dimensions is not so difficult. Here the scale of the second variable could pose a problem as small local minima might be overlooked easily. Have taken a (quick) second look at the problem, I agree that scaling and centering are more likely to be useful solutions than stochastic global optimization stuff. Even using 'parscale' (see the optim help page) may clear up the problem. 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] Title for a group of plots?
Or even simpler (unbelievable!), see ?title: par(mfcol=c(2,2), oma=c(0,0,1,0)) replicate(4, plot(1)) title(Hello World!, outer=TRUE) Best, Uwe Ligges On 10.11.2011 17:13, Sarah Goslee wrote: Try ?mtext for information on how to plot into the outer margin of a group of plots. Sarah On Thu, Nov 10, 2011 at 11:07 AM, Kevin Burtonrkevinbur...@charter.net wrote: I can get multiple plots on a page like: op- par(mfcol = c(3, 1)) What I was wondering is if there is a way to have a title for the whole page? I can specify the title for each individual plot like: plot(xxx, main=.) But I would like a 'title' for the group of plots. Is this possible? Thank you. Kevin __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] optim seems to be finding a local minimum
I won't requote all the other msgs, but the latest (and possibly a bit glitchy) version of optimx on R-forge 1) finds that some methods wander into domains where the user function fails try() (new optimx runs try() around all function calls). This includes L-BFGS-B 2) reports that the scaling is such that you really might not expect to get a good solution then 3) Actually gets a better result than the xlf-myfunc(c(0.888452533990788,94812732.0897449)) xlf [1] 334.607 with Kelley's variant of Nelder Mead (from dfoptim package), with myoptx methodpar fvalues fns grs itns conv KKT1 4 LBFGSB NA, NA 8.988466e+307 NA NULL NULL NA 2 Rvmmin 0.1, 200186870.6 25593.83 201 NULL0 FALSE 3 bobyqa 6.987875e-01, 2.001869e+08 1933.229 44 NA NULL0 FALSE 1 nmkb 8.897590e-01, 9.470163e+07 334.1901 204 NA NULL0 FALSE KKT2 xtimes meths 4NA 0.01 LBFGSB 2 FALSE 0.11 Rvmmin 3 FALSE 0.24 bobyqa 1 FALSE 1.08 nmkb But do note the terrible scaling. Hardly surprising that this function does not work. I'll have to delve deeper to see what the scaling setup should be because of the nature of the function setup involving some of the data. (optimx includes parscale on all methods). However, original poster DID include code, so it was easy to do a quick check. Good for him. JN ## Comparing this solution to Excel Solver solution: myfunc(c(0.888452533990788,94812732.0897449)) -- Dimitri Liakhovitski marketfusionanalytics.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] Sum of the deviance explained by each term in a gam model does not equal to the deviance explained by the full model.
It's the same problem as with any regression model where the predictors are not strictly orthogonal. You can make the explained deviances per term sum to the whole model explained deviance by dropping terms sequentially, but then the explained deviance per term depends on the order in which you drop terms. The alternative in the code you pasted in below is at least independent of the ordering. I suppose you could make things add up in the way you want by symmetrising the computation like this... ## reduce model one way dev.1- (deviance(b1)-deviance(b))/deviance(b0) ## prop expl by s(x2) dev.2- (deviance(b0)-deviance(b1))/deviance(b0) ## prop expl by s(x1) ## reduce it the other way... dev.2[2]- (deviance(b2)-deviance(b))/deviance(b0) ## prop expl by s(x2) dev.1[2]- (deviance(b0)-deviance(b2))/deviance(b0) # prop expl by s(x1) ## average the alternatives dev.1 - mean(dev.1) dev.2 - mean(dev.2) ## so now explained deviance adds up... dev.1+dev.2 summary(b)$dev.expl ... but I'm not convinced that this really adds anything useful, and it gets very tedious as the number of model terms increases... -- Simon Wood, Mathematical Science, University of Bath BA2 7AY UK +44 (0)1225 386603 http://people.bath.ac.uk/sw283 On 10/11/11 14:02, Huidong TIAN wrote: Dear R users, I read your methods of extracting the variance explained by each predictor in different places. My question is: using the method you suggested, the sum of the deviance explained by all terms is not equal to the deviance explained by the full model. Could you tell me what caused such problem? set.seed(0) n-400 x1- runif(n, 0, 1) ## to see problem with not fixing smoothing parameters ## remove the `##' from the next line, and the `sp' ## arguments from the `gam' calls generating b1 and b2. x2- runif(n, 0, 1) ## *.1 + x1 f1- function(x) exp(2 * x) f2- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10 f- f1(x1) + f2(x2) e- rnorm(n, 0, 2) y- f + e ## fit full and reduced models... b- gam(y~s(x1)+s(x2)) b1- gam(y~s(x1),sp=b$sp[1]) b2- gam(y~s(x2),sp=b$sp[2]) b0- gam(y~1) ## calculate proportions deviance explained... dev.1- (deviance(b1)-deviance(b))/deviance(b0) ## prop explained by s(x2) dev.2- (deviance(b2)-deviance(b))/deviance(b0) ## prop explained by s(x1) dev.1 + dev.2 [1] 0.6974949 summary(b)$dev.expl [1] 0.7298136 I checked the two models (b1 b2), found the model coefficients are different with model b, so I feel it could be the problem. wish to hear your comments. Huidong Tian [[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] proportional area venn diagrams?
In http://finzi.psych.upenn.edu/R/Rhelp02a/archive/14637.html some rudimentary R functions were given for drawing proportional area venn diagrams with area of each intersection ~ the count in a 2 x 2 x 2 table. I'm interested in this, for another application: showing the correlations among Y, X1, X2 using area ~ r^2 of each pair (sometimes called a Ballantine diagram). Before I attempt to hack the code given in that post, are there any packages/functions that do venn diagrams with proportional areas? Rseek turned up quite a view venn-like functions, but I couldn't find any that drew them with proportional areas. -Michael -- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 4700 Keele StreetWeb: http://www.datavis.ca Toronto, ONT M3J 1P3 CANADA __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] A question on Programming
Hi, On Fri, Nov 11, 2011 at 7:26 AM, Oliver oli...@first.in-berlin.de wrote: Christofer Bogaso bogaso.christofer at gmail.com writes: [...] Here my question is, is there any speed reduction if I put them within a function (I think there may be some speed reduction at least within for-loop, because that loop needs to call that function many times), relative to if I used that group of codes as-it-is in many places? [...] You did not asked for it, you may know it, or may not know it: if you use apply functions and other vector oriented functions, this can bring you a huge speedup, compared to a for-loop. I don't think this is exactly true. I believe the majority *apply functions are really more-or-less just sugar that boil down to having a for loop buried within them. This is just to say that I'm not sure how much of a speedup anybody should expect switching from a `for` loop to `*apply` function (assuming you aren't growing an array w/in each iteration of a for loop, for example), but if you can replace the for/*apply loop with a few lines of vectorized functions, then yeah -- you can expect big gains. -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] optim seems to be finding a local minimum
Hi Dimitri, Your problem has little to do with local versus global optimum. You can convince yourself that the solution you got is not even a local optimum by checking the gradient at the solution. The main issue is that your objective function is not differentiable everywhere. So, you have 2 options: either you use a smooth objective function (e.g. squared residuals) or you use an optimization algorithm than can handle non-smooth objective function. Here I show that your problem is well solved by the `nmkb' function (a bound-constraints version of Nelder-Mead simplex method) from the dfoptim package. library(dfoptim) myopt2 - nmkb(fn=myfunc, par=c(0.1,max(IV)), lower=0) myopt2 $par [1] 8.897590e-01 9.470163e+07 $value [1] 334.1901 $feval [1] 204 $restarts [1] 0 $convergence [1] 0 $message [1] Successful convergence Then, there is also the issue of properly scaling your function, because it is poorly scaled. Look how different the 2 parameters are - they are 7 orders of magnitude apart. You are really asking for trouble here. Hope this is helpful, Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edumailto:rvarad...@jhmi.edu [[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] Predicting x from y
Dear list members, I have just a quick question: I fitted a non-linear model y=a/x+b to describe my data (x=temperature and y=damage in %) and it works really nicely (see example below). I have 7 different species and 8 individuals per species. I measured damage for each individual per species at 4 different temperatures (e.g. -5, -10, -20, -40). Using the individuals per species, I fitted one model per species. Now I'd like to use the fitted model to go back and predict the temperature that causes 50% damage (and it's error). Basically, it pretty easy by just rearranging the formula to x=a/(y-b). But that way I don't get a measure of that temperature's error, do I? Can I use the residual standard error that R gave me for the non-linear model fit? Or do I have to fit 8 lines (each individual) per species, calculate x based on the 8 individuals and then take the mean? Unfortunately, dose.p from the MASS package doesn't work for non-linear models. When I take the log(abs(x)) the relationship becomes not satisfactory linear either. Any suggestions are highly appreciated! Thank you! Stefan EXAMPLE for species #1: y.damage-c(5.7388985,1.7813519,3.7321461,2.9671031, 0.3223196,0.3207941,-1.4197658,-5.3472160, 41.1826677,29.3115243,31.3208841,35.3934115, 58.5848778,31.1541049,42.2983479,27.0615648, 64.1037728,54.7003353,66.7317044,65.4725881, 72.5755056,67.2683495,64.8717942,65.9603073, 75.0762273,56.7041960,60.0049429,70.0286506, 73.2801947,72.7015642,75.0944694,81.0361280) x.temp-c(-5,-5,-5,-5,-5,-5,-5,-5,-10,-10,-10,-10,-10,-10,-10, -10,-20,-20,-20,-20,-20,-20,-20,-20,-40,-40,-40,-40,-40, -40,-40,-40) nls(y.damage~a/x.temp+b,start=list(a=400,b=80)) plot(y.damage~x.temp,xlab='Temperature',ylab='Damage [%]') curve(409.61/x+81.84,from=min(x.temp),to=max(x.temp),add=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.
Re: [R] barplot names.arg
That line in the documentation refers to groups (i.e. columns) on each row of data. C.1 and C.2 are your groups and if you use the beside=TRUE option C.1 and C.2 will appear side-by-side with a single labels on the axis (instead of stacked one atop the other. You could rearrange your data so that columns C.1 and C.2 in the duplicate weekdays became C3. And C.4 in the grouped bar graph and the total number of rows in the matrix would be 4. Using beside=TRUE would give you four bars for Mon/Tue and two bars for Thur/Fri. But you would not have C.1/C.2 and C.3/C.4 stacked. -- 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 Diviya Smith Sent: Thursday, November 10, 2011 11:10 PM To: R. Michael Weylandt Cc: r-help@r-project.org Subject: Re: [R] barplot names.arg The documentation for bar plots includes - names.arga vector of names to be plotted below each bar or *group of bars*. If this argument is omitted, then the names are taken from the names attribute of height if this is a vector, or the column names if it is a matrix Thanks for your suggestion. However, I am still not sure how to group all the entries for Mon and include one label. Your solution sort of does that but it actually does not print the individual values . I would like to show that there is a lot of variation in C.1 for Mon. On Thu, Nov 10, 2011 at 11:32 PM, R. Michael Weylandt michael.weyla...@gmail.com wrote: I'm not sure you can do that data aggregation in barplot directly (a quick skim doesn't reveal anything in the documentation that suggests it to me, thought I might have missed it) though I think this does what you are talking about: barplot(sapply(unique(rownames(mdat)), function(n) colSums(mdat[n,,drop=F])), col = rainbow(2)) Michael On Thu, Nov 10, 2011 at 9:21 PM, Diviya Smith diviya.sm...@gmail.com wrote: Hello there, I have a question regarding bar plots. I am trying to plot the data from the following matrix as a barplot - # input data mdat - matrix(c(0.1,0.9,0.9,0.1,0.5,0.5,0.45,1-0.45,0.6,0.4,0.8,0.2), nrow = 6, ncol=2, byrow=TRUE, +dimnames = list(c(Mon, Mon, Tues, Tues, Thurs, Friday), +c(C.1, C.2))) # plot barplot(t(as.matrix(mdat)), col=rainbow(2), ylab=Sales, names.arg = rownames(mdat), border=NA, cex.names=0.5) I am using names.arg to print the label for the bars. However, I would like to group all the entries for the Mon and print the label only once. Is there a way to do this? The documentation for barplots suggests that this can be done but I was not able to figure it out. Please help. Thanks in advance, Diviya [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Upgrade R?
On 10.11.2011 17:53, Kevin Burton wrote: The problem with this documentation is two-fold. One it seems to concentrate on building from source which I don't need. Two it doesn't address the upgade. I have a number of packages and so I need to do what has been suggested and install the latest version *first*. Then copy the libraries (packages). Then uninstall the previous version. It is on this last step that I am stuck on right now. The last link on uninstalling R manually was what I needed. Thank you. Just delete the directory. Uwe Ligges Kevin From: jose Bartolomei [mailto:surfpr...@hotmail.com] Sent: Thursday, November 10, 2011 10:19 AM To: rkevinbur...@charter.net; R Help Subject: RE: [R] Upgrade R? Hi, Don't know if this will help you but... In my short experience and following the guidelines you should first uninstall R. http://cran.r-project.org/doc/manuals/R-admin.html#Installing-R-under-Window s Unistall it from the Windows control panel. The old R version libraries file will be kept on machine. For example : C:\Program Files\R\R-2.13.0\library Then install the new version via: http://cran.r-project.org/doc/manuals/R-admin.html#Installing-R-under-Window s You can copy/paste libraries from the old version R library file to the new one. C:\Program Files\R\R-2.14.0\library There is too an function named: ?update.packagee If above was what you did, then there is a post on Uninstalling R manually: http://learnserver.csd.univie.ac.at/rcomwiki/doku.php?id=wiki:uninstalling_r _manually Regards, Jose From: rkevinbur...@charter.net To: r-help@r-project.org Date: Thu, 10 Nov 2011 09:07:20 -0600 Subject: Re: [R] Upgrade R? Since apparently there is no one familiar with this error message let me rephrase the question. Is there a 'manual' process to fully remove a version of 'R' from my machine? This is a Window PC running Windows 7. Thank you. Kevin From: Kevin Burton [mailto:rkevinbur...@charter.net] Sent: Monday, November 07, 2011 2:23 PM To: 'r-help@r-project.org' Subject: Upgrade R? I am trying to upgrade to R 2.14 from R 2.13.1 I have compied all the libraries from the 'library' directory in my existing installation (2.13.1) to the installed R 2.14. Now I want to uninstall the old installation (R 2.13.1) and I get the error: Internal Error: Cannot find utCompiledCode record for this version of the uninstaller. Any ideas? Kevin [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] optim seems to be finding a local minimum
Thank you very much to everyone who replied! As I mentioned - I am not a mathematician, so sorry for stupid comments/questions. I intuitively understand what you mean by scaling. While the solution space for the first parameter (.alpha) is relatively compact (probably between 0 and 2), the second one (.beta) is all over the place - because it is a function of IV (input vector). And that's, probably, my main challenge - that I am trying to write a routine for different possible IVs that I might be facing (they may be in hundreds, in thousands, in millions). Should I be rescaling the IV somehow (e.g., by dividing it by its max) - or should I do something with the parameter .beta inside my function? So far, I've written a loop over many different starting points for both parameters. Then, I take the betas around the best solution so far, split it into smaller steps for beta (as starting points) and optimize again for those starting points. What disappoints me is that even when I found a decent solution (the minimized value of 336) it was still worse than the Solver solution! And I am trying to prove to everyone here that we should do R, not Excel :-) Thanks again for your help, guys! Dimitri On Fri, Nov 11, 2011 at 9:10 AM, John C Nash nas...@uottawa.ca wrote: I won't requote all the other msgs, but the latest (and possibly a bit glitchy) version of optimx on R-forge 1) finds that some methods wander into domains where the user function fails try() (new optimx runs try() around all function calls). This includes L-BFGS-B 2) reports that the scaling is such that you really might not expect to get a good solution then 3) Actually gets a better result than the xlf-myfunc(c(0.888452533990788,94812732.0897449)) xlf [1] 334.607 with Kelley's variant of Nelder Mead (from dfoptim package), with myoptx method par fvalues fns grs itns conv KKT1 4 LBFGSB NA, NA 8.988466e+307 NA NULL NULL NA 2 Rvmmin 0.1, 200186870.6 25593.83 20 1 NULL 0 FALSE 3 bobyqa 6.987875e-01, 2.001869e+08 1933.229 44 NA NULL 0 FALSE 1 nmkb 8.897590e-01, 9.470163e+07 334.1901 204 NA NULL 0 FALSE KKT2 xtimes meths 4 NA 0.01 LBFGSB 2 FALSE 0.11 Rvmmin 3 FALSE 0.24 bobyqa 1 FALSE 1.08 nmkb But do note the terrible scaling. Hardly surprising that this function does not work. I'll have to delve deeper to see what the scaling setup should be because of the nature of the function setup involving some of the data. (optimx includes parscale on all methods). However, original poster DID include code, so it was easy to do a quick check. Good for him. JN ## Comparing this solution to Excel Solver solution: myfunc(c(0.888452533990788,94812732.0897449)) -- Dimitri Liakhovitski marketfusionanalytics.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. -- Dimitri Liakhovitski marketfusionanalytics.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] optim seems to be finding a local minimum
Some tips: 1) Excel did not, as far as I can determine, find a solution. No point seems to satisfy the KKT conditions (there is a function kktc in optfntools on R-forge project optimizer. It is called by optimx). 2) Scaling of the input vector is a good idea given the seeming wide range of values. That is, assuming this can be done. If the function depends on the relative values in the input vector rather than magnitude, this may explain the trouble with your function. That is, if the function depends on the relative change in the input vector and not its scale, then optimizers will have a lot of trouble if the scale factor for this vector is implicitly one of the optimization parameters. 3) If you can get the gradient function you will almost certainly be able to do better, especially in finding whether you have a minimum i.e., null gradient, positive definite Hessian. When you have gradient function, kktc uses Jacobian(gradient) to get the Hessian, avoiding one level of digit cancellation. JN On 11/11/2011 10:20 AM, Dimitri Liakhovitski wrote: Thank you very much to everyone who replied! As I mentioned - I am not a mathematician, so sorry for stupid comments/questions. I intuitively understand what you mean by scaling. While the solution space for the first parameter (.alpha) is relatively compact (probably between 0 and 2), the second one (.beta) is all over the place - because it is a function of IV (input vector). And that's, probably, my main challenge - that I am trying to write a routine for different possible IVs that I might be facing (they may be in hundreds, in thousands, in millions). Should I be rescaling the IV somehow (e.g., by dividing it by its max) - or should I do something with the parameter .beta inside my function? So far, I've written a loop over many different starting points for both parameters. Then, I take the betas around the best solution so far, split it into smaller steps for beta (as starting points) and optimize again for those starting points. What disappoints me is that even when I found a decent solution (the minimized value of 336) it was still worse than the Solver solution! And I am trying to prove to everyone here that we should do R, not Excel :-) Thanks again for your help, guys! Dimitri On Fri, Nov 11, 2011 at 9:10 AM, John C Nash nas...@uottawa.ca wrote: I won't requote all the other msgs, but the latest (and possibly a bit glitchy) version of optimx on R-forge 1) finds that some methods wander into domains where the user function fails try() (new optimx runs try() around all function calls). This includes L-BFGS-B 2) reports that the scaling is such that you really might not expect to get a good solution then 3) Actually gets a better result than the xlf-myfunc(c(0.888452533990788,94812732.0897449)) xlf [1] 334.607 with Kelley's variant of Nelder Mead (from dfoptim package), with myoptx methodpar fvalues fns grs itns conv KKT1 4 LBFGSB NA, NA 8.988466e+307 NA NULL NULL NA 2 Rvmmin 0.1, 200186870.6 25593.83 201 NULL0 FALSE 3 bobyqa 6.987875e-01, 2.001869e+08 1933.229 44 NA NULL0 FALSE 1 nmkb 8.897590e-01, 9.470163e+07 334.1901 204 NA NULL0 FALSE KKT2 xtimes meths 4NA 0.01 LBFGSB 2 FALSE 0.11 Rvmmin 3 FALSE 0.24 bobyqa 1 FALSE 1.08 nmkb But do note the terrible scaling. Hardly surprising that this function does not work. I'll have to delve deeper to see what the scaling setup should be because of the nature of the function setup involving some of the data. (optimx includes parscale on all methods). However, original poster DID include code, so it was easy to do a quick check. Good for him. JN ## Comparing this solution to Excel Solver solution: myfunc(c(0.888452533990788,94812732.0897449)) -- Dimitri Liakhovitski marketfusionanalytics.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] Fwd: Use of R for VECM
- Forwarded Message - From: vrama...@neo.tamu.edu To: bernhard pfaff bernhard.pf...@pfaffikus.de Sent: Friday, November 11, 2011 9:03:11 AM GMT -06:00 US/Canada Central Subject: Use of R for VECM Hello Fellow R'ers I am a new user of R and I am applying it for solving Bi-Variate (Consumption and Output) VECM with Co-Integration (I(1)) with three lags on Consumption and Output expressed aa lags of differences. R Code and one segment fo the output (other parts of the output are repeatitive) are as follows. us=read.table(usdata.1995-2009.txt,header=T) sjd-us[,c(Y,C)] sjd.vecm1 - ca.jo(sjd, ecdet='const', type=eigen, K=3, spec=longrun, + season=4) sjd.vecm2 - ca.jo(sjd, ecdet='const', type=eigen, K=3, spec=transitory, + season=4) sjd.vecm.ols1 - cajools(sjd.vecm1) sjd.vecm.ols2 - cajools(sjd.vecm2) summary(sjd.vecm.ols1) Response Y.d : Call: lm(formula = substitute(Y.d), data = data.mat) Residuals: Min 1Q Median 3QMax -0.0049787 -0.0012948 0.703 0.0009653 0.0063192 Coefficients: Estimate Std. Error t value Pr(|t|) sd1 -0.0002012 0.0007653 -0.263 0.793820 sd2 0.0013339 0.0007616 1.752 0.086378 . sd3 0.0007372 0.0007947 0.928 0.358348 Y.dl1-0.2215246 0.1600532 -1.384 0.172875 C.dl1 0.9220846 0.1646314 5.601 1.08e-06 *** Y.dl2-0.1600219 0.1273838 -1.256 0.215245 C.dl2 0.4712112 0.2124175 2.218 0.031401 * Y.l3 -0.3227708 0.0881079 -3.663 0.000631 *** C.l3 0.2376579 0.0688854 3.450 0.001194 ** constant 0.4707624 0.1182284 3.982 0.000236 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.002013 on 47 degrees of freedom Multiple R-squared: 0.7053, Adjusted R-squared: 0.6426 F-statistic: 11.25 on 10 and 47 DF, p-value: 1.707e-09 I am trying to decipher the output and have the following questions. 1. If you have published a book with explanations of the output, I would like to acquire one. Please advise me where I can get one. It might me and you lot of time. In the meantime 2. Are Y.dl1, C.dl1, Y.dl2, C.dl2, Y.l3 and C.l3 are the coefficients I am looking for? 3. My equition requires Y.l1 and C.l1. Is there a way I could change from Y.l3 and C.l3 4. I need to use these coefficients for forecasting. What command I use to extract the coefficients from this program? 5. What are sd1, sd2, sd3? I did not input any seasonal dummies and yet the output seems to have included them. Does it affect the value of the coefficients? Is there a better way to reframe the command to avoid them? I hope it is not too much to ask for your help. I thank you in advance. Ram Ramaiah __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] performance of adaptIntegrate vs. integrate
The integrand is highly peaked. It is approximately an impulse function where much of the mass is concentrated at a very small interval. Plot the function and see for yourself. This is the likely cause of the problem. Other types of integrands where you could experience problems are: integrands with singularity at either limit and slowly decaying oscillatory integrands. As to why integrate performs better than adaptIntegrate in this situation, I don't know. You have to study the two implementations. Wynn's epsilon algorithm is an extrapolation method for improving the convergence of a sequence. This could be an explanation for the better performance, but I cannot say for sure. Hope this is helpful, Ravi --- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edumailto:rvarad...@jhmi.edu [[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] drawing ellipses in R
Those formulas are the standard way to convert from polar coordinates to Euclidean coordinates. The polar coordinates are 'r' which is the radius or distance from the center point and 'theta' which is the angle (0 is pointing in the positive x direction). If r is constant and theta coveres a full cycle then you will get a circle. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of mms...@comcast.net Sent: Monday, October 31, 2011 10:50 PM To: r-h...@stat.math.ethz.ch Subject: [R] drawing ellipses in R Hello, I have been following the thread dated Monday, October 9, 2006 when Kamila Naxerova asked a question about plotting elliptical shapes. Can you explain the equations for X and Y. I believe they used the parametric form of x and y (x=r cos(theta), y=r sin(theta). I don't know what r is here ? Can you explain 1)the origin of these equations and 2) what is r? Sincerely, Mary A. Marion [[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] multivariate random variable
Your question is a bit too general to give a useful answer. One possible answer to your question is: mrv - matrix( runif(1000), ncol=10 ) Which generates multivariate random observations, but is unlikely to be what you are really trying to accomplish. There are many tools for generating multivariate random data including Metropolis-Hastigns, Gibbs sampling, rejection sampling, conditional generation, copulas, and many others, which one will be best (or which combination will be best) depends on what you are actually trying to accomplish. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Anera Salucci Sent: Tuesday, November 01, 2011 5:22 AM To: r-help@r-project.org Subject: [R] multivariate random variable Dear All, How can I generate multivariate random variable (not multivariate normal ) I am in urgent [[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] Building a statically-linked extension?
On Fri, 11 Nov 2011, Tyler Pirtle wrote: On Fri, Nov 11, 2011 at 2:09 AM, Prof Brian Ripley rip...@stats.ox.ac.uk wrote: On Thu, 10 Nov 2011, Tyler Pirtle wrote: Hi there, I'm writing an R extension that has a C component that relies on two third party libraries that I'm bundling with the extension. I'd like to statically link my resulting extension so that I can rely on the bundled versions of the libraries I'm distributing, so there's two questions - 1, does R allow statically linked C extensions to be used at runtime? Yes 2, are there any standard ways of having R build my extension statically? Yes. If you want to know more, follow the posting guide: this is not an R-help topic, the missing 'at a minimum' information is vital Apologies for spamming R-help, I suppose this should have gone to R-devel. I can't ever remember the difference between the two. So, Professor, if I re-post to R-devel would you be willing to reveal any more detail for me? ;) I also apologize for a lack of additional information, I'm happy to provide you with any, but I'm not exactly sure what else I can offer..I just want to build a statically linked R extension, and it doesn't seem to be covered in Writing R Extensions. You seem to know the answer, any pointers would be greatly appreciated. You will need to tell us your OS, for a start. Tyler Thanks, Tyler [[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. -- 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, UK Fax: +44 1865 272595 -- 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.
Re: [R] Building a statically-linked extension?
On Fri, Nov 11, 2011 at 2:09 AM, Prof Brian Ripley rip...@stats.ox.ac.ukwrote: On Thu, 10 Nov 2011, Tyler Pirtle wrote: Hi there, I'm writing an R extension that has a C component that relies on two third party libraries that I'm bundling with the extension. I'd like to statically link my resulting extension so that I can rely on the bundled versions of the libraries I'm distributing, so there's two questions - 1, does R allow statically linked C extensions to be used at runtime? Yes 2, are there any standard ways of having R build my extension statically? Yes. If you want to know more, follow the posting guide: this is not an R-help topic, the missing 'at a minimum' information is vital Apologies for spamming R-help, I suppose this should have gone to R-devel. I can't ever remember the difference between the two. So, Professor, if I re-post to R-devel would you be willing to reveal any more detail for me? ;) I also apologize for a lack of additional information, I'm happy to provide you with any, but I'm not exactly sure what else I can offer..I just want to build a statically linked R extension, and it doesn't seem to be covered in Writing R Extensions. You seem to know the answer, any pointers would be greatly appreciated. Tyler Thanks, Tyler [[alternative HTML version deleted]] __** R-help@r-project.org mailing list https://stat.ethz.ch/mailman/**listinfo/r-helphttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/** posting-guide.html http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~**ripley/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 [[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] Jordan Form of a matrix
Hello. Is it possible to find the Jordan Form of a matrix with R? Arnau. Arnau Mir Torres Edifici A. Turmeda Campus UIB Ctra. Valldemossa, km. 7,5 07122 Palma de Mca. tel: (+34) 971172987 fax: (+34) 971173003 email: arnau@uib.es URL: http://dmi.uib.es/~arnau __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Fwd: Use of R for VECM
I suspect that part of the hesitancy in replying to your query is that your academic e-mail and subsequent google-ability suggest you are a student in the economics department and that this might be homework. In the meanwhile, you might want to look into Pfaff's book on cointegration and time series analysis in R. Michael On Nov 11, 2011, at 10:06 AM, vrama...@neo.tamu.edu wrote: - Forwarded Message - From: vrama...@neo.tamu.edu To: bernhard pfaff bernhard.pf...@pfaffikus.de Sent: Friday, November 11, 2011 9:03:11 AM GMT -06:00 US/Canada Central Subject: Use of R for VECM Hello Fellow R'ers I am a new user of R and I am applying it for solving Bi-Variate (Consumption and Output) VECM with Co-Integration (I(1)) with three lags on Consumption and Output expressed aa lags of differences. R Code and one segment fo the output (other parts of the output are repeatitive) are as follows. us=read.table(usdata.1995-2009.txt,header=T) sjd-us[,c(Y,C)] sjd.vecm1 - ca.jo(sjd, ecdet='const', type=eigen, K=3, spec=longrun, + season=4) sjd.vecm2 - ca.jo(sjd, ecdet='const', type=eigen, K=3, spec=transitory, + season=4) sjd.vecm.ols1 - cajools(sjd.vecm1) sjd.vecm.ols2 - cajools(sjd.vecm2) summary(sjd.vecm.ols1) Response Y.d : Call: lm(formula = substitute(Y.d), data = data.mat) Residuals: Min 1Q Median 3QMax -0.0049787 -0.0012948 0.703 0.0009653 0.0063192 Coefficients: Estimate Std. Error t value Pr(|t|) sd1 -0.0002012 0.0007653 -0.263 0.793820 sd2 0.0013339 0.0007616 1.752 0.086378 . sd3 0.0007372 0.0007947 0.928 0.358348 Y.dl1-0.2215246 0.1600532 -1.384 0.172875 C.dl1 0.9220846 0.1646314 5.601 1.08e-06 *** Y.dl2-0.1600219 0.1273838 -1.256 0.215245 C.dl2 0.4712112 0.2124175 2.218 0.031401 * Y.l3 -0.3227708 0.0881079 -3.663 0.000631 *** C.l3 0.2376579 0.0688854 3.450 0.001194 ** constant 0.4707624 0.1182284 3.982 0.000236 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.002013 on 47 degrees of freedom Multiple R-squared: 0.7053,Adjusted R-squared: 0.6426 F-statistic: 11.25 on 10 and 47 DF, p-value: 1.707e-09 I am trying to decipher the output and have the following questions. 1. If you have published a book with explanations of the output, I would like to acquire one. Please advise me where I can get one. It might me and you lot of time. In the meantime 2. Are Y.dl1, C.dl1, Y.dl2, C.dl2, Y.l3 and C.l3 are the coefficients I am looking for? 3. My equition requires Y.l1 and C.l1. Is there a way I could change from Y.l3 and C.l3 4. I need to use these coefficients for forecasting. What command I use to extract the coefficients from this program? 5. What are sd1, sd2, sd3? I did not input any seasonal dummies and yet the output seems to have included them. Does it affect the value of the coefficients? Is there a better way to reframe the command to avoid them? I hope it is not too much to ask for your help. I thank you in advance. Ram Ramaiah __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] rbind.data.frame drops attributes for factor variables
As the doctor says, if it hurts don't do that. A factor is a sequence of integers with a corresponding list of character strings. Factors in two separate vectors can and usually do map the same integer to different strings, and R cannot tell how you want that resolved. Convert these columns to character before combining them, and only convert to factor when you have all of your possibilities present (or you specify them in the creation of the factor vector). --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. 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. Sammy Zee szee2...@gmail.com wrote: Hi all, When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. I see the following post with same problem. However i did not see any reply to the following posting offering a solution. Could someone please help. http://r.789695.n4.nabble.com/rbind-data-frame-drops-attributes-for-factor-variables-td919575.html Thanks, Sanjay [[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] mc.cores and computer settings on osx and linux
for the googleable r-help archives, I thought I would post what I wrote into my .Rprofile to automatically set some system information. the most relevant aspect is the determination of mc.cores. this is useful when users want to use the parallel package options(uname= system(uname, intern=TRUE)) options(os= if (getOption(uname)==Darwin) osx else linux) if ((getOption(os) != osx) (getOption(os) != linux)) stop(You need to set options yourself. I only grok linux and osx\n) options(mc.cores= as.numeric(if (getOption(os)==osx) system(sysctl -n hw.ncpu, intern=TRUE) else system(grep 'core id' /proc/cpuinfo | sort | uniq | wc -l, intern=TRUE))) options(hostname= system(hostname, intern=TRUE)) [below, I am also posting my current wrapper for the parallel library. I know it is amateurish, but it may be useful for novices exploring parallel calculations. it is a friendlier face.] [gripes: R is powerful, and the team that maintains it are saints. But R is not friendly. it lacks the ability to turn off recycling for enhanced error detection. it does not throw clear errors when one accesses a non-existing column in a data frame. it does not print out the user program line number where an error occurred. it lacks an end-user documentation system [like POD], though it does have good package documentation. it does have some unexpected behavior: mymatrix[1:2,] is a matrix, but mymatrix[1:1,] is a numeric. huh? data.table is necessary for reasonably fast data manipulation, but data.table giveth and taketh. it has some really strange unexpected behavior---mydatatable[,1] is not the second column, as one would expect it to be. yes, it is documented, but syntax should be as expected.] /iaw Ivo Welch (ivo.we...@gmail.com) ### ### these R functions are very type-limited wrappers for ### by()-like operations, using the multicore library. this ### means effort-less multi-CPU calculations. ### ### the user functions MUST return a numeric scalar, a vector, a matrix, or ### a data frame. to enhance speed, internally the user function is ### wrapped, too. ### ### the output is ONE matrix, whose row-names are the categories. check.output - function( mc.rv ) { ## check that we have a list of matrices, and that each matrix has the same number of columns numofcols - (-1) for (i in 1:length(mc.rv)) { if (is.null(mc.rv[[i]])) next; if (! (is.matrix(mc.rv[[i]])|is.data.frame(mc.rv[[i]])) ) abort(iaw-mc.R:check.output: Element, i, is not a matrix/dataframe, but a , whatis(mc.rv[[i]])) if (numofcols0) numofcols - ncol(mc.rv[[i]]) if (numofcols0) next if (ncol(mc.rv[[i]]) != numofcols) { print(head(mc.rv[[i]])) abort(iaw-mc.R:check.output: Element, i, should have, numofcols, columns, but has, ncol(mc.rv[[i]]), columns instead.) } } } add.by.names - function( mc.rv ) { for (i in 1:length(mc.rv)) if (!is.null(mc.rv[[i]])) row.names(mc.rv[[i]]) - rep( names(mc.rv)[i], nrow(mc.rv[[i]]) ) mc.rv } .mc.by - function(lcapplyversion, data, INDICES, FUN, ...) { si - split(1:nrow(data), INDICES) ## input = set of row indexes ; output = one row in a matrix or data frame, that can be stacked up FUN.ON.ROWS - function(.index, ...) { rv - FUN(data[.index,], ...); if (is.null(rv)) rv else if (is.vector(rv)) matrix(rv, nrow=1) else rv } soln - lcapplyversion( si, FUN.ON.ROWS, ... ) check.output(soln) rv - do.call(rbind, add.by.names(soln)) if (is.null(rv)) { print(head(soln)); abort(Sorry, but in .mc.by, the rv is null!\n) } if (ncol(rv)==1) { nm - rownames(rv) rv - as.vector(rv) names(rv) - nm } rv } mc.by - function(data, INDICES, FUN, ...) { .mc.by(mclapply, data, INDICES, FUN, ...) } oc.by - function(data, INDICES, FUN, ...) { .mc.by(lapply, data, INDICES, FUN, ...) } mc.byallrows - function(data, FUN, ...) { si - as.list(1:nrow(data)) ## a little faster than the split for large data sets FUN.ON.ROWS - function(.index, ...) { rv - FUN(data[.index,], ...); if (is.null(rv)) rv else if (is.vector(rv)) matrix(rv, nrow=1) else rv } soln - mclapply( si, FUN.ON.ROWS, ..., mc.cores= 4 ) check.output(soln) rv - do.call(rbind, soln) ## omits naming. if (ncol(rv)==1) rv - as.vector(rv) rv } if (0) { function.sample - function(d) cbind(d$x+d$y, d$x, d$y) function.sample.simpler - function(d) (d$x+d$y) d - data.frame( i=c( rep(1,2), rep(2,3), rep(3,4) ), x=rnorm(9), y=rnorm(9) ) report - function( text2print, f.output ) { cat(\n\n, text2print, :\n); print(f.output); cat(\n\n) } report( the original R by() function, by( d, d$i, function.sample )) report( wrappled multicore by mc.by with user function returning scalar, mc.by( d, d$i, function.sample.simpler )) report( wrappled multicore by mc.by with user function returning vector, mc.by( d, d$i, function.sample
[R] Random-walk Metropolis-Hasting
Following is my code, can some one help on the error at the bottom? mh-function(iterations,alpha,beta){ + data-read.table(epidemic.txt,header = TRUE) + attach(data, warn.conflicts = F) + k-97 + d - (sqrt((x-x[k])^2 + (y-y[k])^2)) + p - 1-exp(-alpha*d^(-beta)) + p.alpha-1 - exp(-3*d^(-beta)) + p.beta - 1 - exp(alpha*d^(-2)) + iterations-1000 + mu.lambda - c(0,0);s.lambda - c(100,100) + prop.s - c(0.1,0.1) + lambda - matrix(nrow=iterations, ncol=2) + acc.prob-0 + current.lambda - c(0,0) + for(t in 1:iterations){ + prop.lambda - rnorm(2,current.lambda,prop.s) + a - p.beta/p.alpha *(dnorm(prop.lambda,mu.lambda,s.lambda))*dnorm(current.lambda,mu.lambda,s.lambda) + accept - min(1,a) + u-runif(1) + if(u[t]=accept[t]){ + current.lambda - prop.lambda + acc.prob - acc.prob +1 + } + lambda[t,] - current.lambda + } + lambda + } mh(1000,0,0) Error in if (u[t] = accept[t]) { : missing value where TRUE/FALSE needed [[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 variable help
I have an R script with the following applicable lines: xshort - window(s, start=st, end=ed) . . . xshort - ts(xshort, frequency=1, start=1) . . . m1 - m2 - m3 - m4 - m5 - m6 - NULL m1 - tslm(xshort ~ trend) I get an error: Error in get(dataname) : object 'xshort' not found When I do traceback() I get: 3: get(dataname) 2: tslm(xshort ~ trend) at #19 1: model.cross.validation(l[[MEN]]$series) Which points to the call to tslm above. Since I am not supply 'data' to the tslm call (in the forecast package),, I am assuming that the code is dying here (in tslm): if (missing(data)) { dataname - as.character(formula)[2] x - get(dataname) data - data.frame(x) colnames(data) - dataname } Any ideas what is failing? Kevin [[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] Formula variable help
Seems like it could work---can you save your script as a .txt file and send it as an attachment, upload it online, or make a reproducible example? Cheers, Josh On Fri, Nov 11, 2011 at 9:52 AM, rkevinbur...@charter.net wrote: I have an R script with the following applicable lines: xshort - window(s, start=st, end=ed) . . . xshort - ts(xshort, frequency=1, start=1) . . . m1 - m2 - m3 - m4 - m5 - m6 - NULL m1 - tslm(xshort ~ trend) I get an error: Error in get(dataname) : object 'xshort' not found When I do traceback() I get: 3: get(dataname) 2: tslm(xshort ~ trend) at #19 1: model.cross.validation(l[[MEN]]$series) Which points to the call to tslm above. Since I am not supply 'data' to the tslm call (in the forecast package),, I am assuming that the code is dying here (in tslm): if (missing(data)) { dataname - as.character(formula)[2] x - get(dataname) data - data.frame(x) colnames(data) - dataname } Any ideas what is failing? Kevin [[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. -- 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] Random-walk Metropolis-Hasting
On 11.11.2011 18:48, Gyanendra Pokharel wrote: Following is my code, can some one help on the error at the bottom? mh-function(iterations,alpha,beta){ + data-read.table(epidemic.txt,header = TRUE) + attach(data, warn.conflicts = F) + k-97 + d- (sqrt((x-x[k])^2 + (y-y[k])^2)) + p- 1-exp(-alpha*d^(-beta)) + p.alpha-1 - exp(-3*d^(-beta)) + p.beta- 1 - exp(alpha*d^(-2)) + iterations-1000 + mu.lambda- c(0,0);s.lambda- c(100,100) + prop.s- c(0.1,0.1) + lambda- matrix(nrow=iterations, ncol=2) + acc.prob-0 + current.lambda- c(0,0) + for(t in 1:iterations){ + prop.lambda- rnorm(2,current.lambda,prop.s) + a- p.beta/p.alpha *(dnorm(prop.lambda,mu.lambda,s.lambda))*dnorm(current.lambda,mu.lambda,s.lambda) + accept- min(1,a) + u-runif(1) + if(u[t]=accept[t]){ + current.lambda- prop.lambda + acc.prob- acc.prob +1 + } + lambda[t,]- current.lambda + } + lambda + } mh(1000,0,0) Error in if (u[t]= accept[t]) { : missing value where TRUE/FALSE needed As the error message tells you: accept[t] is NA t 1 since accept- min(1,a) Uwe Ligges [[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] symbols and legend: how to harmonize point size ?
Hi, I was wondering if it is possible to harmonize the ouput of symbols() and legend() both from the graphics package. Let us take this example: x-runif(10) y-runif(10) z-runif(10) leg-round(seq(min(z),max(z),l=4),2) # 4 values rounded up to 2 decimals for the legend symbols(x,y,circles=z,inches=0.2) legend(topright,legend=leg,pch=1,pt.cex=leg/max(leg)*2) # multiplied by 2 arbitrarily just to make it visible Actually, what I want to do is to pass to pt.cex a value which would make the biggest circle in the legend (leg/max(leg) = 1) exactly the same size as the one specified in symbols (here 0.2 inches). I suppose this is possible using par(cin) but I cannot figure out how to do it properly. Any hint appreciated, Best, Patrick __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Encryption package of R
Dear all: Is there any encryption package of R? For instance, 3DES or AES algorithm package or whatever... Thanks a lot in advance. Best regards, Jason [[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 Startup options (MDI and SDI)
Here's my setup: - I'm on a Windows machine (I don't have full admin rights) - I have a folder with an *.RData file and an .RProfile file - I want the user to be able to start R by double clicking on the *.RData file Can I specify the application start up options (like --no-save --mdi) in the local RProfile file? 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.
Re: [R] Jordan Form of a matrix
Hi Arnau, Not aware of direct implementation. It was discussed in octave project as well (http://octave.1599824.n4.nabble.com/Jordan-canonical-form-td2216965.html) It is numerically ill-conditioned to compute that. See this http://dl.acm.org/citation.cfm?id=355912 Best, Mehmet Süzen, PhD Mango Solutions data analysis that delivers t: +44 (0) 1249 767700 m: +44 (0) 7517 639318 f: +44 (0) 1249 767707 e: msu...@mango-solutions.com w: www.mango-solutions.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Arnau Mir Torres Sent: 11 November 2011 17:19 To: r-help@r-project.org Subject: [R] Jordan Form of a matrix Hello. Is it possible to find the Jordan Form of a matrix with R? Arnau. Arnau Mir Torres Edifici A. Turmeda Campus UIB Ctra. Valldemossa, km. 7,5 07122 Palma de Mca. tel: (+34) 971172987 fax: (+34) 971173003 email: arnau@uib.es URL: http://dmi.uib.es/~arnau __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. LEGAL NOTICE This message is intended for the use o...{{dropped:10}} __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Formula variable help
It seems that there is a bug in the forecast::tslm function. I have forwarded what I think is the bug (the call to get() should supply the argument 'envir=parent.frame()'). Thank you. On Nov 11, 2011, at 12:11 PM, Joshua Wiley jwiley.ps...@gmail.com wrote: Seems like it could work---can you save your script as a .txt file and send it as an attachment, upload it online, or make a reproducible example? Cheers, Josh On Fri, Nov 11, 2011 at 9:52 AM, rkevinbur...@charter.net wrote: I have an R script with the following applicable lines: xshort - window(s, start=st, end=ed) . . . xshort - ts(xshort, frequency=1, start=1) . . . m1 - m2 - m3 - m4 - m5 - m6 - NULL m1 - tslm(xshort ~ trend) I get an error: Error in get(dataname) : object 'xshort' not found When I do traceback() I get: 3: get(dataname) 2: tslm(xshort ~ trend) at #19 1: model.cross.validation(l[[MEN]]$series) Which points to the call to tslm above. Since I am not supply 'data' to the tslm call (in the forecast package),, I am assuming that the code is dying here (in tslm): if (missing(data)) { dataname - as.character(formula)[2] x - get(dataname) data - data.frame(x) colnames(data) - dataname } Any ideas what is failing? Kevin [[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. -- 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.
[R] multivariate modeling codes
HI, I am relatively new to R and would appreciate some help or directions for this. I am trying to model 3 longitudinal outcomes jointly and to identify some predictors for these 3 joint outcomes (all continuous). I am trying to find some codes that I may modify to do this but cannot seem to find anything. -- View this message in context: http://r.789695.n4.nabble.com/multivariate-modeling-codes-tp4032624p4032624.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
Dear Contributors I would like to perform this operation using a loop, instead of repeating the same operation many times. The numbers from 1 to 4 related to different groups that are in the database and for which I have the same data. x-c(1,3,7) datiP1 - datiP[datiP$city ==1,x]; datiP2 - datiP[datiP$city ==2,x]; datiP3 - datiP[datiP$city ==3,x] datiP4 - datiP[datiP$city ==4,x]; -- Thank you for any help you can provide. Francesca [[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] multivariate modeling codes
Not sure if this helps, but did you try Google? http://www.jstatsoft.org/v35/i09/paper http://cran.r-project.org/web/packages/lcmm/lcmm.pdf http://www.warwick.ac.uk/statsdept/useR-2011/abstracts/010411-liquetbenoit.pdf yurirouge wrote: HI, I am relatively new to R and would appreciate some help or directions for this. I am trying to model 3 longitudinal outcomes jointly and to identify some predictors for these 3 joint outcomes (all continuous). I am trying to find some codes that I may modify to do this but cannot seem to find anything. -- View this message in context: http://r.789695.n4.nabble.com/multivariate-modeling-codes-tp4032624p4032643.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] multivariate modeling codes
Hi, It is not at all clear to me what type of model you are trying to fit. You could consider ?manova for a multivariate analysis of variance require(lme4) ?lmer for longitudinal regression type models or perhaps OpenMx: http://openmx.psyc.virginia.edu/ for structural equation modelling in R. Cheers, Josh On Fri, Nov 11, 2011 at 12:25 PM, yurirouge lilysh...@msn.com wrote: HI, I am relatively new to R and would appreciate some help or directions for this. I am trying to model 3 longitudinal outcomes jointly and to identify some predictors for these 3 joint outcomes (all continuous). I am trying to find some codes that I may modify to do this but cannot seem to find anything. -- View this message in context: http://r.789695.n4.nabble.com/multivariate-modeling-codes-tp4032624p4032624.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] Help
Hi Francesca, Try something like this: x - c(1, 3, 7) dati - lapply(1:4, function(i) {datiP[datiP$city == i, x]}) dati[[1]] # datiP1 dati[[2]] # datiP2 dati[[3]] # datiP3 dati[[4]] # datiP4 if the *only* groups are 1, 2, 3, 4 (i.e., 1:4 is exhaustive), this can be simplified: dati - by(datiP, datiP$city, `[`, x) For documentation, see: ?lapply ?by ?[ # to see how I use the extraction operator as a function Hope this helps, Josh On Fri, Nov 11, 2011 at 12:24 PM, Francesca francesca.panco...@gmail.com wrote: Dear Contributors I would like to perform this operation using a loop, instead of repeating the same operation many times. The numbers from 1 to 4 related to different groups that are in the database and for which I have the same data. x-c(1,3,7) datiP1 - datiP[datiP$city ==1,x]; datiP2 - datiP[datiP$city ==2,x]; datiP3 - datiP[datiP$city ==3,x] datiP4 - datiP[datiP$city ==4,x]; -- Thank you for any help you can provide. Francesca [[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. -- 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] Predicting x from y
On 11-Nov-11 14:51:19, Schreiber, Stefan wrote: Dear list members, I have just a quick question: I fitted a non-linear model y=a/x+b to describe my data (x=temperature and y=damage in %) and it works really nicely (see example below). I have 7 different species and 8 individuals per species. I measured damage for each individual per species at 4 different temperatures (e.g. -5, -10, -20, -40). Using the individuals per species, I fitted one model per species. Now I'd like to use the fitted model to go back and predict the temperature that causes 50% damage (and it's error). Basically, it pretty easy by just rearranging the formula to x=a/(y-b). But that way I don't get a measure of that temperature's error, do I? Can I use the residual standard error that R gave me for the non-linear model fit? Or do I have to fit 8 lines (each individual) per species, calculate x based on the 8 individuals and then take the mean? Unfortunately, dose.p from the MASS package doesn't work for non-linear models. When I take the log(abs(x)) the relationship becomes not satisfactory linear either. Any suggestions are highly appreciated! Thank you! Stefan EXAMPLE for species #1: y.damage-c(5.7388985,1.7813519,3.7321461,2.9671031, 0.3223196,0.3207941,-1.4197658,-5.3472160, 41.1826677,29.3115243,31.3208841,35.3934115, 58.5848778,31.1541049,42.2983479,27.0615648, 64.1037728,54.7003353,66.7317044,65.4725881, 72.5755056,67.2683495,64.8717942,65.9603073, 75.0762273,56.7041960,60.0049429,70.0286506, 73.2801947,72.7015642,75.0944694,81.0361280) x.temp-c(-5,-5,-5,-5,-5,-5,-5,-5,-10,-10,-10,-10,-10,-10,-10, -10,-20,-20,-20,-20,-20,-20,-20,-20,-40,-40,-40,-40,-40, -40,-40,-40) nls(y.damage~a/x.temp+b,start=list(a=400,b=80)) plot(y.damage~x.temp,xlab='Temperature',ylab='Damage [%]') curve(409.61/x+81.84,from=min(x.temp),to=max(x.temp),add=T) A couple of comments. First, in general it is not straightforward to estimate the value of a covariate (here temperature) by inverting the regression of a response (here damage) on that covariate. This the inverse regression or calibration problem, (and it is problematic)! For instance, in linear regression the estimate obtained by inversion has (theoretically) no expectation, and has infinite variance. For an outline, and a few references, see the Wikipedia article: http://en.wikipedia.org/wiki/Calibration_(statistics) Second, I would be inclined to try nls() on a reformulated version of the problem. Let T50 denote the temperature for 50% damage, and introduce this as a parameter (displacing your parameter a): y = 50*(b + T50)/(b + x) where T50 = a/50 - b in terms of your original parameters a and b. With this formula for the non-linear dependence of damage on temperature, it is no longer necessAry to invert the regression equation, since the parameter you want is already there and will be estimated directly. Hoping this helps, Ted. E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 11-Nov-11 Time: 21:15:57 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] (no subject)
Dear colleagues, I'm trying to fit a multinomial logistic regression for an ordinal variable. I see in the help pages for multinom in nnet that one should scale the predictors from 0-1. Is that really necessary? Also: can anyone clarify what the difference between alternative-specific and individual specific variables are? Yours, Simon Kiss l * Simon J. Kiss, PhD Assistant Professor, Wilfrid Laurier University 73 George Street Brantford, Ontario, Canada N3T 2C9 Cell: +1 905 746 7606 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] performance of adaptIntegrate vs. integrate
Dear Ravi, Thank you for your answer. The integrand I proposed was a dummy example for demonstration purposes. I experienced a similar slowdown in a real problem, where knowing in advance the shape of the integrand would not be so easy. Your advice is sound; I would have to study the underlying code of the two implementations to know where the difference lies. Delving into the source code and the algorithms gets quite technical though, so I was hoping someone already familiar with integrate's internals might shed some light. Thanks, baptiste On 12 November 2011 03:55, Ravi Varadhan rvarad...@jhmi.edu wrote: The integrand is highly peaked. It is approximately an impulse function where much of the mass is concentrated at a very small interval. Plot the function and see for yourself. This is the likely cause of the problem. Other types of integrands where you could experience problems are: integrands with singularity at either limit and slowly decaying oscillatory integrands. As to why integrate performs better than adaptIntegrate in this situation, I don’t know. You have to study the two implementations. Wynn’s epsilon algorithm is an extrapolation method for improving the convergence of a sequence. This could be an explanation for the better performance, but I cannot say for sure. Hope this is helpful, Ravi --- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] performance of adaptIntegrate vs. integrate
Dear Hans, [see inline below] On 11 November 2011 22:44, Hans W Borchers hwborch...@googlemail.com wrote: baptiste auguie baptiste.auguie at googlemail.com writes: Dear list, [cross-posting from Stack Overflow where this question has remained unanswered for two weeks] I'd like to perform a numerical integration in one dimension, I = int_a^b f(x) dx where the integrand f: x in IR - f(x) in IR^p is vector-valued. integrate() only allows scalar integrands, thus I would need to call it many (p=200 typically) times, which sounds suboptimal. The cubature package seems well suited, as illustrated below, library(cubature) Nmax - 1e3 tolerance - 1e-4 integrand - function(x, a=0.01) c(exp(-x^2/a^2), cos(x)) adaptIntegrate(integrand, -1, 1, tolerance, 2, max=Nmax)$integral [1] 0.01772454 1.68294197 However, adaptIntegrate appears to perform quite poorly when compared to integrate. Consider the following example (one-dimensional integrand), library(cubature) integrand - function(x, a=0.01) exp(-x^2/a^2)*cos(x) Nmax - 1e3 tolerance - 1e-4 # using cubature's adaptIntegrate time1 - system.time(replicate(1e3, { a - adaptIntegrate(integrand, -1, 1, tolerance, 1, max=Nmax) }) ) # using integrate time2 - system.time(replicate(1e3, { b - integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax) }) ) time1 user system elapsed 2.398 0.004 2.403 time2 user system elapsed 0.204 0.004 0.208 a$integral [1] 0.0177241 b$value [1] 0.0177241 a$functionEvaluations [1] 345 b$subdivisions [1] 10 Somehow, adaptIntegrate was using many more function evaluations for a similar precision. Both methods apparently use Gauss-Kronrod quadrature, though ?integrate adds a Wynn's Epsilon algorithm. Could that explain the large timing difference? Cubature is astonishingly slow here though it was robust and accurate in most cases I used it. You may write to the maintainer for more information. Indeed, I have been a happy user of cubature too, for higher-dimensional problems. And I've used other codes from Steve Johnson before which were all of high quality. I might send an email to both him and the package developer to inquire about this poor performance. Even a pure R implementation of adaptive Gauss-Kronrod as in pracma::quadgk is faster than cubature: library(pracma) time3 - system.time(replicate(1e3, { d - quadgk(integrand, -1, 1) } # 0.0177240954011303 ) ) # user system elapsed # 0.899 0.002 0.897 If your functions are somehow similar and you are willing to dispense with the adaptive part, then compute Gaussian quadrature nodes xi and weights wi for the fixed interval and compute sum(wi * fj(xi)) for each function fj. This avoids recomputing nodes and weights for every call or function. I've used such a technique before, in C++ code intefaced with Rcpp. however, I do like the adaptative part, and implementing it seems less trivial. I don't really want to reinvent the wheel if I can help it find a faster track. Package 'gaussquad' provides such nodes and weights. Or copy the relevant calculations from quadgk (the usual (7,15)-rule). I've used the nodes and weights from statmod::gauss.quad. Thanks for the tips. Best regards, baptiste Regards, Hans Werner I'm open to suggestions of alternative ways of dealing with vector-valued integrands. Thanks. baptiste __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Encryption package of R
Jun Ji junj at stanford.edu writes: Is there any encryption package of R? For instance, 3DES or AES algorithm package or whatever... Thanks a lot in advance. library(sos) findFn(encryption) finds a partial AES implementation, but not much else. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Win upgrade pb (virus)
I would suggestto Mario that, if he wishes to use colorspace and keep AVG, he should install R 2.13.2. The colorspace library built with that version of R does not give a false positive with AVG. To update paackages he must remove not only colorspace but also its reverse dependencies, the reverse dependencies of its dependencies etc. I have done this and am running 2.13.2 and 2.14 in parallel. If you upgrade R to a new version it is always worthwhile to keep the older version until you are sure that the new version can run all the packages that you need. Apart from problems with anti-virus software there is a posibility that there may be a delay in one of your packages or a dependency becoming available. I have reported the false positive to AVG and have received an acknowledgment from them. The AVG program contains a facility for reporting false positives. Please use it and put a little more pressure on them to update their software. The web site http://virusscan.jotti.org allows you to submit a problem file to 20 different anti-virus programs. Only 1 program (AVG) finds a problem in colorspace.dll. The other 19 pass the file. Best Regards John On 11 November 2011 10:54, Prof Brian Ripley rip...@stats.ox.ac.uk wrote: On Fri, 11 Nov 2011, Mario Valle wrote: I just upgraded my Win7 32bits installation to 2.14.0 after deinstalling 2.12.x First thing I moved the win-library from 2.12 to 2.14 and executed a update.packages(ask='graphics',checkBuilt=TRUE) (Swiss mirror). This aborts with the console message: Error in if (any(diff)) { : missing value where TRUE/FALSE needed And the antivirus (AVG) pops up complaining that colorspace.dll contains the virus Win32/Heur I cannot ignore the thread because update.packages has already crashed when I can reach the antivirus dialogbox to push the ignore button. To continue I had to remove the colorspace and ggplot2 packages. Is it a known problem? What else can I do (except removing the above packages or changing antivirus)? Yes (a known problem in AVG), and the standard advice is to change antivirus (or set exceptions, which I gather AVG does not allow). And please report the false positive to AVG: the more they hear about the more attention they will give it. The binary distribution has been checked by e.g. Sophos (which both winbuilder and ox.ac.uk use). Thanks! mario -- Ing. Mario Valle Data Analysis and Visualization Group | http://www.cscs.ch/~mvalle Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60 v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82 -- 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, UK Fax: +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. -- 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.
[R] Combining Overlapping Data
I've scoured the archives but have found no concrete answer to my question. Problem: Two data sets 1st data set(x) = 20,000 rows 2nd data set(y) = 5,000 rows Both have the same column names, the column of interest to me is a variable called strain. For example, a strain named Chab1405 appears in x 150 times and in y 25 times... strain Chab1999 only appears 200 times in x and none in y (so i dont want that retained). I want to create a new data frame that has all 175 measurements for Chab1405 and any other 'strain' that appears in both the two data sets.. but not strains that appear in only one data set...So i want the intersection of two data sets (maybe?). I've tried x %in% y, but that only gives TRUE/FALSE -- View this message in context: http://r.789695.n4.nabble.com/Combining-Overlapping-Data-tp4032719p4032719.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] need help
hello all R experts, how do I calculate the reliability between the two groups using the ICCs? I'll appreciate your reply, Thanks Sincerely, Supreet kaur, Biomedical research engineer, Nationwide Childrens Hospital, Columbus, OH (614)355-3509 [[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] list.dir() function
Hi I have an organism directory that contains two folders galGal3 and hg19 and many other files. orgDir = '/home/mary/org' When I try to use list.dir() function, it gives me the same answer, no matter what is the value of full.names argument. list.dirs(path = indexDir, full.names = FALSE)[1] /home/mary/org [2] /home/mary/org/galGal3 [3] /home/mary/org/hg19 list.dirs(path = indexDir, full.names = TRUE) [1] /home/mary/org [2] /home/mary/org/galGal3 [3] /home/mary/org/hg19 Also, It prints the directory itself which I don't want to be printed. Why it is so? Any workaround for this problem? Thanks -- - Mary Kindall Yorktown Heights, NY 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.
Re: [R] Predicting x from y
Follow-up: See at end. On 11-Nov-11 21:16:02, Ted Harding wrote: On 11-Nov-11 14:51:19, Schreiber, Stefan wrote: Dear list members, I have just a quick question: I fitted a non-linear model y=a/x+b to describe my data (x=temperature and y=damage in %) and it works really nicely (see example below). I have 7 different species and 8 individuals per species. I measured damage for each individual per species at 4 different temperatures (e.g. -5, -10, -20, -40). Using the individuals per species, I fitted one model per species. Now I'd like to use the fitted model to go back and predict the temperature that causes 50% damage (and it's error). Basically, it pretty easy by just rearranging the formula to x=a/(y-b). But that way I don't get a measure of that temperature's error, do I? Can I use the residual standard error that R gave me for the non-linear model fit? Or do I have to fit 8 lines (each individual) per species, calculate x based on the 8 individuals and then take the mean? Unfortunately, dose.p from the MASS package doesn't work for non-linear models. When I take the log(abs(x)) the relationship becomes not satisfactory linear either. Any suggestions are highly appreciated! Thank you! Stefan EXAMPLE for species #1: y.damage-c(5.7388985,1.7813519,3.7321461,2.9671031, 0.3223196,0.3207941,-1.4197658,-5.3472160, 41.1826677,29.3115243,31.3208841,35.3934115, 58.5848778,31.1541049,42.2983479,27.0615648, 64.1037728,54.7003353,66.7317044,65.4725881, 72.5755056,67.2683495,64.8717942,65.9603073, 75.0762273,56.7041960,60.0049429,70.0286506, 73.2801947,72.7015642,75.0944694,81.0361280) x.temp-c(-5,-5,-5,-5,-5,-5,-5,-5,-10,-10,-10,-10,-10,-10,-10, -10,-20,-20,-20,-20,-20,-20,-20,-20,-40,-40,-40,-40,-40, -40,-40,-40) nls(y.damage~a/x.temp+b,start=list(a=400,b=80)) plot(y.damage~x.temp,xlab='Temperature',ylab='Damage [%]') curve(409.61/x+81.84,from=min(x.temp),to=max(x.temp),add=T) A couple of comments. First, in general it is not straightforward to estimate the value of a covariate (here temperature) by inverting the regression of a response (here damage) on that covariate. This the inverse regression or calibration problem, (and it is problematic)! For instance, in linear regression the estimate obtained by inversion has (theoretically) no expectation, and has infinite variance. For an outline, and a few references, see the Wikipedia article: http://en.wikipedia.org/wiki/Calibration_(statistics) Second, I would be inclined to try nls() on a reformulated version of the problem. Let T50 denote the temperature for 50% damage, and introduce this as a parameter (displacing your parameter a): y = 50*(b + T50)/(b + x) where T50 = a/50 - b in terms of your original parameters a and b. With this formula for the non-linear dependence of damage on temperature, it is no longer necessAry to invert the regression equation, since the parameter you want is already there and will be estimated directly. Hoping this helps, Ted. I think I have mis-read your model: I read it as y = a/(x+b) whereas you wrote y/x+b and your model formula in nls() is y.damage~a/x.temp+b, i.e. y.damage ~ (a/x.temp) + b which confirms it. In that case, you may be able to get a satisfactory result by using a linear regression with y.damage = a*z.temp + b where z.temp = 1/x.temp so the model formula would be y.damage ~ z.temp You then have the straightforward inverse regression problem (aka calibration problem). The solution to this takes a bit of explanation, which I do not have the time for right now. I will write further about it in the morning. Ted. E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 11-Nov-11 Time: 23:07:35 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Combining Overlapping Data
What about merge() with all=FALSE? x - data.frame(a=letters[1:6], b=1:6) y - data.frame(a=letters[4:9], b=11:16) x a b 1 a 1 2 b 2 3 c 3 4 d 4 5 e 5 6 f 6 y a b 1 d 11 2 e 12 3 f 13 4 g 14 5 h 15 6 i 16 merge(x, y, by=a, all=FALSE) a b.x b.y 1 d 4 11 2 e 5 12 3 f 6 13 If that doesn't work, some sample data would be useful. Sarah On Fri, Nov 11, 2011 at 4:07 PM, kickout kyle.ko...@gmail.com wrote: I've scoured the archives but have found no concrete answer to my question. Problem: Two data sets 1st data set(x) = 20,000 rows 2nd data set(y) = 5,000 rows Both have the same column names, the column of interest to me is a variable called strain. For example, a strain named Chab1405 appears in x 150 times and in y 25 times... strain Chab1999 only appears 200 times in x and none in y (so i dont want that retained). I want to create a new data frame that has all 175 measurements for Chab1405 and any other 'strain' that appears in both the two data sets.. but not strains that appear in only one data set...So i want the intersection of two data sets (maybe?). I've tried x %in% y, but that only gives TRUE/FALSE -- 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.
Re: [R] need help
On Fri, Nov 11, 2011 at 4:41 PM, Supreet kaur ksupre...@gmail.com wrote: hello all R experts, how do I calculate the reliability between the two groups using the ICCs? Possibly by using GmeanRel() from the multilevel package. Searching for reliability group ICC at http://www.rseek.org offers several other options, one of which may be more appropriate for your unspecified application. Sarah -- 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.
Re: [R] Combining Overlapping Data
Use 'intersect' to get the items common in both dataframes and then use that to extract the data in common. On Friday, November 11, 2011, kickout kyle.ko...@gmail.com wrote: I've scoured the archives but have found no concrete answer to my question. Problem: Two data sets 1st data set(x) = 20,000 rows 2nd data set(y) = 5,000 rows Both have the same column names, the column of interest to me is a variable called strain. For example, a strain named Chab1405 appears in x 150 times and in y 25 times... strain Chab1999 only appears 200 times in x and none in y (so i dont want that retained). I want to create a new data frame that has all 175 measurements for Chab1405 and any other 'strain' that appears in both the two data sets.. but not strains that appear in only one data set...So i want the intersection of two data sets (maybe?). I've tried x %in% y, but that only gives TRUE/FALSE -- View this message in context: http://r.789695.n4.nabble.com/Combining-Overlapping-Data-tp4032719p4032719.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. [[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] Predicting x from y
Thanks Ted! I really appreciate your time! Thanks for the link about the 'problem of calibration', and your suggestion to reformulate my model. I had no idea about it before. I certainly learnt something today. I will try your suggestions later today and let you know how it works out. Stefan -Original Message- From: r-help-boun...@r-project.org on behalf of ted.hard...@wlandres.net Sent: Fri 11/11/2011 4:07 PM To: r-help@r-project.org Subject: Re: [R] Predicting x from y Follow-up: See at end. On 11-Nov-11 21:16:02, Ted Harding wrote: On 11-Nov-11 14:51:19, Schreiber, Stefan wrote: Dear list members, I have just a quick question: I fitted a non-linear model y=a/x+b to describe my data (x=temperature and y=damage in %) and it works really nicely (see example below). I have 7 different species and 8 individuals per species. I measured damage for each individual per species at 4 different temperatures (e.g. -5, -10, -20, -40). Using the individuals per species, I fitted one model per species. Now I'd like to use the fitted model to go back and predict the temperature that causes 50% damage (and it's error). Basically, it pretty easy by just rearranging the formula to x=a/(y-b). But that way I don't get a measure of that temperature's error, do I? Can I use the residual standard error that R gave me for the non-linear model fit? Or do I have to fit 8 lines (each individual) per species, calculate x based on the 8 individuals and then take the mean? Unfortunately, dose.p from the MASS package doesn't work for non-linear models. When I take the log(abs(x)) the relationship becomes not satisfactory linear either. Any suggestions are highly appreciated! Thank you! Stefan EXAMPLE for species #1: y.damage-c(5.7388985,1.7813519,3.7321461,2.9671031, 0.3223196,0.3207941,-1.4197658,-5.3472160, 41.1826677,29.3115243,31.3208841,35.3934115, 58.5848778,31.1541049,42.2983479,27.0615648, 64.1037728,54.7003353,66.7317044,65.4725881, 72.5755056,67.2683495,64.8717942,65.9603073, 75.0762273,56.7041960,60.0049429,70.0286506, 73.2801947,72.7015642,75.0944694,81.0361280) x.temp-c(-5,-5,-5,-5,-5,-5,-5,-5,-10,-10,-10,-10,-10,-10,-10, -10,-20,-20,-20,-20,-20,-20,-20,-20,-40,-40,-40,-40,-40, -40,-40,-40) nls(y.damage~a/x.temp+b,start=list(a=400,b=80)) plot(y.damage~x.temp,xlab='Temperature',ylab='Damage [%]') curve(409.61/x+81.84,from=min(x.temp),to=max(x.temp),add=T) A couple of comments. First, in general it is not straightforward to estimate the value of a covariate (here temperature) by inverting the regression of a response (here damage) on that covariate. This the inverse regression or calibration problem, (and it is problematic)! For instance, in linear regression the estimate obtained by inversion has (theoretically) no expectation, and has infinite variance. For an outline, and a few references, see the Wikipedia article: http://en.wikipedia.org/wiki/Calibration_(statistics) Second, I would be inclined to try nls() on a reformulated version of the problem. Let T50 denote the temperature for 50% damage, and introduce this as a parameter (displacing your parameter a): y = 50*(b + T50)/(b + x) where T50 = a/50 - b in terms of your original parameters a and b. With this formula for the non-linear dependence of damage on temperature, it is no longer necessAry to invert the regression equation, since the parameter you want is already there and will be estimated directly. Hoping this helps, Ted. I think I have mis-read your model: I read it as y = a/(x+b) whereas you wrote y/x+b and your model formula in nls() is y.damage~a/x.temp+b, i.e. y.damage ~ (a/x.temp) + b which confirms it. In that case, you may be able to get a satisfactory result by using a linear regression with y.damage = a*z.temp + b where z.temp = 1/x.temp so the model formula would be y.damage ~ z.temp You then have the straightforward inverse regression problem (aka calibration problem). The solution to this takes a bit of explanation, which I do not have the time for right now. I will write further about it in the morning. Ted. E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 11-Nov-11 Time: 23:07:35 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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
Re: [R] multivariate modeling codes
Hi: R-Bloggers picked up on this web site that contains several interesting posts re usage of R, including this one on using a Bayesian approach to multivariate mixed effects models using the MCMCglmm package: http://www.quantumforest.com/2011/11/coming-out-of-the-bayesian-closet-multivariate-version/ With multiple responses in a longitudinal study, this might be of interest to you. A similar question re multivariate response was asked a few days ago on R-sig-mixed-models, where Kevin Wright posted the above link. I'd suggest that future questions on this subject be moved to that group. HTH, Dennis On Fri, Nov 11, 2011 at 12:25 PM, yurirouge lilysh...@msn.com wrote: HI, I am relatively new to R and would appreciate some help or directions for this. I am trying to model 3 longitudinal outcomes jointly and to identify some predictors for these 3 joint outcomes (all continuous). I am trying to find some codes that I may modify to do this but cannot seem to find anything. -- View this message in context: http://r.789695.n4.nabble.com/multivariate-modeling-codes-tp4032624p4032624.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] Combining Overlapping Data
Hi: This doesn't sort the data by strain level, but I think it does what you're after. It helps if strain is either a factor or character vector in each data frame. h - function(x, y) { tbx - table(x$strain) tby - table(y$strain) # Select the strains who have more than one member # in each data frame mgrps - intersect(names(tbx[tbx 0]), names(tby[tby 0])) # concatenate the data with common strains rbind(subset(x, gp %in% mgrps), subset(y, gp %in% mgrps)) } # Result: dc - h(x, y) HTH, Dennis On Fri, Nov 11, 2011 at 1:07 PM, kickout kyle.ko...@gmail.com wrote: I've scoured the archives but have found no concrete answer to my question. Problem: Two data sets 1st data set(x) = 20,000 rows 2nd data set(y) = 5,000 rows Both have the same column names, the column of interest to me is a variable called strain. For example, a strain named Chab1405 appears in x 150 times and in y 25 times... strain Chab1999 only appears 200 times in x and none in y (so i dont want that retained). I want to create a new data frame that has all 175 measurements for Chab1405 and any other 'strain' that appears in both the two data sets.. but not strains that appear in only one data set...So i want the intersection of two data sets (maybe?). I've tried x %in% y, but that only gives TRUE/FALSE -- View this message in context: http://r.789695.n4.nabble.com/Combining-Overlapping-Data-tp4032719p4032719.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] Time series analysis with random effects
Hi folks, i have some problems with my evaluation. We have collect tons of data from 23 testpersons for our new road study. I have now a time series for each person and all the logs when he accelerates or hits the break trying to solve five different tasks. The dataset lools like: #Unixtime time accelerate break UserID task #10372 1312358742 10:05:41.600 3.054 0.000 2 lane #10373 1312358742 10:05:41.700 3.056 1.000 2 emergency The break variable is binary and i want to test if there is a significant difference between the tasks using the break variable. I have not really a clue what i have to look for. My idea was to use a 2 Sample t-test for the accelerate variable but the conditions are not fulfilled for this test. I am not really experienced with time series analysis so thanks a lot for your advice! Bernd -- View this message in context: http://r.789695.n4.nabble.com/Time-series-analysis-with-random-effects-tp4033263p4033263.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] predictive apriori
Dear list members, I know that there is the arules package with the implementation of the apriori algorithm. However i want to use the predictive apriori instead. These algorithm can mine as rules as i want and there is an implementation on weka. There is some implementation on R? -- Att, Flávio Barros [[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] Encryption package of R
On Nov 11, 2011, at 4:04 PM, Ben Bolker wrote: Jun Ji junj at stanford.edu writes: Is there any encryption package of R? For instance, 3DES or AES algorithm package or whatever... Thanks a lot in advance. library(sos) findFn(encryption) finds a partial AES implementation, but not much else. Hi, When it comes to encryption, you are best to use a library that has a good history behind it. I would suggest looking at http://www.gnupg.org/ and perhaps specifically, libgcrypt (http://directory.fsf.org/wiki/Libgcrypt). HTH, Marc Schwartz __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Odds ratios from lrm plot
The code library(Design) f - lrm(y~x1+x2+x1*x2, data=data) plot(f) produces a plot of log odds vs x2 with 0.95 confidence intervals. How do I get a plot of odds ratios vs x2 instead? Thanks -- View this message in context: http://r.789695.n4.nabble.com/Odds-ratios-from-lrm-plot-tp4033340p4033340.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] Odds ratios from lrm plot
The Design package is obsolete. Use its replacement rms - see http://biostat.mc.vanderbilt.edu/Rrms. Use something like the following to reproduce what you already have: To plot odds ratios against a specific value of x2 and at a specific x1, type ?contrast.rms to see examples. Frank RD235 wrote: The code library(Design) f - lrm(y~x1+x2+x1*x2, data=data) plot(f) produces a plot of log odds vs x2 with 0.95 confidence intervals. How do I get a plot of odds ratios vs x2 instead? Thanks - Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/Odds-ratios-from-lrm-plot-tp4033340p4033414.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] Odds ratios from lrm plot
On Nov 11, 2011, at 9:08 PM, RD235 wrote: The code library(Design) f - lrm(y~x1+x2+x1*x2, data=data) plot(f) produces a plot of log odds vs x2 with 0.95 confidence intervals. How do I get a plot of odds ratios vs x2 instead? You would construct a dataset that had selected combinations of x1 and x2 and then plot the response. It's been a couple of years since I used Design. Harrell migrated to lattice plots when he upgraded to 'rms' and changed the Predict/plot interface a bit. In the new version it might something like plot(predict(f, x1=seq( ...), x2=c(1,2), type= fitted), deopending on teh data arrangement. Why don't you look more closely at the help page for `lrm` and work through the examples. At least in the current 'rms' version there is more than one worked example, and I'm pretty sure there were one in Design before that. -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] fitting a distribtution to event data
Hi all I am trying to fit a distribution (i.e. gamma) to some data and I understand how to use the fitdistr from the MASS package. However, what should you do when the data has a probability associated with it that are not all equal. e.g. Pr(Occurrence) Size 0.00460 0.02 30 0.00670 0.06 40 regards -- View this message in context: http://r.789695.n4.nabble.com/fitting-a-distribtution-to-event-data-tp4033760p4033760.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] Odp: Error in matrix, not ordered vectors or numerical value, and SIAR.
Hello Petr, The demo's don't run either, with the same errors. Thanks for your help. Best wishes. -- View this message in context: http://r.789695.n4.nabble.com/Error-in-matrix-not-ordered-vectors-or-numerical-value-and-SIAR-tp4024578p4033682.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] lm: how are polynomial functions interpreted?
A http://polynomial-trimonial-binomial.blogspot.com/2011/10/how-to-deal-with-polynomials.html polynomial function is evaluated by foiling out the equation so you can solve for 'x'. -- View this message in context: http://r.789695.n4.nabble.com/lm-how-are-polynomial-functions-interpreted-tp875663p4033696.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.