Perhaps you would be interested in the nlsList() function in the nlme
package.
HTH,
Andy
> From: Jose A. Hernandez
>
> R-community,
>
> Assuming a file with 3 columns: site, nrate, yield. In SAS I
> can easily
> run the analysis BY SITE using the following code:
>
> *--------------------------------------------------*;
> proc nlin method=MARQUARDT;
> by farm;
> parms b0=100 b1=0.33 b2=-0.00111 x0=120;
> model yield = (b0 + b1 * nrate + b2 * nrate * nrate) * (nrate <= x0) +
> (b0 + b1 * x0 + b2 * x0 * x0) * (nrate>x0);
> output out=out predicted=pred ess=rss parms=b0 b1 b2 x0;
> *--------------------------------------------------*;
>
> In R I tried using nls in a loop, but it appears the analysis is very
> sensitive to the initial values and I keep getting the "singular
> gradient matrix at initial parameter estimates" error.
>
> Any ideas on how to work around this issue would be greatly
> appreciated.
>
> Thanks,
>
> # Example 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)
> bar_1 <- NULL
> for (i in 1:length(levels(site))) {
> j <- site == levels(site)[i]
> x <- nrate[j]
> y <- yield[j]
> nls.st <- list(b0=140, b1=0.5, b2=-0.002, x0=125)
> out<- try(nls(y ~ (b0 + b1*x + b2*I(x^2))*(x <= x0)
> + (b0 + b1*x0 + b2*I(x0^2))*(x > x0),
> start = nls.st))
> fit1 <- summary(out)$parameters
> qp.resi <- summary(out)$residuals
> rss <- sum(qp.resi^2)
> bar_1 <- rbind(bar_1, c(i, fit1[1, 1:2], fit1[2, 1:2], fit1[3,
> 1:2],
> fit1[4, 1:2], rss))
> }
> dimnames(bar_1) <- list(NULL, c("Site", "b0", "s.e.", "b1", "s.e.",
> "b2", "s.e.", "x0", "s.e.", "RSS"))
> print(bar_1)
>
> ______________________________________________
> [EMAIL PROTECTED] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
>
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
