Load packages

library(readxl) # lee excel
library(plyr)
library(tidyverse)# manipula datos eficientemente
library(openxlsx) # escribe excel
library(sp)
library(magrittr)
# library(Hmsc)
library(raster)
library(rasterVis)
library(sf) # new spatial package
library(mapview) # plotting maps in html
library(knitr)
library(colorspace)
library(DT)
library (sars) # Species area relationship

Read data

Read excel files from GAICA and paste all in aves_full

###################################
###
### AVES
###
###################################

 aves_atlantico <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "ATLANTICO"))
#
 aves_quindio <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "QUINDIO"))
#
 aves_santander <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "SANTANDER"))
#
 aves_meta <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "META"))


# aves_atlantico <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "ATLANTICO"))
# 
# aves_quindio <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "QUINDIO"))
# 
# aves_santander <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "SANTANDER"))
# 
# aves_meta <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Tablas_GCS/data/GAICA_Aves_2013.xlsx", sheet = "META"))


aves_full_2013 <- rbind(aves_atlantico, aves_quindio, aves_santander, aves_meta)

aves_full_2013$LOC <- as.numeric(aves_full_2013$LOC)
aves_full_2013$Fecha <- as.POSIXct(paste(
      aves_full_2013$YEAR, 
      aves_full_2013$MES,
      aves_full_2013$DIA,
      sep="-"),format = c("%Y-%m-%d"), tz="UTC")

rm(aves_atlantico)
rm(aves_quindio)
rm(aves_santander)
rm(aves_meta)

aves_full_2016 <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/AVES_GAICA_FINAL_2016_2017.xlsx", sheet = "INVENTARIO"))

# aves_full_2016 <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Tablas_GCS/data/AVES_GAICA_FINAL_2016_2017.xlsx", sheet = "INVENTARIO"))


aves_full_2016$YEAR <- lubridate::year(aves_full_2016$Fecha)
aves_full_2016$MES <- as.character(lubridate::month(aves_full_2016$Fecha))
aves_full_2016$DIA <- as.character(lubridate::day(aves_full_2016$Fecha))
aves_full <- bind_rows (aves_full_2016, aves_full_2013)

Species Area Models

The species–area relationship (SAR) describes the near universally observed pattern whereby the number of species increases with the area sampled, and it has been described as one of ecology’s few laws.

acording to:

Matthews, T.J., Triantis, K., Whittaker, R.J. and Guilhaumon, F. (2019) sars: an R package for fitting, evaluating and comparing species–area relationship models. Ecography, In press.

Species Area Models for Guajira

Bosque

Model: Cumulative Weibull 3 par

Formula: S == d(1 - exp(-c * A^z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_guajira <- puntos %>% filter(DEPARTAMENTO == "Guajira", Cobertura==c("Bosque", "Bosque ripario"))


count_guajira <- puntos_guajira %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_guajira <- 
count_guajira %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_guajira$area <- 7854 #en m2

sp_guajira <- sp_guajira %>%
mutate(a = cumsum(area ))

guajira <- cbind(sp_guajira[,5], sp_guajira[,3])
guajira[17,2] <- 235
guajira[18,2] <- 236
guajira[19,2] <- 237

#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= guajira, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
                                    normaTest = "none", homoTest = "none", 
                              neg_check = FALSE, confInt = TRUE, ciN = 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 10  models successfully fitted
## 
## The following model could not be fitted or was removed due to model checks:
## chapman 
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight    AICc    R2   R2a     Shape Asymptote
## 1  weibull3  0.607  65.179 0.996 0.996 convex up     FALSE
## 2   negexpo  0.180  67.610 0.995 0.995 convex up     FALSE
## 3       mmf  0.102  68.755 0.996 0.995 convex up     FALSE
## 4     asymp  0.055  69.969 0.995 0.995 convex up     FALSE
## 5     monod  0.044  70.445 0.994 0.994 convex up     FALSE
## 6      koba  0.013  72.874 0.994 0.993 convex up     FALSE
## 7     power  0.000  84.300 0.989 0.987 convex up     FALSE
## 8     heleg  0.000  87.557 0.989 0.986 convex up     FALSE
## 9    linear  0.000  98.874 0.975 0.972    linear     FALSE
## 10     loga  0.000 119.928 0.925 0.916 convex up     FALSE
fit2 <- sar_multi(data= guajira, obj =c("asymp", "negexpo"))
## 
##  Now attempting to fit the 2 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > asymp   : v
## > negexpo : v
fit3 <- sar_weibull3(guajira)# best AIC
print(summary(fit3))#, scipen = 999))
## 
## Model:
## Cumulative Weibull 3 par.
## 
## Call: 
## S == d(1 - exp(-c * A^z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##    0%   25%   50%   75%  100% 
## -8.10 -2.90  0.30  2.45  9.30 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%    97.5%
## d  3.1408e+02  2.3005e+01  1.3653e+01  3.1055e-10  2.6807e+02 360.0933
## c  1.8112e-06  1.0180e-06  1.7792e+00  9.4207e-02 -2.2477e-07   0.0000
## z  1.1457e+00  5.9085e-02  1.9391e+01  1.5406e-12  1.0275e+00   1.2639
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: 62.32, AICc: 65.18, BIC: 66.1
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Valle_del_Rio_Cesar_bosque_secundario <- fit3

# Save an object to a file
saveRDS(Valle_del_Rio_Cesar_bosque_secundario, file = "Valle_del_Rio_Cesar_bosque_secundario.rds")

Valle_del_Rio_Cesar_bosque_secundario <- readRDS(file = "Valle_del_Rio_Cesar_bosque_secundario.rds")

sar_pred(Valle_del_Rio_Cesar_bosque_secundario, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Silvopastoril

Model: Cumulative Weibull 3 par.

Formula: S == d(1 - exp(-c * A^z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")#, YEAR==2017)
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")#,  YEAR==2017)
puntos <- rbind(puntos_a, puntos_b)

puntos_guajira <- puntos %>% filter(DEPARTAMENTO == "Guajira", Cobertura==c("Bosque/ potrero con arboles dispersos", "Bosque/cultivo productivo", "Bosque/potrero con arboles dispersos"))

count_guajira <- puntos_guajira %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_guajira <- 
count_guajira %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_guajira$area <- 7854 #en m2

sp_guajira <- sp_guajira %>%
mutate(a = cumsum(area ))

guajira <- cbind(sp_guajira[,5], sp_guajira[,3])


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= guajira, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1  weibull3  0.469 11.336 0.999 0.999 convex up     FALSE
## 2   negexpo  0.334 12.015 0.999 0.998 convex up     FALSE
## 3       mmf  0.091 14.619 0.999 0.999 convex up     FALSE
## 4     monod  0.075 15.015 0.998 0.998 convex up     FALSE
## 5      koba  0.022 17.459 0.998 0.997 convex up     FALSE
## 6   chapman  0.009 19.215 0.999 0.998 convex up     FALSE
## 7     power  0.000 27.657 0.993 0.990 convex up     FALSE
## 8     asymp  0.000 31.786 0.995 0.992 convex up     FALSE
## 9     heleg  0.000 34.625 0.993 0.989 convex up     FALSE
## 10   linear  0.000 36.527 0.981 0.974    linear     FALSE
## 11     loga  0.000 41.789 0.965 0.954 convex up     FALSE
fit2 <- sar_multi(data= guajira, obj =c("negexpo", "weibull3"))
## 
##  Now attempting to fit the 2 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo  : v
## > weibull3 : v
fit3 <- sar_weibull3(guajira)# best AIC
print(summary(fit3))
## 
## Model:
## Cumulative Weibull 3 par.
## 
## Call: 
## S == d(1 - exp(-c * A^z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##   0%  25%  50%  75% 100% 
## -0.9 -0.7 -0.1  0.2  1.1 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%    97.5%
## d 1.4875e+02 8.0830e+00 1.8403e+01 1.6595e-06 1.3258e+02 164.9163
## c 5.6354e-06 1.7430e-06 3.2332e+00 1.7839e-02 2.1494e-06   0.0000
## z 1.1075e+00 3.6338e-02 3.0477e+01 8.2812e-08 1.0348e+00   1.1802
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: 1.34, AICc: 11.34, BIC: 2.13
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)

sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Valle_del_Rio_Cesar_silvopastoriles <- fit3

# Save an object to a file
saveRDS(Valle_del_Rio_Cesar_silvopastoriles, file = "Valle_del_Rio_Cesar_silvopastoriles.rds")

Valle_del_Rio_Cesar_silvopastoriles <- readRDS(file = "Valle_del_Rio_Cesar_silvopastoriles.rds")

sar_pred(Valle_del_Rio_Cesar_silvopastoriles, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Pasturas

Model: Monod

Formula: S = d/(1 + c * A^(-1))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")#, YEAR==2017)
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")#,  YEAR==2017)
puntos <- rbind(puntos_a, puntos_b)

puntos_guajira <- puntos %>% filter(DEPARTAMENTO == "Guajira", Cobertura==c("Bosque/pastos limpios", "Bosque/potrero con arboles dispersos"))

count_guajira <- puntos_guajira %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_guajira <- 
count_guajira %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_guajira$area <- 7854 #en m2

sp_guajira <- sp_guajira %>%
mutate(a = cumsum(area ))

guajira <- cbind(sp_guajira[,5], sp_guajira[,3])


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= guajira, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1     monod  0.333 15.261 0.999 0.999 convex up     FALSE
## 2   negexpo  0.330 15.283 0.999 0.999 convex up     FALSE
## 3      koba  0.213 16.157 0.999 0.999 convex up     FALSE
## 4     power  0.119 17.323 0.999 0.998 convex up     FALSE
## 5    linear  0.005 23.689 0.997 0.995    linear     FALSE
## 6      loga  0.000 40.591 0.954 0.923 convex up     FALSE
## 7  weibull3  0.000 44.541 0.999 0.998 convex up     FALSE
## 8       mmf  0.000 44.690 0.999 0.998 convex up     FALSE
## 9   chapman  0.000 45.283 0.999 0.998 convex up     FALSE
## 10    heleg  0.000 47.323 0.999 0.998 convex up     FALSE
## 11    asymp  0.000 61.551 0.990 0.974 convex up     FALSE
fit2 <- sar_multi(data= guajira, obj =c("negexpo", "monod", "koba", "power"))
## 
##  Now attempting to fit the 4 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > monod   : v
## > koba    : v
## > power   : v
fit3 <- sar_monod(guajira)# best AIC
print(summary(fit3))
## 
## Model:
## Monod
## 
## Call: 
## S == d/(1 + c * A^(-1))
## 
## Did the model converge: TRUE
## 
## Residuals:
##     0%    25%    50%    75%   100% 
## -1.200 -0.675  0.100  0.575  0.900 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%    97.5%
## d 5.3265e+02 5.6122e+01 9.4909e+00 6.8777e-04 4.2041e+02    644.9
## c 1.7603e+05 2.2379e+04 7.8659e+00 1.4118e-03 1.3127e+05 220786.3
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: 3.26, AICc: 15.26, BIC: 2.64
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Species Area Model for Meta

Bosque

Model: Negative exponential

Formula: S = d * (1 - exp(-z * A))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_meta <- puntos %>% filter(DEPARTAMENTO == "Meta", Cobertura==c("Bosque", "bosque", "Bosque/ rastrojo", "Bosque ripario", "Rastrojo/Arboles dispersos" ))

count_meta <- puntos_meta %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_meta <- 
count_meta %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_meta$area <- 7854 #en m2

sp_meta <- sp_meta %>%
mutate(a = cumsum(area ))

meta <- cbind(sp_meta[,5], sp_meta[,3])


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= meta, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE)#a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 9  models successfully fitted
## 
## The following models could not be fitted or were  removed due to model checks:
## chapman, asymp 
## 
## Ranked models based on AICc  weights:
## 
##      Model Weight    AICc    R2   R2a     Shape Asymptote
## 1  negexpo  0.555   0.525 1.000 1.000 convex up     FALSE
## 2 weibull3  0.445   0.969 1.000 1.000 convex up     FALSE
## 3      mmf  0.000  14.877 1.000 1.000 convex up     FALSE
## 4    monod  0.000  27.572 1.000 1.000 convex up     FALSE
## 5     koba  0.000  70.834 0.998 0.998 convex up     FALSE
## 6    power  0.000 102.687 0.994 0.993 convex up     FALSE
## 7   linear  0.000 143.803 0.972 0.970    linear     FALSE
## 8     loga  0.000 167.113 0.936 0.931 convex up     FALSE
## 9    heleg  0.000 242.733 0.135 0.027 convex up     FALSE
fit2 <- sar_multi(data= meta, obj =c("negexpo", "weibull3"))
## 
##  Now attempting to fit the 2 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo  : v
## > weibull3 : v
fit3 <- sar_negexpo(meta)# best AIC
print(summary(fit3))
## 
## Model:
## Negative exponential
## 
## Call: 
## S == d * (1 - exp(-z * A))
## 
## Did the model converge: TRUE
## 
## Residuals:
##    0%   25%   50%   75%  100% 
## -1.30 -0.50 -0.15  0.40  2.70 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%  97.5%
## d 3.4050e+02 2.1405e+00 1.5907e+02 2.1765e-40 3.3622e+02 344.78
## z 5.9101e-06 5.9875e-08 9.8708e+01 5.2158e-35 5.7904e-06   0.00
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: -0.47, AICc: 0.53, BIC: 3.52
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Piedemonte_del_Meta_bosque_secundario <- fit3

# Save an object to a file
saveRDS(Piedemonte_del_Meta_bosque_secundario, file = "Piedemonte_del_Meta_bosque_secundario.rds")

Piedemonte_del_Meta_bosque_secundario <- readRDS(file = "Piedemonte_del_Meta_bosque_secundario.rds")

sar_pred(Piedemonte_del_Meta_bosque_secundario, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Silvopastoril

Model: Cumulative Weibull 3 par.

Formula: S = d(1 - exp(-c * A^z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_meta <- puntos %>% filter(DEPARTAMENTO == "Meta", Cobertura==c("Potrero con arboles dispersos", "Bosque/potrero con arboles dispersos",  "Cerca viva", "Bosque / Sistema silvopastoril" , "Bosque/ potrero con arboles dispersos", "cerca viva", "Area de cultivo/ Bosque", "Bosque/ Sistema silvopastroril (Yopo)"))

count_meta <- puntos_meta %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_meta <- 
count_meta %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_meta$area <- 7854 #en m2

sp_meta <- sp_meta %>%
mutate(a = cumsum(area ))

meta <- cbind(sp_meta[,5], sp_meta[,3])


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= meta, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 10  models successfully fitted
## 
## The following model could not be fitted or was removed due to model checks:
## chapman 
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1  weibull3  0.455  9.562 0.998 0.998 convex up     FALSE
## 2       mmf  0.343 10.125 0.998 0.998 convex up     FALSE
## 3     asymp  0.156 11.708 0.998 0.998 convex up     FALSE
## 4      koba  0.040 14.411 0.997 0.997 convex up     FALSE
## 5     power  0.004 18.978 0.996 0.996 convex up     FALSE
## 6     monod  0.001 21.684 0.996 0.995 convex up     FALSE
## 7     heleg  0.001 22.229 0.997 0.996 convex up     FALSE
## 8   negexpo  0.000 28.573 0.994 0.993 convex up     FALSE
## 9    linear  0.000 51.873 0.972 0.968    linear     FALSE
## 10     loga  0.000 58.233 0.959 0.953 convex up     FALSE
fit2 <- sar_multi(data= meta, obj =c("weibull3", "mmf", "asymp"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > mmf      : v
## > asymp    : v
fit3 <- sar_weibull3(meta)# best AIC
print(summary(fit3))
## 
## Model:
## Cumulative Weibull 3 par.
## 
## Call: 
## S == d(1 - exp(-c * A^z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##     0%    25%    50%    75%   100% 
## -2.100 -0.650  0.250  0.625  1.300 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%    97.5%
## d 1.7647e+02 2.7226e+01 6.4817e+00 2.0617e-05 1.2202e+02 230.9183
## c 8.4434e-05 1.9237e-05 4.3892e+00 7.3210e-04 4.5960e-05   0.0001
## z 7.7596e-01 3.7050e-02 2.0944e+01 2.1212e-11 7.0186e-01   0.8501
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: 5.93, AICc: 9.56, BIC: 9.02
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Piedemonte_del_Meta_silvopastoriles <- fit3

# Save an object to a file
saveRDS(Piedemonte_del_Meta_silvopastoriles, file = "Piedemonte_del_Meta_silvopastoriles.rds")

Piedemonte_del_Meta_silvopastoriles <- readRDS(file = "Piedemonte_del_Meta_silvopastoriles.rds")

sar_pred(Piedemonte_del_Meta_silvopastoriles, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Pasturas

Model: Logarithmic

Formula: S = c + z * log(A)

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_meta <- puntos %>% filter(DEPARTAMENTO == "Meta", Cobertura==c("Rastrojo/Arboles dispersos", "Bosque secundario/pastos limpios", NA))

count_meta <- puntos_meta %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_meta <- 
count_meta %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(nn) %>%
  #tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(nn))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_meta$area <- as.numeric(7854) #en m2

sp_meta <- sp_meta %>%
mutate(a = cumsum(area))

meta <- cbind(sp_meta[,6], sp_meta[,4])
meta[3,1] <- 7854 +7854+7854
meta[4,1] <- 7854 +7854+7854+7854
meta[5,1] <- 7854 +7854+7854+7854+7854

#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= meta, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a       Shape Asymptote
## 1      loga  0.327 40.290 0.835 0.670   convex up     FALSE
## 2    linear  0.158 41.738 0.780 0.559      linear     FALSE
## 3     power  0.130 42.132 0.761 0.523 convex down     FALSE
## 4   negexpo  0.128 42.157 0.760 0.520   convex up     FALSE
## 5     monod  0.128 42.160 0.760 0.520   convex up     FALSE
## 6      koba  0.128 42.162 0.760 0.520   convex up     FALSE
## 7       mmf  0.000    Inf 0.956 0.824     sigmoid      TRUE
## 8     heleg  0.000    Inf 0.761 0.046 convex down     FALSE
## 9   chapman  0.000    Inf 0.760 0.041   convex up     FALSE
## 10 weibull3  0.000    Inf 0.963 0.852     sigmoid      TRUE
## 11    asymp  0.000    Inf 0.854 0.416   convex up     FALSE
fit2 <- sar_multi(data= meta, obj =c( "loga", "linear", "power"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > loga   : v
## > linear : v
## > power  : v
fit3 <- sar_loga(meta)# best AIC
print(summary(fit3))
## 
## Model:
## Logarithmic
## 
## Call: 
## S == c + z * log(A)
## 
## Did the model converge: TRUE
## 
## Residuals:
##   0%  25%  50%  75% 100% 
## -3.7 -1.8 -1.0  2.8  3.7 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%   97.5%
## c  -99.305060   28.259943   -3.513987    0.039086 -155.824946 -42.785
## z   11.072134    2.842327    3.895447    0.030012    5.387480  16.757
## 
## R-squared: 0.83, Adjusted R-squared: 0.67
## AIC: 16.29, AICc: 40.29, BIC: 15.12
## Observed shape: convex up, Asymptote: FALSE
## 
## 
## Warning: The fitted values of the model contain negative values (i.e. negative species richness values)
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = F)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 50000, 100000))

Species Area Model for Quindio

Bosque

Model: MMF

formula: S = d/(1 + c * A^(-z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_quindio <- puntos %>% filter(DEPARTAMENTO == "Quindio", Cobertura==c("Bosque ripario", "Bosque"))

count_quindio <- puntos_quindio %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_quindio <- 
count_quindio %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>%
  tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_quindio$area <- 7854 #en m2

sp_quindio <- sp_quindio %>%
mutate(a = cumsum(area ))

quindio <- cbind(sp_quindio[,5], sp_quindio[,3])

#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= quindio, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE)#a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 9  models successfully fitted
## 
## The following models could not be fitted or were  removed due to model checks:
## chapman, asymp 
## 
## Ranked models based on AICc  weights:
## 
##      Model Weight    AICc     R2    R2a     Shape Asymptote
## 1      mmf  0.976 -31.415  1.000  1.000 convex up     FALSE
## 2 weibull3  0.020 -23.611  1.000  1.000 convex up     FALSE
## 3  negexpo  0.003 -19.713  1.000  1.000 convex up     FALSE
## 4    monod  0.001 -18.196  1.000  1.000 convex up     FALSE
## 5    power  0.000  86.333  0.997  0.996 convex up     FALSE
## 6    heleg  0.000  89.238  0.997  0.996 convex up     FALSE
## 7   linear  0.000 125.729  0.982  0.980    linear     FALSE
## 8     loga  0.000 158.598  0.930  0.923 convex up     FALSE
## 9     koba  0.000 264.295 -4.746 -5.293    linear     FALSE
fit2 <- sar_multi(data= quindio, obj =c("power", "linear", "monod"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > power  : v
## > linear : v
## > monod  : v
fit3 <- sar_mmf(quindio)# best AIC
print(summary(fit3))
## 
## Model:
## MMF
## 
## Call: 
## S == d/(1 + c * A^(-z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##     0%    25%    50%    75%   100% 
## -0.600 -0.400  0.000  0.325  0.700 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%      97.5%
## d 7.4575e+02 9.9353e+00 7.5061e+01 5.2365e-27 7.2588e+02 7.6562e+02
## c 3.3691e+05 1.4126e+04 2.3849e+01 1.0809e-16 3.0865e+05 3.6516e+05
## z 1.0232e+00 5.1951e-03 1.9696e+02 8.6034e-36 1.0128e+00 1.0336e+00
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: -33.52, AICc: -31.41, BIC: -28.81
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Eje_Cafetero_bosque_secundario <- fit3

# Save an object to a file
saveRDS(Eje_Cafetero_bosque_secundario, file = "Eje_Cafetero_bosque_secundario.rds")

readRDS(file = "Eje_Cafetero_bosque_secundario.rds")
## 
## Model:
## MMF
## 
## Call:
## S == d/(1 + c * A^(-z))
## 
## Coefficients:
##            d            c            z 
## 7.457507e+02 3.369069e+05 1.023233e+00
sar_pred(Eje_Cafetero_bosque_secundario, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Silvopastoril

Model: Logarithmic

formula: S == c + z * log(A)

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_quindio <- puntos %>% filter(DEPARTAMENTO == "Quindio", Cobertura==c("Cerca viva", "Potrero con arboles dispersos", "Area de cultivo", "Sistema Silvopastoril", "cerca viva" , "Potrero con arboles dispersos/ Cerca viva natural"))

count_quindio <- puntos_quindio %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_quindio <- 
count_quindio %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(nn) %>%
  # tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(nn))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_quindio$area <- as.numeric(7854) #en m2

sp_quindio <- sp_quindio %>%
mutate(a = cumsum(area))

quindio <- cbind(sp_quindio[,6], sp_quindio[,4])
#quindio[3,1] <- 7854 +7854+7854
#quindio[4,1] <- 7854 +7854+7854+7854
#quindio[5,1] <- 7854 +7854+7854+7854+7854
quindio[6,1] <- 7854 +7854
quindio[7,1] <- 7854 +7854

#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= quindio, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2    R2a       Shape Asymptote
## 1      loga  0.175 24.510 0.389  0.084   convex up     FALSE
## 2    linear  0.175 24.510 0.389  0.084      linear     FALSE
## 3     power  0.175 24.510 0.389  0.084 convex down     FALSE
## 4     monod  0.161 24.673 0.375  0.062      linear     FALSE
## 5   negexpo  0.161 24.673 0.375  0.062      linear     FALSE
## 6      koba  0.151 24.805 0.363  0.044   convex up     FALSE
## 7       mmf  0.000 38.510 0.389 -0.222 convex down     FALSE
## 8     asymp  0.000 38.510 0.389 -0.222   convex up     FALSE
## 9  weibull3  0.000 38.510 0.389 -0.222   convex up     FALSE
## 10    heleg  0.000 38.510 0.389 -0.222 convex down     FALSE
## 11  chapman  0.000 38.674 0.375 -0.251      linear     FALSE
fit2 <- sar_multi(data= quindio, obj =c("power", "linear", "koba"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > power  : v
## > linear : v
## > koba   : v
fit3 <- sar_loga(quindio)# best AIC
print(summary(fit3))
## 
## Model:
## Logarithmic
## 
## Call: 
## S == c + z * log(A)
## 
## Did the model converge: TRUE
## 
## Residuals:
##   0%  25%  50%  75% 100% 
## -3.3 -1.5  0.3  1.0  3.7 
## 
## Parameters:
##    Estimate Std. Error   t value  Pr(>|t|)      2.5%   97.5%
## c -41.87564   25.88335  -1.61786   0.16662 -93.64234  9.8911
## z   4.92921    2.76203   1.78463   0.13439  -0.59485 10.4533
## 
## R-squared: 0.39, Adjusted R-squared: 0.08
## AIC: 16.51, AICc: 24.51, BIC: 16.35
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = TRUE)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Eje_Cafetero_silvopastoriles <- fit3

# Save an object to a file
saveRDS(Eje_Cafetero_silvopastoriles, file = "Eje_Cafetero_silvopastoriles.rds")

Eje_Cafetero_silvopastoriles <- readRDS(file = "Eje_Cafetero_silvopastoriles.rds")

sar_pred(Eje_Cafetero_silvopastoriles, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Pasturas

Model: Cumulative Weibull 3 par.

formula: S = d(1 - exp(-c * A^z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_quindio <- puntos %>% filter(DEPARTAMENTO == "Quindio", Cobertura==c("Cerca viva/Potrero", "Area de cultivo"))

count_quindio <- puntos_quindio %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_quindio <- 
count_quindio %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(nn) %>%
  #tally(nn, sort = TRUE) %>%
  mutate(s = cumsum(nn))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_quindio$area <- 7854 #en m2

sp_quindio <- sp_quindio %>%
mutate(a = cumsum(area ))

quindio <- cbind(sp_quindio[,6], sp_quindio[,4])

quindio[1,1] <- 7854
quindio[2,1] <- 7854+7854
quindio[3,1] <- 7854+7854+7854
quindio[4,1] <- 7854+7854+7854+7854
quindio[5,1] <- 7854+7854+7854+7854+7854
quindio[6,1] <- 7854+7854+7854+7854+7854+7854
quindio[7,1] <- 7854+7854+7854+7854+7854+7854+7854
quindio[6,2] <- 23
quindio[7,2] <- 25


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= quindio, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a       Shape Asymptote
## 1  weibull3  0.374 22.704 0.995 0.989     sigmoid      TRUE
## 2    linear  0.320 23.020 0.958 0.937      linear     FALSE
## 3     power  0.127 24.869 0.945 0.918 convex down     FALSE
## 4       mmf  0.120 24.983 0.992 0.985     sigmoid     FALSE
## 5   negexpo  0.017 28.852 0.903 0.854      linear     FALSE
## 6     monod  0.017 28.852 0.903 0.854      linear     FALSE
## 7      koba  0.017 28.884 0.902 0.854   convex up     FALSE
## 8      loga  0.007 30.639 0.875 0.812   convex up     FALSE
## 9     asymp  0.000 36.888 0.959 0.917   convex up     FALSE
## 10    heleg  0.000 38.869 0.945 0.890 convex down     FALSE
## 11  chapman  0.000 42.852 0.903 0.806      linear     FALSE
fit2 <- sar_multi(data= quindio, obj =c( "weibull3", "power", "linear"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > power    : v
## > linear   : v
fit3 <- sar_weibull3(quindio)# best AIC
print(summary(fit3))
## 
## Model:
## Cumulative Weibull 3 par.
## 
## Call: 
## S == d(1 - exp(-c * A^z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##    0%   25%   50%   75%  100% 
## -1.00 -0.70 -0.30  0.35  0.80 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%   97.5%
## d  2.4814e+01  9.0500e-01  2.7419e+01  1.0523e-05  2.3004e+01 26.6239
## c  1.4006e-15  5.0787e-15  2.7579e-01  7.9637e-01 -8.7567e-15  0.0000
## z  3.2847e+00  3.5338e-01  9.2950e+00  7.4534e-04  2.5779e+00  3.9914
## 
## R-squared: 0.99, Adjusted R-squared: 0.99
## AIC: 2.7, AICc: 22.7, BIC: 2.49
## Observed shape: sigmoid - observed shape algorithm failed: observed shape set to theoretical shape (sigmoid), Asymptote: TRUE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = F)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 50000, 100000))

Species Area Models for Atlantico

Bosque

Model: Cumulative Weibull 3 par

formula: S = d(1 - exp(-c * A^z))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_atlantico <- puntos %>% filter(DEPARTAMENTO == "Atlantico", Cobertura==c("Bosque"))

count_atlantico <- puntos_atlantico %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_atlantico <- 
count_atlantico %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(nn) %>% tally(nn, sort = TRUE) %>% mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_atlantico$area <- 7854 #en m2

sp_atlantico <- sp_atlantico %>%
mutate(a = cumsum(area ))

atlantico <- cbind(sp_atlantico[,5], sp_atlantico[,3])

# for(i in 1:52){
#   atlantico[i,1]<- 7854*i
# }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= atlantico, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 9  models successfully fitted
## 
## The following models could not be fitted or were  removed due to model checks:
## chapman, asymp 
## 
## Ranked models based on AICc  weights:
## 
##      Model Weight    AICc     R2    R2a     Shape Asymptote
## 1 weibull3      1  48.935  1.000  0.999 convex up     FALSE
## 2      mmf      0  66.587  0.999  0.999 convex up     FALSE
## 3  negexpo      0  78.526  0.998  0.998 convex up     FALSE
## 4    monod      0  91.129  0.997  0.996 convex up     FALSE
## 5    power      0 126.577  0.986  0.984 convex up     FALSE
## 6    heleg      0 129.473  0.986  0.984 convex up     FALSE
## 7   linear      0 152.137  0.958  0.954    linear     FALSE
## 8     loga      0 158.373  0.946  0.941 convex up     FALSE
## 9     koba      0 271.854 -5.096 -5.676    linear     FALSE
fit2 <- sar_multi(data= atlantico, obj =c( "negexpo", "mmf", "weibull3"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo  : v
## > mmf      : v
## > weibull3 : v
fit3 <- sar_weibull3(atlantico)# best AIC
 print(summary(fit3))
## 
## Model:
## Cumulative Weibull 3 par.
## 
## Call: 
## S == d(1 - exp(-c * A^z))
## 
## Did the model converge: TRUE
## 
## Residuals:
##     0%    25%    50%    75%   100% 
## -3.800 -1.675  0.300  1.525  4.400 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%    97.5%
## d 4.0960e+02 6.3369e+00 6.4637e+01 1.1938e-25 3.9693e+02 422.2742
## c 1.9286e-06 3.4533e-07 5.5849e+00 1.5278e-05 1.2380e-06   0.0000
## z 1.1398e+00 1.7772e-02 6.4133e+01 1.4062e-25 1.1042e+00   1.1753
## 
## R-squared: 1, Adjusted R-squared: 1
## AIC: 46.83, AICc: 48.94, BIC: 51.54
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Bajo_Magdalena_bosque_secundario <- fit3

# Save an object to a file
saveRDS(Bajo_Magdalena_bosque_secundario, file = "Bajo_Magdalena_bosque_secundario.rds")

Bajo_Magdalena_bosque_secundario <- readRDS(file = "Bajo_Magdalena_bosque_secundario.rds")

sar_pred(Bajo_Magdalena_bosque_secundario, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Silvopastoril

Model: Negative exponential

formula: S = d * (1 - exp(-z * A))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_atlantico <- puntos %>% filter(DEPARTAMENTO == "Atlantico", Cobertura==c("Bosque/potrero con arboles dispersos", "Bosque/ potrero con arboles dispersos", "Sistema Silvopastoril", "Bosque/cultivo productivo"))

count_atlantico <- puntos_atlantico %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_atlantico <- 
count_atlantico %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(n) %>% tally(n, sort = TRUE) %>% mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_atlantico$area <- 7854 #en m2

sp_atlantico <- sp_atlantico %>%
mutate(a = cumsum(area ))

atlantico <- cbind(sp_atlantico[,5], sp_atlantico[,3])

# for(i in 1:52){
#   atlantico[i,1]<- 7854*i
# }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= atlantico, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1   negexpo  0.500 -0.409 0.993 0.990 convex up     FALSE
## 2  weibull3  0.185  1.578 0.997 0.995 convex up     FALSE
## 3     monod  0.134  2.221 0.990 0.986 convex up     FALSE
## 4       mmf  0.049  4.244 0.996 0.993 convex up     FALSE
## 5      koba  0.046  4.349 0.987 0.981 convex up     FALSE
## 6     asymp  0.044  4.456 0.996 0.993 convex up     FALSE
## 7      loga  0.034  4.994 0.986 0.980 convex up     FALSE
## 8   chapman  0.005  8.924 0.993 0.987 convex up     FALSE
## 9     power  0.003 10.184 0.972 0.961 convex up     FALSE
## 10   linear  0.000 16.319 0.941 0.917    linear     FALSE
## 11    heleg  0.000 19.517 0.972 0.952 convex up     FALSE
fit2 <- sar_multi(data= atlantico, obj =c( "negexpo", "monod", "weibull3"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo  : v
## > monod    : v
## > weibull3 : v
fit3 <- sar_negexpo(atlantico)# best AIC
 print(summary(fit3))
## 
## Model:
## Negative exponential
## 
## Call: 
## S == d * (1 - exp(-z * A))
## 
## Did the model converge: TRUE
## 
## Residuals:
##     0%    25%    50%    75%   100% 
## -1.000 -0.125  0.250  0.400  0.500 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%  97.5%
## d 2.6117e+01 1.5932e+00 1.6392e+01 3.2830e-06 2.2930e+01 29.303
## z 2.4207e-05 2.6191e-06 9.2425e+00 9.0598e-05 1.8969e-05  0.000
## 
## R-squared: 0.99, Adjusted R-squared: 0.99
## AIC: -6.41, AICc: -0.41, BIC: -6.17
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Bajo_Magdalena_silvopastoriles <- fit3

# Save an object to a file
saveRDS(Bajo_Magdalena_silvopastoriles, file = "Bajo_Magdalena_silvopastoriles.rds")

Bajo_Magdalena_silvopastoriles <- readRDS(file = "Bajo_Magdalena_silvopastoriles.rds")

sar_pred(Bajo_Magdalena_silvopastoriles, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Pasturas

Model: Kobayashi

formula: S = c * log(1 + A/z)

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_atlantico <- puntos %>% filter(DEPARTAMENTO == "Atlantico", Cobertura==c("Potrero con arboles dispersos", "Bosque/ Potrero con arboles dispersos"))

count_atlantico <- puntos_atlantico %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_atlantico <- 
count_atlantico %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(n) %>% 
  #tally(n, sort = TRUE) %>% 
  mutate(s = cumsum(nn))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_atlantico$area <- 7854 #en m2

sp_atlantico <- sp_atlantico %>%
mutate(a = cumsum(area ))

atlantico <- cbind(sp_atlantico[,6], sp_atlantico[,4])
atlantico[5,2] <- 15
atlantico[6,2] <- 15
atlantico[7,2] <- 15

 for(i in 1:7){
   atlantico[i,1]<- 7854*i
 }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= atlantico, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1     monod  0.237 14.417 0.790 0.685 convex up     FALSE
## 2      loga  0.210 14.663 0.783 0.674 convex up     FALSE
## 3      koba  0.206 14.701 0.782 0.672 convex up     FALSE
## 4     power  0.164 15.161 0.767 0.650 convex up     FALSE
## 5   negexpo  0.128 15.652 0.750 0.625 convex up      TRUE
## 6    linear  0.055 17.343 0.681 0.522    linear     FALSE
## 7       mmf  0.000 28.393 0.791 0.582 convex up     FALSE
## 8  weibull3  0.000 28.408 0.791 0.581 convex up     FALSE
## 9     heleg  0.000 29.151 0.767 0.534 convex up     FALSE
## 10    asymp  0.000 29.606 0.751 0.503 convex up     FALSE
## 11  chapman  0.000 29.741 0.747 0.493 convex up      TRUE
fit2 <- sar_multi(data= atlantico, obj =c( "monod", "koba", "koba"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > monod : v
## > koba  : v
## > koba  : v
fit3 <- sar_koba(atlantico)# best AIC
 print(summary(fit3))
## 
## Model:
## Kobayashi
## 
## Call: 
## S == c * log(1 + A/z)
## 
## Did the model converge: TRUE
## 
## Residuals:
##    0%   25%   50%   75%  100% 
## -1.40 -0.85  0.10  0.55  1.90 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%     97.5%
## c    3.188204    0.836232    3.812582    0.012467    1.515739    4.8607
## z  419.807307  472.847468    0.887828    0.415294 -525.887630 1365.5022
## 
## R-squared: 0.78, Adjusted R-squared: 0.67
## AIC: 6.7, AICc: 14.7, BIC: 6.54
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 50000, 100000))

Species Area Models for Santander

Bosque

Model: Logarithmic

formula: S == c + z * log(A)

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_santander <- puntos %>% filter(DEPARTAMENTO == "Santander")

puntos_santander$Cobertura[1:150] <- "Bosque"
puntos_santander <- puntos_santander[1:150,]

count_santander <- puntos_santander %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_santander <- 
count_santander %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(n) %>% tally(nn) %>% mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_santander$area <- 7854 #en m2

sp_santander <- sp_santander %>%
mutate(a = cumsum(area ))

santander <- cbind(sp_santander[,5], sp_santander[,3])

# for(i in 1:52){
#   santander[i,1]<- 7854*i
# }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= santander, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1      koba  0.541 34.614 0.997 0.996 convex up     FALSE
## 2     monod  0.253 36.137 0.996 0.996 convex up     FALSE
## 3       mmf  0.075 38.576 0.997 0.996 convex up     FALSE
## 4   negexpo  0.065 38.849 0.996 0.996 convex up     FALSE
## 5  weibull3  0.051 39.328 0.996 0.996 convex up     FALSE
## 6   chapman  0.014 41.869 0.996 0.995 convex up     FALSE
## 7     power  0.001 47.111 0.994 0.994 convex up     FALSE
## 8     heleg  0.000 50.128 0.994 0.993 convex up     FALSE
## 9     asymp  0.000 54.121 0.993 0.992 convex up     FALSE
## 10   linear  0.000 75.962 0.978 0.976    linear     FALSE
## 11     loga  0.000 98.587 0.939 0.933 convex up     FALSE
fit2 <- sar_multi(data= santander, obj =c( "koba", "monod", "loga"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > koba  : v
## > monod : v
## > loga  : v
fit3 <- sar_loga(santander)# best AIC
 print(summary(fit3))
## 
## Model:
## Logarithmic
## 
## Call: 
## S == c + z * log(A)
## 
## Did the model converge: TRUE
## 
## Residuals:
##      0%     25%     50%     75%    100% 
## -21.400  -4.950   2.900   5.975  10.600 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%    97.5%
## c -3.6373e+02  2.4807e+01 -1.4662e+01  3.6596e-12 -4.1335e+02 -314.117
## z  3.8953e+01  2.2148e+00  1.7588e+01  1.2366e-13  3.4524e+01   43.383
## 
## R-squared: 0.94, Adjusted R-squared: 0.93
## AIC: 97.25, AICc: 98.59, BIC: 100.53
## Observed shape: convex up, Asymptote: FALSE
## 
## 
## Warning: The fitted values of the model contain negative values (i.e. negative species richness values)
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Boyaca_y_Santander_bosque_secundario <- fit3

# Save an object to a file
saveRDS(Boyaca_y_Santander_bosque_secundario, file = "Boyaca_y_Santander_bosque_secundario.rds")

Boyaca_y_Santander_bosque_secundario <- readRDS(file = "Boyaca_y_Santander_bosque_secundario.rds")

sar_pred(Boyaca_y_Santander_bosque_secundario, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Silvopastoril

Model: Negative exponential

formula: S = d * (1 - exp(-z * A))

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_santander <- puntos %>% filter(DEPARTAMENTO == "Santander")

puntos_santander$Cobertura[151:275] <- "SPP"
puntos_santander <- puntos_santander[151:275,]

count_santander <- puntos_santander %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_santander <- 
count_santander %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(n) %>% tally(nn) %>% mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_santander$area <- 7854 #en m2

sp_santander <- sp_santander %>%
mutate(a = cumsum(area ))

santander <- cbind(sp_santander[,5], sp_santander[,3])
santander[10,2] <- 80
santander[11,2] <- 80
# for(i in 1:52){
#   santander[i,1]<- 7854*i
# }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= santander, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 10  models successfully fitted
## 
## The following model could not be fitted or was removed due to model checks:
## chapman 
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1   negexpo  0.405 28.080 0.989 0.987 convex up     FALSE
## 2     monod  0.278 28.832 0.989 0.986 convex up     FALSE
## 3      koba  0.189 29.604 0.988 0.985 convex up     FALSE
## 4  weibull3  0.031 33.232 0.989 0.985 convex up     FALSE
## 5     asymp  0.030 33.311 0.989 0.985 convex up     FALSE
## 6       mmf  0.023 33.828 0.989 0.984 convex up     FALSE
## 7     heleg  0.023 33.828 0.989 0.984 convex up     FALSE
## 8     power  0.020 34.051 0.982 0.977 convex up     FALSE
## 9    linear  0.000 42.105 0.962 0.952    linear     FALSE
## 10     loga  0.000 42.885 0.959 0.949 convex up     FALSE
fit2 <- sar_multi(data= santander, obj =c( "negexpo", "monod", "koba"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > monod   : v
## > koba    : v
fit3 <- sar_negexpo(santander)# best AIC
 print(summary(fit3))
## 
## Model:
## Negative exponential
## 
## Call: 
## S == d * (1 - exp(-z * A))
## 
## Did the model converge: TRUE
## 
## Residuals:
##   0%  25%  50%  75% 100% 
## -4.6 -1.3 -0.2  1.8  3.8 
## 
## Parameters:
##     Estimate Std. Error    t value   Pr(>|t|)       2.5%  97.5%
## d 1.1579e+02 9.1695e+00 1.2627e+01 4.9830e-07 9.7447e+01 134.12
## z 1.4889e-05 1.9135e-06 7.7810e+00 2.7612e-05 1.1062e-05   0.00
## 
## R-squared: 0.99, Adjusted R-squared: 0.99
## AIC: 24.65, AICc: 28.08, BIC: 25.84
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha
Boyaca_y_Santander_silvopastoriles <- fit3

# Save an object to a file
saveRDS(Boyaca_y_Santander_silvopastoriles, file = "Boyaca_y_Santander_silvopastoriles.rds")

Boyaca_y_Santander_silvopastoriles <- readRDS(file = "Boyaca_y_Santander_silvopastoriles.rds")

sar_pred(Boyaca_y_Santander_silvopastoriles, area = c(10000, 100000, 1000000, 10000000)) # predict 1ha, 10ha, 100ha, 1000ha

Pasturas

Model: Logarithmic

formula: S = c + z * log(A)

puntos_a <- aves_full %>% filter(METODO == "Puntos de Muestreo")
puntos_b <- aves_full %>% filter(METODO == "PUNTOS")
puntos <- rbind(puntos_a, puntos_b)

puntos_santander <- puntos %>% filter(DEPARTAMENTO == "Santander")

puntos_santander$Cobertura[276:334] <- "SPP"
puntos_santander <- puntos_santander[276:334,]

count_santander <- puntos_santander %>% dplyr::select(NOMBRE_CIENTIFICO, PUNTO_MUESTREO, Fecha) %>% 
                      # filter(YEAR == unique(site_species_a$YEAR)[1]) %>%  # yr 2013
                      dplyr::count(NOMBRE_CIENTIFICO, PUNTO_MUESTREO)#, Fecha) 



sp_santander <- 
count_santander %>% group_by(PUNTO_MUESTREO ) %>% count(n) %>% arrange(n) %>% tally(nn) %>% mutate(s = cumsum(n))

# 3.1416 * (50 *50) es area del sitio de conteo

sp_santander$area <- 7854 #en m2

sp_santander <- sp_santander %>%
mutate(a = cumsum(area ))

santander <- cbind(sp_santander[,5], sp_santander[,3])
santander[1,2] <- 5
santander[5,2] <- 35
santander[6,2] <- 38
santander[7,2] <- 37
# for(i in 1:52){
#   santander[i,1]<- 7854*i
# }


#run the ‘sar_average’ function using a vector of model names 
fit <- sar_average(data= santander, obj =c("power","loga","koba","mmf","monod", "heleg",
                                   "negexpo","chapman","weibull3","asymp", "linear"),
normaTest = "none", homoTest = "none", neg_check = FALSE, confInt = TRUE, ciN
= 50, verb = FALSE) #a message is provided indicating that one model (koba) could not be
## 
## Calculating sar_multi confidence intervals - this may take  some time:
##fitted

summary(fit) # see models and select lower AIC
## 
## Sar_average object summary:
## 
## 11  models successfully fitted
## 
## All models were fitted successfully
## 
## Ranked models based on AICc  weights:
## 
##       Model Weight   AICc    R2   R2a     Shape Asymptote
## 1      loga  0.908 17.857 0.986 0.980 convex up     FALSE
## 2  weibull3  0.022 25.336 0.995 0.989 convex up     FALSE
## 3   negexpo  0.020 25.535 0.959 0.939 convex up     FALSE
## 4       mmf  0.014 26.201 0.994 0.988 convex up     FALSE
## 5     asymp  0.013 26.338 0.994 0.988 convex up     FALSE
## 6     monod  0.012 26.427 0.954 0.931 convex up     FALSE
## 7      koba  0.009 27.176 0.949 0.923 convex up     FALSE
## 8     power  0.002 29.795 0.925 0.888 convex up     FALSE
## 9    linear  0.001 32.652 0.887 0.831    linear     FALSE
## 10  chapman  0.000 39.535 0.959 0.919 convex up     FALSE
## 11    heleg  0.000 43.795 0.925 0.850 convex up     FALSE
fit2 <- sar_multi(data= santander, obj =c( "weibull3", "negexpo", "loga"))
## 
##  Now attempting to fit the 3 SAR models: 
## 
## --  multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > negexpo  : v
## > loga     : v
fit3 <- sar_loga(santander)# best AIC
 print(summary(fit3))
## 
## Model:
## Logarithmic
## 
## Call: 
## S == c + z * log(A)
## 
## Did the model converge: TRUE
## 
## Residuals:
##   0%  25%  50%  75% 100% 
## -1.4 -1.2  0.2  0.5  2.6 
## 
## Parameters:
##      Estimate  Std. Error     t value    Pr(>|t|)        2.5%    97.5%
## c -1.5356e+02  9.4825e+00 -1.6194e+01  1.6365e-05 -1.7253e+02 -134.597
## z  1.7697e+01  9.2907e-01  1.9048e+01  7.3498e-06  1.5839e+01   19.555
## 
## R-squared: 0.99, Adjusted R-squared: 0.98
## AIC: 9.86, AICc: 17.86, BIC: 9.69
## Observed shape: convex up, Asymptote: FALSE
# par(mfrow = c(2,1)) #plot all model fits with the multimodel SAR curve
plot(fit, ModTitle = "a) Multimodel SAR")

#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own 
plot(fit, allCurves = FALSE, ModTitle =
      "c) Multimodel SAR with confidence intervals", confInt = T)

#Barplot of the information criterion weights of each model 
# plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)


sar_pred(fit3, area = c(10000, 50000, 100000))