#Packages library(AICcmodavg) library(sjPlot) #Input land = read.table("aurita_matrizcat_vec.txt", header = T, dec = ",") View(land) summary(land) #Correlation test - Pearson cor(land[,c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)]) #Doing the balanced models library (MuMIn) library(stats) #All of the combinations full<- glm(presence~area+altitude+enn+d_min_pen+d_min_jac+mat_3+mat_4+mat_9+mat_15+mat_21+mat_a_n_veg+mat_a_alag+mat_agricu, family=binomial(link="logit"), data=land2) options(na.action = "na.fail") fm1<-dredge(full, m.lim = c(0,4)) selec <- dredge(full, rank="AICc", extra = "adjR^2", m.lim = c(0,4)) fix(selec) ##check c-hat for global model #with simulation - 10000 install.packages("devtools") library("devtools") install.packages("DHARMa") library("DHARMa") require(DHARMa) sim_mod1 <- simulateResiduals(full, refit=F, n=10000) testDispersion(sim_mod1) plot(sim_mod1) ##compute importance (cumulative weight - w+) importance(fm1) #Moran's test res = simulateResiduals(full) testSpatialAutocorrelation(res, x = land$latitude, y = land$longitude) #CI of the most parcimonious models containing the variables #Variables: d_min_jac; mat_21; mat_4; mat_a_n_veg modelo1 = glm(presence~d_min_jac+ mat_21+ mat_4+ mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo1) cbind(coef(modelo1), confint(modelo1)) #Variable: mat3 modelo2 = glm(presence~d_min_jac+ mat_21+ mat_3+ mat_4, family = binomial (link="logit"), data = land2) summary(modelo2) cbind(coef(modelo2), confint(modelo2)) #Variable: average altitude modelo3 = glm(presence~d_min_jac+ mat_21+ mat_4+ altitude, family = binomial (link="logit"), data = land2) summary(modelo3) cbind(coef(modelo3), confint(modelo3)) #Variable: enn modelo7 = glm(presence~d_min_jac+ enn+ mat_4+ altitude, family = binomial (link="logit"), data = land2) summary(modelo7) cbind(coef(modelo7), confint(modelo7)) #Variable: pasture modelo8 = glm(presence~d_min_jac+ mat_15+ mat_4+ mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo8) cbind(coef(modelo8), confint(modelo8)) #Variable: C. penicillata modelo9 = glm(presence~d_min_jac+ d_min_pen+mat_21+mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo9) cbind(coef(modelo9), confint(modelo9)) #Variable: fragment size modelo17 = glm(presence~ area+ d_min_jac+ mat_21+mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo17) cbind(coef(modelo17), confint(modelo17)) #Variable: planted matrix modelo21 = glm(presence~ d_min_jac+ mat_21+ mat_9+ mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo21) cbind(coef(modelo21), confint(modelo21)) #Variable: flooded areas modelo24 = glm(presence~ d_min_jac+ mat_21+ mat_a_alag+ mat_a_n_veg, family = binomial (link="logit"), data = land2) summary(modelo24) cbind(coef(modelo24), confint(modelo24)) #Variable: matrix agriculture modelo25 = glm(presence~ altitude+ d_min_jac+ mat_4+ mat_agricu, family = binomial (link="logit"), data = land2) summary(modelo25) cbind(coef(modelo25), confint(modelo25)) ##Figure's code example ggplot(land, aes (x=mat_3, y=presence, color = as.factor(presence))) + geom_point(size = 2) + geom_smooth(method = "glm", se = T, method.args = list(family = "binomial"), linetype = 1, col = 1, size = 0.7) + labs(x= expression("Forest formation matrix (%)"), y= NULL) + theme(axis.text = element_text(color="black", size = 10), axis.text.y = element_text(angle = 90, hjust = 0.5), panel.background = element_rect(fill = F), axis.line = element_line(size = 0.5, colour = "black", linetype=1), axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0), size = 10), axis.title.x = element_text( size = 10), legend.title = element_blank(), legend.position = "none") + scale_y_continuous(labels=function(x) format(x, decimal.mark = ".", scientific = F)) + scale_x_continuous(labels= function(x) paste0(x*100)) + scale_color_manual(values=c("dimgray","black"), labels = c("Controle", expression(paste(italic("C. aurita")))))