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.

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.

Draft Idea

draft

Idea Evolution

The initial idea changed to sample just ants and dung bettles in a 80 m grid. The Carbon data came in September 2018.

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

Read Data

For Dung beetle and ants

path <- "C:/Users/diego.lizcano/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

For carbono

carbono_data <- read_excel(paste(path,"/CARBONO-06-09-2018.xlsx", sep=""), 
    sheet = "Densidad-carbono")

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_data$Carbono_percent,  na.rm = T)
mediana <-median(carbono_data$Carbono_percent, na.rm = T)

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

c1 <- ggplot (data=carbono_data, aes(Carbono_percent) ) +
  geom_histogram (color = "white", bins = 100) +
  geom_density(color = "sky blue") +
  geom_rug (aes(Carbono_percent) ) +
  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

As boxplot

c2 <- ggplot (data=carbono_data, aes(x=factor(Profundidad), y=Carbono_percent)) +
  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 <- lm(Carbono_percent ~ factor(Profundidad), data=carbono_data)
anova (glm1)

summary(glm1)

## 
## Call:
## lm(formula = Carbono_percent ~ factor(Profundidad), data = carbono_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.06146 -0.02169 -0.00579  0.01593  0.40659 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.061464   0.001412  43.536  < 2e-16 ***
## factor(Profundidad)20 -0.013511   0.001990  -6.790 1.53e-11 ***
## factor(Profundidad)30 -0.017307   0.001991  -8.694  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03435 on 1788 degrees of freedom
## Multiple R-squared:  0.04456,    Adjusted R-squared:  0.04349 
## F-statistic:  41.7 on 2 and 1788 DF,  p-value: < 2.2e-16

Si hay diferencias!

Carbon content - Soil compactation

is there any relation?

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

glm2 <- glm(Carbono_percent ~ Densidad_g + Profundidad, data=carbono_data)

ggplot(glm2, aes(x=Densidad_g, y=Carbono_percent)) + 
  geom_point (alpha = 1/5) + 
  geom_smooth (method = lm) #+ 

  #facet_grid (Profundidad ~ .) 

# wiew as sjPlot  
plot_model(glm2, type = "pred", terms = c("Densidad_g", "Profundidad"))

summary(glm2)

## 
## Call:
## glm(formula = Carbono_percent ~ Densidad_g + Profundidad, data = carbono_data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.08362  -0.02058  -0.00627   0.01314   0.40912  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.090359   0.004273  21.145  < 2e-16 ***
## Densidad_g    -0.023707   0.003315  -7.152 1.24e-12 ***
## Profundidad20 -0.012132   0.001972  -6.152 9.43e-10 ***
## Profundidad30 -0.015292   0.001984  -7.709 2.08e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.001147769)
## 
##     Null deviance: 2.2082  on 1790  degrees of freedom
## Residual deviance: 2.0511  on 1787  degrees of freedom
## AIC: -7036.3
## 
## 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

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

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

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


# average carbono and compactation data 
mean_carbono_data <- carbono_data %>% 
  group_by (Punto) %>% 
  summarise (Carbono = mean(Carbono_percent),
             SoilCompactation= mean(Densidad_g)) 

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


glm3 <- glm(Richness ~  Carbono + SoilCompactation * Carbono* SoilCompactation, 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(glm3, glm4, glm5, glm6)

summary(glm3)

## 
## Call:
## glm(formula = Richness ~ Carbono + SoilCompactation * Carbono * 
##     SoilCompactation, family = poisson(), data = Full_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9629  -1.0228  -0.2231   0.5415   4.0674  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                2.9658     0.4082   7.266 3.70e-13 ***
## Carbono                  -12.5609     7.1739  -1.751    0.080 .  
## SoilCompactation          -1.4406     0.3135  -4.596 4.31e-06 ***
## Carbono:SoilCompactation   8.6428     5.7409   1.505    0.132    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 720.56  on 593  degrees of freedom
## Residual deviance: 665.80  on 590  degrees of freedom
## AIC: 2316
## 
## Number of Fisher Scoring iterations: 5

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

plot_model(glm3, type = "pred", terms = c("SoilCompactation"))

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)  

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

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

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


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

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

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 16299)
## 
## Matrix products: default
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2  fields_9.6      maps_3.3.0      spam_2.2-0     
##  [5] dotCall64_1.0-0 sqldf_0.4-11    RSQLite_2.1.1   gsubfn_0.7     
##  [9] proto_1.0.0     sjmisc_2.7.4    sjPlot_2.6.0    mapview_2.5.1  
## [13] sf_0.6-3        readxl_1.1.0    lubridate_1.7.4 ggmap_2.6.1    
## [17] rgdal_1.3-4     rgeos_0.3-28    raster_2.6-7    sp_1.3-1       
## [21] forcats_0.3.0   stringr_1.3.1   dplyr_0.7.6     purrr_0.2.5    
## [25] readr_1.1.1     tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0  
## [29] 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.16     
##  [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.5       lme4_1.1-18-1     
##  [25] bindr_0.1.1        survival_2.41-3    zoo_1.8-3         
##  [28] glue_1.3.0         gtable_0.2.0       emmeans_1.2.3     
##  [31] webshot_0.5.0      sjstats_0.17.0     scales_1.0.0      
##  [34] mvtnorm_1.0-8      DBI_1.0.0          ggeffects_0.5.0   
##  [37] Rcpp_0.12.18       viridisLite_0.3.0  xtable_1.8-3      
##  [40] spData_0.2.9.3     units_0.6-0        foreign_0.8-67    
##  [43] bit_1.1-14         mapproj_1.2.6      stats4_3.4.0      
##  [46] prediction_0.3.6   survey_3.33-2      htmlwidgets_1.2   
##  [49] httr_1.3.1         RColorBrewer_1.1-2 geosphere_1.5-7   
##  [52] modeltools_0.2-22  pkgconfig_2.0.2    nnet_7.3-12       
##  [55] labeling_0.3       tidyselect_0.2.4   rlang_0.2.2       
##  [58] reshape2_1.4.3     later_0.7.4        munsell_0.5.0     
##  [61] cellranger_1.1.0   tools_3.4.0        cli_1.0.0         
##  [64] sjlabelled_1.0.13  broom_0.5.0        ggridges_0.5.0    
##  [67] evaluate_0.11      yaml_2.2.0         knitr_1.20        
##  [70] bit64_0.9-7        satellite_1.0.1    RgoogleMaps_1.4.2 
##  [73] coin_1.2-2         nlme_3.1-137       mime_0.5          
##  [76] xml2_1.2.0         compiler_3.4.0     bayesplot_1.6.0   
##  [79] rstudioapi_0.7     png_0.1-7          e1071_1.7-0       
##  [82] stringi_1.2.4      lattice_0.20-35    Matrix_1.2-14     
##  [85] classInt_0.2-3     psych_1.8.4        nloptr_1.0.4      
##  [88] effects_4.0-2      stringdist_0.9.5.1 pillar_1.3.0      
##  [91] pwr_1.2-2          estimability_1.3   data.table_1.11.4 
##  [94] httpuv_1.4.5       R6_2.2.2           promises_1.0.1    
##  [97] codetools_0.2-15   MASS_7.3-47        assertthat_0.2.0  
## [100] chron_2.3-52       rjson_0.2.20       withr_2.1.2       
## [103] mnormt_1.5-5       multcomp_1.4-8     parallel_3.4.0    
## [106] hms_0.4.2          coda_0.19-1        glmmTMB_0.2.2.0   
## [109] class_7.3-14       minqa_1.2.4        rmarkdown_1.10.12 
## [112] snakecase_0.9.2    carData_3.0-1      shiny_1.1.0       
## [115] base64enc_0.1-3

Cited Literature