Exploring Biodiversity, Carbon and Soil compactation. A spatially explicit aproach
Diego J. Lizcano, Diego Nararrete, Alejandro Lopera, et al.
September 10, 2018
knitr::opts_chunk$set(echo = TRUE)
# put on output yaml
# link-citations: yes
# csl: C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/ecology.csl
# bibliography: C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/AAA.bib
# toc: yesTo 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. The initial idea changed to sample just ants and dung bettles in a 80 m grid. The Carbon data came in September 2018.
Data analysis sugestions y Diego N.
Question: To what extent soil compaction influences biodiversity and carbon content in soils of natural and human-modified ecosystems?
Results
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) # GLMMRead Data
For Dung beetle and ants
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(path2,"/TNC_Escarabajo_Carbono_June_2018_no_cultivo.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(path2,"/TNC_Escarabajo_Carbono_June_2018_no_cultivo.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"## [1] 722## [1] "min"## [1] 16For carbono
carbono_data <- read_excel(paste(path2,"/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
Is the carbon content diferent per deept?
Test for diferences in deept in an ANOVA
##
## 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## 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##
## 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## 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##
## 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## 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.1530252Si 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) 
##
## 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: 2YES!!!! ***
Make sf
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")
Plot Ant Data
Ant Richness
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_data2 %>%
group_by (Punto) %>%
summarise (Carbono = mean(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 ~ Carbono * 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)##
## Call:
## glm(formula = Richness ~ Carbono + SoilCompactation + factor(Cobertura),
## family = poisson(), data = Full_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0539 -0.8812 -0.1935 0.5388 3.7815
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.35653 0.23069 10.215 < 2e-16 ***
## Carbono -0.03056 0.01285 -2.378 0.0174 *
## SoilCompactation -0.73591 0.16395 -4.489 7.17e-06 ***
## factor(Cobertura)Pastizal -0.37026 0.05724 -6.468 9.91e-11 ***
## factor(Cobertura)Rastrojo -0.65462 0.14843 -4.410 1.03e-05 ***
## factor(Cobertura)Silvopastoril -0.52559 0.08895 -5.909 3.45e-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.19 on 583 degrees of freedom
## AIC: 2228.7
##
## Number of Fisher Scoring iterations: 5
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=4)
# 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){# length cobertura
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
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## 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.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.0 gsubfn_0.7
## [13] proto_1.0.0 sjmisc_2.7.5 sjPlot_2.6.1
## [16] mapview_2.3.0 sf_0.6-3 readxl_1.1.0
## [19] lubridate_1.7.4 ggmap_2.6.1 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.0 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.3.3 crosstalk_1.0.0
## [7] leaflet_2.0.2 TH.data_1.0-9 digest_0.6.18
## [10] foreach_1.4.4 htmltools_0.3.6 memoise_1.1.0
## [13] magrittr_1.5 modelr_0.1.2 R.utils_2.7.0
## [16] sandwich_2.5-0 jpeg_0.1-8 colorspace_1.3-2
## [19] blob_1.1.1 rvest_0.3.2 haven_1.1.2
## [22] tcltk_3.3.3 crayon_1.3.4 jsonlite_1.5
## [25] bindr_0.1.1 survival_2.40-1 zoo_1.8-4
## [28] iterators_1.0.10 glue_1.3.0 gtable_0.2.0
## [31] emmeans_1.2.4 webshot_0.5.1 sjstats_0.17.1
## [34] scales_1.0.0 mvtnorm_1.0-8 DBI_1.0.0
## [37] ggeffects_0.6.0 Rcpp_0.12.19 viridisLite_0.3.0
## [40] xtable_1.8-3 spData_0.2.9.4 units_0.6-1
## [43] foreign_0.8-67 bit_1.1-12 mapproj_1.2.6
## [46] stats4_3.3.3 prediction_0.3.6 htmlwidgets_1.3
## [49] httr_1.3.1 RColorBrewer_1.1-2 geosphere_1.5-7
## [52] modeltools_0.2-22 pkgconfig_2.0.2 R.methodsS3_1.7.1
## [55] labeling_0.3 tidyselect_0.2.5 rlang_0.2.2
## [58] reshape2_1.4.3 later_0.7.5 munsell_0.5.0
## [61] cellranger_1.1.0 tools_3.3.3 cli_1.0.1
## [64] sjlabelled_1.0.14 broom_0.5.0 ggridges_0.5.1
## [67] evaluate_0.12 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-131 mime_0.6
## [76] R.oo_1.22.0 xml2_1.2.0 bayesplot_1.6.0
## [79] rstudioapi_0.8 png_0.1-7 e1071_1.7-0
## [82] stringi_1.1.6 lattice_0.20-34 classInt_0.2-3
## [85] psych_1.8.4 nloptr_1.0.4 stringdist_0.9.5.1
## [88] pillar_1.3.0 pwr_1.2-2 estimability_1.3
## [91] data.table_1.11.8 httpuv_1.4.5 R6_2.3.0
## [94] promises_1.0.1 codetools_0.2-15 gdalUtils_2.0.1.14
## [97] MASS_7.3-45 assertthat_0.2.0 chron_2.3-52
## [100] rprojroot_1.3-2 rjson_0.2.20 withr_2.1.2
## [103] mnormt_1.5-5 multcomp_1.4-8 parallel_3.3.3
## [106] hms_0.4.2 coda_0.19-2 glmmTMB_0.2.2.0
## [109] class_7.3-14 minqa_1.2.4 rmarkdown_1.10
## [112] snakecase_0.9.2 shiny_1.1.0 base64enc_0.1-3