Re: [R] help: cannot allocate vector of length 828310236
Does it make any statistical sense to do polr or probit regression (not the same thing) on `really huge data'? There are few regression-like problems in which model inadequacy does not swamp estimation uncertainty for as few as a 1000 cases. If you want to do that sort of thing, by all means use SAS to do it. But if you are not prepared to spend a few $$ on adequate RAM, don't expect free technical consultancy, especially not from those whose work you are using and not crediting. - The uncredited author of polr(). On Mon, 14 Aug 2006, T Mu wrote: Hi all, I was trying a probit regression using polr() and got this message, polr is a strange choice of tool for 'probit regression' as the term is usually used. It does 'ordered probit regression'. Error in model.matrix.default(Terms, m, contrasts) : cannot allocate vector of length 828310236 The data is about 20M (a few days ago I asked a question about large file, thank you for responses, then I use MS Access to select those columns I would use). R is 2.3.1, Windows XP, 512M Ram. I am going to read some help on memory use in R, but hope anybody can give me some quick hints. Quick hint: read and follow the posting guide BEFORE posting. Is it because iphysical memory runs out, or some other things could be wrong with data or polr()? Does R use virtual memory? If so, what options can I set? If not, can R deal with really huge data (except adding RAM according to data size)? If this is the case, it is too bad that I have to tell my boss to go back to SAS. Now it is not a speed issue yet. Thank you. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help with workaround for: Function '`[`' is not in thederivatives table
Earl F. Glynn asks: -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Earl F. Glynn Sent: Tuesday, 15 August 2006 8:44 AM To: r-help@stat.math.ethz.ch Subject: [R] Help with workaround for: Function '`[`' is not in thederivatives table # This works fine: a - 1 b - 2 c - 3 E - expression(a * exp(b*X) + c) X - c(0.5, 1.0, 2.0) eval(E) [1] 5.718282 10.389056 57.598150 D(E, b) a * (exp(b * X) * X) eval(D(E, b)) [1] 1.359141 7.389056 109.196300 # But if (a,b,c) are replaced with (A[1], A[2], A[3]), how can I get a derivative using D? It's well to note what D can differentiate with respect to. The second argument is, to quote the help page, a character string giving the name of a variable... 'A[1]' is a character string, but it is not the name of a variable. When parsed it becomes a call to the function, literally, `[`. A - c(1, 2, 3) E - expression(A[1] * exp(A[2]*X) + A[3]) X - c(0.5, 1.0, 2.0) eval(E) [1] 5.718282 10.389056 57.598150 No problem, because the evaluator does know how to evaluate `[`, along with everything else in this expr. # Why doesn't this work? Any workarounds? D(E, A[2]) Error in D(E, A[2]) : Function '`[`' is not in the derivatives table It fails because 'A[2]' is not a (character string giving the) name of a variable. I think the error message could be a bit more informative, but ... all you need to do is read he help page, really. If I want to have a long vector of coefficients, A, (perhaps dozens) how can I use D to compute partial derivatives? Essentially you need to turn calls to '[' into names of variables, do the derivative business and then turn them back again. This is not easy to do in complete generality, but if you are only talking about singly subscripted arrays, you can get somewhere. Here is an outline of what I mean. Ident - ([A-z][A-z0-9_.]*) Subsc - ([A-z][A-z0-9_.]*|[0-9]+) patn - paste(Ident, \\[, Subsc, \\], sep = ) repl - \\1__\\2 E - expression(A[1] * exp(A[2]*X) + A[3]) Es - deparse(E[[1]]) Es [1] A[1] * exp(A[2] * X) + A[3] Ess - gsub(patn, repl, Es) Ess [1] A__1 * exp(A__2 * X) + A__3 Ex - parse(text = Ess)[[1]] Ex A__1 * exp(A__2 * X) + A__3 OK, the calls to `[` have been replaced by variables with two underscores in the middle. We hope this works - there is a strong assumption here on just how complicated your indices are, for example. We are assuming they are either identifiers (not using two successive underscores in their names) or numbers. If they are not, you need to get even craftier. Ex1 - D(Ex, A__2) Ex1 A__1 * (exp(A__2 * X) * X) Ex1s - deparse(Ex1) Ex1s [1] A__1 * (exp(A__2 * X) * X) pat1 - paste(Ident, __, Subsc, sep = ) rep1 - \\1\\[\\2\\] Ex1ss - gsub(pat1, rep1, Ex1s) Ex1ss [1] A[1] * (exp(A[2] * X) * X) Ex2 - parse(text = Ex1ss)[[1]] Ex2 A[1] * (exp(A[2] * X) * X) Which is the required result. This is messy and gets messier if you are looking for some kind of generality, but you need to remember, R is not, and never will be, a replacement for Maple. Thanks for any help with this. Best of luck in automating this, but the tools are there. Bill Venables. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to call forth a class definition buried in a package
Dear list, I would like to use this C call: tmp - .Call(R_MPinv, covariance, tol, svdmem) but it gives me the following error: Error in getClass(LinStatExpectCovarMPinv) : LinStatExpectCovarMPinv is not a defined class I think the reason for this is that the C code starts with this call: PROTECT(ans = NEW_OBJECT(MAKE_CLASS(LinStatExpectCovarMPinv))); and the Class definition for LinStatExpectCovarMPinv hides deep buried in package party. Were it not for the C call, i could avoid this error by specifying a where= argument to GetClass(). But what can I do now? Thank you! Bálint -- Czúcz Bálint PhD hallgató BCE KTK Talajtan és Vízgazdálkodás Tanszék 1118 Budapest, Villányi út 29-43. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Question on .Options$max.print - print/cat extremely long strings on a screen
Dear list members, sorry for my incompleteness! My problem is the following: (R 2.3.1 on Mac OS X 10.4.7 RAM 1GByte using Mac GUI) I have a function like this: foo1 - function() { out - NULL for(i in 1:10010) out - paste(out, i, . line\n, sep=) return(out) } a - foo1() Now I want to display 'a' on the screen: cat(a) This doesn't work. 'a' is displayed only until '830. line' print(a) The same, 'a' is displayed only until '\n754. line\n755'. cat(a, file='test.txt') # OK works fine. That means, internally 'a' is fine. Then I tried this way: foo2 - function() { out - NULL for(i in 1:10010) out - c(out, paste(i, . line, sep=)) return(out) } a - foo2() With cat(a,sep=\n) I see the complete content of 'a'. paste(a, collapse = \n) I only see 'a' until '\n754. line\n755'. cat(paste(a, collapse =\n)) I only see 'a' until '830. line'. cat(paste(a, collapse =\n),file='test.txt') This is OK. My question now is whether there is an option to specifiy the maximum size of a string which is displayed on the screen (running R in a GUI)? Or is this fixed? I read the help page about '.Options' and I found a variable 'max.print' with the comment that is not yet used in basic R. I don't know whether this variable is responsible for that. I increased 'max.print' but nothing changed. ## If I try this code on a Windows XP machine with 756 MByte RAM R-GUI says after executing foo1 - function() { out - NULL for(i in 1:10010) out - paste(out, i, . line\n, sep=) return(out) } a - foo1() cat(a) ... 7548. lineWarning message: printing of extremely long output is truncated At least Windows writes a warning message. ### If I start a R session without the R-GUI via Mac Terminal typing 'R' a-foo1() cat(a) everything works perfectly!!! I know the issue of outputting long strings and the way to display long strings via foo2() would be ok for me, but I spent some time to figure out why my function foo1() didn't work. R, running in GUI on Mac, don't give you a warning. Now I know that if I print a long string at the Mac-GUI-console the missing final quote character is an indicator for a truncated output on a Mac. Maybe it would be nice to output a warning(?) Cheers, Hans __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Protection stack overflow
G'day all. I'm a new user of R -- but an arms-length user, as I'm running it from Octave via the ROctave interface that Duncan Temple Lang wrote some years ago and Stephan Van Der Walt recently updated for use with Octave 2.1.71. I'm using R version 2.1.1. ROctave uses libR.so to provide the interface. My system is Ubuntu linux 5.10 and I'm using the packages that come with this distro. I'm getting a protection stack overflow error when recursively calculating AR(p) time series models using the arima() function. The recursion is involved because I calculate a new model with each new time series data point. (I'm trying to reproduce some results in a research paper and this is what they're doing.) I've tried setting the expressions parameter to a higher number using options(expressions=50) but I'm still getting this problem with stack overflow. I can get about 400-odd iterations and then the overflow error appears. I yet haven't tried running the recursion natively in R, and I realise I should do that. Your advice would be much appreciated. Cheers, Paul. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Query the kind of calling a funcion
Dear all, My question is concerned to the kind how a function is called. Example A: foo(1) Example B: a - foo(1) Is there any way for the function foo() to recognise whether the returned value of foo() is stored in a variable or not, i.e. to distinguish between Example A and B? Any comments are welcomed. Many thanks in advance Hans __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Protection stack overflow
R has a command-line option to set the ppstack size, --max-ppsize=NSet max size of protect stack to N Looks like you need to supply this (and it can be done with embedded R) if the problem persists with current R. You could also try arima0 or even ar to do the fitting. On Tue, 15 Aug 2006, Paul Koufalas wrote: G'day all. I'm a new user of R -- but an arms-length user, as I'm running it from Octave via the ROctave interface that Duncan Temple Lang wrote some years ago and Stephan Van Der Walt recently updated for use with Octave 2.1.71. I'm using R version 2.1.1. ROctave uses libR.so to provide the interface. My system is Ubuntu linux 5.10 and I'm using the packages that come with this distro. Note that your version of R is well outdated, and the posting guide did ask you to update it BEFORE posting. I'm getting a protection stack overflow error when recursively calculating AR(p) time series models using the arima() function. The recursion is involved because I calculate a new model with each new time series data point. (I'm trying to reproduce some results in a research paper and this is what they're doing.) I've tried setting the expressions parameter to a higher number using options(expressions=50) but I'm still getting this problem with stack overflow. I can get about 400-odd iterations and then the overflow error appears. I yet haven't tried running the recursion natively in R, and I realise I should do that. Your advice would be much appreciated. Cheers, Paul. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to reply to a thread if .. R-help mails .. in digest
Thanks a lot, Ted, for the good extensive explanation. Let me try summarize, confirm and add a bit w/o repeating too much: - If you get regular mailing list e-mails, and reply to postings, any decent mail software will take the 'Message-Id:' header of the message you reply to, and produce 'In-Reply-To:' and 'References:' headers from it, and will add them to your own e-mail. [ Most e-mail softwares keep these headers hidden, and quite a few pieces of e-mail crapware don't even allow you to see these headers.] { And as Ted mentioned; unfortunately, we still have too many r-help posters who **misuse** their e-mail-software's `reply' and then inadvertently hijack threads... } - Yes, the threads are *not* built from 'Subject:'s but rather using the e-mail headers 'References:' and/or 'In-Reply-To:' This is all based on international e-mail standards (RFC/..) mostly going back into the age where most R-help readers have not ever used e-mail.. - With mailman, there are 2 ways to receive digests: (1) plain text --- which is default, since some dumb e-mail programs can only deal with those (2) MIME --- which uses the MIME standard to send the digest. With a good e-mail software, this means that you get ``each message as attachment'' {that's how it typically looks to the user} which you then can open - in your mail software(!) - and it will behave like a regular e-mail; it has a (typically hidden) 'Message-ID:' etc == when replying to such a message you automatically get both: 1) correct Subject 2) correct In-Reply-To + References === correct thread *So* : What we (and many) recommend to all digest subscribers is to activate the MIME option - on their mailing list membership info page to where you get from https://stat.ethz.ch/mailman/listinfo/r-help after entering your e-mail address at the very bottom of the page (and then your list password). -- but, as Michael Dewey just mentioned, unfortunately there are (too many) pieces of e-mail crapware around which even don't support the above correctly... Martin Maechler, ETH Zurich, provider of (most of) the mailing list infrastructure for R. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] split a y-axis to show data on different scales
The pro's and con's of using scale breaks were discussed by Cleveland (1985) The Elements of Graphing Data (Wadsworth, pp. 85-91, 149). I don't know what Cleveland said about this is the second edition Spencer Graves: but I believe there are times when scale breaks are appropriate, but the display should make this nonstandard transition very clear; ... in which case you are close to having two graphs sharing an x-axis and therefore saving on ink (yay!). __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] A plot with a bisector
Dear R Users, How to make a plot with a bisector is my question. Indeed, I want to plot a QQplot. But it doesn't have such a bisector. Thanks a lot. Amir - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Aliases for arguments in a function
Hi all. I have a function that I would like to use either the argument name as originally defined or another name. Within the function (and other functions) use the argument name as originally written, so I don't want to simply remove the old argument name for the new one, but simply allow the function to treat both argument names as equivalent. Here is an example: foo - function(arg1, this) { if(this 0) stop(this must be positive) return(arg1/this) } foo(arg1=5, this=10) But, I also want foo() to work equivalently with the following (where 'this' and 'that 'are treated as if they were the same): foo(arg1=5, that=10) Any thoughts would be appreciated. Thanks, Ken - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Aliases for arguments in a function
foo - function(arg1, this, that) { if (missing(this) !missing(that)) this - that if(this 0) stop(this must be positive) return(arg1/this) } foo(arg1=5, this=10) foo(arg1=5, that=10) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] nls
Is there anyway to change any y[i] value (i=2,...6) to make following NLS workable? x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) lm(1/y ~~ x) nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE) #0.0920573 : 1.16122 0.01565 1.0 #Error in numericDeriv(form[[3]], names(ind), env) : #Missing value or an infinity produced when evaluating the model - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls
You problem is x^c for x = 0. If you intended only c 1, try a starting value meeting that condition (but it seems that the optimal c is about 0.27 is you increase x slightly). Why have you used ~~ ? (Maybe because despite being asked not to, you sent HTML mail?) On Tue, 15 Aug 2006, Xiaodong Jin wrote: Is there anyway to change any y[i] value (i=2,...6) to make following NLS workable? x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) lm(1/y ~~ x) nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE) #0.0920573 : 1.16122 0.01565 1.0 #Error in numericDeriv(form[[3]], names(ind), env) : #Missing value or an infinity produced when evaluating the model - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Aliases for arguments in a function
foo - function(arg1, this=that, that) ... works. On Tue, 15 Aug 2006, [EMAIL PROTECTED] wrote: Hi all. I have a function that I would like to use either the argument name as originally defined or another name. Within the function (and other functions) use the argument name as originally written, so I don't want to simply remove the old argument name for the new one, but simply allow the function to treat both argument names as equivalent. Here is an example: foo - function(arg1, this) { if(this 0) stop(this must be positive) return(arg1/this) } foo(arg1=5, this=10) But, I also want foo() to work equivalently with the following (where 'this' and 'that 'are treated as if they were the same): foo(arg1=5, that=10) Any thoughts would be appreciated. Thanks, Ken - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] ARCH Jump diffusion models
Hi I was wondering if you could tell me how would I go abt fitting the following model to a data series. It is an ARCH jump diffusion model. X(t)=k+ ó (t)z(t)+J x Y, and ó(t)= sqrt[a +b [X(t-1)-k]^2] Here k is mean of the series. z(t) is iid standard normal error J is a jump variable (either 1 or 0) Y is normally distributed with mean m and variance v Regards Avishek If you are not an intended recipient of this e-mail, please notify the sender, delete it and do not read, act upon, print, disclose, copy, retain or redistribute it. Click here for important additional terms relating to this e-mail. http://www.ml.com/email_terms/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls
Prof Brian Ripley [EMAIL PROTECTED] writes: You problem is x^c for x = 0. If you intended only c 1, try a starting value meeting that condition (but it seems that the optimal c is about 0.27 is you increase x slightly). Surely you mean c 0. nls(1/y ~ a+b*x^exp(c), start=list(a=1.16122,b=0.01565,c=0)) Nonlinear regression model model: 1/y ~ a + b * x^exp(c) data: parent.frame() a b c 0.9944025 0.1953168 -1.1495206 residual sum-of-squares: 0.03303464 nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=exp(-1.1))) Nonlinear regression model model: 1/y ~ a + b * x^c data: parent.frame() a b c 0.9944026 0.1953165 0.3167891 residual sum-of-squares: 0.03303464 (but even setting c=exp(-1) triggers the error; there could be cause for robustifying the nls algorithm) Why have you used ~~ ? (Maybe because despite being asked not to, you sent HTML mail?) On Tue, 15 Aug 2006, Xiaodong Jin wrote: Is there anyway to change any y[i] value (i=2,...6) to make following NLS workable? x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) lm(1/y ~~ x) nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE) #0.0920573 : 1.16122 0.01565 1.0 #Error in numericDeriv(form[[3]], names(ind), env) : #Missing value or an infinity produced when evaluating the model - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls
Hi Why do you want to change your variable values? It smells a rat to me. If you just change your a,b,c values nls arrives to some finite result (e.g. c=1.5 or c=0.3) . BTW by what magic you obtained such precise and wrong estimates for a,b,c? HTH Petr On 15 Aug 2006 at 5:54, Xiaodong Jin wrote: Date sent: Tue, 15 Aug 2006 05:54:51 -0700 (PDT) From: Xiaodong Jin [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Subject:[R] nls Is there anyway to change any y[i] value (i=2,...6) to make following NLS workable? x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) lm(1/y ~~ x) nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE) #0.0920573 : 1.16122 0.01565 1.0 #Error in numericDeriv(form[[3]], names(ind), env) : #Missing value or an infinity produced when evaluating the #model - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls
On Tue, 15 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: You problem is x^c for x = 0. If you intended only c 1, try a starting value meeting that condition (but it seems that the optimal c is about 0.27 is you increase x slightly). Surely you mean c 0. I did. nls(1/y ~ a+b*x^exp(c), start=list(a=1.16122,b=0.01565,c=0)) Nonlinear regression model model: 1/y ~ a + b * x^exp(c) data: parent.frame() a b c 0.9944025 0.1953168 -1.1495206 residual sum-of-squares: 0.03303464 nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=exp(-1.1))) Nonlinear regression model model: 1/y ~ a + b * x^c data: parent.frame() a b c 0.9944026 0.1953165 0.3167891 residual sum-of-squares: 0.03303464 (but even setting c=exp(-1) triggers the error; there could be cause for robustifying the nls algorithm) Well, there is an option to use a bounded-region algorithm, e.g. x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE, algorithm=port, lower=c(-Inf, -Inf, 0), upper=rep(Inf, 3)) works. Why have you used ~~ ? (Maybe because despite being asked not to, you sent HTML mail?) On Tue, 15 Aug 2006, Xiaodong Jin wrote: Is there anyway to change any y[i] value (i=2,...6) to make following NLS workable? x - c(0,5,10,15,20,25,30) y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000) lm(1/y ~~ x) nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE) #0.0920573 : 1.16122 0.01565 1.0 #Error in numericDeriv(form[[3]], names(ind), env) : #Missing value or an infinity produced when evaluating the model - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] question re: summarry.lm and NA values
Is there a way to get the following code to include NA values where the coefficients are NA? ((summary(reg))$coefficients) explanation: Using a loop, I am running regressions on several subsets of data1. reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) My regression has 10 independent variables, and I therefore expect 11 coefficients. After each regression, I wish to save the coefficients and standard errors of the coefficients in a table with 22 columns. I successfully extract the coefficients using the following code: reg$coefficients I attempt to extract the standard errors using : aperm((summary(reg))$coefficients)[2,] ((summary(reg))$coefficients) My problem: For some of my subsets, I am missing data for one or more of the independent variables. This of course causes the coefficients and standard erros for this variable to be NA. Is there a way to include the NA standard errors, so that I have the same number of standard erros and coefficients for each regression, and can then store the coefficients and standard erros in my table of 22 columns? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Presentation of multiple models in one table using xtable
Consider this situation: x1 - runif(100); x2 - runif(100); y - 2 + 3*x1 - 4*x2 + rnorm(100) m1 - summary(lm(y ~ x1)) m2 - summary(lm(y ~ x2)) m3 - summary(lm(y ~ x1 + x2)) Now you have estimated 3 different competing models, and suppose you want to present the set of models in one table. xtable(m1) is cool, but doing that thrice would give us 3 different tables. What I want is this one table: --- M1 M2 M3 --- Intercept 0.0816 3.6292 2.2272 (0.5533) (0.2316)***(0.2385)*** x12.81512.7606 (0.5533)*** (0.3193)*** x2 -4.2899-4.2580 (0.401)*** (0.3031)*** $\sigma_e$1.538 1.175 0.8873 $R^2$ 0.2089 0.5385 0.7393 --- How would one set about doing this? I am hoping that it's possible to write a function xtable.multi.lm where one would say xtable.multi.lm(m1,m2,m3) and get the above table. My sense is there are three challenges: 1. How to write a general R function which eats a unpredictable number of summary(lm()) objects, and fill out a matrix with results such as the above. 2. How to get a good xtable(), with decimal alignment and with the *** stuff (actually, $^{***}$). Will there be any catch in dropping into mathmode for $R^2$? After each pair of lines, I'd like to have a \\[2mm] so as to get a nice spacing in the table. 3. This style of presentation seems relevant to a whole host of models - whether ARCH models or survival models - not just OLS regressions. It would be very nice if one supported all manner of model objects and not just what comes out of lm(). I'm happy to take a crack at writing this function, which should ideally go back into the xtable library. But I don't have an idea on how to go about these two questions. If you will guide me, I am happy to work on it. :-) -- Ajay Shah http://www.mayin.org/ajayshah [EMAIL PROTECTED] http://ajayshahblog.blogspot.com *(:-? - wizard who doesn't know the answer. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] fMultivar OLS - how to do dynamic regression?
Hi folks! Does anybody know how to use the OLS function in fMultivar to do dynamic regression? I've tried specifying lags in OLS using a data series created in fSeries and it doesn't seem to work. I've done dynamic regression using dyn$lm and I was wondering how to accomplish the same thing using the OLS function from fMultivar. Thanks! John [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] question re: summarry.lm and NA values
Hi On 15 Aug 2006 at 7:01, r user wrote: Date sent: Tue, 15 Aug 2006 07:01:13 -0700 (PDT) From: r user [EMAIL PROTECTED] To: rhelp r-help@stat.math.ethz.ch Subject:[R] question re: summarry.lm and NA values Is there a way to get the following code to include NA values where the coefficients are NA ? ((summary(reg))$coefficients) better coef(reg) explanation: Using a loop, I am running regressions on several subsets of data1 . reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) My regression has 10 independent variables, and I therefore expect 11 coefficients. After each regression, I wish to save the coefficients and standard errors of the coefficients in a table with 22 columns. I successfully extract the coefficients using the following code: reg$coefficients I attempt to extract the standard errors using : aperm((summary(reg))$coefficients)[2,] ((summary(reg))$coefficients) My problem: For some of my subsets, I am missing data for one or more of the independent variables. This of course causes the coefficients and standard erros for this variable to be NA . ??%^*^?? What version? My lm behaves in accordance with na.action and it throws an error in case na.fail, computes a value in case of na.omit or na.exclude and again throws an error if the variable consist exclusively from NA values. The only way how to get NA in coeficient is when a variable is either constant or linear combination of other variable(s). Then coef(reg) will give you correctly NA in the variable which appears constant and in this case you could use it for setting standard error also as NA let say by using ifelse statement and matching of names. HTH Petr Is there a way to include the NA standard errors, so that I have the same number of standard erros and coefficients for each regression, and can then store the coefficients and standard erros in my table of 22 columns? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] question re: summarry.lm and NA values
Hi just as a quick workaround you probably can use aliased value from summary fff-rep(summary(reg)$aliased,4) dim(fff)-c(no.of.your.variables,4) fff[which(fff)]-NA fff[which(!fff)]-coef(summary(reg)) to get coef matrix with NA values HTH Petr On 15 Aug 2006 at 17:15, r user [EMAIL PROTECTED], wrote: From: Petr Pikal [EMAIL PROTECTED] To: r user [EMAIL PROTECTED], rhelp r-help@stat.math.ethz.ch Subject:Re: [R] question re: summarry.lm and NA values Date sent: Tue, 15 Aug 2006 17:15:01 +0200 Hi On 15 Aug 2006 at 7:01, r user wrote: Date sent:Tue, 15 Aug 2006 07:01:13 -0700 (PDT) From: r user [EMAIL PROTECTED] To: rhelp r-help@stat.math.ethz.ch Subject: [R] question re: summarry.lm and NA values Is there a way to get the following code to include NA values where the coefficients are NA ? ((summary(reg))$coefficients) better coef(reg) explanation: Using a loop, I am running regressions on several subsets of data1 . reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) My regression has 10 independent variables, and I therefore expect 11 coefficients. After each regression, I wish to save the coefficients and standard errors of the coefficients in a table with 22 columns. I successfully extract the coefficients using the following code: reg$coefficients I attempt to extract the standard errors using : aperm((summary(reg))$coefficients)[2,] ((summary(reg))$coefficients) My problem: For some of my subsets, I am missing data for one or more of the independent variables. This of course causes the coefficients and standard erros for this variable to be NA . ??%^*^?? What version? My lm behaves in accordance with na.action and it throws an error in case na.fail, computes a value in case of na.omit or na.exclude and again throws an error if the variable consist exclusively from NA values. The only way how to get NA in coeficient is when a variable is either constant or linear combination of other variable(s). Then coef(reg) will give you correctly NA in the variable which appears constant and in this case you could use it for setting standard error also as NA let say by using ifelse statement and matching of names. HTH Petr Is there a way to include the NA standard errors, so that I have the same number of standard erros and coefficients for each regression, and can then store the coefficients and standard erros in my table of 22 columns? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Looking for info on the Regression Modeling Strategies in R course in DC area
Hello list, A colleague of mine mentioned a great course on Regression Modeling Strategies in R. Anyone knows if this course is offered as public course in DC area? Thanks a bunch - TM - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] REML with random slopes and random intercepts giving strange results
Hi everyone, I have been using REML to derive intercepts and coeficients for each individual in a growth study. So the code is m2 - lmer(change.wt ~ newwt+(newwt|id), data = grow) Calling coef(model.lmer) gives a matrix with this information which is what I want. However, as a test I looked at each individual on its own and used a simple linear regression to obtain the same information, then I compared the results. It looks like the REML method doesnt seem to approximate the two parameters as well as using the simple linear regression on each individual separately, as judged by looking at graphs. Indeed, why do the results differ at all? Excuse my naivety if this is a silly question. Thanks to everyone for replying to my previous questions, very much appreciated. Simon Pickett PhD student Centre For Ecology and Conservation Tremough Campus University of Exeter in Cornwall TR109EZ Tel 01326371852 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] question re: summarry.lm and NA values
On Tue, 15 Aug 2006, Petr Pikal wrote: Hi On 15 Aug 2006 at 7:01, r user wrote: Date sent:Tue, 15 Aug 2006 07:01:13 -0700 (PDT) From: r user [EMAIL PROTECTED] To: rhelp r-help@stat.math.ethz.ch Subject: [R] question re: summarry.lm and NA values Is there a way to get the following code to include NA values where the coefficients are NA ? ((summary(reg))$coefficients) better coef(reg) coef(summary(reg)), perhaps. explanation: Using a loop, I am running regressions on several subsets of data1 . reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) My regression has 10 independent variables, and I therefore expect 11 coefficients. After each regression, I wish to save the coefficients and standard errors of the coefficients in a table with 22 columns. I successfully extract the coefficients using the following code: reg$coefficients I attempt to extract the standard errors using : aperm((summary(reg))$coefficients)[2,] ((summary(reg))$coefficients) My problem: For some of my subsets, I am missing data for one or more of the independent variables. This of course causes the coefficients and standard erros for this variable to be NA . ??%^*^?? What version? My lm behaves in accordance with na.action and it throws an error in case na.fail, computes a value in case of na.omit or na.exclude and again throws an error if the variable consist exclusively from NA values. The only way how to get NA in coeficient is when a variable is either constant or linear combination of other variable(s). Then coef(reg) will give you correctly NA in the variable which appears constant and in this case you could use it for setting standard error also as NA let say by using ifelse statement and matching of names. That happens in the print method, stats:::print.summary.lm contains coefs - x$coefficients if (!is.null(aliased - x$aliased) any(aliased)) { cn - names(aliased) coefs - matrix(NA, length(aliased), 4, dimnames = list(cn, colnames(coefs))) coefs[!aliased, ] - x$coefficients } so the code is already available Is there a way to include the NA standard errors, so that I have the same number of standard erros and coefficients for each regression, and can then store the coefficients and standard erros in my table of 22 columns? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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, [EMAIL PROTECTED] 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] question re: summarry.lm and NA values
Is there a way to... always has the answer yes in R (or C or any language for that matter). The question is: Is there a GOOD way...? where good depends on the specifics of the situation. So after that polemic, below is an effort to answer, (adding to what Petr Pikal already said): -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of r user Sent: Tuesday, August 15, 2006 7:01 AM To: rhelp Subject: [R] question re: summarry.lm and NA values Is there a way to get the following code to include NA values where the coefficients are NA? ((summary(reg))$coefficients) BAAAD! Don't so this. Use the extractor on the object: coef(reg) This suggests that you haven't read the documentation carefully, which tends to arouse the ire of would-be helpers. explanation: Using a loop, I am running regressions on several subsets of data1. reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) ??? There's an error here I think. Do you mean update()? Do you have your subscripting correct? My regression has 10 independent variables, and I therefore expect 11 coefficients. After each regression, I wish to save the coefficients and standard errors of the coefficients in a table with 22 columns. I successfully extract the coefficients using the following code: reg$coefficients Use the extractor, coef() I attempt to extract the standard errors using : aperm((summary(reg))$coefficients)[2,] BAAAD! Use the extractor vcov(): sqrt(diag(vcov(reg))) ((summary(reg))$coefficients) My problem: For some of my subsets, I am missing data for one or more of the independent variables. This of course causes the coefficients and standard erros for this variable to be NA. Not it doesn't, as Petr said. One possible approach: Assuming that a variable is actually missing (all NA's), note that coef(reg) is a named vector, so that the character string names of the regressors actually used are available. You can thus check for what's missing and add them as NA's at each return. Though I confess that I see no reason to put things ina matrix rather than just using a list. But that's a matter of personal taste I suppose. Is there a way to include the NA standard errors, so that I have the same number of standard erros and coefficients for each regression, and can then store the coefficients and standard erros in my table of 22 columns? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] REML with random slopes and random intercepts giving strange results
I don't this is because you are using REML. The BLUPs from a mixed model experience some shrinkage whereas the OLS estimates would not. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Simon Pickett Sent: Tuesday, August 15, 2006 11:34 AM To: r-help@stat.math.ethz.ch Subject: [R] REML with random slopes and random intercepts giving strange results Hi everyone, I have been using REML to derive intercepts and coeficients for each individual in a growth study. So the code is m2 - lmer(change.wt ~ newwt+(newwt|id), data = grow) Calling coef(model.lmer) gives a matrix with this information which is what I want. However, as a test I looked at each individual on its own and used a simple linear regression to obtain the same information, then I compared the results. It looks like the REML method doesnt seem to approximate the two parameters as well as using the simple linear regression on each individual separately, as judged by looking at graphs. Indeed, why do the results differ at all? Excuse my naivety if this is a silly question. Thanks to everyone for replying to my previous questions, very much appreciated. Simon Pickett PhD student Centre For Ecology and Conservation Tremough Campus University of Exeter in Cornwall TR109EZ Tel 01326371852 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] nls convergence problem
I'm having problems getting nls to agree that convergence has occurred in a toy problem. nls.out never gets defined when there is an error in nls. Reaching the maximum number of iterations is alway an error, so nls.out never gets defined when the maximum number of iterations is reched. From ?nls.control: tol: A positive numeric value specifying the tolerance level for the relative offset convergence criterion. From some S-Plus documentation I found online via Google: http://www.uni-muenster.de/ZIV/Mitarbeiter/BennoSueselbeck/s-html/helpfiles/nls.control.html tolerance: tolerance for the relative offset convergence criterion in the algorithm. Default 0.001. Note that the convergence test used in nls() is strictly relative. Therefore if the solution to a problem turned out to be a perfect fit (unlikely except in artificial examples), convergence is not guaranteed to be recognized by the algorithm. Here's my toy problem: ?nls.control ?nls # Method 2 X - 0:15 Y - 9.452 * exp(-0.109*X) + 5.111 # Toy problem nls.out - nls(Y ~ a*exp(b*X)+c, +start=list(a=6,b=-0.5,c=1), +control=nls.control(maxiter=15, tol=0.01), # nothing makes sense here +trace=TRUE) 1016.507 : 6.0 -0.5 1.0 143.5807 : 6.1680290 -0.1506021 4.4013020 7.306365 : 9.10093164 -0.09114858 5.44620298 0.0342703 : 9.2801070 -0.1109063 5.2795803 3.001985e-05 : 9.4506654 -0.1089749 5.1122982 1.918531e-14 : 9.452 -0.109 5.111 6.894644e-28 : 9.452 -0.109 5.111 2.208811e-29 : 9.452 -0.109 5.111 7.888609e-30 : 9.452 -0.109 5.111 7.888609e-30 : 9.452 -0.109 5.111 7.099748e-30 : 9.452 -0.109 5.111 3.155444e-30 : 9.452 -0.109 5.111 3.155444e-30 : 9.452 -0.109 5.111 3.155444e-30 : 9.452 -0.109 5.111 3.155444e-30 : 9.452 -0.109 5.111 3.155444e-30 : 9.452 -0.109 5.111 Error in nls(Y ~ a * exp(b * X) + c, start = list(a = 6, b = -0.5, c = 1), : number of iterations exceeded maximum of 15 There is near-perfect convergence after 12 iterations, but I cannot figure out a way for R to recognize it. What does relative offset convergence criterion mean? How can nls.control be used to say this problem has converged, and have nls exit without an error? Should the R documentation be modified to explain what relative offset convergence criterion means? Should the R documentation be expanded to include the comment from the S-Plus: Therefore if the solution to a problem turned out to be a perfect fit (unlikely except in artificial examples), convergence is not guaranteed to be recognized by the algorithm. If this is true, this seems like a suboptimal design. Thanks for any insight about this. efg Earl F. Glynn Scientific Programmer Stowers Institute for Medical Research __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to show classes of all columns of a data frame?
Hi all, Suppose I have a data frame myDF, col A is factor, col B is numeric, col C is character. I can get their classes by class(myDF$A) but is there a quick way to show what classes of all columns are? Thank you. Tian [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to show classes of all columns of a data frame?
On Tue, 2006-08-15 at 13:10 -0400, T Mu wrote: Hi all, Suppose I have a data frame myDF, col A is factor, col B is numeric, col C is character. I can get their classes by class(myDF$A) but is there a quick way to show what classes of all columns are? Thank you. Tian Depending upon the output format you desire: lapply(iris, class) $Sepal.Length [1] numeric $Sepal.Width [1] numeric $Petal.Length [1] numeric $Petal.Width [1] numeric $Species [1] factor or sapply(iris, class) Sepal.Length Sepal.Width Petal.Length Petal.Width Species numericnumericnumericnumeric factor See ?lapply and ?sapply HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] zlim not working in persp3d
Hi list, I'm trying to limit the z axis plotted using persp3d {rgl} to a given range. I'm trying the following statement persp3d(g0,pa,D, xlim=range(g0),ylim=range(pa), zlim=range (zbounds),col=lightblue) where g0 and pa are both a range of 20 numbers from 0.05 to 1, and D is a 20 x 20 matrix with values ranging from 0 to 3. I'm trying to set a cap on the Z displayed in the surface to 1000 so I can see the lower levels in greater detail, but the persp3d statement doesn't appear to respond to any zlim parameter that I have tried. Can anyone suggest a method? I'm using R Gui for OSX ver 1.16, R version 2.3.1, and rgl is version 0.67-2 Thanks Ken [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] rexp question
I am using rexp to generate several exponential distributions. I am passing rexp a vector of rates , r. I am wanting to simulate a sample of size 200 for each rate so the code looks like: rexp(n=200*length(r),rate=r) this gives me a vector of the random exponential variables, but they are all disjointed b/c rexp goes through and simulates an exponential variable for each rate and it does that 200 times, is there any way to get it to do it the oppisite way,i.e., 200 times in a row for each rate. I have already tried using a for loop but it is slow. thanks, Spencer [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Grasper model error
I tried this over a the grasp users yahoo group and got no responseSo I wonder if anyone here knows about grasper I keep getting this error when trying to run a model. Error in smooth.construct.tp.smooth.spec(object, data, knots) : Too many knots for t.p.r.s term: see `gam.control' to increase limit, or use a different basis, or see large data set help for `gam'. I'm using R for OSX Gui Version 1.16 R-version 2.3.1, running grasper version 0.4-4 My dataset has ~7500 rows of input, and 21 cols of input. Does anyone have an idea how to troubleshoot this? Thanks Ken [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] rexp question
Try: rexp(n=200*length(r), rate=rep(r, each=200) ) Hope this helps, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Spencer Jones Sent: Tuesday, August 15, 2006 12:15 PM To: R-Help Subject: [R] rexp question I am using rexp to generate several exponential distributions. I am passing rexp a vector of rates , r. I am wanting to simulate a sample of size 200 for each rate so the code looks like: rexp(n=200*length(r),rate=r) this gives me a vector of the random exponential variables, but they are all disjointed b/c rexp goes through and simulates an exponential variable for each rate and it does that 200 times, is there any way to get it to do it the oppisite way,i.e., 200 times in a row for each rate. I have already tried using a for loop but it is slow. thanks, Spencer [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls convergence problem
Earl F. Glynn efg at stowers-institute.org writes: Here's my toy problem: ?nls.control ?nls # Method 2 X - 0:15 Y - 9.452 * exp(-0.109*X) + 5.111 # Toy problem nls.out - nls(Y ~ a*exp(b*X)+c, +start=list(a=6,b=-0.5,c=1), +control=nls.control(maxiter=15, tol=0.01), # nothing makes sense here +trace=TRUE) This toy problem is exactly what the warning is for: Warning Do not use nls on artificial zero-residual data. Add some noise and try again. Dieter Menne __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Hierarchical clustering
Hi, I'm new to R; this is my second email to this forum. I have a dataset that looks like this: GeneA B C D A 511 6116 15151515 B 164 115 1515494 C 21646 D 1616 E 1461 F G H There are about 4 genes with 3-10 patients per dataset. Gene expression values fill the matrix. I want to do a simple hierarchical clustering along BOTH rows and columns and get a tree dendrogram output for both in one panel. Is there a simple code for this, leaving out fancy options? Thank you very much. -DS. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] model = F causing error in polr()
Hi all, I got an error message if I set model =F in polr(), like polr(y ~ x1 + x2, data1, model = F, method = probit) Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : variable lengths differ (found for '(model)') but polr(y ~ x1 + x2, data1, method = probit) would work. Why? Thank you, Tian [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] model = F causing error in polr()
Try ``model = FALSE'' rather than ``model = F'' and see if it makes a difference. You make have an unwanted variable named ``F'' lurking somewhere. (In general it is a *bad* idea to use ``F'' when ``FALSE'' is intended.) cheers, Rolf Turner [EMAIL PROTECTED] ===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+=== Original message: Hi all, I got an error message if I set model =F in polr(), like polr(y ~ x1 + x2, data1, model = F, method = probit) Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : variable lengths differ (found for '(model)') but polr(y ~ x1 + x2, data1, method = probit) would work. Why? Thank you, Tian __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls convergence problem
Dieter Menne [EMAIL PROTECTED] wrote in message news:[EMAIL PROTECTED] Earl F. Glynn efg at stowers-institute.org writes: This toy problem is exactly what the warning is for: Warning Do not use nls on artificial zero-residual data. Add some noise and try again. Thank you! I had adapted some code and must confess I had read ?nls.control thoroughly, but not ?nls. I had even used debug on nls, traced it through line by line to the .Call statement, trying to figure out why nls.out never got defined. The source code has no comments at all. IMHO, the warning should be in the Description at the top of the ?nls page, not at the bottom of the page. The warning should also appear on the ?nls.control page. But, a better way would be to have a software design that eliminated the warning. It's not clear to me why this problem cannot be fixed somehow. You shouldn't need to add noise to a problem to solve it. (It's a bit like saying addition works, but not for integers without adding some noise.) If there can be arbitrary defaults of maxiter=50, and (relative) tol=1e-5 in nls.control, there could be another arbitrary (absolute) convergence criterion. Or, maybe there's something I don't understand about the algorithm being used. Just my $0.02 and minority opinion, efg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Hierarchical clustering
Hi, I'm new to R; this is my second email to this forum. I have a dataset that looks like this: Gene AB 12536.25 2532.2 22527.35 4583.3 There are about 4 genes with 3-10 patients per dataset. Gene expression values fill the matrix. I want to do a simple hierarchical clustering along BOTH rows and columns and get a tree dendrogram output for both in one panel. Is there a simple code for this, leaving out fancy options? If a heat map can be generated, that'll be even better. Thank you very much. -DS. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls convergence problem
Or, maybe there's something I don't understand about the algorithm being used. Indeed! So before making such comments, why don't you try to learn about it? Doug Bates is a pretty smart guy, and I think you do him a disservice when you assume that he somehow overlooked something that he explicitly warned you about. I am fairly confident that if he could have made the problem go away, he would have. So I think your vent was a bit inconsiderate and perhaps even intemperate. The R Core folks have produced a minor miracle IMO, and we should all be careful before assuming that they have overlooked easily fixable problems. They're certainly not infallible -- but they're a lot less fallible than most of the rest of us when it comes to R. Just my $0.02 and minority opinion, efg ... and mine. -- Bert Gunter __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] A model for possibly periodic data with varying amplitude [repost, much edited]
Hi dear R community, I have up to 12 measures of a protein for each of 6 patients, taken every two or three days. The pattern of the protein looks periodic, but the height of the peaks is highly variable. I'm testing for periodicity using a Monte Carlo simulation envelope approach applied to a cumulative periodogram. Now I want to predict the location of the peaks in time. Of course, the peaks might be occurring on unmeasured days. Sadly, an NDA prohibits me from sharing the real data. The data look something like this: ## patient - data.frame( day = c(1, 3, 5, 8, 10, 12, 15, 17, 19, 22, 24, 26), protein = c(5, 3, 10, 7, 2, 8, 25, 22, 7, 10, 12, 5) ) plot(patient$day, patient$protein, type=b) # This is my model: wave.form - deriv3( ~ sin(2*pi*((day-offset)/period + 0.25)) * amplitude + mean, c(period, offset, amplitude, mean), function(day, period, offset, amplitude, mean){}) curve(wave.form(x, period=7, offset=2, mean=5, amplitude=4), from=1, to=30) require(nlme) wave.1 - gnls(protein ~ wave.form(day, period, offset, amplitude, mean), start=list(period=7, offset=0, amplitude=10, mean=5), weights=varPower(), data=patient) wave.1 # I emphasize that I don't care about the mean or amplitude, just the offset and the period. The problem for my model is that spikes in the data, such as at day 15, seem to make fitting the model quite unstable, (although it is fine in this artificial case). The best I can think of right now is to use the weights=varPower() model in gnls(), which I expect will down-weight the extreme spikes. Often that model fails to converge. Then, all I can do is fall back on modelling the log response. If I do both, then when the weights model converges the estimates sometimes differ, and sometimes considerably. I know that I need more data, but more data are not available. So, instead I need to do the best I can :) I have a half-baked idea of conditionally shrinking the peaks, so that they all look like the same size - perhaps by discretizing the data into, say, three possible values defined by quantiles of the response variable or the estimated slope - but I'm concerned that doing so risks creating the periodicity that I want to detect. Also the series is not necessarily stationary, so I wouild have to smooth it somehow first. On the other hand, maybe I can test for periodicity in the original data, and then fit on some transformed data. If I reject the null hypothesis of no periodic oscillation, then is it reasonable to condition subsequent analysis on the belief that periodicity exists, and try to determine its characteristics? I wonder if anyone can suggest an analysis technique or model that more directly accommodates the varying heights of the peaks? Or is this pretty much the best I can do? Any suggestions will be gratefully received. Cheers, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls convergence problem
Berton Gunter [EMAIL PROTECTED] wrote in message news:[EMAIL PROTECTED] Or, maybe there's something I don't understand about the algorithm being used. Indeed! So before making such comments, why don't you try to learn about it? Doug Bates is a pretty smart guy, and I think you do him a disservice when you assume that he somehow overlooked something that he explicitly warned you about. I am fairly confident that if he could have made the problem go away, he would have. So I think your vent was a bit inconsiderate and perhaps even intemperate. The R Core folks have produced a minor miracle IMO, and we should all be careful before assuming that they have overlooked easily fixable problems. They're certainly not infallible -- but they're a lot less fallible than most of the rest of us when it comes to R. I meant no disrespect to Doug Bates or any of the R Core folks. I thought what I wrote had a neutral tone and was respectful. I am sorry if anyone was offended by my comments and suggestions. I am certainly thankful for all the hard work that has gone into developing R. efg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] coefficients' order in polr()?
Hi all, I am using polr(). The resulting coefficients of first levels are always 0. What to do if I wnat to get the coefficients of the last level 0. For example, suppose x has 3 levels, 1, 2, 3 probit - plor(y ~ x, data1, method='probit') will get coefficients of level 2, 3 of x, but I want coefficients of level 1, 2 Thank you, Tian [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Looking for info on the Regression Modeling Strategies in R course in DC area
Tim McDonald wrote: Hello list, A colleague of mine mentioned a great course on Regression Modeling Strategies in R. Anyone knows if this course is offered as public course in DC area? Thanks a bunch - TM I'll be teaching this course in DC Sept 28-29 for XL Solutions. Information about the course may be obtained from biostat.mc.vanderbilt.edu/rms Frank Harrell - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] coefficients' order in polr()?
you can use _relevel_ to Reorder Levels of Factor. use ?relevel to get more information. 2006/8/16, T Mu [EMAIL PROTECTED]: Hi all, I am using polr(). The resulting coefficients of first levels are always 0. What to do if I wnat to get the coefficients of the last level 0. For example, suppose x has 3 levels, 1, 2, 3 probit - plor(y ~ x, data1, method='probit') will get coefficients of level 2, 3 of x, but I want coefficients of level 1, 2 Thank you, Tian [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- 黄荣贵 Department of Sociology Fudan University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Specifying Path Model in SEM for CFA
I'm using specify.model for the sem package. I can't figure out how to represent the residual errors for the observed variables for a CFA model. (Once I get this working I need to add some further constraints.) Here is what I've tried: model.sa - specify.model() F1 - X1,l11, NA F1 - X2,l21, NA F1 - X3,l31, NA F1 - X4,l41, NA F1 - X5, NA, 0.20 F2 - X1,l12, NA F2 - X2,l22, NA F2 - X3,l32, NA F2 - X4,l42, NA F2 - X6, NA, 0.25 F1- F2,g12, 1 F1- F1,g11, 1 F2- F2,g22, 1 X1- X1, NA, 1 X2- X2, NA, 1 X3- X3, NA, 1 X4- X4, NA, 1 X5- X5, NA, 1 X6- X6, NA, 1 This at least converges: summary(fit.sem) Model Chisquare = 2147 Df = 10 Pr(Chisq) = 0 Chisquare (null model) = 2934 Df = 15 Goodness-of-fit index = 0.4822 Adjusted goodness-of-fit index = -0.087387 RMSEA index = 0.66107 90 % CI: (NA, NA) Bentler-Bonnett NFI = 0.26823 Tucker-Lewis NNFI = -0.098156 Bentler CFI = 0.26790 BIC = 2085.1 Normalized Residuals Min. 1st Qu. MedianMean 3rd Qu.Max. -5.990 -0.618 0.192 0.165 1.700 3.950 Parameter Estimates Estimate Std Error z value Pr(|z|) l11 -0.245981 0.21863 -1.12510 0.26054748 X1 --- F1 l21 -0.308249 0.22573 -1.36555 0.17207875 X2 --- F1 l31 0.202590 0.079102.56118 0.01043175 X3 --- F1 l41 -0.235156 0.21980 -1.06985 0.28468885 X4 --- F1 l12 0.839985 0.219623.82476 0.00013090 X1 --- F2 l22 0.828460 0.225483.67418 0.00023862 X2 --- F2 l32 0.066722 0.083690.79725 0.42530606 X3 --- F2 l42 0.832037 0.218403.80963 0.00013917 X4 --- F2 g12 0.936719 0.643311.45609 0.14536647 F2 -- F1 g11 2.567669 1.256082.04418 0.04093528 F1 -- F1 g22 1.208497 0.550402.19567 0.02811527 F2 -- F2 Iterations = 59 And it produces the following path diagram: path.diagram(fit.sem) digraph fit.sem { rankdir=LR; size=8,8; node [fontname=Helvetica fontsize=14 shape=box]; edge [fontname=Helvetica fontsize=10]; center=1; F2 [shape=ellipse] F1 [shape=ellipse] F1 - X1 [label=l11]; F1 - X2 [label=l21]; F1 - X3 [label=l31]; F1 - X4 [label=l41]; F1 - X5 [label=]; F2 - X1 [label=l12]; F2 - X2 [label=l22]; F2 - X3 [label=l32]; F2 - X4 [label=l42]; F2 - X6 [label=]; } But I don't see the residual error terms that go into each of the observed variables X1 - X6. I've tried: model.sa - specify.model() E1 - X1, e1, 1 E2 - X2, e2, 1 E3 - X3, e3, 1 E4 - X4, e4, 1 E5 - X5, e5, 1 E6 - X6, e6, 1 E1- E1, s1, NA E2- E2, s2, NA E3- E3, s3, NA E4- E4, s4, NA E5- E5, s5, NA E6- E6, s6, NA F1 - X1,l11, NA F1 - X2,l21, NA F1 - X3,l31, NA F1 - X4,l41, NA F1 - X5, NA, 1 F2 - X1,l12, NA F2 - X2,l22, NA F2 - X3,l32, NA F2 - X4,l42, NA F2 - X6, NA, 1 F1- F2, NA, 1 F1- F1, NA, 1 F2- F2,g22, NA X1- X1, NA, 1 X2- X2, NA, 1 X3- X3, NA, 1 X4- X4, NA, 1 X5- X5, NA, 1 X6- X6, NA, 1 I'm trying to use E1 - E6 as the residual error terms. But I get warning messages about no variances for X1-X6 and it doesn't converge. Also, the associated path diagram: digraph fit.sem { rankdir=LR; size=8,8; node [fontname=Helvetica fontsize=14 shape=box]; edge [fontname=Helvetica fontsize=10]; center=1; E1 [shape=ellipse] E2 [shape=ellipse] E3 [shape=ellipse] E4 [shape=ellipse] E5 [shape=ellipse] E6 [shape=ellipse] F2 [shape=ellipse] F1 [shape=ellipse] E1 - X1 [label=]; E2 - X2 [label=]; E3 - X3 [label=]; E4 - X4 [label=]; E5 - X5 [label=]; E6 - X6 [label=]; F1 - X1 [label=l11]; F1 - X2 [label=l21]; F1 - X3 [label=l31]; F1 - X4 [label=l41]; F1 - X5 [label=]; F2 - X1 [label=l12]; F2 - X2 [label=l22]; F2 - X3 [label=l32]; F2 - X4 [label=l42]; F2 - X6 [label=]; } Has ellipses around the E1-E6 which I believe indicates they are latent factors and not residual errors. If anyone could point in the right direction I would appreciate it. Rick B. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.