[R] Quadratic curve on a loop
Liaw, Andy
andy_liaw at merck.com
Tue Nov 23 22:40:09 CET 2004
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
>
>
More information about the R-help
mailing list