Dinâmica florestal com R

Géraldine Derroire & Sylvain Schmitt

Cirad

2025-11-07

Pacotes R

Usaremos os pacotes tidyverse e GGally para manipular tabelas e criar gráficos, vegan para fazer análises da vegetação e brms e tidybayes para fazer modelização bayesiana.

library(tidyverse)
library(GGally)
library(vegan)
library(brms)
library(tidybayes)

Dados

Dados de parcelas florestais permanentes

  • Dados de uma parcela permanente (parcela 6) na estação de pesquisa florestal de Paracou na Guiana Francesa

  • Todas as árvores com um diâmetro à altura do peito (DAP) ≥ 10 cm são mapeadas, identificadas botanicamente e medidas em cada censo

Dados de parcelas florestais permanentes

Começamos carregando os dados na forma de objetos no formato Rdata.

dt_P6 <- read_csv(file = "data/2024-08-29_ParacouP6AllYears.csv")
dt_P6 <- dt_P6 %>% 
  select(Forest, Plot, PlotArea, SubPlot, 
         idTree, Xutm, Yutm, 
         Family, Genus, Species,
         CensusYear, CodeAlive, CircCorr)

⚠️ Esses dados NÃO são de acesso livre, por favor, não os distribua nem utilize para qualquer outro fim.

Dados de parcelas florestais permanentes

Esses dados contêm:

  • Forest: Nome da floresta

  • Plot: Número da parcela

  • PlotArea: Área da parcela (hectares)

  • SubPlot: Número da subparcela

  • idTree: Identificador único da árvore no banco de dados

  • Xutm e Yutm: Coordenadas da árvore (UTM 22N (EPSG: 32 622))

  • Family

  • Genus

  • Species

  • CensusYear: Ano da medição

  • CodeAlive: A árvore está viva (1) ou morta (0)?

  • CircCorr: Circunferência (cm) da árvore no ponto de medição

Dados de características funcionais

Vamos usar os dados das características funcionais das árvores da Guiana Francesa obtidos na base de dados TRY e compilados no pacote rcontroll.

dt_traits <- read_tsv("https://raw.githubusercontent.com/sylvainschmitt/rcontroll/refs/heads/main/inst/extdata/TROLLv3_species.txt") %>% 
  select(-s_seedmass, -s_regionalfreq, -s_drymass) %>% 
  rename(species = s_name, LMA = s_LMA, Nmass = s_Nmass,
         Pmass = s_Pmass, wsg = s_wsg, dbhmax = s_dbhmax, 
         hmax = s_hmax, ah = s_ah, tlp = s_tlp)
dt_traits
# A tibble: 45 × 9
   species                     LMA  Nmass   Pmass   wsg dbhmax  hmax    ah   tlp
   <chr>                     <dbl>  <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
 1 Abarema_jupunba           110.  0.0232 5.7 e-4 0.576   0.66  50.9 0.308 -1.60
 2 Bocoa_prouacensis          95.6 0.0255 5.05e-4 0.798   0.4   47.9 0.251 -2.08
 3 Brosimum_rubescens         85.8 0.0188 6.88e-4 0.619   0.52  50.9 0.308 -2.28
 4 Carapa_procera            130.  0.0177 6.36e-4 0.55    0.49  41.4 0.187 -1.93
 5 Caryocar_glabrum           84.2 0.0217 8.13e-4 0.612   1.15  48.6 0.275 -1.88
 6 Cecropia_obtusa            69.5 0.0276 1.38e-3 0.388   0.35  50.9 0.308 -1.43
 7 Chrysophyllum_sanguinole… 151.  0.0154 4.69e-4 0.566   0.67  50.9 0.308 -1.45
 8 Conceveiba_guianensis      86.7 0.0235 7.29e-4 0.537   0.23  33.3 0.16  -1.70
 9 Cordia_sagotii            134.  0.0253 5.92e-4 0.443   0.47  36.4 0.139 -1.91
10 Couepia_bracteosa         123.  0.015  4.82e-4 0.784   0.45  50.9 0.308 -2.26
# ℹ 35 more rows

Dados de características funcionais

Esses dados contêm os valores de 11 características funcionais para 45 espécies. Vamos focar nos seguintes características:

  • massa foliar por área (LMA, em \(g~m^{-2}\))

  • massa de nitrogênio foliar (Nmass, em \(g~g^{-1}\))

  • massa de fósforo foliar (Pmass, em \(g~g^{-1}\))

  • gravidade específica da madeira (wsg, em \(g~cm^{-3}\))

  • potencial hídrico foliar no ponto de perda de turgor (tlp, em \(MPa\))

dt_traits <- dt_traits %>% 
  select(species, LMA, Nmass, Pmass, wsg, tlp)

dt_traits %>% slice_head(n = 5)
# A tibble: 5 × 6
  species              LMA  Nmass    Pmass   wsg   tlp
  <chr>              <dbl>  <dbl>    <dbl> <dbl> <dbl>
1 Abarema_jupunba    110.  0.0232 0.00057  0.576 -1.60
2 Bocoa_prouacensis   95.6 0.0255 0.000505 0.798 -2.08
3 Brosimum_rubescens  85.8 0.0188 0.000688 0.619 -2.28
4 Carapa_procera     130.  0.0177 0.000636 0.55  -1.93
5 Caryocar_glabrum    84.2 0.0217 0.000813 0.612 -1.88

Exploração de dados de características funcionais

Distribuição das valores das características

Podemos começar examinando a distribuição de valores para cada um dos recursos com um histograma usando o pacote ggplot (parte de tidyverse).

dt_traits %>% 
  pivot_longer(LMA : tlp, 
               names_to = "traits", values_to = "values") %>% 
  ggplot(aes(x = values)) + geom_histogram() +
  facet_wrap(~ traits, scales = "free") + theme_bw()

Coeficientes de variação das características

Podemos ver que as características são mais ou menos variáveis, mas para quantificar isso adequadamente, podemos usar o coeficiente de variação definido como:

\[CV = \frac{desvio\ padrão}{média}\]

dt_traits %>% summarise(
  across(LMA : tlp,
         function(v) {round(sd(v)/abs(mean(v)), 3)}))
# A tibble: 1 × 5
    LMA Nmass Pmass   wsg   tlp
  <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.295 0.253 0.455 0.161 0.127

Correlações entre pares de características

As características não variam de forma independente, portanto, podemos explorar as relações entre cada par de características usando os coeficientes de correlação.

ggpairs(dt_traits, columns = names(dt_traits)[-1])

Axes de covariância de características funcionais

Para entender melhor a organização geral de todos os recursos em “estratégias funcionais”, podemos construir uma análise de componentes principais (ACP) com o pacote vegan.

trait_pca <- dt_traits %>% select(-species) %>% 
  pca(scale = TRUE)
summary(trait_pca)

Call:
pca(X = ., scale = TRUE) 

Partitioning of correlations:
              Inertia Proportion
Total               5          1
Unconstrained       5          1

Eigenvalues, and their contribution to the correlations 

Importance of components:
                         PC1    PC2    PC3     PC4     PC5
Eigenvalue            2.2584 1.3822 0.5916 0.43239 0.33541
Proportion Explained  0.4517 0.2764 0.1183 0.08648 0.06708
Cumulative Proportion 0.4517 0.7281 0.8464 0.93292 1.00000

Axes de covariância de características funcionais

Podemos representar a PCA em um biplot de eixos 1-2 com as características como variáveis representadas por setas e as espécies por pontos posicionados de acordo com os valores de suas características.

# display as características
biplot(trait_pca, scaling = 2, 
       display = "species") 
# display as espécies
points(trait_pca, scaling = 2, 
       display = "sites", type ="points") 

  • Axe 1: espectro econômico das folhas: estratégias de aquisição de recursos conservadoras (lentas) versus aquisitivas (rápidas) Wright et al 2004

  • As características da madeira Baraloto et al 2010 e as características relacionadas com a água Maréchaux et al 2019 são independentes do espectro econômico das folhas.

Exploração dos dados de crescimento

Preparação dos dados

Para usar os dados de crescimento, começamos filtrando os dados da parcela 6.

  • Manter apenas as árvores totalmente determinadas
dt_P6 <- dt_P6 %>% 
  filter(Species != "Indet.")
  • Calcular o DAP
dt_P6 <- dt_P6 %>% 
  mutate(DBH = CircCorr / pi)
  • Manter apenas a observação de árvores vivas com \(DAP ≥ 10 cm\)
dt_P6 <- dt_P6 %>% 
  filter(DBH >= 10 & CodeAlive ==1) 
  • Adicionar os nomes completos das espécies
dt_P6 <- dt_P6 %>% 
  unite(GenSp,
        Genus, Species,
        sep = "_", 
        remove = FALSE)

Número de indivíduos

  • Número total de indivíduos
n_distinct(dt_P6$idTree)
[1] 4262
  • Número de indivíduos por censo e duração do monitoramento
dt_P6 %>% count(CensusYear)  %>% 
  ggplot(aes(x = CensusYear, y = n)) + 
  geom_point() + geom_line() + 
  ylim(0, 4000) +
  labs(x = "Ano da medição", y = "Número de indivíduos") +
  theme_bw()

Cálculo do crescimento anual

dt_P6 <- dt_P6 %>% 
  group_by(idTree) %>% 
  arrange(CensusYear) %>% 
  mutate(Growth = (DBH - lag(DBH)) / (CensusYear - lag(CensusYear))) %>% 
  filter(Growth >= 0) %>% 
  ungroup()

dt_P6 %>% select(idTree, CensusYear, DBH, Growth) %>% 
  filter(idTree == 100621) %>% 
  head(n = 15)
# A tibble: 15 × 4
   idTree CensusYear   DBH Growth
    <dbl>      <dbl> <dbl>  <dbl>
 1 100621       1985  12.7 0     
 2 100621       1986  12.7 0     
 3 100621       1987  12.9 0.159 
 4 100621       1988  12.9 0     
 5 100621       1989  13.2 0.318 
 6 100621       1991  12.9 0     
 7 100621       1992  12.9 0     
 8 100621       1993  13.1 0.159 
 9 100621       1994  13.2 0.159 
10 100621       1997  13.1 0     
11 100621       1999  13.2 0.0796
12 100621       2001  13.2 0     
13 100621       2003  13.2 0     
14 100621       2005  13.2 0     
15 100621       2007  13.2 0     

Trajetórias de crescimento

plot_traj_DBH <- function(id, dt) {
    sp <- dt %>% filter(idTree == id) %>% distinct(GenSp) # especie
    n <- dt %>% filter(idTree == id) %>% count() # medições
    filter(dt, idTree == id) %>% 
      ggplot(aes(x=CensusYear, y = DBH)) + 
      geom_point() + 
      geom_smooth(method = "loess", se = FALSE) + 
      labs(x = "Ano da medição",
           y = "Diâmetro à altura do peito (em cm)",
           title = paste0("Trajetória do diâmetro - Indivíduo ", id, 
                          " - ", sp, " (", n, " medições)")) + theme_bw() 
}
plot_traj_DBH(id = "103768", dt = dt_P6)

Trajetórias de crescimento

plot_traj_Growth <- function(id, dt) {
    sp <- dt %>% filter(idTree == id) %>% distinct(GenSp) # especie
    n <- dt %>% filter(idTree == id) %>% count() # medições
    filter(dt, idTree == id) %>% 
      ggplot(aes(x=DBH, y = Growth)) + 
      geom_point() + 
      geom_smooth(method = "loess", se = FALSE) + 
      labs(x = "Ano da medição",
           y = "Crescimento anual (em cm/ano)",
           title = paste0("Trajetória de crecimento - Indivíduo ", id, 
                          " - ", sp, " (", n, " medições)")) + theme_bw() 
}
plot_traj_Growth(id = "103768", dt = dt_P6)

Trajetórias de crescimento

plot_traj_sp <- function(sp, dt) {
    n <- dt %>% filter(GenSp == sp) %>% count() # medições
    filter(dt, GenSp == sp) %>% 
      ggplot(aes(x=DBH, y = Growth)) + 
      geom_point() + 
      geom_smooth(method = "loess", se = FALSE) + 
      labs(x = "Ano da medição",
           y = "Crescimento anual (em cm/ano)",
           title = paste0("Trajetória de crecimento - ",
                          sp, " (", n, " medições)")) + theme_bw() 
}
plot_traj_sp(sp = "Cecropia_obtusa", dt = dt_P6)

Modelagem do crescimento

Modelo de Canham

\[log(AGR+1)=G \times exp[-0.5 \times( \frac{log(DBH/D)}{K}^2)]\]

data_frame(DBH = 10:100) %>% mutate(log_agr = 0.5*exp(-0.5*(log(DBH/50)^2/0.2))) %>% 
  ggplot(aes(DBH, exp(log_agr)-1)) + geom_line() + theme_bw() +
  labs(x = "Ano da medição",
       y = "Crescimento anual") + theme_bw() +
  geom_vline(xintercept = 50, colour = "red", linetype = "dashed") + 
  geom_hline(yintercept = exp(0.5)-1, colour = "red", linetype = "dashed")

Uma inferência com Cecropia obtusa

fit <- brm(bf(log(Growth+1) ~ G*exp(-0.5*(log(DBH/D)^2/K)),
              G ~ 1, D ~ 1, K ~ 1,
              nl = TRUE),
           prior = c(
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "G"),
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "K"),
                 prior("uniform(1,30)", lb = 1, ub = 30, nlpar = "D")
          ),
           data = filter(dt_P6, GenSp  == "Cecropia_obtusa"),
           cores = 4)
save(fit, file = "data/fit_cecro.Rdata")
load(file = "data/fit_cecro.Rdata") ; fit
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: log(Growth + 1) ~ G * exp(-0.5 * (log(DBH/D)^2/K)) 
         G ~ 1
         D ~ 1
         K ~ 1
   Data: filter(dt_P6, GenSp == "Cecropia_obtusa") (Number of observations: 94) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
G_Intercept     0.44      0.09     0.31     0.65 1.05      114      373
D_Intercept    18.81      4.60     8.33    28.89 1.01     1127      793
K_Intercept     0.94      0.96     0.04     2.87 1.05      107      602

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.33      0.03     0.29     0.39 1.01      585      760

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Uma inferência com Cecropia obtusa

plot_fit <- function(fit, Species) {
  G <- posterior_summary(fit)[1, "Estimate"]
  D <- posterior_summary(fit)[2, "Estimate"]
  K <- posterior_summary(fit)[3, "Estimate"]
  filter(dt_P6, GenSp  == !!Species) %>% 
    ggplot(aes(x=DBH, y = Growth)) + geom_point() + 
    stat_function(fun = function(x) exp(G*exp(-0.5*(log(x/D)^2/K))) - 1,
                  colour = "red", size = 2) +
        labs(x = "Ano da medição",
             y = "Crescimento anual") + theme_bw() +
    geom_vline(xintercept = D, colour = "red", linetype = "dashed") + 
    geom_hline(yintercept = exp(G)-1, colour = "red", linetype = "dashed") +
    ylim(0, NA) + xlim(10, NA)
}

Uma inferência com Cecropia obtusa

plot_fit(fit = fit, Species = "Cecropia_obtusa")

Uma inferência com Dicorynia guianensis

fit <- brm(bf(log(Growth+1) ~ G*exp(-0.5*(log(DBH/D)^2/K)),
              G ~ 1, D ~ 1, K ~ 1,
              nl = TRUE),
           prior = c(
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "G"),
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "K"),
                 prior("uniform(1,300)", lb = 1, ub = 300, nlpar = "D")
          ),
           data = filter(dt_P6, GenSp  == "Dicorynia_guianensis"),
           cores = 4)
save(fit, file = "data/fit_dico.Rdata")
load(file = "data/fit_dico.Rdata") ; fit
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: log(Growth + 1) ~ G * exp(-0.5 * (log(DBH/D)^2/K)) 
         G ~ 1
         D ~ 1
         K ~ 1
   Data: filter(dt_P6, GenSp == "Dicorynia_guianensis") (Number of observations: 1268) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
G_Intercept     0.29      0.01     0.27     0.30 1.00     2015     2518
D_Intercept    33.76      1.39    31.47    36.86 1.00     2212     2101
K_Intercept     0.62      0.09     0.48     0.83 1.00     1682     1955

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.18      0.00     0.17     0.18 1.00     3586     2904

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Uma inferência com Dicorynia guianensis

plot_fit(fit = fit, Species = "Dicorynia_guianensis")

Uma inferência com Eschweilera sagotiana

fit <- brm(bf(log(Growth+1) ~ G*exp(-0.5*(log(DBH/D)^2/K)),
              G ~ 1, D ~ 1, K ~ 1,
              nl = TRUE),
           prior = c(
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "G"),
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "K"),
                 prior("uniform(1,300)", lb = 1, ub = 300, nlpar = "D")
          ),
           data = filter(dt_P6, GenSp  == "Eschweilera_sagotiana"),
           cores = 4)
save(fit, file = "data/fit_esch.Rdata")
load(file = "data/fit_esch.Rdata") ; fit
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: log(Growth + 1) ~ G * exp(-0.5 * (log(DBH/D)^2/K)) 
         G ~ 1
         D ~ 1
         K ~ 1
   Data: filter(dt_P6, GenSp == "Eschweilera_sagotiana") (Number of observations: 4529) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
G_Intercept     0.10      0.00     0.10     0.11 1.00     4335     2730
D_Intercept    27.63      3.92    20.82    36.10 1.00     3677     3077
K_Intercept     2.77      0.20     2.25     2.99 1.00     2754     1740

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.13      0.00     0.13     0.13 1.00     4047     2565

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Uma inferência com Eschweilera sagotiana

plot_fit(fit = fit, Species = "Eschweilera_sagotiana")

Modelo de crescimento individual

\[log(AGR+1)= G_{tree} \times exp[-0.5 \times( \frac{log(DBH/D)^2}{K})] ~~|~~G_{tree} = G_{species}+\epsilon_{tree}\]

expand_grid(DBH = 10:100, G_tree = rnorm(10, 0.5, 0.1)) %>% 
  mutate(log_agr = G_tree*exp(-0.5*(log(DBH/50)^2/0.2))) %>% 
  ggplot(aes(DBH, exp(log_agr)-1,
             group = as.character(G_tree))) + geom_line() + theme_bw() +
  labs(x = "Ano da medição",
       y = "Crescimento anual") + theme_bw() +
  geom_vline(xintercept = 50, colour = "red", linetype = "dashed") + 
  geom_hline(yintercept = exp(0.5)-1, colour = "red", linetype = "dashed")

Modelo de crescimento individual

fit <- brm(bf(log(Growth+1) ~ G*exp(-0.5*(log(DBH/D)^2/K)),
              G ~ 1|idTree, D ~ 1, K ~ 1,
              nl = TRUE),
           prior = c(
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "G"),
                 prior("uniform(0.001,3)", lb = 0.001, ub = 3, nlpar = "K"),
                 prior("uniform(1,300)", lb = 1, ub = 300, nlpar = "D")
          ),
           data = filter(dt_P6, GenSp  == "Dicorynia_guianensis"),
           cores = 4)
save(fit, file = "data/fit_dico_ind.Rdata")
load(file = "data/fit_dico_ind.Rdata") ; fit
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: log(Growth + 1) ~ G * exp(-0.5 * (log(DBH/D)^2/K)) 
         G ~ 1 | idTree
         D ~ 1
         K ~ 1
   Data: filter(dt_P6, GenSp == "Dicorynia_guianensis") (Number of observations: 1268) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~idTree (Number of levels: 51) 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(G_Intercept)     0.10      0.01     0.07     0.12 1.00      843     1708

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
G_Intercept     0.24      0.02     0.20     0.28 1.00      775     1235
D_Intercept    44.30     11.20    29.64    72.26 1.01     1085     1531
K_Intercept     2.04      0.57     1.02     2.95 1.00     1274     1447

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.16      0.00     0.16     0.17 1.00     3887     2804

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Modelo de crescimento individual

fit %>% 
  spread_draws(b_G_Intercept, r_idTree__G[idTree,]) %>% 
  mutate(G_tree = exp(b_G_Intercept + r_idTree__G) - 1) %>% 
  group_by(idTree) %>% 
  summarise(G_tree = median(G_tree)) %>% 
  ggplot(aes(G_tree)) + geom_histogram() + theme_bw() + 
  labs(x = "Crescimento individual (em cm/ano)")

Usar características funcionais para explicar o crescimento

Distribuição das crescimentos

Começamos recuperando todos os crescimentos das espécies de Paracou disponíveis em Schmitt, Hérault, Derroire (2023), e selecionamos as espécies para as quais temos valores de características

dt <- read_csv("data/growth_raw.csv", locale = locale(decimal_mark = ",")) %>% 
  select(Species, `Median Gmax`) %>% 
  rename(species = Species, Gmax =  `Median Gmax`) %>% 
  mutate(species = gsub(" ", "_", species)) %>% 
  right_join(dt_traits) %>% 
  filter(Gmax < 3)
dt %>% ggplot(aes(Gmax)) + geom_histogram() + theme_bw() + scale_x_log10()

Correlações entre características e crescimento

Primeiro, podemos analisar a correlação com o crescimento. Como as características estão correlacionadas, usamos apenas wsg e sla para representar as diferentes estratégias.

ggpairs(dt, columns = c("Gmax", "LMA", "wsg"))

Regressão linear

Também podemos usar a regressão linear para identificar estratégias relacionadas ao crescimento das espécies.

lm(Gmax ~ LMA + wsg, dt) %>% 
  summary()

Call:
lm(formula = Gmax ~ LMA + wsg, data = dt)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.239838 -0.082860 -0.002604  0.090345  0.208692 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.2191956  0.2169689   5.619 2.48e-05 ***
LMA          0.0001114  0.0011878   0.094  0.92633    
wsg         -0.9041136  0.2920908  -3.095  0.00624 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1309 on 18 degrees of freedom
Multiple R-squared:  0.3476,    Adjusted R-squared:  0.2752 
F-statistic: 4.796 on 2 and 18 DF,  p-value: 0.0214

Conclusão

Parabéns 👏 , agora vocês são profissionais do modelagem do crescimento!

Falando mais seriamente, ficaremos felizes em responder a quaisquer perguntas que você possa ter.