Avaliação final dos modelos
Arthur Barros, Rodrigo Carmo e FeLipe Melo
08/12/2022
Dado
## tibble [934 x 12] (S3: tbl_df/tbl/data.frame)
## $ frag_cat : Factor w/ 4 levels "Control","Large",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ area : num [1:934] 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 ...
## $ site : Factor w/ 8 levels "C1","C2","G1",..: 7 7 7 7 7 7 7 7 7 7 ...
## $ plot : chr [1:934] "Q1" "Q1" "Q1" "Q1" ...
## $ disp_mode : Factor w/ 2 levels "abio","zoo": 2 2 1 2 1 2 2 2 2 2 ...
## $ seed_size_cat: chr [1:934] "M" "P" "P" "G" ...
## $ habit : chr [1:934] "tree" "und-tree" "tree" "tree" ...
## $ origen : Factor w/ 3 levels "d","nd","x": 2 2 2 2 1 2 2 2 2 1 ...
## $ id_real : Factor w/ 43 levels "1","5","6","7",..: 32 30 26 10 15 7 18 37 39 22 ...
## $ abund : num [1:934] 40 1 38 5 20 0 5 1 114 0 ...
## $ yr : Factor w/ 2 levels "2006","2007": 1 1 1 1 1 1 1 1 1 1 ...
## $ site_plot : chr [1:934] "P1 Q1" "P1 Q1" "P1 Q1" "P1 Q1" ...Vamos calcular a diversidade interpolada usando o inext. Bom, esse é o nosso dado:
## # A tibble: 95 x 12
## frag_cat area site plot disp_mode seed_size_cat habit origen id_real abund
## <fct> <dbl> <fct> <chr> <fct> <chr> <chr> <fct> <fct> <dbl>
## 1 Large 640 G1 Q1 zoo M tree nd 71 12
## 2 Large 640 G1 Q1 zoo G tree nd 18 45
## 3 Large 640 G1 Q1 zoo M tree nd 13 25
## 4 Large 640 G1 Q1 zoo G tree d 38 2
## 5 Large 640 G1 Q1 zoo G tree d 8 2
## 6 Large 640 G1 Q1 zoo M tree nd 91 3
## 7 Large 640 G1 Q1 zoo G tree d 105 2
## 8 Large 640 G1 Q1 zoo P tree nd 100 200
## 9 Large 640 G1 Q1 zoo M tree d 73 0
## 10 Large 640 G1 Q1 zoo G tree nd 55 0
## # ... with 85 more rows, and 2 more variables: yr <fct>, site_plot <chr>Vou tratar o fragmento como variável contínua. E para facilitar os gráficos, vou usar sempre Log10(area)
div_bruto %>%
filter(yr =="2006") %>%
rename("div0_06"=div0, "div1_06"=div1, "div2_06"=div2) %>%
select(area, disp_mode, origen, div0_06, div1_06, div2_06)-> div_cont_area_dummy_06
div_bruto %>%
filter(yr =="2007") %>%
rename("div0_07"=div0, "div1_07"=div1, "div2_07"=div2) %>%
select(area, disp_mode, origen, div0_07, div1_07, div2_07) %>%
mutate(div0_06 = div_cont_area_dummy_06$div0_06,
div1_06 = div_cont_area_dummy_06$div1_06,
div2_06 = div_cont_area_dummy_06$div2_06,
delta_div0= div0_07 - div0_06,
delta_div1= div1_07 - div1_06,
delta_div2= div2_07 - div2_06)->div_cont_area_geral
histogram(~delta_div0, data = div_cont_area_geral) #gaussian
histogram(~delta_div1, data = div_cont_area_geral) #gaussian
histogram(~delta_div2, data = div_cont_area_geral) #gaussian
# div 0
div_cont_area_geral %>%
ggplot(aes(x=log10(area), y=delta_div0))+
geom_jitter(alpha=0.5)+
geom_smooth(method =lm, se=T)+
facet_grid(disp_mode ~ origen, scales="free")
d0_area_m1<-lm(delta_div0~area*disp_mode+area*origen, data = div_cont_area_geral,
na.action = "na.fail")
summary(d0_area_m1)
##
## Call:
## lm(formula = delta_div0 ~ area * disp_mode + area * origen, data = div_cont_area_geral,
## na.action = "na.fail")
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1020 -0.7088 0.1472 0.5680 4.0122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.910e+00 1.348e-01 -14.166 < 2e-16 ***
## area -2.457e-04 6.499e-05 -3.781 0.000178 ***
## disp_modezoo -2.103e+00 1.589e-01 -13.236 < 2e-16 ***
## origennd 1.616e+00 1.589e-01 10.175 < 2e-16 ***
## area:disp_modezoo 5.001e-04 7.537e-05 6.635 9.77e-11 ***
## area:origennd 3.301e-04 7.537e-05 4.379 1.50e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.46 on 431 degrees of freedom
## Multiple R-squared: 0.4805, Adjusted R-squared: 0.4745
## F-statistic: 79.74 on 5 and 431 DF, p-value: < 2.2e-16
check_model(d0_area_m1)
check_distribution(d0_area_m1)
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 69%
## tweedie 19%
## weibull 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## normal 62%
## tweedie 25%
## weibull 6%Tenho dificuldades para interpretar variáveis contínuas de maneira biológica. Se fosse categórica, eu conseguiria dizer qual categoria tem maior variação da diversidade interanualmente, mas aqui eu não sei o que o p significativo de área significa. Quanto maior a área, menor é a variação internual? ou o contrário, quanto menor a área, menor variação anual da diversidade sensível as espécies raras? Confuso :(
Sobre os modelos, é interessante notar que o modelo tem um R de 45%, F = 72 e outliers que podem influenciar o modelo (gráfico de influential observations)
Agora para a diverisidade sensível as espécies comuns, temos:
# div1
div_cont_area_geral %>%
ggplot(aes(x=log10(area), y=delta_div1))+
geom_jitter(alpha=0.5)+
geom_smooth(method =lm, se=T)+
facet_grid(disp_mode ~ origen, scales="free")
d1_area_m2<-lm(delta_div1~area*disp_mode+area*origen, data = div_cont_area_geral,
na.action = "na.fail")
summary(d1_area_m2)
##
## Call:
## lm(formula = delta_div1 ~ area * disp_mode + area * origen, data = div_cont_area_geral,
## na.action = "na.fail")
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0696 -0.7090 -0.1481 0.4222 3.6708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.371e+00 1.201e-01 -11.422 < 2e-16 ***
## area -8.026e-05 5.787e-05 -1.387 0.166
## disp_modezoo -1.163e+00 1.414e-01 -8.219 2.42e-15 ***
## origennd 2.112e+00 1.414e-01 14.929 < 2e-16 ***
## area:disp_modezoo 6.639e-05 6.711e-05 0.989 0.323
## area:origennd 3.015e-05 6.711e-05 0.449 0.653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.3 on 431 degrees of freedom
## Multiple R-squared: 0.4553, Adjusted R-squared: 0.449
## F-statistic: 72.05 on 5 and 431 DF, p-value: < 2.2e-16
check_model(d1_area_m2)
check_distribution(d1_area_m2)
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 69%
## tweedie 31%
##
## Predicted Distribution of Response
##
## Distribution Probability
## normal 69%
## tweedie 22%
## weibull 9%A variação na produção de sementes é equivalente em todos os fragmentos (area - p value > 0,1) e atinge os padrões de dispersão de diferentes maneiras. De modo geral, encontramos uma menor variação anual nas nas sementes zoocóricas que abióticas (p value < 0,001) e, uma maior variação nas sementes não dispersadas quando comparado com as dispersadas (p value < 0,001). No entanto, somente para as sementes zoocóricas não dispersadas é que obtivemos uma diferença na variação da diversidade entre os diferentes tamanhos de fragmentos (Figura X). Fragmentos até 30 hectares apresentam variação negativa (menos diversidade em 2007 do que em 2006) enquanto os maiores apresentam uma variação positiva entre anos (maior diversidade em 2007 do que em 2006).
Aqui o modelo parece um pouco melhor já que não apresenta outliers com capacidade de modificar o modelo. R e F statitic continuam semelhantes.
Agora para a diverisidade sensível as espécies dominantes, temos:
div_cont_area_geral %>%
ggplot(aes(x=log10(area), y=delta_div2))+
geom_jitter(alpha=0.5)+
geom_smooth(method =lm, se=T)+
facet_grid(disp_mode ~ origen, scales="free")
d2_area_m2<-lm(delta_div2~area*disp_mode+area*origen, data = div_cont_area_geral,
na.action = "na.fail")
summary(d2_area_m2)
##
## Call:
## lm(formula = delta_div2 ~ area * disp_mode + area * origen, data = div_cont_area_geral,
## na.action = "na.fail")
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3793 -0.5474 -0.1205 0.3802 3.0599
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.169e+00 1.020e-01 -11.457 < 2e-16 ***
## area -2.323e-05 4.918e-05 -0.472 0.637
## disp_modezoo -8.724e-01 1.202e-01 -7.258 1.85e-12 ***
## origennd 1.863e+00 1.202e-01 15.495 < 2e-16 ***
## area:disp_modezoo -2.180e-05 5.703e-05 -0.382 0.703
## area:origennd -1.609e-05 5.703e-05 -0.282 0.778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.105 on 431 degrees of freedom
## Multiple R-squared: 0.4576, Adjusted R-squared: 0.4513
## F-statistic: 72.73 on 5 and 431 DF, p-value: < 2.2e-16
check_model(d2_area_m2)
check_distribution(d2_area_m2)
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 69%
## tweedie 31%
##
## Predicted Distribution of Response
##
## Distribution Probability
## normal 69%
## tweedie 22%
## weibull 9%E minha interpretação é a mesma do que as espécies comuns, os gráficos e modelos apresentam resultados bem semelhantes. Acredito que a única diferença é que nesse modelo, outliers parecem puxar a curva para fora dos pontos (influential observations) . Então, ressaltando, a área só é significativa no modelo de espécies raras.
Anteriormente, o que estava fazendo era isso
A título de comparação, caso queiram, deixei o que tinha feito antes aqui embaixo. Os modelos apresentam maiores valores de R e F, de modo geral, apesar das variáveis serem colineares.
Esse aqui é o dado com as diversidades interpoladas. Diferente do dado lá de cima, aqui eu usei apenas o fragmento categórico ao invés da área. Isso gera menos colunas (afinal temos 4 classes de fragmentos e 6 áreas).
Buscando responder as perguntas iniciais (Qual o efeito individual e composto da área do fragmento e da produção de sementes na diversidade e nos padrões de dispersão) eu pensei nesses modelos:
div_0~ yr+frag_cat+disp_mod+origen- Modelo só com efeitos individuaisdiv_0~ frag_cat+disp_mod+origen+(1|yr)- Modelo interpretando a produção de sementes como variável aleatoriadiv_0~ yr*frag_cat+disp_mod+origen- Modelo com efeito compostodiv_0~ yr*frag_cat + yr*frag_cat*disp_mod + yr*frag_cat*origen- Modelo ideal, mas impraticável de interpretar
Usei o Munim (função dredge) para vê qual séria o melhor modelo e para cada uma das opções (modelo só com efeitos individuais, modelo com efeito composto, ano como variável), os modelos estípuladas inicialmente (estes que escrevi acima para vocês) eram os melhores (AIC value).
Aqui está a comparação entre cada modelo
library(performance)
library(see)
library(MuMIn)
#histogram(~div0, data = div_data_geral)
#histogram(~div1, data = div_data_geral)
#histogram(~div2, data = div_data_geral)
# model without interaction. Here, the effect of patch are and yr is isolated
mod_patch_div0<-lm(div0~yr+frag_cat+origen+disp_mode, data = div_data_geral,
na.action = "na.fail")
#check_model(mod_patch_div0) # the model is ok
#dredge(mod_patch_div0, evaluate = TRUE, rank = "AICc", fixed = NULL, trace = TRUE)
# the result show that the best model is this
# here, year is considered a random effect in model. It's a very common form to analyze the year
model_lmer_div0<-lmer(div0~frag_cat+origen+disp_mode+(1|yr), data = div_data_geral,
na.action = "na.fail")
#check_model(model_lmer_div0) # the model is ok too
#dredge(model_lmer_div0, evaluate = TRUE, rank = "AICc", fixed = NULL, trace = TRUE)
# Again, the result show that the best model is this
# and Here, we have the compound effect between year(seed production) and each variable
model_compound_div0<-lm(div0~yr*frag_cat+origen+disp_mode, data = div_data_geral,
na.action = "na.fail")
check_model(model_compound_div0)
#dredge(model_compound_div0, evaluate = TRUE, rank = "AICc", fixed = NULL, trace = TRUE)
# The result show that the model with e without the collinearity lose few AIC, but the weight
# the role of collinearity. I think that, in this case, colinearity explain the compound effect.
model_dreams_div0<-lm(div0~yr*frag_cat+yr*frag_cat*origen+yr*frag_cat*disp_mode, data = div_data_geral,
na.action = "na.fail")
check_model(model_dreams_div0)
compare_performance(mod_patch_div0,
model_lmer_div0,
model_compound_div0,
model_dreams_div0,
metrics = "all", rank = T)
## # Comparison of Model Performance Indices
##
## Name | Model | RMSE | Sigma | AIC weights | BIC weights | Performance-Score
## ---------------------------------------------------------------------------------------------
## model_dreams_div0 | lm | 2.234 | 2.279 | 1.00 | 1.000 | 100.00%
## model_compound_div0 | lm | 2.493 | 2.513 | 1.32e-23 | 3.37e-10 | 7.43%
## mod_patch_div0 | lm | 2.539 | 2.553 | 3.84e-27 | 7.33e-11 | 1.15e-03%
## model_lmer_div0 | lmerMod | 2.539 | 2.553 | 1.54e-29 | 2.94e-13 | 4.28e-08%
compare_performance(mod_patch_div0,
#model_lmer_div0,
model_compound_div0,
model_dreams_div0,
metrics = "all", rank = T)
## # Comparison of Model Performance Indices
##
## Name | Model | R2 | R2 (adj.) | RMSE | Sigma | AIC weights | BIC weights | Performance-Score
## ---------------------------------------------------------------------------------------------------------------
## model_dreams_div0 | lm | 0.757 | 0.748 | 2.234 | 2.279 | 1.00 | 1.000 | 100.00%
## model_compound_div0 | lm | 0.698 | 0.693 | 2.493 | 2.513 | 1.32e-23 | 3.37e-10 | 10.16%
## mod_patch_div0 | lm | 0.686 | 0.683 | 2.539 | 2.553 | 3.84e-27 | 7.33e-11 | 0.00%Os melhores modelos foram o modelo com efeito composto e o modelo ideal. O interressante aqui é a colinearidade. No modelo de efeito composto, ano e fragmento são colineares e no modelo ideal qualquer variável com interação com o fragmento colinearidade. Aqui temos duas coisas importantes:
Primeiro, colinearidade significa que duas variáveis apresentam uma relação linear entre si e “Most commonly, collinearity is intrinsic, meaning that collinear variables are different manifestations of the same underlying, and in some cases, immeasurable process (or latent variable).” (Dorrnam et al., 2013).
Aqui, a colinearidade entre yr*frag_cat não é um problema pois alegamos desde o começo um efeito composto entre essa duas variáveis (i.e. a colinearidade é o que confirma que existe um efeito composto aqui). No entanto, é inviável interpretar yr*frag_cat*origen ou yr*frag_cat*disp_mode. O que podemos fazer para entender melhor esse efeito é usar apenas uma das variáveis (nessa caso o ano - o fragmento é colinear com todas as outras variáveis) dado que ambas tem uma relação linear e possuem parcialmente a explicação uma da outra.
E note que é possível interpretar a colinearidade entre fragmento e as variáveis de dispersão (disp_mode e origen) pelos resultados de defaunação que Melo et al apresentam na tabela 1. Assim, uma possível explicação para a colinearidade que vemos entre as variáveis é devido ao efeito ecológico da defaunação que conecta menores fragmentos a menores diversidade de sementes zoocóricas dispersadas.
Assim, sugeri esse modelo a posteriori
div0~yr*frag_cat+yr*origen+yr*disp_mode
model_final_div0<-lm(div0~yr*frag_cat+yr*origen+yr*disp_mode, data = div_data_geral,
na.action = "na.fail")
check_model(model_final_div0)
compare_performance(model_final_div0,
model_compound_div0,
#model_dreams_div0,
metrics = "all", rank = T)
## # Comparison of Model Performance Indices
##
## Name | Model | R2 | R2 (adj.) | RMSE | Sigma | AIC weights | BIC weights | Performance-Score
## ---------------------------------------------------------------------------------------------------------------
## model_final_div0 | lm | 0.727 | 0.722 | 2.368 | 2.392 | 1.000 | 1.000 | 100.00%
## model_compound_div0 | lm | 0.698 | 0.693 | 2.493 | 2.513 | 2.09e-13 | 1.72e-11 | 0.00%
summary(model_final_div0)
##
## Call:
## lm(formula = div0 ~ yr * frag_cat + yr * origen + yr * disp_mode,
## data = div_data_geral, na.action = "na.fail")
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.4163 -0.7086 0.2031 1.2873 6.5677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3561 0.3360 12.965 < 2e-16 ***
## yr2007 -1.3853 0.4752 -2.915 0.00369 **
## frag_catlarge -0.6535 0.3880 -1.684 0.09263 .
## frag_catmedium 2.0521 0.3880 5.289 1.73e-07 ***
## frag_catsmall 0.6045 0.3880 1.558 0.11973
## origennd -1.3931 0.2743 -5.078 5.11e-07 ***
## disp_modezoo 8.0082 0.2743 29.191 < 2e-16 ***
## yr2007:frag_catlarge -0.6714 0.5487 -1.224 0.22159
## yr2007:frag_catmedium -2.4989 0.5487 -4.554 6.37e-06 ***
## yr2007:frag_catsmall -1.7552 0.5487 -3.199 0.00145 **
## yr2007:origennd 2.2756 0.3880 5.865 7.44e-09 ***
## yr2007:disp_modezoo -2.1253 0.3880 -5.478 6.36e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.392 on 596 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7221
## F-statistic: 144.4 on 11 and 596 DF, p-value: < 2.2e-16
div_data_geral %>%
ggplot(aes(x=frag_cat , y=div0, fill=yr, color=yr))+
geom_boxplot(alpha=0.5)+
facet_grid(disp_mode ~ origen, scales="free")+
scale_color_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
scale_fill_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
theme_classic()+
labs(y="Alpha diversity", x="Patch Area", caption = "Q=0")
div_data_geral %>%
ggplot(aes(x=frag_cat , y=div1, fill=yr, color=yr))+
geom_boxplot(alpha=0.5)+
facet_grid(disp_mode ~ origen, scales="free")+
scale_color_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
scale_fill_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
theme_classic()+
labs(y="Alpha diversity", x="Patch Area", caption = "Q=1")
div_data_geral %>%
ggplot(aes(x=frag_cat , y=div2, fill=yr, color=yr))+
geom_boxplot(alpha=0.5)+
facet_grid(disp_mode ~ origen, scales="free")+
scale_color_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
scale_fill_manual(name="Year",
labels=c("2006", "2007"),
values = c("#FFCCCC", "#FF9966")) +
theme_classic()+
labs(y="Alpha diversity", x="Patch Area", caption = "Q=2")