Ajuste de equações lineares e não lineares por grupo"

In forestmangr: Forest Mensuration and Management

knitr::opts_chunk$set(collapse = T, comment = "#>")
knitr::opts_chunk$set(fig.width=7, fig.height=5)
options(tibble.print_min = 6L, tibble.print_max = 6L)
library(forestmangr)
library(dplyr)
library(tidyr)

Vamos ajustar alguns modelos para estimar a altura dominante, lineares e não-lineares, e compará-los em seguida. Vamos utilizar os primeiros 10 talhões do dado de exmplo exfm16.

library(forestmangr)
library(dplyr)
library(tidyr)

data(exfm14)
dados <- exfm14 %>% filter(strata%in%1:10)
dados

Para ajustar o modelo de Schumacher para altura dominante em função da idade, podemos utilizar lm_table. Observe que, não há a necessidade de criar novas variáveis para realizar o ajuste, graças às funções log e inv:

mod1 <- lm_table(dados, log(dh) ~ inv(age))
mod1

Para ajustar um modelo não linear, como o de Chapman-Richards, podemos utilizar a função nls_table. Esta utiliza o algoritimo Levenberg-Marquardt por padrão, garantindo um ótimo ajuste. Por ser um ajuste não linear, valores iniciais para os coeficientes devem ser informados:

mod2 <- nls_table(dados, dh ~ b0 * (1 - exp( -b1 * age )  )^b2, 
          mod_start = c( b0=23, b1=0.03, b2 = 1.3  ) )
mod2

Se quisermos fazer esses ajustes por estrato, basta utilizar o argumento .groups:

mod1 <- lm_table(dados, log(dh) ~ inv(age), .groups = "strata")
mod1

mod2 <- nls_table(dados, dh ~ b0 * (1 - exp( -b1 * age )  )^b2, 
          mod_start = c( b0=23, b1=0.03, b2 = 1.3  ),
          .groups = "strata" )
mod2

Se o ajuste não ficar adequado, é possível utilizar um dataframe com chutes específicos para cada grupo, e utilizá-lo como start:

tab_start <- data.frame(strata = c(1:10), 
              rbind(
              data.frame(b0=rep(23, 5),b1=rep(0.03,5),b2=rep(1.3,5) ), 
              data.frame(b0=rep(23, 5),b1=rep(0.03,5),b2=rep(.5,5) )))
tab_start

mod2 <- nls_table(dados, dh ~ b0 * (1 - exp( -b1 * age )  )^b2, 
          mod_start = tab_start,
          .groups = "strata" )
mod2

Agora vamos ajustar mais alguns modelos. Os modelos utilizados serão:

Schumacher: $$ Ln(DH) = \beta_0 + \beta_1 * \frac{1}{age} $$

Chapman-Richards: $$ DH = \beta_0 * (1 - exp^{-\beta_1 * age})^{\beta_2} $$

Bayley-Clutter: $$ Ln(DH) = \beta_0 + \beta_1 * \begin{pmatrix} \frac{1}{age} \end{pmatrix} ^{\beta_2} $$

Curtis: $$ DH = \beta_0 + \beta_1 * \frac{1}{age} $$

E iremos adicionar os valores estimados aos dados originais, utilizando o output merge_est, nomeando as variáveis estimadas utilizando est.name:

dados_est <- dados %>% 
  lm_table(log(dh) ~ inv(age), .groups = "strata",
           output = "merge_est", est.name = "Schumacher") %>% 

  nls_table(dh ~ b0 * (1 - exp( -b1 * age )  )^b2, 
          mod_start = c( b0=23, b1=0.03, b2 = 1.3  ),.groups="strata",
          output ="merge_est",est.name="Chapman-Richards") %>% 

  nls_table(log(dh) ~ b0 + b1 * ( inv(age)^b2 ) , 
          mod_start = c( b0=3, b1=-130, b2 = 1.5),.groups = "strata",
          output ="merge_est",est.name = "Bailey-Clutter") %>% 

  lm_table(dh ~ inv(age), .groups = "strata",
           output = "merge_est", est.name = "Curtis") 

head(dados_est)  

Obs: as funções lm_table e nls_table verificam se o modelo possui log na variável y, e caso possua, ele o retira automaticamente. Por isso, não há a necessidade de calcular a exponencial dos dados estimados.

Para comparar os modelos, podemos calcular a raiz qudrada do erro médio, e o bias de todos os modelos. Para isso, vamos unir as variáveis estimadas em uma única coluna com tidyr::gather, agrupar por variável, e utilizar as funções rmse_per e bias_per:

dados_est %>% 
  gather(Model, Value, 
         Schumacher, `Chapman-Richards`, `Bailey-Clutter`, Curtis) %>% 
  group_by(Model) %>% 
  summarise(
    RMSE = rmse_per(y = dh, yhat = Value),
    BIAS = bias_per(y = dh, yhat = Value) )

Outra forma de avaliar estes modelos é utilizando gráficos de resíduos. Para isso, podemos utilizar a função resid_plot:

resid_plot(dados_est, "dh", "Schumacher", "Chapman-Richards", "Bailey-Clutter", "Curtis")

Podemos utlizar outros tipos de gráficos, como histogramas:

resid_plot(dados_est, "dh", "Schumacher","Chapman-Richards","Bailey-Clutter", "Curtis",
           type = "histogram_curve")

E gráfico de versus:

resid_plot(dados_est, "dh", "Schumacher", "Chapman-Richards", "Bailey-Clutter", "Curtis",
           type = "versus")

forestmangr documentation built on Nov. 24, 2023, 1:07 a.m.