Ligia Análise
Arthur Barros
29/09/2023
Organizando os dados
paged_table(arvores_individuos)paged_table(vars_by_parcela)paged_table(biomassa_herbacea)# Biomassa total das herbáceas
biomassa_herbacea %>%
group_by(Espécie) %>%
summarise(quant_biomass_esp=sum(Biomassa_Kg_ha)) # Isso aqui é só pra notar que existe espécies com 2 kg e espécies com 2000## # A tibble: 128 × 2
## Espécie quant_biomass_esp
## <fct> <dbl>
## 1 Aeschynomene viscidula 370.
## 2 alternantera brasiliana 402.
## 3 Angelonia campestris 28.8
## 4 Anthephora hermaphrodita 230.
## 5 Anthurium affinr 2.46
## 6 aristida 533.
## 7 Aristida 48.1
## 8 Arvore de Natal 2.36
## 9 asteraceae 2107.
## 10 asteraceae 2 56.6
## # ℹ 118 more rowsbiomassa_herbacea %>%
group_by(Parcela) %>%
summarise(quant_biomass_esp=sum(Biomassa_Kg_ha)) %>%
glimpse-> biomass_herb_total_by_parcela## Rows: 20
## Columns: 2
## $ Parcela <fct> p11, p14, p15, p16, p17, p2, p23, p25, p26, p27, r12…
## $ quant_biomass_esp <dbl> 12550.376, 5804.918, 10135.692, 8218.796, 5884.866, …# Riqueza de herbáceas (contagem de espécies herbáceas)
biomassa_herbacea %>%
group_by(Parcela) %>%
summarise(contagem_spp_herbacea = sum(length(unique(Espécie)))) %>%
glimpse -> contagem_spp_herbacea## Rows: 20
## Columns: 2
## $ Parcela <fct> p11, p14, p15, p16, p17, p2, p23, p25, p26, p27,…
## $ contagem_spp_herbacea <int> 26, 17, 40, 15, 19, 38, 38, 31, 24, 33, 22, 25, …# proporção de monocotiledônias
biomassa_herbacea %>%
select(2:3, 5) %>%
group_by(Parcela) %>%
summarise(prop_mono = mean(Monocotiledonea == 'sim')) %>%
glimpse->proporcao_de_mono # Quanto mais próximo de 1,00, maior a quantidade de monocotiledônias## Rows: 20
## Columns: 2
## $ Parcela <fct> p11, p14, p15, p16, p17, p2, p23, p25, p26, p27, r12, r18, r…
## $ prop_mono <dbl> 0.07031250, 0.17021277, 0.17605634, 0.23076923, 0.83950617, …# juntar tudo num dado só
vars_by_parcela %>%
left_join(arvores_individuos, by = "Parcela") %>%
left_join(biomass_herb_total_by_parcela, by = "Parcela") %>%
left_join(contagem_spp_herbacea, by = "Parcela") %>%
left_join(proporcao_de_mono, by = "Parcela") %>%
mutate(new_regeneration_col=case_when(
is.na(regeneracao_idade) ~ "sem_regeneração",
!is.na(regeneracao_idade) ~ "com_regeneração",
),
new_regeneration_col=as.factor(new_regeneration_col)) %>%
select(-regeneracao_idade) %>%
glimpse -> data_model## Rows: 20
## Columns: 8
## $ Parcela <fct> p11, p14, p15, p16, p17, p2, p23, p25, p26, p27,…
## $ mm_precipitacao <dbl> 673, 540, 510, 555, 940, 647, 785, 588, 645, 903…
## $ CAD <dbl> 0.22194715, 0.23711884, 0.06373218, 0.26504915, …
## $ quant_indiv_arvores <int> 74, 138, 86, 152, 167, 170, 151, 216, 80, 181, 1…
## $ quant_biomass_esp <dbl> 12550.376, 5804.918, 10135.692, 8218.796, 5884.8…
## $ contagem_spp_herbacea <int> 26, 17, 40, 15, 19, 38, 38, 31, 24, 33, 22, 25, …
## $ prop_mono <dbl> 0.07031250, 0.17021277, 0.17605634, 0.23076923, …
## $ new_regeneration_col <fct> sem_regeneração, sem_regeneração, sem_regeneraçã…Modelos
Pressupostos
# vendo as variáveis
summary(data_model)## Parcela mm_precipitacao CAD quant_indiv_arvores
## p11 : 1 Min. :510.0 Min. :0.0000 Min. : 63.0
## p14 : 1 1st Qu.:581.2 1st Qu.:0.2581 1st Qu.:102.5
## p15 : 1 Median :773.5 Median :0.3919 Median :159.5
## p16 : 1 Mean :747.7 Mean :0.4158 Mean :166.7
## p17 : 1 3rd Qu.:900.8 3rd Qu.:0.5858 3rd Qu.:209.2
## p2 : 1 Max. :940.0 Max. :0.9813 Max. :337.0
## (Other):14
## quant_biomass_esp contagem_spp_herbacea prop_mono
## Min. : 3280 Min. :15.00 Min. :0.02439
## 1st Qu.: 6381 1st Qu.:24.00 1st Qu.:0.17299
## Median : 7571 Median :27.00 Median :0.30558
## Mean : 8505 Mean :28.30 Mean :0.34111
## 3rd Qu.:10153 3rd Qu.:33.25 3rd Qu.:0.44679
## Max. :22423 Max. :40.00 Max. :0.83951
##
## new_regeneration_col
## com_regeneração:10
## sem_regeneração:10
##
##
##
##
## plot_histogram(data_model$quant_biomass_esp)
full_model<-glm(quant_biomass_esp ~ mm_precipitacao + CAD + quant_indiv_arvores + contagem_spp_herbacea + prop_mono + new_regeneration_col, data = data_model, family = "gaussian")
shapiro.test(residuals(full_model)) # dados não são normais pq tem o range da biomassa é gigante##
## Shapiro-Wilk normality test
##
## data: residuals(full_model)
## W = 0.79605, p-value = 0.0007549#check_homogeneity(full_model) # tá dando erro. mas dá pra vê que na função abaixo
check_model(full_model)
Basicamente o modelo tá razoável, mas não cumpre o pressuposto da normalidade. Isso se deve aos valores da biomassa que estão muito distantes da reta. Afinal, o menor valor é na casa dos 3000 e o maior valor é na casa dos 22000. A partir daqui a gente tem algumas opções.
A primeira opção é transformar a nossa variável dependente. A primeira tentativa é com o a raiz quadrada:
full_model_sqrt<-glm(sqrt(quant_biomass_esp) ~ mm_precipitacao + CAD + quant_indiv_arvores + contagem_spp_herbacea + prop_mono + new_regeneration_col, data = data_model, family = "gaussian")
shapiro.test(residuals(full_model_sqrt))##
## Shapiro-Wilk normality test
##
## data: residuals(full_model_sqrt)
## W = 0.8761, p-value = 0.01507Continua sem atender a normalidade. Então vou transformar com o log.
full_model_log10<-glm(log10(quant_biomass_esp) ~ mm_precipitacao + CAD + quant_indiv_arvores + contagem_spp_herbacea + prop_mono + new_regeneration_col, data = data_model, family = "gaussian", na.action = "na.fail")
shapiro.test(residuals(full_model_log10)) # residuos normais##
## Shapiro-Wilk normality test
##
## data: residuals(full_model_log10)
## W = 0.93371, p-value = 0.1819check_model(full_model_log10) # Nota que basicamente o log achatou os teus dados e deixou três outcomes ainda
dredge(full_model_log10) -> first_try
subset(first_try, delta < 4) # coloquei 4, mas se usa 3 ou 2## Global model call: glm(formula = log10(quant_biomass_esp) ~ mm_precipitacao + CAD +
## quant_indiv_arvores + contagem_spp_herbacea + prop_mono +
## new_regeneration_col, family = "gaussian", data = data_model,
## na.action = "na.fail")
## ---
## Model selection table
## (Int) CAD cnt_spp_hrb mm_prc new_rgn_col prp_mon qnt_ind_arv df
## 33 4.098 -0.001246 3
## 35 3.878 0.007180 -0.001145 4
## 34 4.029 0.1608 -0.001236 4
## 37 3.950 0.0002210 -0.001354 4
## 53 3.794 0.0006011 -0.38980 -0.001322 5
## 49 4.122 -0.09808 -0.001191 4
## 1 3.890 2
## 41 4.089 + -0.001224 4
## 3 3.637 0.008921 3
## 36 3.853 0.1229 0.006288 -0.001150 5
## 51 3.896 0.007569 -0.12400 -0.001071 5
## logLik AICc delta weight
## 33 8.468 -9.4 0.00 0.268
## 35 9.554 -8.4 0.99 0.163
## 34 9.109 -7.6 1.88 0.104
## 37 8.922 -7.2 2.26 0.087
## 53 10.583 -6.9 2.56 0.075
## 49 8.645 -6.6 2.81 0.066
## 1 5.512 -6.3 3.12 0.056
## 41 8.476 -6.3 3.15 0.055
## 3 6.800 -6.1 3.34 0.050
## 36 9.948 -5.6 3.83 0.040
## 51 9.869 -5.5 3.98 0.037
## Models ranked by AICc(x)summary(model.avg(first_try))##
## Call:
## model.avg(object = first_try)
##
## Component model call:
## glm(formula = log10(quant_biomass_esp) ~ <64 unique rhs>, family =
## gaussian, data = data_model, na.action = na.fail)
##
## Component models:
## df logLik AICc delta weight
## 6 3 8.47 -9.44 0.00 0.18
## 26 4 9.55 -8.44 0.99 0.11
## 16 4 9.11 -7.55 1.88 0.07
## 36 4 8.92 -7.18 2.26 0.06
## 356 5 10.58 -6.88 2.56 0.05
## 56 4 8.65 -6.62 2.81 0.04
## (Null) 2 5.51 -6.32 3.12 0.04
## 46 4 8.48 -6.29 3.15 0.04
## 2 3 6.80 -6.10 3.34 0.03
## 126 5 9.95 -5.61 3.83 0.03
## 256 5 9.87 -5.45 3.98 0.02
## 236 5 9.78 -5.28 4.16 0.02
## 246 5 9.59 -4.89 4.54 0.02
## 136 5 9.44 -4.59 4.85 0.02
## 146 5 9.42 -4.56 4.88 0.02
## 1 3 6.03 -4.56 4.88 0.02
## 5 3 5.96 -4.42 5.02 0.01
## 2356 6 11.41 -4.36 5.08 0.01
## 4 3 5.90 -4.29 5.14 0.01
## 25 4 7.45 -4.23 5.21 0.01
## 156 5 9.23 -4.17 5.27 0.01
## 24 4 7.28 -3.90 5.54 0.01
## 346 5 9.08 -3.88 5.56 0.01
## 14 4 7.24 -3.82 5.62 0.01
## 3 3 5.55 -3.59 5.84 0.01
## 12 4 7.08 -3.49 5.95 0.01
## 35 4 6.86 -3.05 6.39 0.01
## 456 5 8.65 -3.01 6.43 0.01
## 1356 6 10.72 -2.98 6.45 0.01
## 3456 6 10.71 -2.97 6.47 0.01
## 23 4 6.80 -2.93 6.50 0.01
## 124 5 8.24 -2.20 7.23 0.00
## 15 4 6.39 -2.11 7.32 0.00
## 1246 6 10.27 -2.07 7.36 0.00
## 1256 6 10.18 -1.91 7.53 0.00
## 1346 6 10.18 -1.89 7.55 0.00
## 235 5 8.07 -1.86 7.58 0.00
## 1236 6 10.12 -1.78 7.65 0.00
## 34 4 6.17 -1.67 7.77 0.00
## 45 4 6.13 -1.60 7.84 0.00
## 2346 6 9.95 -1.44 8.00 0.00
## 13 4 6.04 -1.40 8.03 0.00
## 2456 6 9.87 -1.28 8.16 0.00
## 134 5 7.65 -1.02 8.41 0.00
## 245 5 7.65 -1.02 8.42 0.00
## 125 5 7.64 -0.98 8.45 0.00
## 345 5 7.45 -0.61 8.83 0.00
## 234 5 7.39 -0.50 8.94 0.00
## 1456 6 9.44 -0.41 9.02 0.00
## 145 5 7.26 -0.23 9.20 0.00
## 123 5 7.08 0.12 9.56 0.00
## 23456 7 11.55 0.23 9.67 0.00
## 135 5 7.01 0.27 9.71 0.00
## 12356 7 11.47 0.40 9.84 0.00
## 13456 7 11.14 1.05 10.48 0.00
## 2345 6 8.64 1.18 10.62 0.00
## 1234 6 8.46 1.54 10.98 0.00
## 1245 6 8.33 1.80 11.23 0.00
## 12346 7 10.76 1.81 11.24 0.00
## 1345 6 8.25 1.97 11.40 0.00
## 1235 6 8.12 2.22 11.66 0.00
## 12456 7 10.34 2.65 12.09 0.00
## 12345 7 9.13 5.07 14.50 0.00
## 123456 8 11.80 5.48 14.92 0.00
##
## Term codes:
## CAD contagem_spp_herbacea mm_precipitacao
## 1 2 3
## new_regeneration_col prop_mono quant_indiv_arvores
## 4 5 6
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value
## (Intercept) 3.911e+00 2.268e-01 2.360e-01 16.568
## quant_indiv_arvores -9.375e-04 6.848e-04 7.081e-04 1.324
## contagem_spp_herbacea 2.504e-03 4.739e-03 4.906e-03 0.510
## CAD 3.827e-02 1.080e-01 1.127e-01 0.340
## mm_precipitacao 8.406e-05 2.301e-04 2.380e-04 0.353
## prop_mono -5.599e-02 1.567e-01 1.623e-01 0.345
## new_regeneration_colsem_regeneração 9.808e-03 4.643e-02 4.887e-02 0.201
## Pr(>|z|)
## (Intercept) <2e-16 ***
## quant_indiv_arvores 0.186
## contagem_spp_herbacea 0.610
## CAD 0.734
## mm_precipitacao 0.724
## prop_mono 0.730
## new_regeneration_colsem_regeneração 0.841
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value
## (Intercept) 3.9107101 0.2267925 0.2360402 16.568
## quant_indiv_arvores -0.0012190 0.0005163 0.0005556 2.194
## contagem_spp_herbacea 0.0074313 0.0054822 0.0059018 1.259
## CAD 0.1638150 0.1712210 0.1837857 0.891
## mm_precipitacao 0.0003335 0.0003561 0.0003763 0.886
## prop_mono -0.2295834 0.2466756 0.2611279 0.879
## new_regeneration_colsem_regeneração 0.0540458 0.0974124 0.1037727 0.521
## Pr(>|z|)
## (Intercept) <2e-16 ***
## quant_indiv_arvores 0.0282 *
## contagem_spp_herbacea 0.2080
## CAD 0.3727
## mm_precipitacao 0.3755
## prop_mono 0.3793
## new_regeneration_colsem_regeneração 0.6025
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Resultados: 6, 26 e 16 CAD, contagem_spp_herbacea Oa melhores modelos foram, em ordem, o modelo apenas com a variável quant_indiv_arvores, depois o modelo com as variáveis independetes: contagem_spp_herbacea + quant_indiv_arvores, e por fim o modelo com as seguintes variáveis independentes: CAD + quant_indiv_arvores. Mas note que a diferença entre eles é de menos de 2.0 AIC. O que pode indicar que eles são praticamente iguais.
- Observação: Existe uma galera que fala que até 3.0 AIC é tudo igual.
Dentre os três melhores modelos, apenas o primeiro modelo apresentou alguma diferença significativa. Veja:
result_dredge_m1<-glm(log10(quant_biomass_esp) ~ quant_indiv_arvores, data = data_model,
family = "gaussian", na.action = "na.fail")
summary(result_dredge_m1)##
## Call:
## glm(formula = log10(quant_biomass_esp) ~ quant_indiv_arvores,
## family = "gaussian", data = data_model, na.action = "na.fail")
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0975024 0.0914193 44.821 <2e-16 ***
## quant_indiv_arvores -0.0012457 0.0005007 -2.488 0.0229 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.02789541)
##
## Null deviance: 0.67479 on 19 degrees of freedom
## Residual deviance: 0.50212 on 18 degrees of freedom
## AIC: -10.936
##
## Number of Fisher Scoring iterations: 2data_model %>%
ggplot(aes(quant_indiv_arvores, log10(quant_biomass_esp))) +
geom_jitter(alpha=0.5)+
geom_smooth(method = "lm")+
theme_classic()+
labs(x= "Frequência de árvores", y="Log da Biomassa de Herbáceas")
Mas acho que precisamos ser cuidados ao afirmar um efeito da frequência de árvores na biomassa de herbácea. Porque, basicamente, os modelos são muito parecidos e nenhum outro aprensentou significância. Além disso, o estimate e o valor Z são bem baixos, o que indicam um efeito pequeno.
Além disso, existem muitos outliers na biomassa. 3000 a 22000 é muita coisa de diferença. Os valores reais, sem transformar em hectates, apresentam essa range alta também?
Acho que é interessante perceber que vey, quanto mais pontos no gráfico melhor. Então, ao invés de fazer por parcela, a gente pode pegar a precipitação, o CAD, a proporção de monocotiledonias e qualquer var que se repita por parcela e aplicar. Mas o problema da alta range dos dados ainda vai rolar.
E como pode alguma parcela ter 22000 kg de biomassa por 1 hectare? Isso não é muito?