[R] Quadratic curve on a loop

[Liaw, Andy]

AFAICS, there are two errors:

1. curve() is probably not what you want: it plots functions, and that's not
what predict() gives you.

2. predict() needs to be given the correct data frame: The name of the
variable(s) need to match those used in the formula.  The formula has `x' on
the righthand side, so the data frame given to predict() needs to have a
variable named `x', otherwise `x' will be looked up elsewhere, and if it
finds one, it's unlikely to be the one you intended.

Try something like:

lines(predict(out, newdata=data.frame(x=seq(min(x), max(x), length=51)))

HTH,
Andy

> From: Jose A. Hernandez
> 
> Dear R community, once again I request your generous help on 
> an R issue 
> I cannot seem to figure out.
> 
> I am trying to do some analysis for multiple sites, in my 
> example I am 
> using only two, but the same will done on many sites. It's just ONE 
> line, so I hope somebody can give a me a hand by email.
> 
> Look at the loop to make the plot x ~ y. If the curve function is 
> removed from the loop everything works fine, but when I tried 
> to add the 
> curve to the plot, I get the following error:
> 
> Error in xy.coords(x, y) : x and y lengths differ
> In addition: Warning message:
> 'newdata' had 101 rows but variable(s) found have 6 rows
> 
> I know this error should be a hint about how to proceed but I 
> am stuck 
> figuring ways to make the curve to print in the plots.
> 
> Any hints on how to solve this issue would be greatly appreciated.
> 
> Sincerely,
> 
> Jose
> 
> PS.
>  > version
>           _
> platform i386-pc-mingw32
> arch     i386
> os       mingw32
> system   i386, mingw32
> status
> major    2
> minor    0.1
> year     2004
> month    11
> day      15
> language R
> 
> 
> 
> ###################################################################
> # some data
> site <- c(1,1,1,1,1,1,2,2,2,2,2,2)
> nrate <- c(0,60,90,120,150,180,0,60,90,120,150,180)
> yield <- 
> c(161.7,187.1,188.5,196.6,196.0,196.5,160.7,186.1,189.5,194.6,
> 195.0,198.5)
> site <- as.factor(site)
> 
> # file name
> fname <- levels(site)
> i <- nchar(fname) == 1
> fname[i] <- paste("0", fname[i], sep = "")
> fname <- paste("site", fname, ".eps", sep = "")
> 
> ###################################################################
> ###################################################################
> 
> for (i in 1:length(levels(site))) {
>      j <- site == levels(site)[i]
>      x <- nrate[j]
>      y <- yield[j]
>      out <- try(lm(y ~ x + I(x^2)), silent = TRUE)
>      plot(y ~ x,
>           pch = 16,
>           xlab = expression(paste("Nitrogen rate [lbs ac" ^-1,"]")),
>           ylab = expression(paste("Corn yield [bu ac"^-1,"]")))
>           curve(predict(out, data.frame(nrate = x)), add = T)
>           title(paste("Site", i))
>           # dev.copy2eps(file = fname[i], horizontal = TRUE)
> }
> ###################################################################
> ###################################################################
> 
> 
> 
> -- 
> Jose A. Hernandez
> Department of Soil, Water, and Climate
> University of Minnesota
> 1991 Upper Buford Circle
> St. Paul, MN 55108
> 
> Ph. (612) 625-0445, Fax. (612) 625-2208
> 
> ______________________________________________
> R-help at 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
> 
>