Following on my posting of this morning, concerning a
problem that I am having constructing a self-starting function for use with nls
(and eventually with nlsList and nlme), the following is the self-starting
function called NRhyperbola:
> NRhyperbola
function (Irr,theta,Am,alpha,Rd)
{
# Am is the maximum gross photosynthetic rate
# Rd is the dark resiration rate (positive value)
# theta is the shape parameter
# alpha is the initial quantum yeild
#
(1/(2*theta))*
(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2
-4*alpha*theta*Am*Irr))- Rd
}
attr(,"initial")
function (mCall, LHS, data)
{
xy <- sortedXyData(mCall[["Irr"]],
LHS, data)
if (nrow(xy) < 3)
stop("Too
few unique irradiance values")
fit <- smooth.spline(xy[,
"x"], xy[, "y"])
alpha <- predict(fit, x = 0, deriv
= 1)$y
if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)-0)
if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0
# to correct for minimal values greater that zero
Rd <- (predict(fit, 0)$y - delta.for.zero)
# to correct for maximal values
alpha <- predict(fit, x = max(xy[,"x"],na.rm=T),
deriv = 1)$y
if(alpha>0){
Am <- predict(fit, max(xy[,"x"],na.rm=T))$y+alpha*(3000-max(xy[,"x"],na.rm=T))
Am <- Am - Rd
}
theta1 <- rep(NA, 4)
light <- c(100, 200, 500, 800)
for (i in 1:4) {
y.pred
<- predict(fit, light[i])$y
theta1[i]
<- ((y.pred + Rd) * (alpha * light[i] + Am) -
alpha * Am * light[i])/(y.pred + Rd)^2
}
theta <- mean(theta1)
Rd<- -1*Rd
value <- c(theta, Am, alpha,
Rd)
names(value) <- mCall[c("theta",
"Am", "alpha", "Rd")]
value
}
attr(,"class")
[1] "selfStart"
This function was created, following pp. 345-346 of Pinheiro
& Bates (Mixed-effects models in S and S-PLUS) by first creating NRhyperbola:
NRhyperbola<- function (Irr,theta,Am,alpha,Rd)
{
# Am is the maximum gross photosynthetic rate
# Rd is the dark resiration rate (positive value)
# theta is the shape parameter
# alpha is the initial quantum yeild
#
(1/(2*theta))*
(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2
-4*alpha*theta*Am*Irr))- Rd
}
then creating a function
NRhyperbolaInit<-
function (mCall, LHS, data)
{
xy <- sortedXyData(mCall[["Irr"]],
LHS, data)
if (nrow(xy) < 3)
stop("Too
few unique irradiance values")
fit <- smooth.spline(xy[,
"x"], xy[, "y"])
alpha <- predict(fit, x = 0, deriv
= 1)$y
if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)-0)
if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0
# to correct for minimal values greater that zero
Rd <- (predict(fit, 0)$y - delta.for.zero)
# to correct for maximal values
alpha <- predict(fit, x = max(xy[,"x"],na.rm=T),
deriv = 1)$y
if(alpha>0){
Am <- predict(fit, max(xy[,"x"],na.rm=T))$y+alpha*(3000-max(xy[,"x"],na.rm=T))
Am <- Am - Rd
}
theta1 <- rep(NA, 4)
light <- c(100, 200, 500, 800)
for (i in 1:4) {
y.pred
<- predict(fit, light[i])$y
theta1[i]
<- ((y.pred + Rd) * (alpha * light[i] + Am) -
alpha * Am * light[i])/(y.pred + Rd)^2
}
theta <- mean(theta1)
Rd<- -1*Rd
value <- c(theta, Am, alpha,
Rd)
names(value) <- mCall[c("theta",
"Am", "alpha", "Rd")]
value
}
and then using “selfStart”
NRhyperbola<-selfStart(NRhyperbola,initial=NRhyperbolaInit)
Bill Shipley