---
title: "Ajuste de equações lineares e não lineares por grupo"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ajuste de equações lineares e não lineares por grupo}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
  
```{r, echo = FALSE, message = FALSE, warning=FALSE}
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.
```{r}
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`:

```{r}
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:
```{r}
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`:

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


```{r}
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:
```{r}
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
```

```{r}
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`:

```{r}
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`:
```{r}
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`:
```{r, warning=FALSE, message=FALSE}
resid_plot(dados_est, "dh", "Schumacher", "Chapman-Richards", "Bailey-Clutter", "Curtis")
```

Podemos utlizar outros tipos de gráficos, como histogramas:
```{r, warning=FALSE, message=FALSE}
resid_plot(dados_est, "dh", "Schumacher","Chapman-Richards","Bailey-Clutter", "Curtis",
           type = "histogram_curve")
```

E gráfico de versus:
```{r, warning=FALSE, message=FALSE}
resid_plot(dados_est, "dh", "Schumacher", "Chapman-Richards", "Bailey-Clutter", "Curtis",
           type = "versus")
```