[R] Help me to migrate from SAS

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