[R] Help me to migrate from SAS

[Liaw, Andy]

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)
> 
> ______________________________________________
> 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
> 
>