R help - Oct 2005 - self starting function for nonlinear least squares.

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

 


	[[alternative HTML version deleted]]