Fitting taper models
The timbeR package has functions to estimate diameters along the stem, height at which certain diameter values occur and total or partial volumes. For this, it is necessary to fit a taper model that describes the stem profile. This vignette aims to exemplify the regression analysis needed to fit the three models whose the functions mentioned above are implemented in the timbeR package. The possible models to be used are:
5th degree polynomial taper function (Schöepfer, 1966) \[\frac{d_i}{dbh}=\beta_0\frac{h_i}{h}+\beta_1\left(\frac{h_i}{h}\right)^2+\beta_2\left(\frac{h_i}{h}\right)^3+\beta_3\left(\frac{h_i}{h}\right)^4+\beta_4\left(\frac{h_i}{h}\right)^5\]
Kozak (2004) variable-form taper model \[d_i = \beta_0dbh^{\beta_1}\left[\frac{1-\left(\frac{h_i}{ht}\right)^{1/4}}{1-p^{1/4}}\right]^{\beta_2+\beta_3\left(\frac{1}{e^{dbh/ht}}\right)+\beta_4dbh^{\left[\frac{1-\left(\frac{h_i}{ht}\right)^{1/4}}{1-p^{1/4}}\right]}+\beta_5\left[\frac{1-\left(\frac{h_i}{ht}\right)^{1/4}}{1-p^{1/4}}\right]^{dbh/ht}}\]
Bi (2000) trigonometric variable-form taper model \[d_i=dbh\left[ \left( \frac{log\;sin \left( \frac{\pi}{2} \frac{h_i}{ht} \right)} {log\;sin \left( \frac{\pi}{2} \frac{1,3}{ht} \right)} \right) ^{\beta_0+\beta_1sin\left(\frac{\pi}{2}\frac{h_i}{ht}\right)+\beta_2cos\left(\frac{3\pi}{2}\frac{h_i}{ht}\right)+\beta_3sin\left(\frac{\pi}{2}\frac{h_i}{ht}\right)/\frac{h_i}{ht}+\beta_4dbh+\beta_5\frac{h_i}{ht}\sqrt{dbh}+\beta_6\frac{h_i}{ht}\sqrt{ht}} \right]\]
where:
\(\beta_1,\beta_2,...,\beta_n\) = model parameters;
\(h_i\) = height to section i of the stem;
\(d_i\) = diameter in section i of the stem;
\(dbh\) = diameter at breast height;
\(h\) = total height of the tree;
\(p\) = first inflection point calculated in the segmented model of Max and Burkhart (1976).
We will perform a regression analysis on the tree_scaling dataset, using the aforementioned models. The data can be accessed by importing the timbeR package.
# install.packages("devtools")
# options(download.file.method = "libcurl")
# devtools::install_github('sergiocostafh/timbeR')
library(dplyr)
library(timbeR)
glimpse(tree_scaling)
#> Rows: 136
#> Columns: 5
#> $ tree_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,~
#> $ dbh <dbl> 14.96, 14.96, 14.96, 14.96, 14.96, 14.96, 14.96, 14.96, 14.96,~
#> $ h <dbl> 15.61, 15.61, 15.61, 15.61, 15.61, 15.61, 15.61, 15.61, 15.61,~
#> $ hi <dbl> 0.1561, 0.3122, 0.4683, 0.6244, 0.7805, 1.3000, 1.5610, 2.3415~
#> $ di <dbl> 16.55, 15.92, 15.60, 15.18, 14.96, 14.00, 12.73, 12.10, 10.19,~As we can see, there are five columns in the dataset that refer to the tree id (tree_id), diameter at breast height (dbh), tree total height (h), height at section i (hi) and diameter at hi height (di).
A common way to visualize the stem profile from collected data is to plot the relationship between relative diameters and relative heights (di / dbh vs hi / ht), as follows.
library(ggplot2)
tree_scaling <- tree_scaling %>%
mutate(did = di/dbh,
hih = hi/h)
ggplot(tree_scaling, aes(x = hih, y = did, group = tree_id))+
geom_point()+
labs(x = 'hi / h',
y = 'di / dbh')
Now that we understand the dataset, we can start the regression analysis. The first model we will fit is the 5th degree polynomial.
poli5 <- lm(did~hih+I(hih^2)+I(hih^3)+I(hih^4)+I(hih^5),tree_scaling)
summary(poli5)
#>
#> Call:
#> lm(formula = did ~ hih + I(hih^2) + I(hih^3) + I(hih^4) + I(hih^5),
#> data = tree_scaling)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.124049 -0.029700 -0.003642 0.032621 0.112321
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.15657 0.01615 71.596 < 2e-16 ***
#> hih -3.37185 0.46898 -7.190 4.55e-11 ***
#> I(hih^2) 13.57792 3.40969 3.982 0.000113 ***
#> I(hih^3) -29.92092 9.59285 -3.119 0.002235 **
#> I(hih^4) 29.39935 11.38070 2.583 0.010893 *
#> I(hih^5) -10.85327 4.78706 -2.267 0.025028 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.05469 on 130 degrees of freedom
#> Multiple R-squared: 0.9697, Adjusted R-squared: 0.9685
#> F-statistic: 830.8 on 5 and 130 DF, p-value: < 2.2e-16
tree_scaling <- tree_scaling %>%
mutate(di_poli = predict(poli5)*dbh)
poli_rmse <- tree_scaling %>%
summarise(RMSE = sqrt(sum((di_poli-di)^2)/mean(di_poli))) %>%
pull(RMSE) %>%
round(2)
ggplot(tree_scaling,aes(x=hih))+
geom_point(aes(y = (di_poli-di)/di_poli*100))+
geom_hline(aes(yintercept = 0))+
scale_y_continuous(limits=c(-60,60), breaks = seq(-100,100,20))+
scale_x_continuous(limits=c(0,1))+
labs(x = 'hi / h', y = 'Residuals (%)',
title = '5th degree polynomial taper function (Schöepfer, 1966)',
subtitle = 'Dispersion of residuals along the stem',
caption = paste0('Root Mean Squared Error = ', poli_rmse,'%'))+
theme(plot.title.position = 'plot')
The 5th degree polynomial is a fixed-form taper function that represents the average shape of the stem profiles used to fit the model. For this dataset, the Root Mean Square Error of this model was 3.01% and we can see that the residues are heteroskedastic.
Let’s see if we can do better with the Bi model.Due to its non-linear nature, we will use the nlsLM function from the minpack.lm package to estimate the model parameters.
library(minpack.lm)
bi <- nlsLM(di ~ taper_bi(dbh, h, hih, b0, b1, b2, b3, b4, b5, b6),
data=tree_scaling,
start=list(b0=1.8,b1=-0.2,b2=-0.04,b3=-0.9,b4=-0.0006,b5=0.07,b6=-.14))
summary(bi)
#>
#> Formula: di ~ taper_bi(dbh, h, hih, b0, b1, b2, b3, b4, b5, b6)
#>
#> Parameters:
#> Estimate Std. Error t value Pr(>|t|)
#> b0 1.821275 0.455063 4.002 0.000105 ***
#> b1 1.417831 0.287019 4.940 2.38e-06 ***
#> b2 0.142028 0.035637 3.985 0.000112 ***
#> b3 -1.082588 0.301635 -3.589 0.000470 ***
#> b4 -0.004008 0.003088 -1.298 0.196525
#> b5 0.325254 0.053354 6.096 1.17e-08 ***
#> b6 -0.735532 0.106476 -6.908 2.02e-10 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8063 on 129 degrees of freedom
#>
#> Number of iterations to convergence: 5
#> Achieved convergence tolerance: 1.49e-08
tree_scaling <- tree_scaling %>%
mutate(di_bi = predict(bi))
bi_rmse <- tree_scaling %>%
summarise(RMSE = sqrt(sum((di_bi-di)^2)/mean(di_bi))) %>%
pull(RMSE) %>%
round(2)
ggplot(tree_scaling,aes(x=hih))+
geom_point(aes(y = (di_bi-di)/di_bi*100))+
geom_hline(aes(yintercept = 0))+
scale_y_continuous(limits=c(-60,60), breaks = seq(-100,100,20))+
scale_x_continuous(limits=c(0,1))+
labs(x = 'hi / h', y = 'Residuals (%)',
title = 'Bi (2000) trigonometric variable-form taper function',
subtitle = 'Dispersion of residuals along the stem',
caption = paste0('Root Mean Squared Error = ', bi_rmse,'%'))+
theme(plot.title.position = 'plot')
The Bi model performed better than the polynomial function, based on the RMSE value. However, we still have heteroscedasticity in the residues. Let’s see what we get by adjusting the Kozak (2004) model. We will treat the p parameter of this model as one more to be estimated using the nlsLM function.
kozak <- nlsLM(di ~ taper_kozak(dbh, h, hih, b0, b1, b2, b3, b4, b5, b6, b7, b8, p),
start=list(b0=1.00,b1=.97,b2=.03,b3=.49,b4=-
0.87,b5=0.50,b6=3.88,b7=0.03,b8=-0.19, p = .1),
data = tree_scaling,
control = nls.lm.control(maxiter = 1000, maxfev = 2000)
)
summary(kozak)
#>
#> Formula: di ~ taper_kozak(dbh, h, hih, b0, b1, b2, b3, b4, b5, b6, b7,
#> b8, p)
#>
#> Parameters:
#> Estimate Std. Error t value Pr(>|t|)
#> b0 1.057e+00 5.338e-01 1.981 0.0498 *
#> b1 1.088e+00 8.074e-02 13.471 < 2e-16 ***
#> b2 1.120e-01 2.205e-01 0.508 0.6125
#> b3 3.404e-01 5.610e-02 6.067 1.41e-08 ***
#> b4 -2.823e+00 4.948e-01 -5.706 7.83e-08 ***
#> b5 9.832e-01 1.192e-01 8.251 1.79e-13 ***
#> b6 1.197e+01 2.902e+00 4.123 6.73e-05 ***
#> b7 1.082e-01 4.673e-02 2.316 0.0222 *
#> b8 -4.930e-01 2.951e-01 -1.671 0.0972 .
#> p 2.149e-18 6.523e-13 0.000 1.0000
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.6276 on 126 degrees of freedom
#>
#> Number of iterations to convergence: 100
#> Achieved convergence tolerance: 1.49e-08
tree_scaling <- tree_scaling %>%
mutate(di_kozak = predict(kozak))
kozak_rmse <- tree_scaling %>%
summarise(RMSE = sqrt(sum((di_kozak-di)^2)/mean(di_kozak))) %>%
pull(RMSE) %>%
round(2)
ggplot(tree_scaling, aes(x=hih))+
geom_point(aes(y = (di_kozak-di)/di_kozak*100))+
geom_hline(aes(yintercept = 0))+
scale_y_continuous(limits=c(-100,100), breaks = seq(-100,100,20))+
scale_x_continuous(limits=c(0,1))+
labs(x = 'hi / h', y = 'Residuals (%)',
title = 'Kozak (2004) variable-form taper function',
subtitle = 'Dispersion of residuals along the stem',
caption = paste0('Root Mean Squared Error = ', kozak_rmse,'%'))+
theme(plot.title.position = 'plot')
By fitting the Kozak (2004) model, we obtained a lower RMSE and also managed to homogenize the dispersion of the residues.
