Cirad
2025-11-07
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.
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
Começamos carregando os dados na forma de objetos no formato Rdata.
⚠️ Esses dados NÃO são de acesso livre, por favor, não os distribua nem utilize para qualquer outro fim.
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
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
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\))
# 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
Podemos começar examinando a distribuição de valores para cada um dos recursos com um histograma usando o pacote ggplot
(parte de tidyverse
).
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}\]
# 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
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.
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
.
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
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.
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.
Para usar os dados de crescimento, começamos filtrando os dados da parcela 6.
[1] 4262
# 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
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)
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)
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)
\[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")
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")
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).
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)
}
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")
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).
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")
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).
\[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")
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")
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).
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()
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.
Também podemos usar a regressão linear para identificar estratégias relacionadas ao crescimento das espécies.
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
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.
Sylvain Schmitt (sylvain.schmitt@cirad.fr) & Géraldine Derroire (geraldine.derroire@cirad.fr)