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 relationshipRead 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)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.
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:
##
## 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
##
## Now attempting to fit the 2 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > asymp : v
## > negexpo : v
##
## 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, 1000haValle_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 2 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > weibull3 : v
##
## 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, 1000haValle_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 4 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > monod : v
## > koba : v
## > power : v
##
## 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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 2 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > weibull3 : v
##
## 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, 1000haPiedemonte_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > mmf : v
## > asymp : v
##
## 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, 1000haPiedemonte_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > loga : v
## > linear : v
## > power : v
##
## 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))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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > power : v
## > linear : v
## > monod : v
##
## 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, 1000haEje_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > power : v
## > linear : v
## > koba : v
##
## 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, 1000haEje_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > power : v
## > linear : v
##
## 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))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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > mmf : v
## > weibull3 : v
##
## 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, 1000haBajo_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > monod : v
## > weibull3 : v
##
## 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, 1000haBajo_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > monod : v
## > koba : v
## > koba : v
##
## 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))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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > koba : v
## > monod : v
## > loga : v
##
## 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, 1000haBoyaca_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > negexpo : v
## > monod : v
## > koba : v
##
## 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, 1000haBoyaca_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, 1000hapuntos_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:
##
## 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
##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ------------------------------ multi-model SAR --
## > weibull3 : v
## > negexpo : v
## > loga : v
##
## 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))