knitr::opts_chunk$set(echo = TRUE)

To explore the Soil compactation effect on Carbon and Biodiversity

The initial idea:

To sample Carbon, Bettles, Ants and Spiders in a 100-150 m grid in Cubarral (Meta) to model the spatially splicit relatinship between covariates, incorporating remote sensing info. The initial idea changed to sample just ants and dung bettles in a 80 m grid. The Carbon data came in September 2018.

Draft Idea

draft

draft

Data analysis sugestions y Diego Navarrete

Question:

  • To what extent soil compaction influences biodiversity and carbon content in soils of natural and human-modified ecosystems?

Posible Results for Carbon

ia an ANOVA:

  • SOC with depth in each cover and among covers.

  • Bulk density with depth in each cover and among covers.

  • Clay, silt and sand with depth in each cover and among covers.

Correlations (multiple?):

  • SOC and bulk density.
  • All data.
  • In each cover.
  • In each cover and depth.
  • SOC + bulk density + clay, silt, sand
  • All data
  • In each cover;
  • In each cover and depth.

Data management

“Yopo” y “Pastizal arbolado” se convierte en silvopastoril Ripario y bosque se convierten en Bosque

library (tidyverse) # Data handling
library (raster)  # grids, rasters
library (rgeos) # coord transform
library (rgdal) # coord transform
library (ggmap) # visualization
library (lubridate) # date managemgemt
library (readxl) # read excel data
library (sf) # vector maps in R
library (sp) # vector old fashon maps in R
library (mapview) # see maps in web
library (sjPlot) # wiew GLM plots
library (sjmisc) # wiew GLM plots
library (sqldf) # Table with permutations
library (fields) # to do image.plot
library (rsq)
library (easyGgplot2) # multiplot
library (lme4) # GLMM
library (ggeffects)
library (caret)

Read Data

For Dung beetle and ants

#path <- "C:/Users/diego.lizcano/Box Sync/CodigoR/Carbon-Biodiv/Data"
path <- "C:/Users/diego.lizcano/Box Sync/CodigoR/Carbon-Biodiv/Data"
path2 <- "D:/BoxFiles/Box Sync/CodigoR/Carbon-Biodiv/Data"


beetle_data <- read_excel(paste(path,"/TNC_Escarabajo_Carbono_June_2018.xlsx", sep=""), 
    sheet = "Base_Escarabajos")

beetle_data <- beetle_data %>% group_by(Punto) %>% 
  mutate (Richness = n_distinct(Especie))  # New column richness

beetle_data$LongDD <- beetle_data$LonDD * (-1)

ants_data <- read_excel(paste(path,"/TNC_Escarabajo_Carbono_June_2018.xlsx", sep=""), 
    sheet = "Base_Hormigas")

ants_data <- ants_data %>% group_by(Punto) %>% 
  mutate (Richness = n_distinct(Especie))

ants_data$LongDD <- ants_data$LonDD * (-1)
 
#### Check maximum
print("max")
## [1] "max"
max(unique(beetle_data$Punto))
## [1] 722
#### Check minumum
print("min")
## [1] "min"
min(unique(beetle_data$Punto))
## [1] 16
# homogeniza cobertura
ind1 <- which(beetle_data$Cobertura=="Yopo")
beetle_data$Cobertura[ind1]<-"Silvopastoril"

ind2 <- which(beetle_data$Cobertura=="Pastizal arbolado")
beetle_data$Cobertura[ind2]<-"Silvopastoril"

ind3 <- which(beetle_data$Cobertura=="Ripario")
beetle_data$Cobertura[ind3]<-"Bosque"

ind1 <- which(ants_data$Cobertura=="Yopo")
ants_data$Cobertura[ind1]<-"Silvopastoril"

ind2 <- which(ants_data$Cobertura=="Pastizal arbolado")
ants_data$Cobertura[ind2]<-"Silvopastoril"

ind3 <- which(ants_data$Cobertura=="Ripario")
ants_data$Cobertura[ind3]<-"Bosque"

ind4 <- which(ants_data$Cobertura=="Cerca Viva Yopo")
ants_data$Cobertura[ind4]<-"Silvopastoril"

# elimina datos de cultivo
ind5 <- which(beetle_data$Cobertura=="Cultivo")
beetle_data <- beetle_data[-ind5,]

ind6 <- which(ants_data$Cobertura=="Cultivo")
ants_data <- beetle_data[-ind6,]

For carbono

carbono_data <- read_excel(paste(path,"/CARBONO-08-10-2018.xlsx", sep=""), 
    sheet = "Densidad-Carbono")



ind1 <- which(carbono_data$Cobertura=="Yopo")
carbono_data$Cobertura[ind1]<-"Silvopastoril"

ind2 <- which(carbono_data$Cobertura=="Pastizal arbolado")
carbono_data$Cobertura[ind2]<-"Silvopastoril"

ind3 <- which(carbono_data$Cobertura=="Ripario")
carbono_data$Cobertura[ind3]<-"Bosque"

### select < =15
carbono_data2 <- subset(carbono_data, Carbono <= 15)

For carbono check missing points as deep 30 cm in point 711 and deep 10cm in point 718. Take in to account that maximum point number in biodiversity is 722 and minimum 16.

Make training and test data for carbon

set.seed(2000)

#generando datos de prueba y entrenamiento con los datos de carbon. Particion 75-25%
indexes_pres <-  sample(1:nrow(carbono_data2), size=0.25*nrow(carbono_data2))

###
#identifica el 25% de los datos como datos prueba/test
carbono_test <-  carbono_data2[indexes_pres,]
dim(carbono_test)
## [1] 431  16
#identifica el 75% de los datos como datos entrenamiento/train
carbono_train = carbono_data2[-indexes_pres,]
dim(carbono_train)
## [1] 1294   16

Plot Carbon Data

Per depth

As histogram

media <-  mean(carbono_data2$Carbono_percent,  na.rm = T)
mediana <-median(carbono_data2$Carbono_percent, na.rm = T)

media_d <-  mean(carbono_data2$Densidad,  na.rm = T)
mediana_d <-median(carbono_data2$Densidad, na.rm = T)

lineas <- carbono_data2 %>%
  group_by(Profundidad) %>%
  summarise(
    media = mean(Carbono),
    mediana = median(Carbono)
  )

c1 <- ggplot (data=carbono_data2, aes(Carbono) ) +
  geom_histogram (color = "white", bins = 100) +
  geom_density(color = "sky blue") +
  geom_rug (aes(Carbono) ) +
  geom_vline(aes(xintercept=media), 
             data=lineas, linetype="dashed", color = "blue", size=1) +
  geom_vline(aes(xintercept=mediana), 
             data=lineas, linetype="dashed", color = "red", size=1) +
  facet_grid (Profundidad ~ .) +
  ggtitle("Counts of carbon percent by depth, blue=mean, red=median")

c1

lineas_d <- carbono_data2 %>%
  group_by(Profundidad) %>%
  summarise(
    media = mean(Densidad),
    mediana = median(Densidad)
  )


d1 <- ggplot (data=carbono_data2, aes(Densidad) ) +
  geom_histogram (color = "white", bins = 100) +
  geom_density(color = "sky blue") +
  geom_rug (aes(Densidad) ) +
  geom_vline(aes(xintercept=media), 
             data=lineas_d, linetype="dashed", color = "blue", size=1) +
  geom_vline(aes(xintercept=mediana), 
             data=lineas_d, linetype="dashed", color = "red", size=1) +
  facet_grid (Profundidad ~ .) +
  ggtitle("Densidad by depth, blue=mean, red=median")

d1

As boxplot

c2 <- ggplot (data=carbono_data2, aes(x=factor(Profundidad), y=Carbono)) +
  geom_boxplot(width = 0.5, outlier.shape = NA) + 
  geom_jitter(alpha = 1/5, width = 0.2)


c2

Is the carbon content diferent per deept?

Test for diferences in deept in several ANOVAs

glm1 <- glm(Carbono ~ factor(Profundidad), data=carbono_data2)
tab_model (glm1)
  Carbono
Predictors Estimates CI p
(Intercept) 5.70 5.48 – 5.92 <0.001
factor(Profundidad)20 -1.11 -1.42 – -0.80 <0.001
factor(Profundidad)30 -1.35 -1.66 – -1.04 <0.001
Observations 1725
Cox & Snell’s R2 / Nagelkerke’s R2 0.290 / 0.291 
anova (glm1, test = "Chisq")
a1 <- aov(carbono_data2$Carbono ~ factor(carbono_data2$Profundidad))
TukeyHSD(a1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = carbono_data2$Carbono ~ factor(carbono_data2$Profundidad))
## 
## $`factor(carbono_data2$Profundidad)`
##             diff        lwr        upr     p adj
## 20-10 -1.1086376 -1.4781396 -0.7391356 0.0000000
## 30-10 -1.3492339 -1.7184203 -0.9800475 0.0000000
## 30-20 -0.2405962 -0.6080019  0.1268094 0.2742637
# densidad
glm1.2<- glm(Densidad ~ factor(Profundidad), data=carbono_data2)
tab_model (glm1.2)
  Densidad
Predictors Estimates CI p
(Intercept) 1.23 1.21 – 1.25 <0.001
factor(Profundidad)20 0.05 0.02 – 0.07 <0.001
factor(Profundidad)30 0.07 0.05 – 0.10 <0.001
Observations 1725
Cox & Snell’s R2 / Nagelkerke’s R2 0.001 / 0.018 
anova(glm1.2, test = "Chisq")
a3 <- aov(carbono_data2$Densidad ~ factor(carbono_data2$Profundidad))
TukeyHSD(a3)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = carbono_data2$Densidad ~ factor(carbono_data2$Profundidad))
## 
## $`factor(carbono_data2$Profundidad)`
##             diff          lwr        upr     p adj
## 20-10 0.04789779  0.016528537 0.07926705 0.0010220
## 30-10 0.07253243  0.041189966 0.10387489 0.0000002
## 30-20 0.02463463 -0.006556653 0.05582592 0.1530252
# cobertura coarbono
glm1.3<- glm(Carbono ~ factor(Cobertura), data=carbono_data2)
tab_model (glm1.3)
  Carbono
Predictors Estimates CI p
(Intercept) 5.32 5.12 – 5.52 <0.001
Pastizal -1.09 -1.36 – -0.81 <0.001
Rastrojo -1.12 -1.75 – -0.50 <0.001
Silvopastoril 0.54 0.13 – 0.95 0.010
Observations 1725
Cox & Snell’s R2 / Nagelkerke’s R2 0.328 / 0.328 
anova (glm1.3, test = "Chisq")
a1.3 <- aov(carbono_data2$Carbono ~ factor(carbono_data2$Cobertura))
TukeyHSD(a1.3)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = carbono_data2$Carbono ~ factor(carbono_data2$Cobertura))
## 
## $`factor(carbono_data2$Cobertura)`
##                               diff          lwr        upr     p adj
## Pastizal-Bosque        -1.08929882 -1.449955487 -0.7286422 0.0000000
## Rastrojo-Bosque        -1.12442667 -1.939637840 -0.3092155 0.0022612
## Silvopastoril-Bosque    0.53913448  0.003710059  1.0745589 0.0476885
## Rastrojo-Pastizal      -0.03512784 -0.848279398  0.7780237 0.9995097
## Silvopastoril-Pastizal  1.62843330  1.096150009  2.1607166 0.0000000
## Silvopastoril-Rastrojo  1.66356114  0.759226285  2.5678960 0.0000144
# cobertura densidad
glm1.4<- glm(Densidad ~ factor(Cobertura), data=carbono_data2)
tab_model (glm1.4)
  Densidad
Predictors Estimates CI p
(Intercept) 1.19 1.17 – 1.21 <0.001
Pastizal 0.14 0.12 – 0.17 <0.001
Rastrojo 0.15 0.10 – 0.20 <0.001
Silvopastoril 0.07 0.03 – 0.10 <0.001
Observations 1725
Cox & Snell’s R2 / Nagelkerke’s R2 0.005 / 0.089 
anova (glm1.4, test = "Chisq")
a1.4 <- aov(carbono_data2$Densidad ~ factor(carbono_data2$Cobertura))
TukeyHSD(a1.4)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = carbono_data2$Densidad ~ factor(carbono_data2$Cobertura))
## 
## $`factor(carbono_data2$Cobertura)`
##                                diff         lwr          upr     p adj
## Pastizal-Bosque         0.143830455  0.11420808  0.173452830 0.0000000
## Rastrojo-Bosque         0.146939841  0.07998283  0.213896849 0.0000001
## Silvopastoril-Bosque    0.068443105  0.02446626  0.112419952 0.0003804
## Rastrojo-Pastizal       0.003109386 -0.06367846  0.069897229 0.9993867
## Silvopastoril-Pastizal -0.075387350 -0.11910620 -0.031668497 0.0000580
## Silvopastoril-Rastrojo -0.078496736 -0.15277388 -0.004219593 0.0335671

Si hay diferencias!

Plots to support anovas results

carbono_data2$Profundidad <- as.factor(carbono_data2$Profundidad)


# Densidad
BP1 <- ggplot(data=carbono_data2, aes(x=factor(Profundidad), y=Densidad, fill=Cobertura)) + geom_boxplot() +
  theme_classic() + 
  ylab(bquote('Bulk density (g ' ~cm^-3~')')) +
  xlab(bquote('Depth (cm)')) +
  scale_x_discrete(breaks=c("10", "20", "30"),
                   labels=c("0-10", "10-20", "20-30")) +
  theme(axis.text.x = element_text(colour="black",size=14)) +
  theme(axis.text.y = element_text(colour="black",size=14)) +
  theme(axis.title.y = element_text(colour="black",size=16)) +
  theme(axis.title.x = element_text(colour="black",size=16)) +
  theme(strip.text.y = element_text(colour="black",size=14)) +
  theme(legend.text = element_text(colour="black",size=14)) +
  theme(legend.title = element_text(colour="black",size=14)) +
  theme(axis.ticks.length=unit(0.2,"cm"))

BP1

# Carbono
BP2 <- ggplot(data=carbono_data2, aes(x=factor(Profundidad), y=Carbono, fill=Cobertura)) + geom_boxplot() +
  theme_classic() + 
  ylab(bquote('Carbon (%)')) +
  xlab(bquote('Depth (cm)')) +
  scale_x_discrete(breaks=c("10", "20", "30"),
                   labels=c("0-10", "10-20", "20-30")) +
  theme(axis.text.x = element_text(colour="black",size=14)) +
  theme(axis.text.y = element_text(colour="black",size=14)) +
  theme(axis.title.y = element_text(colour="black",size=16)) +
  theme(axis.title.x = element_text(colour="black",size=16)) +
  theme(strip.text.y = element_text(colour="black",size=14)) +
  theme(legend.text = element_text(colour="black",size=14)) +
  theme(legend.title = element_text(colour="black",size=14)) +
  theme(axis.ticks.length=unit(0.2,"cm"))

BP2

# Textura arena
carbono_data2$Arena <- as.numeric(carbono_data2$Arena)
carbono_data2$Limo <- as.numeric(carbono_data2$Limo)
carbono_data2$Arcilla <- as.numeric(carbono_data2$Arcilla)


BP3 <- ggplot(data=carbono_data2, aes(x=factor(Profundidad), y=Arena, fill=Cobertura )) + geom_boxplot() +
 theme_classic() + 
  ylab(bquote('Sand (%)')) +
  xlab(bquote('Depth (cm)')) +
  scale_x_discrete(breaks=c("10", "20", "30")) +
  scale_y_continuous(limits=c(0, 100)) +
  theme(axis.text.x = element_text(colour="white",size=14)) +
  theme(axis.text.y = element_text(colour="black",size=14)) +
  theme(axis.title.y = element_text(colour="black",size=16)) +
  theme(axis.title.x = element_text(colour="white",size=16)) +
  theme(strip.text.y = element_text(colour="black",size=14)) +
  theme(legend.text = element_text(colour="black",size=14)) +
  theme(legend.position = "none") +
  #theme(axis.text.x = element_blank() ) +
  #theme(axis.title.x = element_blank() ) +
  theme(axis.ticks.length=unit(0.2,"cm") )

BP4 <- ggplot(data=carbono_data2, aes(x=factor(Profundidad), y=Limo, fill=Cobertura )) + geom_boxplot() +
 theme_classic() + 
  ylab(bquote('Silt (%)')) +
  xlab(bquote('Depth (cm)')) +
  scale_x_discrete(breaks=c("10", "20", "30")) +
  scale_y_continuous(limits=c(0, 100)) +
  theme(axis.text.x = element_text(colour="white",size=14)) +
  theme(axis.text.y = element_text(colour="black",size=14)) +
  theme(axis.title.y = element_text(colour="black",size=16)) +
  theme(axis.title.x = element_text(colour="white",size=16)) +
  theme(strip.text.y = element_text(colour="black",size=14)) +
  theme(legend.text = element_text(colour="black",size=14)) +
  theme(legend.position = "none") +
  #theme(axis.text.x = element_blank() ) +
  #theme(axis.title.x = element_blank() ) +
  theme(axis.ticks.length=unit(0.2,"cm") )

BP5 <- ggplot(data=carbono_data2, aes(x=factor(Profundidad), y=Arcilla, fill=Cobertura )) + geom_boxplot() +
 theme_classic() + 
  ylab(bquote('Clay (%)')) +
  xlab(bquote('Depth (cm)')) +
  scale_x_discrete(breaks=c("10", "20", "30"),
                   labels=c("0-10", "10-20", "20-30")) +
  scale_y_continuous(limits=c(0, 100)) +
  theme(axis.text.x = element_text(colour="black",size=14)) +
  theme(axis.text.y = element_text(colour="black",size=14)) +
  theme(axis.title.y = element_text(colour="black",size=16)) +
  theme(axis.title.x = element_text(colour="black",size=16)) +
  theme(strip.text.y = element_text(colour="black",size=14)) +
  theme(legend.text = element_text(colour="black",size=14)) +
  theme(legend.title = element_text(colour="black",size=16)) +
  theme(legend.box = "horizontal") +
  theme(legend.position=c(0.08,0.75))  +
  theme(axis.ticks.length=unit(0.2,"cm"))

ggplot2.multiplot(BP3,BP4,BP5, cols=1)

Dung Beattle - Soil compactation

is there any relation?

# make categorical
# carbono_data$Profundidad <- to_factor(carbono_data$Profundidad)


ggplot(carbono_data2, aes(x=Densidad, y=Carbono)) + 
    geom_point (alpha = 1/5) + 
    geom_smooth (method = lm) + geom_smooth (col="red") + facet_grid (Cobertura ~ Profundidad) 

ggplot(carbono_data2, aes(x=Densidad, y=Carbono )) + 
  geom_point (alpha = 1/5) + 
  geom_smooth (method = lm) +
  geom_text(data=carbono_data2, aes(x=Densidad, y=Carbono, label=Id)) +
 facet_grid (Cobertura ~ Profundidad) 

##################   
###### GLM ####### 
##################  

glm2 <- glm(Carbono ~ factor(Profundidad) + Cobertura, data=carbono_data2, family = 'gaussian')
glm3 <- glm(Carbono ~ Densidad + factor(Profundidad) + Cobertura, data=carbono_data2, family = 'gaussian')
glm4 <- glm(Carbono ~ Densidad + factor(Profundidad) + Cobertura + Arena, data=carbono_data2, family = 'gaussian')
glm5 <- glm(Carbono ~ Densidad + factor(Profundidad) + factor(Cobertura) + Arcilla, data=carbono_data2, family = 'gaussian')
glm6 <- glm(Carbono ~ Densidad + factor(Profundidad) + factor(Cobertura) + Limo, data=carbono_data2, family = 'gaussian')
glm7 <- glm(Carbono ~ Arcilla + Limo + Arena, data=carbono_data2, family = 'gaussian')
glm8 <- glm(Carbono ~ Arena + factor(Profundidad) + Densidad, data=carbono_data2, family = 'gaussian')
glm9 <- glm(Carbono ~  Densidad, data=carbono_data2, family = 'gaussian')
glm10 <- glm(Carbono ~  Densidad + factor(Profundidad), data=carbono_data2, family = 'gaussian')


print("Best model is glm3")
## [1] "Best model is glm3"
AIC(glm2, glm3, glm9, glm10)
tab_model(glm3)
  Carbono
Predictors Estimates CI p
(Intercept) 9.33 8.64 – 10.01 <0.001
Densidad -2.75 -3.30 – -2.20 <0.001
factor(Profundidad)20 -1.00 -1.29 – -0.70 <0.001
factor(Profundidad)30 -1.16 -1.45 – -0.86 <0.001
CoberturaPastizal -0.71 -0.98 – -0.43 <0.001
CoberturaRastrojo -0.71 -1.31 – -0.12 0.019
CoberturaSilvopastoril 0.74 0.35 – 1.13 <0.001
Observations 1725
Cox & Snell’s R2 / Nagelkerke’s R2 0.668 / 0.668 
plot_model(glm3, type = "pred", terms = c("Densidad", "Cobertura", "Profundidad" ))

# intercept-only model with Carbono Mass as the response variable and a Gaussian error structure:

GLMM1 <- glmer (Carbono ∼ Densidad * Cobertura + (1|Profundidad), data=carbono_data2, family="gaussian")

GLMM2 <- glmer (Carbono ∼ Densidad * Profundidad + (1|Cobertura), data=carbono_data2, family="gaussian")

GLMM3 <- glmer(Carbono ∼ Densidad * Profundidad + (Densidad|Cobertura), data=carbono_data2, family="gaussian")
 
GLMM4 <- glmer(Carbono ∼ Densidad * Cobertura + (1|Cobertura) +  (1|Profundidad), data=carbono_data2, family="gaussian")
             
GLMM5 <- glmer(Carbono ∼ Densidad * factor(Profundidad) + (1|Cobertura) +  (1|Profundidad), data=carbono_data2, family="gaussian")
             

AIC(GLMM1, GLMM2, GLMM3, GLMM4, GLMM5)
anova(GLMM1, GLMM2, GLMM3, GLMM4, GLMM5)
####################
### Best models ####
####################

tab_model(GLMM1, 
          show.re.var=TRUE, 
          show.aic=TRUE)
  Carbono
Predictors Estimates CI p
(Intercept) 6.29 5.05 – 7.54 <0.001
Densidad -0.81 -1.64 – 0.02 0.056
Pastizal 4.24 2.68 – 5.79 <0.001
Rastrojo 5.94 2.57 – 9.32 0.001
Silvopastoril 2.45 0.37 – 4.53 0.021
Densidad:CoberturaPastizal -3.91 -5.12 – -2.70 <0.001
Densidad:CoberturaRastrojo -5.19 -7.70 – -2.68 <0.001
Densidad:CoberturaSilvopastoril -1.47 -3.11 – 0.18 0.080
Random Effects
σ2 6.19
τ00 Profundidad 0.42
ICC Profundidad 0.06
Observations 1725
Marginal R2 / Conditional R2 0.124 / 0.180
AIC 8064.598 
          #show.r2=TRUE,
          # transform)

tab_model(GLMM4, 
          show.re.var=TRUE, 
          show.aic=TRUE)
  Carbono
Predictors Estimates CI p
(Intercept) 6.29 5.05 – 7.54 <0.001
Densidad -0.81 -1.64 – 0.02 0.056
Pastizal 4.24 2.68 – 5.79 <0.001
Rastrojo 5.94 2.57 – 9.32 0.001
Silvopastoril 2.45 0.37 – 4.53 0.021
Densidad:CoberturaPastizal -3.91 -5.12 – -2.70 <0.001
Densidad:CoberturaRastrojo -5.19 -7.70 – -2.68 <0.001
Densidad:CoberturaSilvopastoril -1.47 -3.11 – 0.18 0.080
Random Effects
σ2 6.19
τ00 Cobertura 0.00
τ00 Profundidad 0.42
ICC Cobertura 0.00
ICC Profundidad 0.06
Observations 1725
Marginal R2 / Conditional R2 0.124 / 0.180
AIC 8066.598 
          #show.r2=TRUE,
          # transform)
#ggplot(glm2, aes(x=Densidad, y=Carbono)) + 
#  geom_point (alpha = 1/5) + 
#  geom_smooth (method = lm) #+
# # facet_grid (Cobertura ~ .) 

plot_model(GLMM1)

plot_model(GLMM1, type = "pred", terms = c("Densidad", "Cobertura", "Profundidad" , title="xxx"))

dat <- ggpredict(GLMM1, terms = c("Densidad", "Cobertura"))
plot(dat, facet = TRUE)

plot(dat, rawdata = TRUE, facet = TRUE)

glm3.3 <- glm3 <- glm(Carbono ~ Densidad + factor(Profundidad) + Cobertura, data=carbono_train, family = 'gaussian')

pred = predict(glm3.3, newdata=carbono_test)
# accuracy <- table(pred, carbono_test$Carbono)
# sum(diag(accuracy))/sum(accuracy)

#confusionMatrix(data=pred, carbono_test$Carbono)

Make sf

beetle_data_sf = st_as_sf(beetle_data, coords = c("LongDD", "LatDD"), crs = "+proj=longlat +datum=WGS84")


 # put coords

ants_data_sf = st_as_sf(ants_data, coords = c("LongDD", "LatDD"), crs = "+proj=longlat +datum=WGS84")

Plot Beetle Data and sampling points

Simple map

# library(tmap)
# qtm(beetle_data_sf)

mapview(beetle_data_sf, map.types = "OpenStreetMap")

Beetle Richness

# 
Beetle_rich_map <- mapview(beetle_data_sf, zcol = "Richness", 
        col.regions = sf.colors(12),
        alpha.regions = 0.2, alpha = 0.2, 
        legend = TRUE, 
        map.types = "Esri.WorldImagery")

Beetle_rich_map

Beetle Abundance

# 
Beetle_abun_map <- mapview(beetle_data_sf, zcol = "Abundancia", 
        col.regions = sf.colors(12),
        alpha.regions = 0.2, alpha = 0.2, 
        legend = TRUE, 
        map.types = "Esri.WorldImagery")

Beetle_abun_map
#sync(Beetle_rich_map, Beetle_abun_map)
# 
p <- ggplot(beetle_data, aes(x=Cobertura, y=Richness)) + 
  geom_boxplot() 
p

Plot Ant Data

Ant Richness

# 
ant_rich_map <- mapview(ants_data_sf, zcol = "Richness", 
        col.regions = sf.colors(12),
        alpha = 0.1, alpha.regions = 0.2,
        legend = TRUE, 
        map.types = "Esri.WorldImagery")

ant_rich_map

Ant Abundance

ant_abun_map <- mapview(ants_data_sf, zcol = "Abundancia", 
         col.regions = sf.colors(15),
         alpha = 0.1, alpha.regions = 0.2,
         legend = TRUE, 
         map.types = "Esri.WorldImagery")

ant_abun_map
# sync(ant_rich_map, ant_abun_map)

As Boxplot

# 
q <- ggplot(ants_data, aes(x=Cobertura, y=Abundancia)) + 
  geom_boxplot()  + scale_y_continuous(trans='log2')

q

Biodiversity - Soil compactation

1st join carbon table and biodiversity table and then standarize

# richness... shirink by punto
rich_beetle_data <- beetle_data %>% 
  group_by (Punto) %>% 
  # select ("Lat", "Lon", "Cobertura")
  summarise (Richness = n_distinct(Especie),
             Abundance = sum (Abundancia),
             Lat = unique (LatDD), 
             Lon = unique (LonDD),
             Cobertura= unique (Cobertura))  # New column richness


# sum carbono and average compactation data 
mean_carbono_data <- carbono_data2 %>% 
  group_by (Punto) %>% 
  summarise (Carbono = sum(Carbono),
             SoilCompactation= mean(Densidad)) 



# merge both data sets
Full_data <- merge(rich_beetle_data ,
                   mean_carbono_data,
                   by="Punto")


# Standarization
Full_data$Carbono_std <-  scale(Full_data[,7])
Full_data$SoilCompactation_std <-  scale(Full_data[,8])

Is diferent the compactation by cover type?

glm10.1 <- lm(Richness ~ Cobertura, data=Full_data, 
             family = 'poisson')
tab_model(glm10.1)
  Richness
Predictors Estimates CI p
(Intercept) 3.76 3.53 – 4.00 <0.001
CoberturaPastizal -1.35 -1.68 – -1.02 <0.001
CoberturaRastrojo -1.95 -2.68 – -1.21 <0.001
CoberturaSilvopastoril -1.68 -2.16 – -1.19 <0.001
Observations 589
R2 / adjusted R2 0.136 / 0.132 
anova(glm10.1)
a2 <- aov(Full_data$Richness ~ Full_data$Cobertura)
TukeyHSD(a2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Full_data$Richness ~ Full_data$Cobertura)
## 
## $`Full_data$Cobertura`
##                              diff        lwr        upr     p adj
## Pastizal-Bosque        -1.3501456 -1.7817626 -0.9185286 0.0000000
## Rastrojo-Bosque        -1.9456894 -2.9144624 -0.9769164 0.0000019
## Silvopastoril-Bosque   -1.6783124 -2.3165755 -1.0400493 0.0000000
## Rastrojo-Pastizal      -0.5955438 -1.5617572  0.3706697 0.3862302
## Silvopastoril-Pastizal -0.3281668 -0.9625383  0.3062047 0.5423433
## Silvopastoril-Rastrojo  0.2673770 -0.8071718  1.3419257 0.9186284

YES!!!! ***

Is there any relation Beattle Richness and compactation?

# make categorical
# carbono_data$Profundidad <- to_factor(carbono_data$Profundidad)


##################   
###### GLM ####### 
##################  

glm21 <- glm(Richness ~ SoilCompactation, data=Full_data, 
             family = 'poisson')

glm22 <- glm(Richness ~ SoilCompactation + Cobertura,
             data=Full_data, family = 'poisson')

glm23 <- glm(Richness ~ SoilCompactation * Cobertura,
             data=Full_data, family = 'poisson')


print("Best model is glm23")
## [1] "Best model is glm23"
AIC(glm21, glm22, glm23)
# summary(glm13)
tab_model(glm23, show.aic=TRUE, transform=NULL)
  Richness
Predictors Log-Mean CI p
(Intercept) 2.38 1.89 – 2.86 <0.001
Soil Compactation -0.89 -1.31 – -0.48 <0.001
CoberturaPastizal -1.19 -2.03 – -0.35 0.005
CoberturaRastrojo -1.02 -2.98 – 0.93 0.306
CoberturaSilvopastoril -1.49 -2.77 – -0.20 0.024
SoilCompactation:CoberturaPastizal 0.66 0.01 – 1.32 0.048
SoilCompactation:CoberturaRastrojo 0.32 -1.16 – 1.80 0.671
SoilCompactation:CoberturaSilvopastoril 0.77 -0.26 – 1.80 0.144
Observations 589
Cox & Snell’s R2 / Nagelkerke’s R2 0.196 / 0.278
AIC 2233.436 
plot_model(glm23, type = "pred", terms = c("SoilCompactation", "Cobertura"))

plot_model(glm23)

# 
# GLMM6 <- glmer(Richness ∼ SoilCompactation * Cobertura + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM7 <- glmer(Richness ∼ SoilCompactation + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM8 <- glmer(Richness ∼ SoilCompactation + Cobertura + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM9 <- glmer(Richness ∼ SoilCompactation * Cobertura + (1|Cobertura) + (1|SoilCompactation) , data=Full_data, family = 'poisson')
# 
# AIC(GLMM6, GLMM7, GLMM8, GLMM9)

dat <- ggpredict(glm23, terms = c("SoilCompactation", "Cobertura"))
plot(dat, facet = TRUE, rawdata = TRUE)

plot(dat, facet = TRUE)

> YES!!!! ***

Is there any relation Beattle Abundance and compactation?

# make categorical
# carbono_data$Profundidad <- to_factor(carbono_data$Profundidad)


##################   
###### GLM ####### 
##################  

glm31 <- glm(Abundance ~ SoilCompactation, data=Full_data, 
             family = 'poisson')

glm32 <- glm(Abundance ~ SoilCompactation + Cobertura,
             data=Full_data, family = 'poisson')

glm33 <- glm(Abundance ~ SoilCompactation * Cobertura,
             data=Full_data, family = 'poisson')


print("Best model is glm33")
## [1] "Best model is glm33"
AIC(glm31, glm32, glm33)
# summary(glm13)
tab_model(glm33, show.aic=TRUE, transform=NULL)
  Abundance
Predictors Log-Mean CI p
(Intercept) 4.01 3.72 – 4.30 <0.001
Soil Compactation -1.44 -1.69 – -1.19 <0.001
CoberturaPastizal -2.23 -2.68 – -1.79 <0.001
CoberturaRastrojo -1.03 -2.44 – 0.38 0.153
CoberturaSilvopastoril -3.19 -4.05 – -2.33 <0.001
SoilCompactation:CoberturaPastizal 1.88 1.53 – 2.23 <0.001
SoilCompactation:CoberturaRastrojo 0.07 -1.02 – 1.16 0.899
SoilCompactation:CoberturaSilvopastoril 2.04 1.36 – 2.71 <0.001
Observations 589
Cox & Snell’s R2 / Nagelkerke’s R2 0.608 / 0.608
AIC 8007.194 
#############
### Riqueza y abundancia
############
tab_model(glm23, glm33, show.aic=TRUE, transform=NULL)
  Richness Abundance
Predictors Log-Mean CI p Log-Mean CI p
(Intercept) 2.38 1.89 – 2.86 <0.001 4.01 3.72 – 4.30 <0.001
Soil Compactation -0.89 -1.31 – -0.48 <0.001 -1.44 -1.69 – -1.19 <0.001
CoberturaPastizal -1.19 -2.03 – -0.35 0.005 -2.23 -2.68 – -1.79 <0.001
CoberturaRastrojo -1.02 -2.98 – 0.93 0.306 -1.03 -2.44 – 0.38 0.153
CoberturaSilvopastoril -1.49 -2.77 – -0.20 0.024 -3.19 -4.05 – -2.33 <0.001
SoilCompactation:CoberturaPastizal 0.66 0.01 – 1.32 0.048 1.88 1.53 – 2.23 <0.001
SoilCompactation:CoberturaRastrojo 0.32 -1.16 – 1.80 0.671 0.07 -1.02 – 1.16 0.899
SoilCompactation:CoberturaSilvopastoril 0.77 -0.26 – 1.80 0.144 2.04 1.36 – 2.71 <0.001
Observations 589 589
Cox & Snell’s R2 / Nagelkerke’s R2 0.196 / 0.278 0.608 / 0.608
AIC 2233.436 8007.194 
plot_model(glm33, type = "pred", terms = c("SoilCompactation", "Cobertura"))

plot_model(glm33)

# 
# GLMM6 <- glmer(Richness ∼ SoilCompactation * Cobertura + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM7 <- glmer(Richness ∼ SoilCompactation + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM8 <- glmer(Richness ∼ SoilCompactation + Cobertura + (1|Cobertura) , data=Full_data, family = 'poisson')
# 
# GLMM9 <- glmer(Richness ∼ SoilCompactation * Cobertura + (1|Cobertura) + (1|SoilCompactation) , data=Full_data, family = 'poisson')
# 
# AIC(GLMM6, GLMM7, GLMM8, GLMM9)

dat3 <- ggpredict(glm33, terms = c("SoilCompactation", "Cobertura"))
plot(dat3, facet = TRUE)#, rawdata = TRUE)

> YES!!!! ***

Tarea Incorporar evaluacion de modelos AUC con 75% 25%

Session Info

Details and pakages used

print(sessionInfo(), locale = FALSE)
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
## 
## Matrix products: default
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] caret_6.0-81           lattice_0.20-38        ggeffects_0.8.0       
##  [4] lme4_1.1-20            Matrix_1.2-15          easyGgplot2_1.0.0.9000
##  [7] rsq_1.1                fields_9.6             maps_3.3.0            
## [10] spam_2.2-1             dotCall64_1.0-0        sqldf_0.4-11          
## [13] RSQLite_2.1.1          gsubfn_0.7             proto_1.0.0           
## [16] sjmisc_2.7.7.9000      sjPlot_2.6.2.9000      mapview_2.6.3         
## [19] sf_0.7-2               readxl_1.3.0           lubridate_1.7.4       
## [22] ggmap_3.0.0            rgdal_1.3-9            rgeos_0.4-2           
## [25] raster_2.8-19          sp_1.3-1               forcats_0.4.0         
## [28] stringr_1.4.0          dplyr_0.8.0.1          purrr_0.3.0           
## [31] readr_1.3.1            tidyr_0.8.3            tibble_2.0.1          
## [34] ggplot2_3.1.0          tidyverse_1.2.1       
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.3    plyr_1.8.4         lazyeval_0.2.1    
##   [4] TMB_1.7.15         splines_3.5.2      crosstalk_1.0.0   
##   [7] leaflet_2.0.2      TH.data_1.0-10     digest_0.6.18     
##  [10] foreach_1.4.4      htmltools_0.3.6    magrittr_1.5      
##  [13] memoise_1.1.0      recipes_0.1.4      gower_0.2.0       
##  [16] modelr_0.1.4       sandwich_2.5-0     jpeg_0.1-8        
##  [19] colorspace_1.4-0   blob_1.1.1         rvest_0.3.2       
##  [22] haven_2.1.0        xfun_0.5           tcltk_3.5.2       
##  [25] crayon_1.3.4       jsonlite_1.6       survival_2.43-3   
##  [28] zoo_1.8-4          iterators_1.0.10   glue_1.3.0        
##  [31] gtable_0.2.0       ipred_0.9-8        emmeans_1.3.2     
##  [34] webshot_0.5.1      sjstats_0.17.3     scales_1.0.0      
##  [37] mvtnorm_1.0-8      DBI_1.0.0          Rcpp_1.0.0        
##  [40] viridisLite_0.3.0  xtable_1.8-3       units_0.6-2       
##  [43] foreign_0.8-71     bit_1.1-14         lava_1.6.5        
##  [46] prodlim_2018.04.18 stats4_3.5.2       prediction_0.3.6.2
##  [49] htmlwidgets_1.3    httr_1.4.0         RColorBrewer_1.1-2
##  [52] modeltools_0.2-22  pkgconfig_2.0.2    nnet_7.3-12       
##  [55] labeling_0.3       reshape2_1.4.3     tidyselect_0.2.5  
##  [58] rlang_0.3.1        later_0.7.5        munsell_0.5.0     
##  [61] cellranger_1.1.0   tools_3.5.2        cli_1.0.1         
##  [64] generics_0.0.2     sjlabelled_1.0.16  broom_0.5.1       
##  [67] ggridges_0.5.1     evaluate_0.13      yaml_2.2.0        
##  [70] ModelMetrics_1.2.2 knitr_1.21         bit64_0.9-7       
##  [73] satellite_1.0.1    RgoogleMaps_1.4.3  coin_1.2-2        
##  [76] nlme_3.1-137       mime_0.6           xml2_1.2.0        
##  [79] compiler_3.5.2     bayesplot_1.6.0    rstudioapi_0.9.0  
##  [82] png_0.1-7          e1071_1.7-0.1      stringi_1.3.1     
##  [85] classInt_0.3-1     psych_1.8.12       nloptr_1.2.1      
##  [88] pillar_1.3.1       pwr_1.2-2          estimability_1.3  
##  [91] data.table_1.12.0  bitops_1.0-6       httpuv_1.4.5.1    
##  [94] R6_2.4.0           promises_1.0.1     codetools_0.2-16  
##  [97] MASS_7.3-51.1      assertthat_0.2.0   chron_2.3-53      
## [100] rjson_0.2.20       withr_2.1.2        mnormt_1.5-5      
## [103] multcomp_1.4-8     parallel_3.5.2     hms_0.4.2         
## [106] rpart_4.1-13       timeDate_3043.102  coda_0.19-2       
## [109] glmmTMB_0.2.3      class_7.3-15       minqa_1.2.4       
## [112] rmarkdown_1.11     snakecase_0.9.2    shiny_1.2.0       
## [115] base64enc_0.1-3

Cited Literature