# R/SelfStartNls.r In OnofriAndreaPG/aomisc: Statistical methods for the agricultural sciences

#### Defines functions polyInv.4InitpolyInv.4meanpolyInv.3InitpolyInv.3meanAvsPFDInitAvsPFDMeanextremeValueInitextremeValueMeangompGrowth3InitgompGrowth3MeangompGrowth2InitgompGrowth2MeangompGrowth1InitgompGrowth1Meanlogistic6Initlogistic6Meanlogistic5Initlogistic5MeanlogiGrowth4InitlogiGrowth4MeanlogiGrowth3InitlogiGrowth3MeanlogiGrowth2InitlogiGrowth2MeanlogiGrowth1InitlogiGrowth1MeanmonoGrowthInitmonoGrowthMean

```#Negative exponential - Monomulecolar growth
monoGrowthMean <- function(predictor, a, b, c) {
x <- predictor
a - (a - b) * exp (- c * x)
}

monoGrowthInit <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]
plateau <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( 1 - (y / plateau ) )
coefs <- coef( lm(pseudoY ~ x) )
temp <- exp(coefs[1])
b <- plateau * (1 - temp)
c <- - coefs[2]
a <- plateau
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSmonoGrowth <- selfStart(monoGrowthMean, monoGrowthInit, parameters=c("a", "b", "c"))

#Logistic growth - 1 ###############################
logiGrowth1Mean <- function(predictor, a, b, c) {
x <- predictor
a / (1 + exp(- b * x + c))
}

logiGrowth1Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

a <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( (a / y ) - 1 )
coefs <- coef( lm(pseudoY ~ x) )
b <- - coefs[2]
c <- coefs[1]
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSlogiGrowth.1 <- selfStart(logiGrowth1Mean, logiGrowth1Init, parameters=c("a", "b", "c"))

#Logistic Growth 2
logiGrowth2Mean <- function(predictor, a, b, c) {
x <- predictor
a / (1 + b * exp(- c * x))
}

logiGrowth2Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

a <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( (a / (y + 0.0001) ) - 1 )
coefs <- coef( lm(pseudoY ~ x) )
c <- - coefs[2]
k <- coefs[1]
b <- exp(k)
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSlogiGrowth.2 <- selfStart(logiGrowth2Mean, logiGrowth2Init, parameters=c("a", "b", "c"))

#Logistic Growth 3
logiGrowth3Mean <- function(predictor, init, m, plateau) {
x <- predictor
init * plateau / (init + (plateau - init) * exp( - m * x))
}

logiGrowth3Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

plateau <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( (plateau / (y + 0.0001) ) - 1 )
coefs <- coef( lm(pseudoY ~ x) )
b <- exp(coefs[1])
init <- plateau / (1 + b)
m <- - coefs[2]
value <- c(init, m, plateau)
names(value) <- mCall[c("init", "m", "plateau")]
value
}

NLSlogiGrowth.3 <- selfStart(logiGrowth3Mean, logiGrowth3Init, parameters=c("init", "m", "plateau"))

#Logistic Growth 4
logiGrowth4Mean <- function(predictor, t50, m, plateau) {
x <- predictor
plateau / (1 + exp(- m * (x - t50)))}

logiGrowth4Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

plateau <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( (plateau / (y + 0.0001) ) - 1 )
coefs <- coef( lm(pseudoY ~ x) )
m <- - coefs[2]
t50 <- coefs[1] / m
value <- c(t50, m, plateau)
names(value) <- mCall[c("t50", "m", "plateau")]
value
}

NLSlogiGrowth.4 <- selfStart(logiGrowth4Mean, logiGrowth4Init, parameters=c("t50", "m", "plateau"))

#Logistic function 5
logistic5Mean <- function(predictor, a, b, c, d) {
x <- predictor
c + (d-c) / (1 + b * exp(a * x))
}

logistic5Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

d <- max(y) * 1.05
c <- min(y) -0.0001
## Linear regression on pseudo y values
pseudoY <- log((d - y) / (y - c ))
coefs <- coef( lm(pseudoY ~ x) )
a <- coefs[2]
k <- coefs[1]
b <- exp(k)
value <- c(a, b, c, d)
names(value) <- mCall[c("a", "b", "c", "d")]
value
}

NLSlogistic.1 <- selfStart(logistic5Mean, logistic5Init, parameters=c("a", "b", "c", "d"))

#Logistic function 6
logistic6Mean <- function(predictor, a, b) {
x <- predictor
1 / (1 + b * exp(- a * x))
}

logistic6Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

## Linear regression on pseudo y values
pseudoY <- log((1 - y) / y)
coefs <- coef( lm(pseudoY ~ x) )
a <- - coefs[2]
k <- coefs[1]
b <- exp(k)
value <- c(a, b)
names(value) <- mCall[c("a", "b")]
value
}

NLSlogistic.2 <- selfStart(logistic6Mean, logistic6Init, parameters=c("a", "b"))

#LOG_LOGISTIC FUNCTIONS##########################################################

#GOMPERTZ MODELS################################################################

#Gompertz growth - 1
gompGrowth1Mean <- function(predictor, a, m, c) {
x <- predictor
a * exp( - (m/c) * exp (-c * x))
}

gompGrowth1Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

plateau <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( - log( y / plateau ) )
coefs <- coef( lm(pseudoY ~ x) )
k <- coefs[1]; c <- - coefs[2]
b <- exp(k)
m <- b * c
a <- plateau
value <- c(a, m, c)
names(value) <- mCall[c("a", "m", "c")]
value
}

NLSgompGrowth.1 <- selfStart(gompGrowth1Mean, gompGrowth1Init, parameters=c("a", "m", "c"))

#Gompertz growth 2
gompGrowth2Mean <- function(predictor, a, b, c) {
x <- predictor
a * exp( - exp (b - c*x))
}

gompGrowth2Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

a <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( - log( y / a ) )
coefs <- coef( lm(pseudoY ~ x) )

k <- coefs[1]
c <- - coefs[2]
b <- k
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSgompGrowth.2 <- selfStart(gompGrowth2Mean, gompGrowth2Init, parameters=c("a", "b", "c"))

#Gompertz growth - 3
gompGrowth3Mean <- function(predictor, a, b, c) {
x <- predictor
a * exp( -b * exp (-c*x))
}

gompGrowth3Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

a <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( - log( y / a ) )
coefs <- coef( lm(pseudoY ~ x) )

k <- coefs[1]
c <- - coefs[2]
b <- exp(k)
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSgompGrowth.3 <- selfStart(gompGrowth3Mean, gompGrowth3Init, parameters=c("a", "b", "c"))

#Extreme value
extremeValueMean <- function(predictor, a, b, c) {
x <- predictor
a * (1 - exp( - exp (b - c*x)))
}

extremeValueInit <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]

a <- max(y) * 1.05

## Linear regression on pseudo y values
pseudoY <- log( - log( (a - y ) / a ) )
coefs <- coef( lm(pseudoY ~ x) )

k <- coefs[1]
c <- - coefs[2]
b <- k
value <- c(a, b, c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSextremeValue <- selfStart(extremeValueMean, extremeValueInit, parameters=c("a", "b", "c"))

#WEIBULL TYPE MODELS
#Weibull-1

#Modified Mitscherlich equation for A vs PFD relationships
AvsPFDMean <- function(predictor, Rd, Amax, Qapp) {
x <- predictor
Rd+Amax*(1-exp((-Qapp/Amax)*x))
}

AvsPFDInit <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]
interc <- min(y)
plateau <- max(y) * 1.05 - interc
## Linear regression on pseudo y values
pseudoY <- log( 1 - (((y - interc) / plateau ) ))
coefs <- coef( lm(pseudoY ~ x - 1) )
Amax <- plateau
Rd <- interc
b <-  coefs[1]
Qapp <- -b*plateau
value <- c(Rd,Amax,Qapp)
names(value) <- mCall[c("Rd", "Amax", "Qapp")]
value
}

NLSAvsPFD <- selfStart(AvsPFDMean, AvsPFDInit, parameters=c("Rd", "Amax", "Qapp"))

#Inverse polynomial
polyInv.3mean <- function(predictor, a, b, c) {
x <- predictor
1/(a + b*x + c*x^2)
}

polyInv.3Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]
## Linear regression on pseudo y values
#pseudoY <- 1/y
coefs <- coef(glm(y ~ x + I(x^2), family=gaussian(link="inverse")))
a <- coefs[1]; b <- coefs[2]; c <- coefs[3]

value <- c(a,b,c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSpolyInv.3 <- selfStart(polyInv.3mean, polyInv.3Init, parameters=c("a", "b", "c"))

#Inverse polynomial 2
polyInv.4mean <- function(predictor, a, b, c) {
x <- predictor
x/(a + b*x + c*x^2)
}

polyInv.4Init <- function(mCall, LHS, data, ...) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
x <-  xy[, "x"]; y <- xy[, "y"]
## Linear regression on pseudo y values
pseudoY <- y/x
coefs <- coef(glm(pseudoY ~ x + I(x^2), family=gaussian(link="inverse")))
a <- coefs[1]; b <- coefs[2]; c <- coefs[3]
print(a);print(b);print(c)
value <- c(a,b,c)
names(value) <- mCall[c("a", "b", "c")]
value
}

NLSpolyInv.4 <- selfStart(polyInv.4mean, polyInv.4Init, parameters=c("a", "b", "c"))
```
OnofriAndreaPG/aomisc documentation built on July 4, 2023, 6:53 p.m.