[R] nls() and lines()

Mon Jul 18 10:55:18 CEST 2011

On 2011-07-17 17:37, Steven Ranney wrote:
> All -
>
> I'm having an issue with trying to plot a model derived from nls()
> onto a simple plot.  I have included a sample data set and the code
> that I've been using.
>
>     year month day       date location mileage  cost gallon      cpg
>    mpg          x
> 2009     1   4   1/4/2009      BZN  124585 19.39  14.37 1.349339
> 10.71677 2009-01-04
> 2009     1  15  1/15/2009      BZN  124888  23.2  16.12 1.439206
> 18.79653 2009-01-15
> 2009     1  27  1/27/2009      BZN  125133 21.51  14.35 1.498955
> 17.07317 2009-01-27
> 2009     2  16  2/16/2009      BZN  125429 27.96  15.54 1.799228
> 19.04762 2009-02-16
> 2009     2  27  2/27/2009      BZN  125702 26.82  14.27 1.879467
> 19.13104 2009-02-27
> 2009     3  19  3/19/2009      BZN  125941 24.38  12.91 1.888459
> 18.51278 2009-03-19
> 2009     4   9   4/9/2009      BZN  126260 32.59  16.30 1.999387
> 19.57055 2009-04-09
> 2009     4  28  4/28/2009      BZN  126587 34.58  16.79 2.059559
> 19.47588 2009-04-28
> 2009     5  17  5/17/2009      BZN  126925 35.78  16.57 2.159324
> 20.39831 2009-05-17
> 2009     5  27  5/27/2009      BZN  127240 35.57  15.01 2.369753
> 20.98601 2009-05-27
> 2009     6   7   6/7/2009      BZN  127590 40.99  16.60 2.469277
> 21.08434 2009-06-07
> 2009     6  21  6/21/2009      BZN  127910 41.52  15.64 2.654731
> 20.46036 2009-06-21
> 2009     7  21  7/21/2009      BZN  128264 43.37  16.67 2.601680
> 21.23575 2009-07-21
> 2009     8  11  8/11/2009      BZN  128618 42.68  16.42 2.599269
> 21.55907 2009-08-11
> 2009     8  27  8/27/2009      BZN  128947 43.12  16.60 2.597590
> 19.81928 2009-08-27
> 2009     9  21  9/21/2009      BZN  129255 40.44  15.56 2.598972
> 19.79434 2009-09-21
> 2009    10   1  10/1/2009      BZN  129541 38.55  14.83 2.599461
> 19.28523 2009-10-01
> 2009    10  11 10/11/2009      BZN  129806 36.65  14.10 2.599291
> 18.79433 2009-10-11
> 2009    10  22 10/22/2009      BZN  130027 30.18  11.61 2.599483
> 19.03531 2009-10-22
> 2009    11   9  11/9/2009      BZN  130358 43.19  16.62 2.598676
> 19.91576 2009-11-09
> 2009    11  22 11/22/2009      BZN  130631 39.23  15.09 2.599735
> 18.09145 2009-11-22
> 2009    12   5  12/5/2009      BZN  130950 44.43  17.09 2.599766
> 18.66589 2009-12-05
> 2009    12  30 12/30/2009      BZN  131239 42.14  16.70 2.523353
> 17.30539 2009-12-30
>
> After converting my dates into R-usable dates:
>
> #convert my dates to R-usable dates
> x<- strptime(date, format="%m/%d/%Y")
> x
> mileage<- cbind(mileage, x)
>
> I plot the data and model mpg as a function of date.  In the nls()
> statement, I convert x back to a numeric value so that I can conduct
> the regression:
>
> plot(mpg~x, data=mileage[year==2009,], ylab="Miles per gallon",
> xlab="2009", yaxs="i", ylim=c(10,30))
> nls.2009<- nls(mpg~(alpha*(as.numeric(x)^2))+(bravo*as.numeric(x))+(charlie),
> data=mileage[year==2009,], start=list(alpha=-2e-14, bravo=5e-5,
> charlie=-31407),
>    trace=T, na.action=na.omit, nls.control(minFactor=0.000000000000000000001))
> plot(mpg~x, data=mileage[year==2009,])
>    modb=seq(min(as.numeric(x)), max(as.numeric(x)), by=10000)
>    lines(modb, predict(nls.2009, lines(as.numeric(x)=modb)))
>
> Unfortunately, when I run the final line of this code, I get the following:
>
> Error: unexpected '=' in "  lines(modb, predict(nls.2009, lines(as.numeric(x)="
>
> In other similar analyses, I've been able to plot an nls() model using
> this exact code--altered of course according to information--but here
> I'm at a loss.  I'm certain it has something to do with the
> lines(...as.numeric(x)) value I'm trying to plot, but I can't figure
> out what I'm doing wrong.

That last line of code doesn't look right to me. The arguments
that you need to supply to predict() are 'object' and 'newdata',
where 'newdata' must have the appropriate form. Unless you have
your own function lines(), I don't think that lines(as.numeric(x)=modb)
would qualify as newdata.

It's usually a bad idea to shove too much stuff into a single command
and a good idea to use str() often.

This 'exact' code worked in the past?

Peter Ehlers

>
> The model is fine, but it's the plotting of the model that escapes me.
>
> I'm running R version 2.12.1 on a Windows 7 machine.
>
> Thanks for your help -
>
> Steven H. Ranney
>
> http://stevenranney.blogspost.com
> http://www.steven-ranney.com
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.