R/DHModels_NRIDbase.R
In gabonNRI/gabontreedata: What the package does (short line)
Defines functions
rmse
AIC_RSS
rse
BAfun
heightfit
####################################################
#### Calculate RMSE (Root mean square error)
rmse <- function(error){
sqrt(mean(error^2))}
####################################################
#### Calculate AIC from residual sum
AIC_RSS <- function(n, k, error){
##AIC = n * ln(RSS/n) + 2k + C, the constant C can be ignored.
AIC_fit <- (n * log(mean(error^2)) + 2*k) }
#### Calculate RSE (Residual standard error): sqrt(SSE/(n-p)), i.e., sqrt(RSS/df), where n = sample size, p = numbers of parameters, df=degree of freedom
rse <- function(error, n, p){
(sum(error^2)/(n-p)) }
####################################################
#### Calculate basal area
BAfun <- function(D){0.00007854*D^2 * 10000} # Multiply by 10000 to express in cm2
####################################################
heightfit <- function(hdat, sub.dat, E_Value){
Hmean <- hdat$Hmean
D <- hdat$D
n <- length(hdat$Hmean)
E_Value <- E_Value
sub.dat <- sub.dat
#Plot_name <- unique(sub.dat$Plot.Code)
##### Calculating RSD for Uncertainity & Organizing Output
results.fun <- function (fit, equation, Hmean, D, error, PredHt, coef){
AIC <- AIC_RSS(n = length(Hmean), length(fit$par), error)
RMSE <- rmse(error)
RSD <- sqrt(sum(error^2)/(length(Hmean)-2))
Xstand <- (D - mean(D))^2
SEM <- RSD * sqrt(1/length(Hmean) + Xstand/sum(Xstand))
res.df <- cbind(D, Hmean, PredHt, Residuals = error, SEM = SEM)
return(list(equation = equation, RSD = RSD, vars = c(coef, AIC = AIC, RMSE = RMSE), res = res.df))
}
plot.fun <- function(res, hdat){
a <- res$vars["a"]
b <- res$vars["b"]
c <- res$vars["c"]
d <- res$vars["d"]
equ.fun <- function(D) eval(parse(text = res$equ))
ggplot(hdat, aes(x = D, y = Hmean)) + geom_point(colour="grey", size = 3) +
xlab("Tree diameter, cm") + ylab("Tree height, m") + theme_bw() +
theme(legend.position = "none") +
stat_function(fun = equ.fun, colour = "darkgreen", size = 1)
#geom_text(data = data_labels, aes(x = time , y = value, label = label), size = 4)
}
##### EQ1: Linear model: H = a + b*D, where D = DBH
startvals1 <- list(a = 2, b = 1)
Fit_nls1 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a + b* D)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq1 <- optim(par=startvals1, fn = Fit_nls1, hessian = FALSE, control = list(parscale = unlist(startvals1)))
a1 = Eq1$par[1]
b1 = Eq1$par[2]
PredHt <- a1 + b1 * D
error <- (Hmean-PredHt)
RMSE_Eq1 <- rmse(error)
AIC_Eq1 <- AIC_RSS(n, length(startvals1), error)
equation <- as.formula(y ~ a + b * D)
res1 <- results.fun(fit = Eq1, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a1, b1))
##### EQ2: Quadratic model, Eq. 27 in Table 2 of page 128 (Huang et al., 2000): H = a + b*D + c*D^2, where D = DBH
startvals2 <- list(a = -5, b = 1, c = 0.01)
Fit_nls2 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + b * D + c * (D^2))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq2 <- optim(par=startvals2, fn = Fit_nls2, hessian = FALSE, control = list(parscale = unlist(startvals2)))
a2 = Eq2$par[1]
b2 = Eq2$par[2]
c2 = Eq2$par[3]
PredHt <- (a2 + b2 * D + c2 * (D^2))
error<- (Hmean-PredHt)
RMSE_Eq2 <- rmse(error)
AIC_Eq2 <- AIC_RSS(n, length(startvals2), error)
equation <- as.formula(y ~ a + b * D + c * (D^2))
res2 <- results.fun(fit = Eq2, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a2, b2, c2))
##### EQ3: Linear models on log scale (Log~Log): H=exp(a+b*ln(D)), where D = DBH (Eq2 of page 8617 (Djomo et al., 2010; Molto et al., 2013))
startvals3 <- list(a = 1, b =1)
Fit_nls3 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- exp(a + b*log(D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq3 <- optim(par=startvals3, fn = Fit_nls3, hessian = FALSE, control = list(parscale = unlist(startvals3)))
a3 = Eq3$par[1]
b3 = Eq3$par[2]
PredHt <- exp(a3 + b3*log(D))
error<- (Hmean-PredHt)
RMSE_Eq3 <- rmse(error)
AIC_Eq3 <- AIC_RSS(n, length(startvals3), error)
equation <- as.formula(y ~ exp(a + b*log(D)))
res3 <- results.fun(fit = Eq3, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a3, b3))
##### EQ4: Quadratic model on log scale: Ln(H) = a + b * ln(D) + c * ln(D^2) (Djomo et al., 2010)
startvals4 <- list(a = -5, b = 1, c = 0.01)
Fit_nls4 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + b * log(D) + c * log(D^2))
t1 <- sum((Hmean-exp(t1))^2)
return(t1)}
Eq4 <- optim(par=startvals4, fn = Fit_nls4, hessian = FALSE, control = list(parscale = unlist(startvals4)))
a4 = Eq4$par[1]
b4 = Eq4$par[2]
c4 = Eq4$par[3]
PredHt <- a4 + b4 * log(D) + c4 * log(D^2)
error<- (Hmean-PredHt)
RMSE_Eq4 <- rmse(error)
AIC_Eq4 <- AIC_RSS(n, length(startvals4), error)
equation <- as.formula(y ~ a + b * log(D) + c * log(D^2))
res4 <- results.fun(fit = Eq4, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a4, b4, c4))
##### EQ5: Non-linear model: H = a + b * ln(D), where D = DBH (Fang and Bailey, 1998; Djomo et al., 2010; Molto et al., 2013)
##### Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
startvals5 <- list(a = 1, b =2)
Fit_nls5 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a + b*log(D)) #Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq5 <- optim(par=startvals5, fn = Fit_nls5, hessian = FALSE, control = list(parscale = unlist(startvals5)))
a5 = Eq5$par[1]
b5 = Eq5$par[2]
PredHt <- a5 + b5*log(D)
error<- (Hmean-PredHt)
RMSE_Eq5 <- rmse(error)
AIC_Eq5 <- AIC_RSS(n, length(startvals5), error)
equation <- as.formula(y ~ a + b*log(D))
res5 <- results.fun(fit = Eq5, equation = equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a5, b5))
##### EQ6: Simplified (two-parameter) Weibull equation: H = a*(1-exp(-D/b)), where D=DBH, page 8617 (Molto et al., 2013)
startvals6 <- list(a = 30, b = 5)
Fit_nls6 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * (1-exp(I(-D/b)))) #Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq6 <- optim(par=startvals6, fn = Fit_nls6, hessian = FALSE, control = list(parscale = unlist(startvals6)))
a6 = Eq6$par[1]
b6 = Eq6$par[2]
PredHt <- a6 * (1-exp(I(-D/b6)))
error<- (Hmean-PredHt)
RMSE_Eq6 <- rmse(error)
AIC_Eq6 <- AIC_RSS(n, length(startvals6), error)
equation <- as.formula(y ~ a * (1-exp(I(-D/b))))
res6 <- results.fun(fit = Eq6, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a6, b6))
##### EQ7: Three-parameter Weibull equation: H = a * (1-exp(-b*D^c)), where D = DBH
##### Equation 13 of Huang et al., 1992; Lewis et al., 2009; Feldpausch et al., 2012; Mitchard et al., 2014)
startvals7 <- list(a = 54.01, b = -0.053, c = 0.0759)
Fit_nls7 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a*(1-exp(-b*D^c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq7 <- optim(par=startvals7, fn = Fit_nls7, hessian = FALSE, control = list(parscale = unlist(startvals7)))
a7 = Eq7$par[1]
b7 = Eq7$par[2]
c7 = Eq7$par[3]
PredHt <- a7*(1-exp(-b7*D^c7))
error<- (Hmean-PredHt)
RMSE_Eq7 <- rmse(error)
AIC_Eq7 <- AIC_RSS(n, length(startvals7), error)
equation <- as.formula(y ~ a*(1-exp(-b*D^c)))
res7 <- results.fun(fit = Eq7, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a7, b7, c7))
##### EQ8: Michaelis-Menten equation: H = a*D /(b+D), in page 8618 (Molto et al., 2013; Equation 3 in Huang et al. (1992))
startvals8 <- list(a = 10, b = 2) #Start points were estimated from Figure 1 in page 8632 (Molto et al., 2013)
Fit_nls8 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * I(D)/(b + I(D)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq8 <- optim(par=startvals8, fn = Fit_nls8, hessian = FALSE, control = list(parscale = unlist(startvals8)))
a8 = Eq8$par[1]
b8 = Eq8$par[2]
PredHt <- a8 * I(D)/(b8 + I(D))
error<- (Hmean-PredHt)
RMSE_Eq8 <- rmse(error)
AIC_Eq8 <- AIC_RSS(n, length(startvals8), error)
equation <- as.formula(y ~ (a * I(D)/(b + I(D))))
res8 <- results.fun(fit = Eq8, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a8, b8))
##### EQ9: Power-law equation: H=a*D^b (Eq 1 of Huang et al., 1992; Chave et al., 2005; Feldpausch et al., 2012; Ngomanda et al., 2014)
##### Start points were estimated from Ngomanda et al.(2014)
startvals9 <- list(a = 26.6049, b = -9.4854)
Fit_nls9 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- ( a* D^b )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq9 <- optim(par=startvals9, fn = Fit_nls9, hessian = FALSE, control = list(parscale = unlist(startvals9)))
a9 = Eq9$par[1]
b9 = Eq9$par[2]
PredHt <- a9* D^b9
error<- (Hmean-PredHt)
RMSE_Eq9 <- rmse(error)
AIC_Eq9 <- AIC_RSS(n, length(startvals9), error)
equation <- as.formula(y ~ a* D^b)
res9 <- results.fun(fit = Eq9, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a9, b9))
##### EQ10:Model 3: Ln(H) = a + b/D or H = a * exp(b/D), where D=DBH (Djomo et al., 2010; Equation 6 in Huang et al. (1992))
startvals10 <- list(a = 26.6049, b = -9.4854)
Fit_nls10 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (exp(a + b/D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq10 <- optim(par=startvals10, fn = Fit_nls10, hessian = FALSE, control = list(parscale = unlist(startvals10)))
a10 = Eq10$par[1]
b10 = Eq10$par[2]
PredHt <- exp(a10 + b10/D)
error<- (Hmean-PredHt)
RMSE_Eq10 <- rmse(error)
AIC_Eq10 <- AIC_RSS(n, length(startvals10), error)
equation <- as.formula(y ~ exp(a + b/D))
res10 <- results.fun(fit = Eq10, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a10, b10))
##### EQ11: H = a + b*(1-exp(-c*(D-Dmin))), where D=DBH (Fang and Bailey, 1998)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals11 <- list(a = 44.7, b =42.3, c = 0.026)
Fit_nls11 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- a + b*(1-exp(-c*(D-min(D))))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq11 <- optim(par=startvals11, fn = Fit_nls11, hessian = FALSE, control = list(parscale = unlist(startvals11)))
a11 = Eq11$par[1]
b11 = Eq11$par[2]
c11 = Eq11$par[3]
PredHt <- a11 + b11 * (1-exp(-c11 * (D - min(D))))
error<- (Hmean-PredHt)
RMSE_Eq11 <- rmse(error)
AIC_Eq11 <- AIC_RSS(n, length(startvals11), error)
equation <- as.formula(y ~ a + b*(1-exp(-c*(D-min(D)))))
res11 <- results.fun(fit = Eq11, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a11, b11, c11))
##### EQ12:H = a + (a-1.3)*b/(D+b), where D=DBH (Fang and Bailey, 1998)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals12 <- list(a = 33.62, b =20.88)
Fit_nls12 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + (a - 1.3) * b/(D + b))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq12 <- optim(par=startvals12, fn = Fit_nls12, hessian = FALSE, control = list(parscale = unlist(startvals12)))
a12 = Eq12$par[1]
b12 = Eq12$par[2]
c12 = Eq12$par[3]
PredHt <- a12 + (a12 - 1.3) * b12/(D + b12)
error<- (Hmean-PredHt)
RMSE_Eq12 <- rmse(error)
AIC_Eq12 <- AIC_RSS(n, length(startvals12), error)
equation <- as.formula(y ~ (a + (a - 1.3) * b/(D + b)))
res12 <- results.fun(fit = Eq12, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a12, b12, c12))
##### EQ13:H = 1.3 + exp(a + b*D^c), where D=DBH (Eq 10 of Huang et al. (1992, 2000); Colbert et al., 2002)
##### Start points were for Black williows in page 173 (Colbert et al., 2002)
startvals13 <- list(a = 4.5535, b =-3.7529, c = 1)
Fit_nls13 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + exp(a + b * D^c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq13 <- optim(par=startvals13, fn = Fit_nls13, hessian = FALSE, control = list(parscale = unlist(startvals13)))
a13 = Eq13$par[1]
b13 = Eq13$par[2]
c13 = Eq13$par[3]
PredHt <- 1.3 + exp(a13 + b13 * D^c13)
error<- (Hmean-PredHt)
RMSE_Eq13 <- rmse(error)
AIC_Eq13 <- AIC_RSS(n, length(startvals13), error)
equation <- as.formula(y ~ 1.3 + exp(a + b * D^c))
res13 <- results.fun(fit = Eq13, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a13, b13, c13))
##### EQ14:H = 1.3 + a*(1-exp(-b*D))^c, where D=DBH (Equation 12 in Huang et al. (1992, 2000); Equation 2 in Sharma and Parton (2007))
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals14 <- list(a = 1.0462, b = 1.0464, c=0.9465)
Fit_nls14 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*(1-exp(-b*D))^c)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq14 <- optim(par=startvals14, fn = Fit_nls14, hessian = FALSE, control = list(parscale = unlist(startvals14)))
a14 = Eq14$par[1]
b14 = Eq14$par[2]
c14 = Eq14$par[3]
PredHt <- 1.3 + a14 *(1-exp(-b14 * D))^c14
error<- (Hmean-PredHt)
RMSE_Eq14 <- rmse(error)
AIC_Eq14 <- AIC_RSS(n, length(startvals14), error)
equation <- as.formula(y ~ 1.3 + a*(1-exp(-b*D))^c)
res14 <- results.fun(fit = Eq14, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a14, b14, c14))
##### Following equations from Table 3 in Huang et al. (1992) and Table 2 in Huang et al. (2000).
##### Parameter estimations of Balsam poplar (broad leaf) in Tables 4 and 5 (Huang et al., 1992)
##### were used as start values in the new fitting for Gabonese rainforests.
##### EQ15: Equation 4 in Huang et al.(1992 and 2000): H = 1.3 + a * (1 - exp(-bD)), where D = DBH
startvals15 <- list(a = 25.3302, b = 0.0512)
Fit_nls15 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * (1 - exp(-b * D)) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq15 <- optim(par=startvals15, fn = Fit_nls15, hessian = FALSE, control = list(parscale = unlist(startvals15)))
a15 = Eq15$par[1]
b15 = Eq15$par[2]
PredHt <- a15 * (1 - exp(-b15 * D))
error<- (Hmean-PredHt)
RMSE_Eq15 <- rmse(error)
AIC_Eq15 <- AIC_RSS(n, length(startvals15), error)
equation <- as.formula(y ~ a * (1 - exp(-b * D)))
res15 <- results.fun(fit = Eq15, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a15, b15))
##### EQ16:Equation 5 in Huang et al. (1992 and 2000): H = 1.3 + D^2 /(a + bD)^2, where D = DBH
startvals16 <- list(a = 1.3209, b = 0.1813)
Fit_nls16 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + D^2 /((a + b * D)^2) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq16 <- optim(par=startvals16, fn = Fit_nls16, hessian = FALSE, control = list(parscale = unlist(startvals16)))
a16 = Eq16$par[1]
b16 = Eq16$par[2]
PredHt <- 1.3 + D^2 /((a16 + b16 * D)^2)
error<- (Hmean-PredHt)
RMSE_Eq16 <- rmse(error)
AIC_Eq16 <- AIC_RSS(n, length(startvals16), error)
equation <- as.formula(y ~ 1.3 + D^2 /((a + b * D)^2))
res16 <- results.fun(fit = Eq16, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a16, b16))
##### EQ17: Equation 7 in Huang et al. (1992 and 2000): H = 1.3 + 10^a*(D^b), where D = DBH
startvals17 <- list(a = 0.4305, b = 0.5871)
Fit_nls17 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + (10^a) * (D^b) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq17 <- optim(par=startvals17, fn = Fit_nls17, hessian = FALSE, control = list(parscale = unlist(startvals17)))
a17 = Eq17$par[1]
b17 = Eq17$par[2]
PredHt <- 1.3 + (10^a17) * (D^b17)
error<- (Hmean-PredHt)
RMSE_Eq17 <- rmse(error)
AIC_Eq17 <- AIC_RSS(n, length(startvals17), error)
equation <- as.formula(y ~ 1.3 + (10^a) * (D^b))
res17 <- results.fun(fit = Eq17, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a17, b17))
##### EQ18:Equation 8 in Huang et al. (1992 and 2000): H = 1.3 + (a*D)/(D+1) + b*D, where D = DBH
startvals18 <- list(a = 6.5487, b = 0.4507)
Fit_nls18 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a*D/(D+1) + b*D )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq18 <- optim(par=startvals18, fn = Fit_nls18, hessian = FALSE, control = list(parscale = unlist(startvals18)))
a18 = Eq18$par[1]
b18 = Eq18$par[2]
PredHt <- 1.3 + a18 * D/(D+1) + b18*D
error<- (Hmean-PredHt)
RMSE_Eq18 <- rmse(error)
AIC_Eq18 <- AIC_RSS(n, length(startvals18), error)
equation <- formula(y ~ 1.3 + a*D/(D+1) + b*D, showEnv = F)
res18 <- results.fun(fit = Eq18, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a18, b18))
##### EQ19:Equation 9 in Huang et al. (1992 and 2000): H = 1.3 + a*(D/(D+1))^b, where D = DBH
startvals19 <- list(a = 27.1752, b = 10.1979)
Fit_nls19 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a*(D/(1 + D))^b )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq19 <- optim(par=startvals19, fn = Fit_nls19, hessian = FALSE, control = list(parscale = unlist(startvals19)))
a19 = Eq19$par[1]
b19 = Eq19$par[2]
PredHt <- 1.3 + a19*(D/(1 + D))^b19
error<- (Hmean-PredHt)
RMSE_Eq19 <- rmse(error)
AIC_Eq19 <- AIC_RSS(n, length(startvals19), error)
equation <- as.formula(y ~ 1.3 + a*(D/(1 + D))^b)
res19 <- results.fun(fit = Eq19, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a19, b19))
##### EQ20:Equation 10 in Huang et al. (1992 and 2000): H = 1.3 + exp(a + b*D^c), where D = DBH
startvals20 <- list(a = 18.804146, b = -16.760486, c=-0.025197)
Fit_nls20 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- 1.3 + exp(a + b * (D^c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq20 <- optim(par=startvals20, fn = Fit_nls20, hessian = FALSE, control = list(parscale = unlist(startvals20)))
a20 = Eq20$par[1]
b20 = Eq20$par[2]
c20 = Eq20$par[3]
PredHt <- 1.3 + exp(a20 + b20 * (D^c20))
error<- (Hmean-PredHt)
RMSE_Eq20 <- rmse(error)
AIC_Eq20 <- AIC_RSS(n, length(startvals20), error)
equation <- as.formula(y ~ 1.3 + exp(a + b * (D^c)))
res20 <- results.fun(fit = Eq20, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a20, b20, c20))
##### EQ21: Equation 11 in Huang et al. (1992 and 2000): H = 1.3 + a/(1+b*exp(-c*D)), where D = DBH
startvals21 <- list(a = 2.5241, b = 0.0012, c=0.1404)
Fit_nls21 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- 1.3 + a/(1 + b*exp(-c*D)) #Equation 11 in Huang et al. (1992)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq21 <- optim(par=startvals21, fn = Fit_nls21, hessian = FALSE, control = list(parscale = unlist(startvals21)))
a21 = Eq21$par[1]
b21 = Eq21$par[2]
c21 = Eq21$par[3]
PredHt <- 1.3 + a21 /(1 + b21 *exp(-c21 * D))
error <- (Hmean - PredHt)
RMSE_Eq21 <- rmse(error)
AIC_Eq21 <- AIC_RSS(n, length(startvals21), error)
equation <- as.formula(y ~ 1.3 + a/(1 + b*exp(-c*D)))
res21 <- results.fun(fit = Eq21, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a21, b21, c21))
##### EQ22: Equation 14 in Huang et al. (1992 and 2000): H = 1.3 + a*exp(-b*exp(-c*D)), where D = DBH (For Eqs. 12 and 13 in Huang et al. (1992 and 2002), please see the above eqs 14 and 7 in the file.)
startvals22 <- list(a = 1.6368, b = 2.1570, c=1.0951)
Fit_nls22 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*exp(-b*exp(-c*D)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq22 <- optim(par=startvals22, fn = Fit_nls22, hessian = FALSE, control = list(parscale = unlist(startvals22)))
a22 = Eq22$par[1]
b22 = Eq22$par[2]
c22 = Eq22$par[3]
PredHt <- 1.3 + a22 * exp(-b22 * exp(-c22 * D))
error <- (Hmean - PredHt)
RMSE_Eq22 <- rmse(error)
AIC_Eq22 <- AIC_RSS(n, length(startvals22), error)
equation <- as.formula(y ~ 1.3 + a*exp(-b*exp(-c*D)))
res22 <- results.fun(fit = Eq22, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a22, b22, c22))
##### EQ23: Equation 15 in Huang et al. (1992) and Eq. 16 in Huang et al. (2000): H = (1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b)
startvals23 <- list(a = 0.0464, b = 1.0716, c=27.0745)
Fit_nls23 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- ((1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq23 <- optim(par=startvals23, fn = Fit_nls23, hessian = FALSE, control = list(parscale = unlist(startvals23)))
a23 = Eq23$par[1]
b23 = Eq23$par[2]
c23 = Eq23$par[3]
PredHt <- (1.3^b23 + (c23^b23-1.3^b23)*(1-exp(-a23*D))/(1-exp(-a23*100)))^(1/b23)
error <- (Hmean - PredHt)
RMSE_Eq23 <- rmse(error)
AIC_Eq23 <- AIC_RSS(n, length(startvals23), error)
equation <- as.formula(y ~ (1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b))
res23 <- results.fun(fit = Eq23, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a23, b23, c23))
##### EQ24: Equation 16 in Huang et al. (1992) and Eq. 17 in Huang et al. (2000): H =1.3 + D^2/(a + b*D + c*D^2)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals24 <- list(a = 0.0038, b = 0.7027, c = 0.0270)
Fit_nls24 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + D^2 /(a + b * D + c * D^2))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq24 <- optim(par=startvals24, fn = Fit_nls24, hessian = FALSE, control = list(parscale = unlist(startvals24)))
a24 = Eq24$par[1]
b24 = Eq24$par[2]
c24 = Eq24$par[3]
PredHt <- 1.3 + D^2 /(a24 + b24 * D + c24 * D^2)
error <- (Hmean - PredHt)
RMSE_Eq24 <- rmse(error)
AIC_Eq24 <- AIC_RSS(n, length(startvals24), error)
equation <- as.formula(y ~ 1.3 + D^2 /(a + b * D + c * D^2))
res24 <- results.fun(fit = Eq24, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a24, b24, c24))
##### EQ25: Equation 17 in Huang et al. (1992) and Eq. 18 in Huang et al. (2000): H = 1.3 + a *D^(b*D^-c), where D = DBH
startvals25 <- list(a = 22.9433, b = -20.9985, c=0.7680)
Fit_nls25 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*D^(b*D^-c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq25 <- optim(par=startvals25, fn = Fit_nls25, hessian = FALSE, control = list(parscale = unlist(startvals25)))
a25 = Eq25$par[1]
b25 = Eq25$par[2]
c25 = Eq25$par[3]
PredHt <- 1.3 + a25 * D^(b25 * D^-c25)
error <- (Hmean - PredHt)
RMSE_Eq25 <- rmse(error)
AIC_Eq25 <- AIC_RSS(n, length(startvals25), error)
equation <- as.formula(y ~ 1.3 + a*D^(b*D^-c))
res25 <- results.fun(fit = Eq25, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a25, b25, c25))
##### EQ26: Equation 18 in Huang et al. (1992) and Eq. 20 in Huang et al. (2000): H = 1.3 + a*exp(b/(D+c)), where D = DBH
startvals26 <- list(a = 3.2971, b = -0.4014, c = 5.5088)
Fit_nls26 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a * exp(b /(D+c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq26 <- optim(par=startvals26, fn = Fit_nls26, hessian = FALSE, control = list(parscale = unlist(startvals26)))
a26 = Eq26$par[1]
b26 = Eq26$par[2]
c26 = Eq26$par[3]
PredHt <- 1.3 + a26 * exp(b26 /(D + c26))
error <- (Hmean - PredHt)
RMSE_Eq26 <- rmse(error)
AIC_Eq26 <- AIC_RSS(n, length(startvals26), error)
equation <- as.formula(y ~ 1.3 + a * exp(b /(D+c)))
res26 <- results.fun(fit = Eq26, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a26, b26, c26))
##### EQ27: Equation 19 in Huang et al.(1992) and Eq.19 in Huang et al. (2000): H = 1.3 + a/(1 + b^-1 * D^-c), where D = DBH
startvals27 <- list(a = 34.4682, b = 0.0369, c = 1.0589)
Fit_nls27 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a/(1+1/(b*D^c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq27 <- optim(par=startvals27, fn = Fit_nls27, hessian = FALSE, control = list(parscale = unlist(startvals27)))
a27 = Eq27$par[1]
b27 = Eq27$par[2]
c27 = Eq27$par[3]
PredHt <- 1.3+a27 /(1+1/(b27 * D^c27))
error <- (Hmean - PredHt)
RMSE_Eq27 <- rmse(error)
AIC_Eq27 <- AIC_RSS(n, length(startvals27), error)
equation <- as.formula(y ~ 1.3 + a/(1+1/(b*D^c)))
res27 <- results.fun(fit = Eq27, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a27, b27, c27))
##### EQ28: Equation 20 in Huang et al. (1992) and Eq. 21 in Huang et al. (2000): H = 1.3 + a*(1 - b*exp(-c*D))^d, where D = DBH
startvals28 <- list(a = 15.2716, b = 0.9574, c = 0.0530, d = 1.1025)
Fit_nls28 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
d <- x[4]
t1 <- (1.3 + a*(1-b*exp(-c*D))^d)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq28 <- optim(par=startvals28, fn = Fit_nls28, hessian = FALSE, control = list(parscale = unlist(startvals28)))
a28 = Eq28$par[1]
b28 = Eq28$par[2]
c28 = Eq28$par[3]
d28 = Eq28$par[4]
PredHt<- 1.3 + a28 *(1- b28 * exp(-c28 * D))^d28
error <- (Hmean - PredHt)
RMSE_Eq28 <- rmse(error)
AIC_Eq28 <- AIC_RSS(n, length(startvals28), error)
equation <- as.formula(y ~ 1.3 + a*(1-b*exp(-c*D))^d)
res28 <- results.fun(fit = Eq28, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a28, b28, c28, d28))
##### EQ29: Equation 23 in Huang et al. (2000): H = 1.3 + a*D*exp(-b*D), where D = DBH.
##### Start values were for white spruce (Table 3 in page 131 of Huang et al. (2000))
startvals29 <- list(a = 1.0264, b = 0.01157)
Fit_nls29 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a * D * exp(-b*D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq29 <- optim(par=startvals29, fn = Fit_nls29, hessian = FALSE, control = list(parscale = unlist(startvals29)))
a29 = Eq29$par[1]
b29 = Eq29$par[2]
PredHt <- 1.3 + a29 * D * exp(-b29 * D)
error <- (Hmean - PredHt)
RMSE_Eq29 <- rmse(error)
AIC_Eq29 <- AIC_RSS(n, length(startvals29), error)
equation <- as.formula(y ~ 1.3 + a * D * exp(-b*D))
res29 <- results.fun(fit = Eq29, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a29, b29))
##### EQ30: Equation 24 in Huang et al. (2000): H = 1.3 + a * D^b * exp(-c*D), where D = DBH. Start values were for white spruce (Table 3 in page 131 of Huang et al. (2000))
startvals30 <- list(a = 0.687824, b = 1.175135, c = 0.017699)
Fit_nls30 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a * D^b * exp(-c*D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq30 <- optim(par=startvals30, fn = Fit_nls30, hessian = FALSE, control = list(parscale = unlist(startvals30)))
a30 = Eq30$par[1]
b30 = Eq30$par[2]
c30 = Eq30$par[3]
PredHt <- 1.3 + a30 * D^b30 * exp(-c30*D)
error <- (Hmean - PredHt)
RMSE_Eq30 <- rmse(error)
AIC_Eq30 <- AIC_RSS(n, length(startvals30), error)
equation <- as.formula(y ~ 1.3 + a * D^b * exp(-c*D))
res30 <- results.fun(fit = Eq30, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a30, b30, c30))
##### EQ31: Pantropical DBH-H model (Chave et al., 2014. GCB 20: 3183): Ln(H) = a - E + b*Ln(D)+ c*(Ln(D))^2 , where D = DBH
startvals31 <- list(a = 0.893, b = 0.760, c = -0.034)
Fit_nls31 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a - E_Value + b * log(D) + c * (log(I(D)))^2)
t1 <- sum((Hmean-t1)^2)
return(t1)}
if(!is.na(E_Value)){
Eq31 <- optim(par=startvals31, fn = Fit_nls31, hessian = FALSE, control = list(parscale = unlist(startvals31)))
a31 = Eq31$par[1]
b31 = Eq31$par[2]
c31 = Eq31$par[3]
PredHt <- a31 - E_Value + b31 * log(D) + c31 * (log(I(D)))^2
error <- (Hmean - PredHt)
RMSE_Eq31 <- rmse(error)
AIC_Eq31 <- AIC_RSS(n, length(startvals31), error)
equation <- as.formula(y ~ a - E_Value + b * log(D) + c * (log(I(D)))^2)
res31 <- results.fun(fit = Eq31, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a31, b31, c31))
}else{AIC_Eq31 <- 9999999}
##### EQ32: The Mitscherlisch model (Ngomanda et al., 2014): H = a - b * exp(-cD), where D = DBH
startvals32 <- list(a = 44.7, b = 42.3, c = 0.026)
Fit_nls32 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a - b * exp(-c * D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq32 <- optim(par=startvals32, fn = Fit_nls32, hessian = FALSE, control = list(parscale = unlist(startvals32)))
a32 = Eq32$par[1]
b32 = Eq32$par[2]
c32 = Eq32$par[3]
PredHt <- a32 - b32 * exp(-c32 * D)
error <- (Hmean - PredHt)
RMSE_Eq32 <- rmse(error)
AIC_Eq32 <- AIC_RSS(n, length(startvals32), error)
equation <- as.formula(y ~ a - b * exp(-c * D))
res32 <- results.fun(fit = Eq32, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a32, b32, c32))
##### Model comparison
AIC_best.fit <- paste("Eq", 1, sep="")
AIC_min_temp = AIC_Eq1
AIC.df <- matrix(nrow = 32, ncol = 2, dimnames = list(c(), c("Equ", "AIC")))
for (k in 1: 32){
AIC_min <- eval(parse(text=paste("AIC_Eq", k, sep="")))
AIC.df[k,] <- c(Equ = paste("Eq", k, sep=""), AIC = round(AIC_min, 1))
AIC_min <- min(AIC_min, AIC_min_temp)
if (! is.nan(AIC_min) && !is.nan(AIC_min_temp) ) {
if ( AIC_min < AIC_min_temp) {AIC_best.fit <- paste("Eq", k, sep="")}
}
AIC_min_temp = AIC_min
}
AIC.df <- data.frame(AIC.df)
AIC.df$AIC <- as.numeric(as.character(AIC.df$AIC))
best.equ <- eval(parse(text=AIC_best.fit))
rnum <- paste("res", gsub("[^[:digit:]]", "", AIC_best.fit), sep="")
res.fin <- get(rnum)
# pltname <- gsub('([[:punct:]])|\\s+','_',Plot_name)
# data_labels <- data.frame(time = c(20), value = c(0.95*max(Hmean)), label = pltname)
## Disable plotting
# gp <- plot.fun(res = res.fin, hdat = hdat)
# print(gp)
# pdf(file = paste("HDgraph", ".pdf", sep=""), width = 4, height = 4)
# print(gp)
# dev.off()
predht.fun <- function(res, sub.dat){
a <- res$vars["a"]
b <- res$vars["b"]
c <- res$vars["c"]
d <- res$vars["d"]
equ.fun <- function(D) eval(parse(text = res$equ))
ph <- equ.fun(sub.dat$D)
}
predht <- predht.fun(res = res.fin, sub.dat)
full.list <- res.fin
full.list$predht <- predht
full.list$Xequ <- as.character(res.fin$equation)[3]
full.list
}
gabonNRI/gabontreedata documentation built on Aug. 26, 2019, 12:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions rmse AIC_RSS rse BAfun heightfit
####################################################
#### Calculate RMSE (Root mean square error)
rmse <- function(error){
sqrt(mean(error^2))}
####################################################
#### Calculate AIC from residual sum
AIC_RSS <- function(n, k, error){
##AIC = n * ln(RSS/n) + 2k + C, the constant C can be ignored.
AIC_fit <- (n * log(mean(error^2)) + 2*k) }
#### Calculate RSE (Residual standard error): sqrt(SSE/(n-p)), i.e., sqrt(RSS/df), where n = sample size, p = numbers of parameters, df=degree of freedom
rse <- function(error, n, p){
(sum(error^2)/(n-p)) }
####################################################
#### Calculate basal area
BAfun <- function(D){0.00007854*D^2 * 10000} # Multiply by 10000 to express in cm2
####################################################
heightfit <- function(hdat, sub.dat, E_Value){
Hmean <- hdat$Hmean
D <- hdat$D
n <- length(hdat$Hmean)
E_Value <- E_Value
sub.dat <- sub.dat
#Plot_name <- unique(sub.dat$Plot.Code)
##### Calculating RSD for Uncertainity & Organizing Output
results.fun <- function (fit, equation, Hmean, D, error, PredHt, coef){
AIC <- AIC_RSS(n = length(Hmean), length(fit$par), error)
RMSE <- rmse(error)
RSD <- sqrt(sum(error^2)/(length(Hmean)-2))
Xstand <- (D - mean(D))^2
SEM <- RSD * sqrt(1/length(Hmean) + Xstand/sum(Xstand))
res.df <- cbind(D, Hmean, PredHt, Residuals = error, SEM = SEM)
return(list(equation = equation, RSD = RSD, vars = c(coef, AIC = AIC, RMSE = RMSE), res = res.df))
}
plot.fun <- function(res, hdat){
a <- res$vars["a"]
b <- res$vars["b"]
c <- res$vars["c"]
d <- res$vars["d"]
equ.fun <- function(D) eval(parse(text = res$equ))
ggplot(hdat, aes(x = D, y = Hmean)) + geom_point(colour="grey", size = 3) +
xlab("Tree diameter, cm") + ylab("Tree height, m") + theme_bw() +
theme(legend.position = "none") +
stat_function(fun = equ.fun, colour = "darkgreen", size = 1)
#geom_text(data = data_labels, aes(x = time , y = value, label = label), size = 4)
}
##### EQ1: Linear model: H = a + b*D, where D = DBH
startvals1 <- list(a = 2, b = 1)
Fit_nls1 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a + b* D)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq1 <- optim(par=startvals1, fn = Fit_nls1, hessian = FALSE, control = list(parscale = unlist(startvals1)))
a1 = Eq1$par[1]
b1 = Eq1$par[2]
PredHt <- a1 + b1 * D
error <- (Hmean-PredHt)
RMSE_Eq1 <- rmse(error)
AIC_Eq1 <- AIC_RSS(n, length(startvals1), error)
equation <- as.formula(y ~ a + b * D)
res1 <- results.fun(fit = Eq1, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a1, b1))
##### EQ2: Quadratic model, Eq. 27 in Table 2 of page 128 (Huang et al., 2000): H = a + b*D + c*D^2, where D = DBH
startvals2 <- list(a = -5, b = 1, c = 0.01)
Fit_nls2 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + b * D + c * (D^2))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq2 <- optim(par=startvals2, fn = Fit_nls2, hessian = FALSE, control = list(parscale = unlist(startvals2)))
a2 = Eq2$par[1]
b2 = Eq2$par[2]
c2 = Eq2$par[3]
PredHt <- (a2 + b2 * D + c2 * (D^2))
error<- (Hmean-PredHt)
RMSE_Eq2 <- rmse(error)
AIC_Eq2 <- AIC_RSS(n, length(startvals2), error)
equation <- as.formula(y ~ a + b * D + c * (D^2))
res2 <- results.fun(fit = Eq2, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a2, b2, c2))
##### EQ3: Linear models on log scale (Log~Log): H=exp(a+b*ln(D)), where D = DBH (Eq2 of page 8617 (Djomo et al., 2010; Molto et al., 2013))
startvals3 <- list(a = 1, b =1)
Fit_nls3 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- exp(a + b*log(D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq3 <- optim(par=startvals3, fn = Fit_nls3, hessian = FALSE, control = list(parscale = unlist(startvals3)))
a3 = Eq3$par[1]
b3 = Eq3$par[2]
PredHt <- exp(a3 + b3*log(D))
error<- (Hmean-PredHt)
RMSE_Eq3 <- rmse(error)
AIC_Eq3 <- AIC_RSS(n, length(startvals3), error)
equation <- as.formula(y ~ exp(a + b*log(D)))
res3 <- results.fun(fit = Eq3, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a3, b3))
##### EQ4: Quadratic model on log scale: Ln(H) = a + b * ln(D) + c * ln(D^2) (Djomo et al., 2010)
startvals4 <- list(a = -5, b = 1, c = 0.01)
Fit_nls4 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + b * log(D) + c * log(D^2))
t1 <- sum((Hmean-exp(t1))^2)
return(t1)}
Eq4 <- optim(par=startvals4, fn = Fit_nls4, hessian = FALSE, control = list(parscale = unlist(startvals4)))
a4 = Eq4$par[1]
b4 = Eq4$par[2]
c4 = Eq4$par[3]
PredHt <- a4 + b4 * log(D) + c4 * log(D^2)
error<- (Hmean-PredHt)
RMSE_Eq4 <- rmse(error)
AIC_Eq4 <- AIC_RSS(n, length(startvals4), error)
equation <- as.formula(y ~ a + b * log(D) + c * log(D^2))
res4 <- results.fun(fit = Eq4, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a4, b4, c4))
##### EQ5: Non-linear model: H = a + b * ln(D), where D = DBH (Fang and Bailey, 1998; Djomo et al., 2010; Molto et al., 2013)
##### Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
startvals5 <- list(a = 1, b =2)
Fit_nls5 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a + b*log(D)) #Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq5 <- optim(par=startvals5, fn = Fit_nls5, hessian = FALSE, control = list(parscale = unlist(startvals5)))
a5 = Eq5$par[1]
b5 = Eq5$par[2]
PredHt <- a5 + b5*log(D)
error<- (Hmean-PredHt)
RMSE_Eq5 <- rmse(error)
AIC_Eq5 <- AIC_RSS(n, length(startvals5), error)
equation <- as.formula(y ~ a + b*log(D))
res5 <- results.fun(fit = Eq5, equation = equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a5, b5))
##### EQ6: Simplified (two-parameter) Weibull equation: H = a*(1-exp(-D/b)), where D=DBH, page 8617 (Molto et al., 2013)
startvals6 <- list(a = 30, b = 5)
Fit_nls6 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * (1-exp(I(-D/b)))) #Log-linear eq.1 in page 8617 (Fang and Bailey, 1998; Molto et al., 2013)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq6 <- optim(par=startvals6, fn = Fit_nls6, hessian = FALSE, control = list(parscale = unlist(startvals6)))
a6 = Eq6$par[1]
b6 = Eq6$par[2]
PredHt <- a6 * (1-exp(I(-D/b6)))
error<- (Hmean-PredHt)
RMSE_Eq6 <- rmse(error)
AIC_Eq6 <- AIC_RSS(n, length(startvals6), error)
equation <- as.formula(y ~ a * (1-exp(I(-D/b))))
res6 <- results.fun(fit = Eq6, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a6, b6))
##### EQ7: Three-parameter Weibull equation: H = a * (1-exp(-b*D^c)), where D = DBH
##### Equation 13 of Huang et al., 1992; Lewis et al., 2009; Feldpausch et al., 2012; Mitchard et al., 2014)
startvals7 <- list(a = 54.01, b = -0.053, c = 0.0759)
Fit_nls7 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a*(1-exp(-b*D^c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq7 <- optim(par=startvals7, fn = Fit_nls7, hessian = FALSE, control = list(parscale = unlist(startvals7)))
a7 = Eq7$par[1]
b7 = Eq7$par[2]
c7 = Eq7$par[3]
PredHt <- a7*(1-exp(-b7*D^c7))
error<- (Hmean-PredHt)
RMSE_Eq7 <- rmse(error)
AIC_Eq7 <- AIC_RSS(n, length(startvals7), error)
equation <- as.formula(y ~ a*(1-exp(-b*D^c)))
res7 <- results.fun(fit = Eq7, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a7, b7, c7))
##### EQ8: Michaelis-Menten equation: H = a*D /(b+D), in page 8618 (Molto et al., 2013; Equation 3 in Huang et al. (1992))
startvals8 <- list(a = 10, b = 2) #Start points were estimated from Figure 1 in page 8632 (Molto et al., 2013)
Fit_nls8 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * I(D)/(b + I(D)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq8 <- optim(par=startvals8, fn = Fit_nls8, hessian = FALSE, control = list(parscale = unlist(startvals8)))
a8 = Eq8$par[1]
b8 = Eq8$par[2]
PredHt <- a8 * I(D)/(b8 + I(D))
error<- (Hmean-PredHt)
RMSE_Eq8 <- rmse(error)
AIC_Eq8 <- AIC_RSS(n, length(startvals8), error)
equation <- as.formula(y ~ (a * I(D)/(b + I(D))))
res8 <- results.fun(fit = Eq8, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a8, b8))
##### EQ9: Power-law equation: H=a*D^b (Eq 1 of Huang et al., 1992; Chave et al., 2005; Feldpausch et al., 2012; Ngomanda et al., 2014)
##### Start points were estimated from Ngomanda et al.(2014)
startvals9 <- list(a = 26.6049, b = -9.4854)
Fit_nls9 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- ( a* D^b )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq9 <- optim(par=startvals9, fn = Fit_nls9, hessian = FALSE, control = list(parscale = unlist(startvals9)))
a9 = Eq9$par[1]
b9 = Eq9$par[2]
PredHt <- a9* D^b9
error<- (Hmean-PredHt)
RMSE_Eq9 <- rmse(error)
AIC_Eq9 <- AIC_RSS(n, length(startvals9), error)
equation <- as.formula(y ~ a* D^b)
res9 <- results.fun(fit = Eq9, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a9, b9))
##### EQ10:Model 3: Ln(H) = a + b/D or H = a * exp(b/D), where D=DBH (Djomo et al., 2010; Equation 6 in Huang et al. (1992))
startvals10 <- list(a = 26.6049, b = -9.4854)
Fit_nls10 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (exp(a + b/D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq10 <- optim(par=startvals10, fn = Fit_nls10, hessian = FALSE, control = list(parscale = unlist(startvals10)))
a10 = Eq10$par[1]
b10 = Eq10$par[2]
PredHt <- exp(a10 + b10/D)
error<- (Hmean-PredHt)
RMSE_Eq10 <- rmse(error)
AIC_Eq10 <- AIC_RSS(n, length(startvals10), error)
equation <- as.formula(y ~ exp(a + b/D))
res10 <- results.fun(fit = Eq10, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a10, b10))
##### EQ11: H = a + b*(1-exp(-c*(D-Dmin))), where D=DBH (Fang and Bailey, 1998)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals11 <- list(a = 44.7, b =42.3, c = 0.026)
Fit_nls11 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- a + b*(1-exp(-c*(D-min(D))))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq11 <- optim(par=startvals11, fn = Fit_nls11, hessian = FALSE, control = list(parscale = unlist(startvals11)))
a11 = Eq11$par[1]
b11 = Eq11$par[2]
c11 = Eq11$par[3]
PredHt <- a11 + b11 * (1-exp(-c11 * (D - min(D))))
error<- (Hmean-PredHt)
RMSE_Eq11 <- rmse(error)
AIC_Eq11 <- AIC_RSS(n, length(startvals11), error)
equation <- as.formula(y ~ a + b*(1-exp(-c*(D-min(D)))))
res11 <- results.fun(fit = Eq11, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a11, b11, c11))
##### EQ12:H = a + (a-1.3)*b/(D+b), where D=DBH (Fang and Bailey, 1998)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals12 <- list(a = 33.62, b =20.88)
Fit_nls12 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a + (a - 1.3) * b/(D + b))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq12 <- optim(par=startvals12, fn = Fit_nls12, hessian = FALSE, control = list(parscale = unlist(startvals12)))
a12 = Eq12$par[1]
b12 = Eq12$par[2]
c12 = Eq12$par[3]
PredHt <- a12 + (a12 - 1.3) * b12/(D + b12)
error<- (Hmean-PredHt)
RMSE_Eq12 <- rmse(error)
AIC_Eq12 <- AIC_RSS(n, length(startvals12), error)
equation <- as.formula(y ~ (a + (a - 1.3) * b/(D + b)))
res12 <- results.fun(fit = Eq12, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a12, b12, c12))
##### EQ13:H = 1.3 + exp(a + b*D^c), where D=DBH (Eq 10 of Huang et al. (1992, 2000); Colbert et al., 2002)
##### Start points were for Black williows in page 173 (Colbert et al., 2002)
startvals13 <- list(a = 4.5535, b =-3.7529, c = 1)
Fit_nls13 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + exp(a + b * D^c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq13 <- optim(par=startvals13, fn = Fit_nls13, hessian = FALSE, control = list(parscale = unlist(startvals13)))
a13 = Eq13$par[1]
b13 = Eq13$par[2]
c13 = Eq13$par[3]
PredHt <- 1.3 + exp(a13 + b13 * D^c13)
error<- (Hmean-PredHt)
RMSE_Eq13 <- rmse(error)
AIC_Eq13 <- AIC_RSS(n, length(startvals13), error)
equation <- as.formula(y ~ 1.3 + exp(a + b * D^c))
res13 <- results.fun(fit = Eq13, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a13, b13, c13))
##### EQ14:H = 1.3 + a*(1-exp(-b*D))^c, where D=DBH (Equation 12 in Huang et al. (1992, 2000); Equation 2 in Sharma and Parton (2007))
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals14 <- list(a = 1.0462, b = 1.0464, c=0.9465)
Fit_nls14 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*(1-exp(-b*D))^c)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq14 <- optim(par=startvals14, fn = Fit_nls14, hessian = FALSE, control = list(parscale = unlist(startvals14)))
a14 = Eq14$par[1]
b14 = Eq14$par[2]
c14 = Eq14$par[3]
PredHt <- 1.3 + a14 *(1-exp(-b14 * D))^c14
error<- (Hmean-PredHt)
RMSE_Eq14 <- rmse(error)
AIC_Eq14 <- AIC_RSS(n, length(startvals14), error)
equation <- as.formula(y ~ 1.3 + a*(1-exp(-b*D))^c)
res14 <- results.fun(fit = Eq14, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a14, b14, c14))
##### Following equations from Table 3 in Huang et al. (1992) and Table 2 in Huang et al. (2000).
##### Parameter estimations of Balsam poplar (broad leaf) in Tables 4 and 5 (Huang et al., 1992)
##### were used as start values in the new fitting for Gabonese rainforests.
##### EQ15: Equation 4 in Huang et al.(1992 and 2000): H = 1.3 + a * (1 - exp(-bD)), where D = DBH
startvals15 <- list(a = 25.3302, b = 0.0512)
Fit_nls15 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (a * (1 - exp(-b * D)) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq15 <- optim(par=startvals15, fn = Fit_nls15, hessian = FALSE, control = list(parscale = unlist(startvals15)))
a15 = Eq15$par[1]
b15 = Eq15$par[2]
PredHt <- a15 * (1 - exp(-b15 * D))
error<- (Hmean-PredHt)
RMSE_Eq15 <- rmse(error)
AIC_Eq15 <- AIC_RSS(n, length(startvals15), error)
equation <- as.formula(y ~ a * (1 - exp(-b * D)))
res15 <- results.fun(fit = Eq15, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a15, b15))
##### EQ16:Equation 5 in Huang et al. (1992 and 2000): H = 1.3 + D^2 /(a + bD)^2, where D = DBH
startvals16 <- list(a = 1.3209, b = 0.1813)
Fit_nls16 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + D^2 /((a + b * D)^2) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq16 <- optim(par=startvals16, fn = Fit_nls16, hessian = FALSE, control = list(parscale = unlist(startvals16)))
a16 = Eq16$par[1]
b16 = Eq16$par[2]
PredHt <- 1.3 + D^2 /((a16 + b16 * D)^2)
error<- (Hmean-PredHt)
RMSE_Eq16 <- rmse(error)
AIC_Eq16 <- AIC_RSS(n, length(startvals16), error)
equation <- as.formula(y ~ 1.3 + D^2 /((a + b * D)^2))
res16 <- results.fun(fit = Eq16, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a16, b16))
##### EQ17: Equation 7 in Huang et al. (1992 and 2000): H = 1.3 + 10^a*(D^b), where D = DBH
startvals17 <- list(a = 0.4305, b = 0.5871)
Fit_nls17 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + (10^a) * (D^b) )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq17 <- optim(par=startvals17, fn = Fit_nls17, hessian = FALSE, control = list(parscale = unlist(startvals17)))
a17 = Eq17$par[1]
b17 = Eq17$par[2]
PredHt <- 1.3 + (10^a17) * (D^b17)
error<- (Hmean-PredHt)
RMSE_Eq17 <- rmse(error)
AIC_Eq17 <- AIC_RSS(n, length(startvals17), error)
equation <- as.formula(y ~ 1.3 + (10^a) * (D^b))
res17 <- results.fun(fit = Eq17, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a17, b17))
##### EQ18:Equation 8 in Huang et al. (1992 and 2000): H = 1.3 + (a*D)/(D+1) + b*D, where D = DBH
startvals18 <- list(a = 6.5487, b = 0.4507)
Fit_nls18 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a*D/(D+1) + b*D )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq18 <- optim(par=startvals18, fn = Fit_nls18, hessian = FALSE, control = list(parscale = unlist(startvals18)))
a18 = Eq18$par[1]
b18 = Eq18$par[2]
PredHt <- 1.3 + a18 * D/(D+1) + b18*D
error<- (Hmean-PredHt)
RMSE_Eq18 <- rmse(error)
AIC_Eq18 <- AIC_RSS(n, length(startvals18), error)
equation <- formula(y ~ 1.3 + a*D/(D+1) + b*D, showEnv = F)
res18 <- results.fun(fit = Eq18, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a18, b18))
##### EQ19:Equation 9 in Huang et al. (1992 and 2000): H = 1.3 + a*(D/(D+1))^b, where D = DBH
startvals19 <- list(a = 27.1752, b = 10.1979)
Fit_nls19 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a*(D/(1 + D))^b )
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq19 <- optim(par=startvals19, fn = Fit_nls19, hessian = FALSE, control = list(parscale = unlist(startvals19)))
a19 = Eq19$par[1]
b19 = Eq19$par[2]
PredHt <- 1.3 + a19*(D/(1 + D))^b19
error<- (Hmean-PredHt)
RMSE_Eq19 <- rmse(error)
AIC_Eq19 <- AIC_RSS(n, length(startvals19), error)
equation <- as.formula(y ~ 1.3 + a*(D/(1 + D))^b)
res19 <- results.fun(fit = Eq19, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a19, b19))
##### EQ20:Equation 10 in Huang et al. (1992 and 2000): H = 1.3 + exp(a + b*D^c), where D = DBH
startvals20 <- list(a = 18.804146, b = -16.760486, c=-0.025197)
Fit_nls20 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- 1.3 + exp(a + b * (D^c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq20 <- optim(par=startvals20, fn = Fit_nls20, hessian = FALSE, control = list(parscale = unlist(startvals20)))
a20 = Eq20$par[1]
b20 = Eq20$par[2]
c20 = Eq20$par[3]
PredHt <- 1.3 + exp(a20 + b20 * (D^c20))
error<- (Hmean-PredHt)
RMSE_Eq20 <- rmse(error)
AIC_Eq20 <- AIC_RSS(n, length(startvals20), error)
equation <- as.formula(y ~ 1.3 + exp(a + b * (D^c)))
res20 <- results.fun(fit = Eq20, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a20, b20, c20))
##### EQ21: Equation 11 in Huang et al. (1992 and 2000): H = 1.3 + a/(1+b*exp(-c*D)), where D = DBH
startvals21 <- list(a = 2.5241, b = 0.0012, c=0.1404)
Fit_nls21 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- 1.3 + a/(1 + b*exp(-c*D)) #Equation 11 in Huang et al. (1992)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq21 <- optim(par=startvals21, fn = Fit_nls21, hessian = FALSE, control = list(parscale = unlist(startvals21)))
a21 = Eq21$par[1]
b21 = Eq21$par[2]
c21 = Eq21$par[3]
PredHt <- 1.3 + a21 /(1 + b21 *exp(-c21 * D))
error <- (Hmean - PredHt)
RMSE_Eq21 <- rmse(error)
AIC_Eq21 <- AIC_RSS(n, length(startvals21), error)
equation <- as.formula(y ~ 1.3 + a/(1 + b*exp(-c*D)))
res21 <- results.fun(fit = Eq21, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a21, b21, c21))
##### EQ22: Equation 14 in Huang et al. (1992 and 2000): H = 1.3 + a*exp(-b*exp(-c*D)), where D = DBH (For Eqs. 12 and 13 in Huang et al. (1992 and 2002), please see the above eqs 14 and 7 in the file.)
startvals22 <- list(a = 1.6368, b = 2.1570, c=1.0951)
Fit_nls22 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*exp(-b*exp(-c*D)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq22 <- optim(par=startvals22, fn = Fit_nls22, hessian = FALSE, control = list(parscale = unlist(startvals22)))
a22 = Eq22$par[1]
b22 = Eq22$par[2]
c22 = Eq22$par[3]
PredHt <- 1.3 + a22 * exp(-b22 * exp(-c22 * D))
error <- (Hmean - PredHt)
RMSE_Eq22 <- rmse(error)
AIC_Eq22 <- AIC_RSS(n, length(startvals22), error)
equation <- as.formula(y ~ 1.3 + a*exp(-b*exp(-c*D)))
res22 <- results.fun(fit = Eq22, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a22, b22, c22))
##### EQ23: Equation 15 in Huang et al. (1992) and Eq. 16 in Huang et al. (2000): H = (1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b)
startvals23 <- list(a = 0.0464, b = 1.0716, c=27.0745)
Fit_nls23 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- ((1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq23 <- optim(par=startvals23, fn = Fit_nls23, hessian = FALSE, control = list(parscale = unlist(startvals23)))
a23 = Eq23$par[1]
b23 = Eq23$par[2]
c23 = Eq23$par[3]
PredHt <- (1.3^b23 + (c23^b23-1.3^b23)*(1-exp(-a23*D))/(1-exp(-a23*100)))^(1/b23)
error <- (Hmean - PredHt)
RMSE_Eq23 <- rmse(error)
AIC_Eq23 <- AIC_RSS(n, length(startvals23), error)
equation <- as.formula(y ~ (1.3^b + (c^b-1.3^b)*(1-exp(-a*D))/(1-exp(-a*100)))^(1/b))
res23 <- results.fun(fit = Eq23, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a23, b23, c23))
##### EQ24: Equation 16 in Huang et al. (1992) and Eq. 17 in Huang et al. (2000): H =1.3 + D^2/(a + b*D + c*D^2)
##### Start points were for Balsam poplar of Equation 12 in page 1289 (Huang et al., 1992)
startvals24 <- list(a = 0.0038, b = 0.7027, c = 0.0270)
Fit_nls24 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + D^2 /(a + b * D + c * D^2))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq24 <- optim(par=startvals24, fn = Fit_nls24, hessian = FALSE, control = list(parscale = unlist(startvals24)))
a24 = Eq24$par[1]
b24 = Eq24$par[2]
c24 = Eq24$par[3]
PredHt <- 1.3 + D^2 /(a24 + b24 * D + c24 * D^2)
error <- (Hmean - PredHt)
RMSE_Eq24 <- rmse(error)
AIC_Eq24 <- AIC_RSS(n, length(startvals24), error)
equation <- as.formula(y ~ 1.3 + D^2 /(a + b * D + c * D^2))
res24 <- results.fun(fit = Eq24, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a24, b24, c24))
##### EQ25: Equation 17 in Huang et al. (1992) and Eq. 18 in Huang et al. (2000): H = 1.3 + a *D^(b*D^-c), where D = DBH
startvals25 <- list(a = 22.9433, b = -20.9985, c=0.7680)
Fit_nls25 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a*D^(b*D^-c))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq25 <- optim(par=startvals25, fn = Fit_nls25, hessian = FALSE, control = list(parscale = unlist(startvals25)))
a25 = Eq25$par[1]
b25 = Eq25$par[2]
c25 = Eq25$par[3]
PredHt <- 1.3 + a25 * D^(b25 * D^-c25)
error <- (Hmean - PredHt)
RMSE_Eq25 <- rmse(error)
AIC_Eq25 <- AIC_RSS(n, length(startvals25), error)
equation <- as.formula(y ~ 1.3 + a*D^(b*D^-c))
res25 <- results.fun(fit = Eq25, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a25, b25, c25))
##### EQ26: Equation 18 in Huang et al. (1992) and Eq. 20 in Huang et al. (2000): H = 1.3 + a*exp(b/(D+c)), where D = DBH
startvals26 <- list(a = 3.2971, b = -0.4014, c = 5.5088)
Fit_nls26 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a * exp(b /(D+c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq26 <- optim(par=startvals26, fn = Fit_nls26, hessian = FALSE, control = list(parscale = unlist(startvals26)))
a26 = Eq26$par[1]
b26 = Eq26$par[2]
c26 = Eq26$par[3]
PredHt <- 1.3 + a26 * exp(b26 /(D + c26))
error <- (Hmean - PredHt)
RMSE_Eq26 <- rmse(error)
AIC_Eq26 <- AIC_RSS(n, length(startvals26), error)
equation <- as.formula(y ~ 1.3 + a * exp(b /(D+c)))
res26 <- results.fun(fit = Eq26, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a26, b26, c26))
##### EQ27: Equation 19 in Huang et al.(1992) and Eq.19 in Huang et al. (2000): H = 1.3 + a/(1 + b^-1 * D^-c), where D = DBH
startvals27 <- list(a = 34.4682, b = 0.0369, c = 1.0589)
Fit_nls27 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a/(1+1/(b*D^c)))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq27 <- optim(par=startvals27, fn = Fit_nls27, hessian = FALSE, control = list(parscale = unlist(startvals27)))
a27 = Eq27$par[1]
b27 = Eq27$par[2]
c27 = Eq27$par[3]
PredHt <- 1.3+a27 /(1+1/(b27 * D^c27))
error <- (Hmean - PredHt)
RMSE_Eq27 <- rmse(error)
AIC_Eq27 <- AIC_RSS(n, length(startvals27), error)
equation <- as.formula(y ~ 1.3 + a/(1+1/(b*D^c)))
res27 <- results.fun(fit = Eq27, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a27, b27, c27))
##### EQ28: Equation 20 in Huang et al. (1992) and Eq. 21 in Huang et al. (2000): H = 1.3 + a*(1 - b*exp(-c*D))^d, where D = DBH
startvals28 <- list(a = 15.2716, b = 0.9574, c = 0.0530, d = 1.1025)
Fit_nls28 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
d <- x[4]
t1 <- (1.3 + a*(1-b*exp(-c*D))^d)
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq28 <- optim(par=startvals28, fn = Fit_nls28, hessian = FALSE, control = list(parscale = unlist(startvals28)))
a28 = Eq28$par[1]
b28 = Eq28$par[2]
c28 = Eq28$par[3]
d28 = Eq28$par[4]
PredHt<- 1.3 + a28 *(1- b28 * exp(-c28 * D))^d28
error <- (Hmean - PredHt)
RMSE_Eq28 <- rmse(error)
AIC_Eq28 <- AIC_RSS(n, length(startvals28), error)
equation <- as.formula(y ~ 1.3 + a*(1-b*exp(-c*D))^d)
res28 <- results.fun(fit = Eq28, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a28, b28, c28, d28))
##### EQ29: Equation 23 in Huang et al. (2000): H = 1.3 + a*D*exp(-b*D), where D = DBH.
##### Start values were for white spruce (Table 3 in page 131 of Huang et al. (2000))
startvals29 <- list(a = 1.0264, b = 0.01157)
Fit_nls29 <- function(x) {
a <- x[1]
b <- x[2]
t1 <- (1.3 + a * D * exp(-b*D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq29 <- optim(par=startvals29, fn = Fit_nls29, hessian = FALSE, control = list(parscale = unlist(startvals29)))
a29 = Eq29$par[1]
b29 = Eq29$par[2]
PredHt <- 1.3 + a29 * D * exp(-b29 * D)
error <- (Hmean - PredHt)
RMSE_Eq29 <- rmse(error)
AIC_Eq29 <- AIC_RSS(n, length(startvals29), error)
equation <- as.formula(y ~ 1.3 + a * D * exp(-b*D))
res29 <- results.fun(fit = Eq29, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a29, b29))
##### EQ30: Equation 24 in Huang et al. (2000): H = 1.3 + a * D^b * exp(-c*D), where D = DBH. Start values were for white spruce (Table 3 in page 131 of Huang et al. (2000))
startvals30 <- list(a = 0.687824, b = 1.175135, c = 0.017699)
Fit_nls30 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (1.3 + a * D^b * exp(-c*D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq30 <- optim(par=startvals30, fn = Fit_nls30, hessian = FALSE, control = list(parscale = unlist(startvals30)))
a30 = Eq30$par[1]
b30 = Eq30$par[2]
c30 = Eq30$par[3]
PredHt <- 1.3 + a30 * D^b30 * exp(-c30*D)
error <- (Hmean - PredHt)
RMSE_Eq30 <- rmse(error)
AIC_Eq30 <- AIC_RSS(n, length(startvals30), error)
equation <- as.formula(y ~ 1.3 + a * D^b * exp(-c*D))
res30 <- results.fun(fit = Eq30, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a30, b30, c30))
##### EQ31: Pantropical DBH-H model (Chave et al., 2014. GCB 20: 3183): Ln(H) = a - E + b*Ln(D)+ c*(Ln(D))^2 , where D = DBH
startvals31 <- list(a = 0.893, b = 0.760, c = -0.034)
Fit_nls31 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a - E_Value + b * log(D) + c * (log(I(D)))^2)
t1 <- sum((Hmean-t1)^2)
return(t1)}
if(!is.na(E_Value)){
Eq31 <- optim(par=startvals31, fn = Fit_nls31, hessian = FALSE, control = list(parscale = unlist(startvals31)))
a31 = Eq31$par[1]
b31 = Eq31$par[2]
c31 = Eq31$par[3]
PredHt <- a31 - E_Value + b31 * log(D) + c31 * (log(I(D)))^2
error <- (Hmean - PredHt)
RMSE_Eq31 <- rmse(error)
AIC_Eq31 <- AIC_RSS(n, length(startvals31), error)
equation <- as.formula(y ~ a - E_Value + b * log(D) + c * (log(I(D)))^2)
res31 <- results.fun(fit = Eq31, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a31, b31, c31))
}else{AIC_Eq31 <- 9999999}
##### EQ32: The Mitscherlisch model (Ngomanda et al., 2014): H = a - b * exp(-cD), where D = DBH
startvals32 <- list(a = 44.7, b = 42.3, c = 0.026)
Fit_nls32 <- function(x) {
a <- x[1]
b <- x[2]
c <- x[3]
t1 <- (a - b * exp(-c * D))
t1 <- sum((Hmean-t1)^2)
return(t1)}
Eq32 <- optim(par=startvals32, fn = Fit_nls32, hessian = FALSE, control = list(parscale = unlist(startvals32)))
a32 = Eq32$par[1]
b32 = Eq32$par[2]
c32 = Eq32$par[3]
PredHt <- a32 - b32 * exp(-c32 * D)
error <- (Hmean - PredHt)
RMSE_Eq32 <- rmse(error)
AIC_Eq32 <- AIC_RSS(n, length(startvals32), error)
equation <- as.formula(y ~ a - b * exp(-c * D))
res32 <- results.fun(fit = Eq32, equation, Hmean = Hmean, D = D, error = error, PredHt = PredHt, coef = c(a32, b32, c32))
##### Model comparison
AIC_best.fit <- paste("Eq", 1, sep="")
AIC_min_temp = AIC_Eq1
AIC.df <- matrix(nrow = 32, ncol = 2, dimnames = list(c(), c("Equ", "AIC")))
for (k in 1: 32){
AIC_min <- eval(parse(text=paste("AIC_Eq", k, sep="")))
AIC.df[k,] <- c(Equ = paste("Eq", k, sep=""), AIC = round(AIC_min, 1))
AIC_min <- min(AIC_min, AIC_min_temp)
if (! is.nan(AIC_min) && !is.nan(AIC_min_temp) ) {
if ( AIC_min < AIC_min_temp) {AIC_best.fit <- paste("Eq", k, sep="")}
}
AIC_min_temp = AIC_min
}
AIC.df <- data.frame(AIC.df)
AIC.df$AIC <- as.numeric(as.character(AIC.df$AIC))
best.equ <- eval(parse(text=AIC_best.fit))
rnum <- paste("res", gsub("[^[:digit:]]", "", AIC_best.fit), sep="")
res.fin <- get(rnum)
# pltname <- gsub('([[:punct:]])|\\s+','_',Plot_name)
# data_labels <- data.frame(time = c(20), value = c(0.95*max(Hmean)), label = pltname)
## Disable plotting
# gp <- plot.fun(res = res.fin, hdat = hdat)
# print(gp)
# pdf(file = paste("HDgraph", ".pdf", sep=""), width = 4, height = 4)
# print(gp)
# dev.off()
predht.fun <- function(res, sub.dat){
a <- res$vars["a"]
b <- res$vars["b"]
c <- res$vars["c"]
d <- res$vars["d"]
equ.fun <- function(D) eval(parse(text = res$equ))
ph <- equ.fun(sub.dat$D)
}
predht <- predht.fun(res = res.fin, sub.dat)
full.list <- res.fin
full.list$predht <- predht
full.list$Xequ <- as.character(res.fin$equation)[3]
full.list
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.