Exploring Biodiversity, Carbon and Soil compactation. A spatially explicit aproach

knitr::opts_chunk$set(echo = TRUE)

To explore the Biodiversity role in the Carbon - Soil compactation correlation found by Diego N. at larger scale.

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

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

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 o All data; o In each cover; o In each cover and depth. • SOC + bulk density + clay, silt, sand o All data o In each cover; o In each cover and depth.

“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

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.

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 an ANOVA

glm1 <- glm(Carbono ~ factor(Profundidad), data=carbono_data2)
anova (glm1)

summary(glm1)

## 
## Call:
## glm(formula = Carbono ~ factor(Profundidad), data = carbono_data2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.6195  -1.9235  -0.3706   1.6909   9.5668  
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             5.6965     0.1119  50.899  < 2e-16 ***
## factor(Profundidad)20  -1.1086     0.1575  -7.038  2.8e-12 ***
## factor(Profundidad)30  -1.3492     0.1574  -8.573  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 7.102021)
## 
##     Null deviance: 12821  on 1724  degrees of freedom
## Residual deviance: 12230  on 1722  degrees of freedom
## AIC: 8282
## 
## Number of Fisher Scoring iterations: 2

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)
summary (glm1.2)

## 
## Call:
## glm(formula = Densidad ~ factor(Profundidad), data = carbono_data2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.77790  -0.15265   0.01518   0.16903   0.58590  
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.226904   0.009501 129.129  < 2e-16 ***
## factor(Profundidad)20 0.047898   0.013373   3.582 0.000351 ***
## factor(Profundidad)30 0.072532   0.013361   5.428 6.49e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.05118666)
## 
##     Null deviance: 89.700  on 1724  degrees of freedom
## Residual deviance: 88.143  on 1722  degrees of freedom
## AIC: -226.84
## 
## Number of Fisher Scoring iterations: 2

a3 <- aov(carbono_data2$Carbono ~ factor(carbono_data2$Cobertura))
TukeyHSD(a3)

##   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
glm1.3<- glm(Carbono ~ factor(Cobertura), data=carbono_data2)
summary (glm1.3)

## 
## Call:
## glm(formula = Carbono ~ factor(Cobertura), data = carbono_data2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.8028  -1.9378  -0.3588   1.5422   9.1843  
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      5.3207     0.1004  52.974  < 2e-16 ***
## factor(Cobertura)Pastizal       -1.0893     0.1402  -7.767 1.37e-14 ***
## factor(Cobertura)Rastrojo       -1.1244     0.3170  -3.547   0.0004 ***
## factor(Cobertura)Silvopastoril   0.5391     0.2082   2.589   0.0097 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 7.051855)
## 
##     Null deviance: 12821  on 1724  degrees of freedom
## Residual deviance: 12136  on 1721  degrees of freedom
## AIC: 8270.8
## 
## Number of Fisher Scoring iterations: 2

a2 <- aov(carbono_data2$Densidad ~ factor(carbono_data2$Profundidad))
TukeyHSD(a2)

##   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

Si hay diferencias!

Grafica para ilustrar resultados de las anovas

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"))

# 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"))


# 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)

Carbon content - Soil compactation

is there any relation?

# make categorical
# carbono_data$Profundidad <- to_factor(carbono_data$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')

# 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 <- lmer(Carbono ∼ Densidad + (1|Cobertura), data=carbono_data2)
              
glm5 <- glm(Carbono ~ Densidad + factor(Profundidad) + factor(Cobertura) + Arcilla, data=carbono_data2, family = "gaussian")


glm2.1 <- glm(Carbono ~ Densidad, data=carbono_data)
rsq(glm2.1)

## [1] 0.07417259

#ggplot(glm2, aes(x=Densidad, y=Carbono)) + 
#  geom_point (alpha = 1/5) + 
#  geom_smooth (method = lm) #+
# # facet_grid (Cobertura ~ .) 


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) 

# wiew as sjPlot  
plot_model(glm2, type = "est", terms = c("Cobertura", "Profundidad"))

summary(glm2)

## 
## Call:
## glm(formula = Carbono ~ factor(Profundidad) + Cobertura, family = "gaussian", 
##     data = carbono_data2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.0462  -1.9021  -0.3685   1.5236   8.5633  
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              6.1635     0.1328  46.420  < 2e-16 ***
## factor(Profundidad)20   -1.1356     0.1531  -7.418 1.86e-13 ***
## factor(Profundidad)30   -1.3615     0.1529  -8.902  < 2e-16 ***
## CoberturaPastizal       -1.1033     0.1368  -8.066 1.35e-15 ***
## CoberturaRastrojo       -1.1174     0.3091  -3.614  0.00031 ***
## CoberturaSilvopastoril   0.5487     0.2031   2.702  0.00695 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.706248)
## 
##     Null deviance: 12821  on 1724  degrees of freedom
## Residual deviance: 11528  on 1719  degrees of freedom
## AIC: 8186.1
## 
## Number of Fisher Scoring iterations: 2

YES!!!! ***

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

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

First model to explain biodiversity in function of carbon and soil compactation. A GLM Poisson model.

Model Algebra for species richness… Simpler model averaging the three deepts.

The N-mixture model (Royle, 2004b) is arguably the Hierarquical model for animal abundance. Here was adapted to species richness.

\(N_i \sim Poisson (\lambda_i)\)

Where: \(N_i\) is the species richness at site \(i\). Which follows a Poisson distribution (count of species). As we are interested in modeling the effect of measurable covariates such as Carbono and Soil compactation, we consider the the log-transformation \(log(\lambda_i)\) to get a linear model with covariates.

\[log(\lambda_i) = \alpha + \beta_1* Carbono + \beta_2* SoilCompactation + \beta_3 *Carbono*SoilCompactation\]
where \(\alpha\) is the intercept and \(\beta_1, \beta_2, \beta_3\) are the coeficients of the covariates.

Coded as GLM in R it has the form:

\[Richness \sim Carbono+SoilCompactation+Carbono*SoilCompactation\] From the Poisson family.

As the data comes from diferent units, it needs to be standarizated.

# 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])

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


glm13 <- glm(Richness ~  Carbono + SoilCompactation * Carbono * SoilCompactation, data = Full_data, family = poisson())

glm14 <- glm(Richness ~  Carbono + SoilCompactation + factor(Cobertura), data = Full_data, family = poisson())

glm15 <- glm(Richness ~  Carbono + SoilCompactation * factor(Cobertura), data = Full_data, family = poisson())

glm16 <- glm(Richness ~  Carbono + SoilCompactation : factor(Cobertura), data = Full_data, family = poisson())

glm17 <- glm(Richness ~  SoilCompactation + Carbono : factor(Cobertura), data = Full_data, family = poisson())

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

# glm3 <- glm(Richness ~  Carbono_std + 
#               SoilCompactation_std * Carbono_std * SoilCompactation_std, data = Full_data, family = poisson())

# glm5 <- glm(Richness ~  I(Carbono_std)^2 + SoilCompactation_std * I(Carbono_std)^2 * SoilCompactation_std, data = Full_data, family = poisson())
# 
# glm6 <- glm(Richness ~  Carbono_std + I(SoilCompactation_std)^2 * Carbono_std * I(SoilCompactation_std)^2, data = Full_data, family = poisson())
# 
 AIC(glm13, glm14, glm15, glm16, glm17, glm18)

summary(glm14)

## 
## Call:
## glm(formula = Richness ~ Carbono + SoilCompactation + factor(Cobertura), 
##     family = poisson(), data = Full_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0359  -0.8773  -0.1953   0.5400   3.7980  
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     2.321469   0.226800  10.236  < 2e-16 ***
## Carbono                        -0.009483   0.004242  -2.235   0.0254 *  
## SoilCompactation               -0.718591   0.162593  -4.420 9.89e-06 ***
## factor(Cobertura)Pastizal      -0.370008   0.057236  -6.465 1.02e-10 ***
## factor(Cobertura)Rastrojo      -0.653946   0.148429  -4.406 1.05e-05 ***
## factor(Cobertura)Silvopastoril -0.530631   0.088784  -5.977 2.28e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 717.53  on 588  degrees of freedom
## Residual deviance: 588.86  on 583  degrees of freedom
## AIC: 2229.3
## 
## Number of Fisher Scoring iterations: 5

# wiew as sjPlot  
# plot_model(glm14, type = "pred", terms = c("Carbono"))

# plot_model(glm14, type = "pred", terms = c("SoilCompactation"))

# plot_model(glm14, type = "pred", terms = c("SoilCompactation", "Cobertura"))

Tarea Incorporar evaluacion de modelos AUC con 75% 25%

It does not have interaction. However we can represent the relation as a surface.

Biodiversity - Richness explained as SoilCompactation and Carbono function.

# Compute expected data 
Carbono <- seq(min(Full_data$Carbono), 
               max(Full_data$Carbono),length.out=100)                # Values of cov 2
SoilCompactation <- seq(min(Full_data$SoilCompactation),  
               max(Full_data$SoilCompactation),length.out=100) 
Cobertura <- rep(unique(Full_data$Cobertura), length.out=100)

# build a mat using library(sqldf)
new_mat <- as.data.frame(expand.grid(Carbono,SoilCompactation, Cobertura)) 
colnames(new_mat) <- c("Carbono", "SoilCompactation", "Cobertura" )

psi.matrix <- array(NA, dim = c(100, 100, 4)) # Prediction matrix, for every

for(i in 1:100){
  for(j in 1:100){
    for(k in 1:4){
       psi.matrix[i, j, k]<-predict(glm14, newdata=data.frame(
          Carbono=Carbono[i],
         # mean=pr.mat$mean[j]), 
        #  range=pr.mat$range[j]), 
          SoilCompactation=SoilCompactation[j],
          Cobertura=Cobertura[k]),
          type="response")
    } 
  }
}


mapPalette <- colorRampPalette(c("blue", "yellow", "orange", "red"))


#Bosque
image.plot(x = Carbono, y = SoilCompactation, z = psi.matrix[,,1], col = mapPalette(100), xlab = "Carbon", 
      ylab = "SoilCompactation", cex.lab = 1.2, legend.width=1, main="Bosque")
contour(x = Carbono, y = SoilCompactation, z = psi.matrix[,,1], add = TRUE, lwd = 1)
matpoints(Full_data$Carbono, Full_data$SoilCompactation, pch="+", cex=0.7, col = "black")

# Pastizal
image.plot(x = Carbono, y = SoilCompactation, z = psi.matrix[,,2], col = mapPalette(100), xlab = "Carbon", 
      ylab = "SoilCompactation", cex.lab = 1.2, legend.width=1, main="Pastizal")
contour(x = Carbono, y = SoilCompactation, z = psi.matrix[,,2], add = TRUE, lwd = 1)
matpoints(Full_data$Carbono, Full_data$SoilCompactation, pch="+", cex=0.7, col = "black")

# Silvopastoril
image.plot(x = Carbono, y = SoilCompactation, z = psi.matrix[,,3], col = mapPalette(100), xlab = "Carbon", 
      ylab = "SoilCompactation", cex.lab = 1.2, legend.width=1, main="Silvopastoril")
contour(x = Carbono, y = SoilCompactation, z = psi.matrix[,,3], add = TRUE, lwd = 1)
matpoints(Full_data$Carbono, Full_data$SoilCompactation, pch="+", cex=0.7, col = "black")

# Rastrojo
image.plot(x = Carbono, y = SoilCompactation, z = psi.matrix[,,4], col = mapPalette(100), xlab = "Carbon", 
      ylab = "SoilCompactation", cex.lab = 1.2, legend.width=1, main="Rastrojo")
contour(x = Carbono, y = SoilCompactation, z = psi.matrix[,,4], add = TRUE, lwd = 1)
matpoints(Full_data$Carbono, Full_data$SoilCompactation, pch="+", cex=0.7, col = "black")

# 
# ########################################################
# ###### Cross-validation for Generalized Linear Models
# #########################################################
# library("boot")
# cost <- function(r, pi = 0) mean(abs(r - pi) > 0.5)  ## cost function necessary for binomial data
# m11.cv <- cv.glm(data = corales, glm14, cost, K = 10)  # use leave-one-out cross validation (can use K-fold cross validation for larger data sets)
# ## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# 
# # Now lets see what our error rate was:
# 
# m11.cv$delta
# ## [1] 0.2381 0.2438
# # That’s not too bad.
# 
# muhat <- fitted(glm14)
# glm14.diag <- glm.diag(glm14)
# (cv.err <- mean((glm14$y - muhat)^2/(1 - glm14.diag$h)^2))
# 
# glm.diag.plots(glm14, glm14.diag)

#####################################################
## Make your reciever-operater curve
#####################################################

# library(pROC)
# 
# 
# m.roc <- multiclass.roc(corales$corales_algunos_sum, predict(glm8, backtransform = TRUE))
# auc(m.roc)
# ci(m.roc)
# plot(m.roc[[1]], m.roc[[2]])
# m.roc[[7]]

Coded as BUGS language

for(i in 1:sites){ 
$Richness[i] ~ pois(lambda[i])$ 
$lambda[i] <- alpha + beta[1] * Carbono[i] + beta[2] * SoilCompactation[i] + beta[3] * Carbono[i] * SoilCompactation[i] $
 }


### Coded as BUGS language and considering deept variation

for(i in 1:sites){ 
$Richness[i] ~ pois(lambda[i])$ 
$lambda[i] <- alpha[deept[i]] + beta[1] * Carbono[i] + beta[2] * SoilCompactation[i] + beta[3] * Carbono[i] * SoilCompactation[i] $
 }

Session Info

Details and pakages used

print(sessionInfo(), locale = FALSE)

## R version 3.4.0 (2017-04-21)
## 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] bindrcpp_0.2.2         lme4_1.1-18-1          Matrix_1.2-14         
##  [4] easyGgplot2_1.0.0.9000 rsq_1.0.1              fields_9.6            
##  [7] maps_3.3.0             spam_2.2-0             dotCall64_1.0-0       
## [10] sqldf_0.4-11           RSQLite_2.1.1          gsubfn_0.7            
## [13] proto_1.0.0            sjmisc_2.7.4           sjPlot_2.6.0          
## [16] mapview_2.6.0          sf_0.6-3               readxl_1.1.0          
## [19] lubridate_1.7.4        ggmap_2.7.904          rgdal_1.3-4           
## [22] rgeos_0.3-28           raster_2.6-7           sp_1.3-1              
## [25] forcats_0.3.0          stringr_1.3.1          dplyr_0.7.6           
## [28] purrr_0.2.5            readr_1.1.1            tidyr_0.8.1           
## [31] tibble_1.4.2           ggplot2_3.0.0          tidyverse_1.2.1       
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.2    plyr_1.8.4         lazyeval_0.2.1    
##   [4] TMB_1.7.14         splines_3.4.0      crosstalk_1.0.0   
##   [7] leaflet_2.0.2      TH.data_1.0-9      digest_0.6.17     
##  [10] htmltools_0.3.6    magrittr_1.5       memoise_1.1.0     
##  [13] modelr_0.1.2       sandwich_2.5-0     htmldeps_0.1.1    
##  [16] jpeg_0.1-8         colorspace_1.3-2   blob_1.1.1        
##  [19] rvest_0.3.2        haven_1.1.2        tcltk_3.4.0       
##  [22] crayon_1.3.4       jsonlite_1.6       bindr_0.1.1       
##  [25] survival_2.41-3    zoo_1.8-3          glue_1.3.0        
##  [28] gtable_0.2.0       emmeans_1.2.3      webshot_0.5.0     
##  [31] sjstats_0.17.0     scales_1.0.0       mvtnorm_1.0-8     
##  [34] DBI_1.0.0          ggeffects_0.5.0    Rcpp_1.0.0        
##  [37] viridisLite_0.3.0  xtable_1.8-3       spData_0.2.9.3    
##  [40] units_0.6-0        foreign_0.8-67     bit_1.1-14        
##  [43] stats4_3.4.0       prediction_0.3.6   survey_3.33-2     
##  [46] htmlwidgets_1.2    httr_1.3.1         modeltools_0.2-22 
##  [49] pkgconfig_2.0.2    nnet_7.3-12        labeling_0.3      
##  [52] reshape2_1.4.3     tidyselect_0.2.4   rlang_0.3.0.1     
##  [55] later_0.7.4        munsell_0.5.0      cellranger_1.1.0  
##  [58] tools_3.4.0        cli_1.0.1          sjlabelled_1.0.13 
##  [61] broom_0.5.0        ggridges_0.5.0     evaluate_0.11     
##  [64] yaml_2.2.0         knitr_1.20         bit64_0.9-7       
##  [67] satellite_1.0.1    RgoogleMaps_1.4.2  coin_1.2-2        
##  [70] nlme_3.1-137       mime_0.6           xml2_1.2.0        
##  [73] compiler_3.4.0     bayesplot_1.6.0    rstudioapi_0.7    
##  [76] png_0.1-7          e1071_1.7-0        stringi_1.2.4     
##  [79] lattice_0.20-35    classInt_0.2-3     psych_1.8.4       
##  [82] nloptr_1.0.4       effects_4.0-2      stringdist_0.9.5.1
##  [85] pillar_1.3.1       pwr_1.2-2          estimability_1.3  
##  [88] data.table_1.11.8  bitops_1.0-6       httpuv_1.4.5      
##  [91] R6_2.3.0           promises_1.0.1     codetools_0.2-15  
##  [94] MASS_7.3-47        assertthat_0.2.0   chron_2.3-52      
##  [97] rjson_0.2.20       withr_2.1.2        mnormt_1.5-5      
## [100] multcomp_1.4-8     parallel_3.4.0     hms_0.4.2         
## [103] coda_0.19-1        glmmTMB_0.2.2.0    class_7.3-14      
## [106] minqa_1.2.4        rmarkdown_1.10.12  snakecase_0.9.2   
## [109] carData_3.0-1      shiny_1.1.0        base64enc_0.1-3

Cited Literature