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)