Spatial Occupancy Model
Alan Ivory
2026-02-19
Set-up
Install necessary packages and import appropriate data
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse, readxl, raster, vegan, tigris, sf, sjPlot, sp, spOccupancy, ggrepel, lme4, lmerTest, MuMIn, brms, MCMCvis, bayesplot, patchwork)
set.seed(97)
# Tree PCQ Data
tree_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Tree_PCQ")
# CWD Data
cwd_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "CWD")
# Soil Data
fuel_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Fuel_Sampling")
# Veg Data
Veg_Cover <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Veg_Cover")
# Shrub Cover Data
shrub_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Shrub_Cover")
# Site Data
CameraData <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraData.xlsx")
CameraLoc <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
sheet = "CameraLocations")
# Add effort data
effort_matrix <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
sheet = "Effort_Matrix_Full") %>%
pivot_longer(cols = matches("^202[4-5]-"), names_to = "week", values_to = "days") %>%
filter(days == "7") %>%
dplyr::select(Plot, week)Coordinate Transformation
coords_sf <- CameraLoc %>%
dplyr::select(Plot, Long, Lat) %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
st_transform(crs = 32616) # UTM zone 16N
coords <- st_coordinates(coords_sf)
rownames(coords) <- CameraLoc$PlotNumber of quadrats sampled per plot
# Count the total number of quadrats per plot
quadrat_count <- Veg_Cover %>%
group_by(Plot) %>%
summarize(total_quadrats = n_distinct(Quadrat), .groups = "drop")Filter All data to only include specified species (Per PLANTS database)
#Filter tree data to only include trees with "tree" in the growth column
tree_data <- dplyr::filter(tree_data, Growth == "Tree")
#Filter Veg Cover to exclude Shrubs and Trees
Veg_Cover <- dplyr::filter(Veg_Cover, Growth != "Shrub" & Growth != "Tree")
#Filter Shrub Cover to only include Shrubs and Trees
shrub_data <- dplyr::filter(shrub_data, Growth == "Shrub" | Growth == "Tree")Provides the option to rarify vegetative species
# Calculate the total number of sites
total_sites <- nrow(CameraLoc)
# Function to filter data by frequency
filter_by_frequency <- function(df) {
# Group data by species and calculate the frequency
freq <- df %>%
group_by(Species) %>%
summarise(Frequency = n_distinct(Plot) / nrow(CameraLoc) * 100) %>%
filter(Frequency >= 0)
# Filter the original data to include only species with frequency >= XX%
filtered_df <- df %>%
filter(Species %in% freq$Species)
return(filtered_df)
}
# Filter tree data by frequency
tree_data <- filter_by_frequency(tree_data)
# Filter Veg Cover data by frequency
Veg_Cover <- filter_by_frequency(Veg_Cover)
# Filter Shrub Cover data by frequency
shrub_data <- filter_by_frequency(shrub_data)Shrub Cover Conversion
# Total length of Shrub cover at a site
shrub_cover <- shrub_data %>%
mutate(Cover = Line_End - Line_Start) %>%
group_by(Species_Name, Plot) %>%
summarise(Shrub_Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
mutate(Shrub_Percent_Cover = Shrub_Total_Cover / 3000 * 100)
# Summed length of shrub over at a site
shrub_cover_summed <- shrub_cover %>%
group_by(Plot) %>%
summarize(total_shrub_cover = sum(Shrub_Total_Cover, na.rm = TRUE), .groups = "drop")Herbacous Cover Conversion
# Combine Plot and Quadrat columns
Veg_Cover <- Veg_Cover %>%
mutate(Plot_Quadrat = paste(Plot, Quadrat, sep = '_'))
# Join with CogonSites to get site information
Veg_Cover <- Veg_Cover %>%
left_join(CameraLoc, by = "Plot")
# Sum species cover across quadrats for each species at each plot
veg_cover_summed <- Veg_Cover %>%
group_by(Plot, Species_Name) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE), .groups = "drop")
# Calculate average herbaceous species cover
avg_species_cover <- veg_cover_summed %>%
left_join(quadrat_count, by = "Plot") %>%
mutate(avg_cover = total_cover / total_quadrats)Merging Herb cover with Shrub
This species matrix includes herbaceous and shrub species
# Merge shrub cover with herbaceous average cover
combined_cover <- avg_species_cover %>%
full_join(
shrub_cover %>%
dplyr::select(Plot, Species_Name, Shrub_Percent_Cover),
by = c("Plot", "Species_Name")
) %>%
mutate(
overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Shrub_Percent_Cover), TRUE, FALSE), # Flag overlaps
final_cover = case_when(
!is.na(avg_cover) & is.na(Shrub_Percent_Cover) ~ avg_cover, # Use herbaceous cover if no shrub data
is.na(avg_cover) & !is.na(Shrub_Percent_Cover) ~ Shrub_Percent_Cover, # Use shrub cover if no herbaceous data
TRUE ~ NA_real_ # Leave as NA where overlaps exist
)
)
# Species Matrix
species_matrix <- combined_cover %>%
dplyr::select(Plot, Species_Name, final_cover) %>%
pivot_wider(
names_from = Species_Name,
values_from = final_cover,
values_fill = 0
)Summarize Cogongrass Cover
avg_cogongrass_cover <- species_matrix %>%
group_by(Plot) %>%
summarize(Avg_Cogongrass_Cover = sum(Imperata_cylindrica, na.rm = TRUE) / n(), .groups = "drop")Herbacous Shannon Diversity Index
# Summarize species cover by site
site_species_cover <- Veg_Cover %>%
group_by(Plot, Species_Name) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'Plot'. You can override using the
## `.groups` argument.## Remove all Imperata_cylindrica_Live and Imperata_cylindrica from species
site_species_cover <- site_species_cover %>%
filter(Species_Name != "Imperata_cylindrica_Live" & Species_Name != "Imperata_cylindrica")
# Calculate Shannon diversity per site
Veg_shannon_diversity <- site_species_cover %>%
group_by(Plot) %>%
mutate(proportion = total_cover / sum(total_cover)) %>%
summarize(Veg_shannon_index = -sum(proportion * log(proportion), na.rm = TRUE))
print(Veg_shannon_diversity)## # A tibble: 206 × 2
## Plot Veg_shannon_index
## <chr> <dbl>
## 1 BI200 2.75
## 2 BI201 2.70
## 3 BI202 2.59
## 4 BI97 1.61
## 5 BI99 2.97
## 6 BN210 2.97
## 7 BN211 2.43
## 8 BN212 2.22
## 9 BN96 3.05
## 10 BN98 2.79
## # ℹ 196 more rowsVegetation Height
if (!is.numeric(fuel_data$Height)) {
fuel_data$Height <- as.numeric(as.character(fuel_data$Height))
}## Warning: NAs introduced by coercion# Calculate average vegetation height per plot
veg_height <- fuel_data %>%
group_by(Plot) %>%
summarize(avg_veg_height = mean(Height, na.rm = TRUE), .groups = "drop")Tree Metrics
# Tree density from point-centered quarter data
if (!is.numeric(tree_data$Distance)) {
tree_data$Distance <- as.numeric(as.character(tree_data$Distance))
}
tree_density_data <- tree_data %>%
group_by(Plot) %>%
summarize(Average_Distance = mean(Distance) / 100, # Convert to meters
Tree_Density = 10000 / (Average_Distance^2)) # Convert to trees per hectare
# Average canopy cover from vegetation quadrats
tree_canopy_data <- Veg_Cover %>%
distinct(Plot, Quadrat, .keep_all = TRUE) %>% # Ensure each quadrat counts once per plot
group_by(Plot) %>%
summarize(Avg_Canopy_Cover = mean(Canopy_Cover, na.rm = TRUE), .groups = "drop") # Calculate the average canopy cover per plot
cor(tree_density_data$Tree_Density, tree_canopy_data$Avg_Canopy_Cover)## [1] 0.2742307CWD Clean
cwd_pa <- cwd_data %>%
group_by(Plot) %>%
summarise(
CWD_presence = as.integer(any(Line_Start > 0 & Line_End > 0)),
.groups = "drop"
)Objective 1: Assessing if cogongrass invasion-related factors influence wildlife presence
Add Covariates
CameraLoc <- CameraLoc %>%
left_join(cwd_pa, by = "Plot") %>%
left_join(Veg_shannon_diversity, by = "Plot") %>%
left_join(avg_cogongrass_cover, by = "Plot") %>%
left_join(shrub_cover_summed %>% dplyr::select(Plot, total_shrub_cover), by = "Plot") %>%
left_join(veg_height, by = "Plot") %>%
left_join(tree_density_data %>% dplyr::select(Plot, Tree_Density), by = "Plot") %>%
left_join(tree_canopy_data %>% dplyr::select(Plot, Avg_Canopy_Cover), by = "Plot") %>%
dplyr::select(-Authority)Distribution of Covariates
par(mfrow = c(3, 3))
hist(CameraLoc$Cogon_Patch_Size, main = "Cogon Patch Size", xlab = "Size (m²)")
hist(CameraLoc$Veg_shannon_index, main = "Vegetation Shannon Index", xlab = "Shannon Index")
hist(CameraLoc$total_shrub_cover, main = "Total Shrub Cover", xlab = "Shrub Cover (%)")
hist(CameraLoc$Avg_Cogongrass_Cover, main = "Average Cogongrass Cover", xlab = "Cogongrass Cover (%)")
hist(CameraLoc$Tree_Density, main = "Tree Density", xlab = "Trees per ha")
hist(CameraLoc$Avg_Canopy_Cover, main = "Average Canopy Cover", xlab = "Canopy Cover (%)")
hist(CameraLoc$avg_veg_height, main = "Average Vegetation Height", xlab = "Height (cm)")
table(CameraLoc$CWD_presence)##
## 0 1
## 25 7par(mfrow = c(1, 1))
Subset to mammals with >10 observations
# Group by Name and count the number of observations
species_counts <- CameraData %>%
filter((Class == "Mammalia" & Name != "Odocoileus_virginianus") | Name == "Meleagris_gallopavo") %>%
group_by(Name) %>%
summarize(count = n(), .groups = "drop")
# Filter for species with count greater than 10
species_subset <- species_counts %>%
filter(count > 10) %>%
pull(Name)
# Filter CameraData to only include species with count > 10
CameraData <- CameraData %>%
filter(Name %in% species_subset)Create Detection Matrix
# Format Data Weekly
observations_weekly <- CameraData %>%
group_by(Plot, week = format(as.Date(Date), "%Y-%U"), Name) %>%
summarise(observations = n(), .groups = 'drop')
# Merge with Effort Matrix to include only valid weeks
observations_weekly <- effort_matrix %>%
left_join(observations_weekly, by = c("Plot" = "Plot", "week")) %>%
replace_na(list(observations = 0))
# Convert to wide format
observations_species <- observations_weekly %>%
pivot_wider(names_from = Name, values_from = observations, values_fill = list(observations = 0)) %>%
dplyr::select(-"NA")
# Create detection array
site_names <- sort(unique(observations_species$Plot))
species_names <- setdiff(colnames(observations_species), c("Plot", "week"))
num_sites <- length(site_names)
num_weeks <- length(unique(observations_species$week))
num_species <- length(species_names)
detection_array <- array(0, dim = c(num_sites, num_weeks, num_species))
dimnames(detection_array) <- list(site_names, unique(observations_species$week), species_names)
for (s in seq_along(species_names)) {
species_col <- species_names[s]
for (i in seq_along(site_names)) {
site <- site_names[i]
for (j in seq_along(unique(observations_species$week))) {
week <- unique(observations_species$week)[j]
detection_array[i, j, s] <- ifelse(
any(observations_species$Plot == site & observations_species$week == week & observations_species[[species_col]] > 0),
1, 0
)
}
}
}
dim(detection_array) # Should be num_sites x num_weeks x num_species## [1] 32 36 8Prepare Site Covariates
# Duplicate CameraLoc to be used in Objective 2
CameraLoc_O2 <- CameraLoc
# Standardize the covariates
CameraLoc <- CameraLoc %>%
dplyr::select(-Plot, -Camera, -Lat, -Long, -Status, - Start_Date)
covariates_matrix <- as.matrix(CameraLoc)
rownames(covariates_matrix) <- site_names
# Standardizing covariates
covariates_matrix <- scale(covariates_matrix)
## Correlation
cor_mat <- cor(covariates_matrix)
# Create week matrix for covariate structure [site x week]
week_vals <- unique(observations_species$week)
week_matrix <- matrix(rep(week_vals, each = num_sites), nrow = num_sites, ncol = num_weeks, byrow = FALSE)
# Create detection covariate list
det.covs <- list(
shrub_cover = covariates_matrix[, "total_shrub_cover"],
veg_height = covariates_matrix[, "avg_veg_height"],
week = week_matrix
)
# Remove dash and convert to numeric
week_numeric <- as.numeric(gsub("-", "", det.covs$week))
## Scale and center week_numeric
week_numeric <- scale(week_numeric)
# Reshape into the original 32x36 matrix
det.covs$week <- matrix(week_numeric, nrow = 32, ncol = 36)
str(det.covs)## List of 3
## $ shrub_cover: Named num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
## ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## $ veg_height : Named num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
## ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## $ week : num [1:32, 1:36] -0.613 -0.613 -0.613 -0.613 -0.613 ...Combine data into a list
This requires combining the presence data and the site covariate data into a single list. This also means that the presence data is in a 3-d format.
# Combine into a named list
data_list <- list(
y = detection_array,
occ.covs = covariates_matrix,
det.covs = det.covs,
coords = coords
)
str(data_list)## List of 4
## $ y : num [1:32, 1:36, 1:8] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## .. ..$ : chr [1:36] "2024-29" "2024-30" "2024-31" "2024-32" ...
## .. ..$ : chr [1:8] "Canis_latrans" "Procyon_lotor" "Dasypus_novemcinctus" "Lynx_rufus" ...
## $ occ.covs: num [1:32, 1:12] -0.237 -0.446 -0.337 -0.308 0.039 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## .. ..$ : chr [1:12] "Cogon_Patch_Size" "VegetationDiversity" "PostTreatmentDensities" "Auth" ...
## ..- attr(*, "scaled:center")= Named num [1:12] 458.388 21.875 0.898 2.844 16.188 ...
## .. ..- attr(*, "names")= chr [1:12] "Cogon_Patch_Size" "VegetationDiversity" "PostTreatmentDensities" "Auth" ...
## ..- attr(*, "scaled:scale")= Named num [1:12] 1027.633 6.871 1.232 0.808 0.397 ...
## .. ..- attr(*, "names")= chr [1:12] "Cogon_Patch_Size" "VegetationDiversity" "PostTreatmentDensities" "Auth" ...
## $ det.covs:List of 3
## ..$ shrub_cover: Named num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
## .. ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## ..$ veg_height : Named num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
## .. ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## ..$ week : num [1:32, 1:36] -0.613 -0.613 -0.613 -0.613 -0.613 ...
## $ coords : num [1:32, 1:2] 513432 512936 494295 491358 494388 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## .. ..$ : chr [1:2] "X" "Y"Convert occupancy and detection covariates to a dataframe (from matrix)
I am unsure why I only had an issue with total shrub cover, but this should fix the “cannot find” issue.
# Convert occupancy and detection covariates to a dataframe
data_list[["occ.covs"]] <- as.data.frame(data_list[["occ.covs"]])
data_list[["occ.covs"]]$total_shrub_cover <- as.numeric(data_list[["occ.covs"]]$total_shrub_cover)
#data_list[["det.covs"]] <- as.data.frame(data_list[["det.covs"]])
#data_list[["det.covs"]]$total_shrub_cover <- as.numeric(data_list[["det.covs"]]$total_shrub_cover)
# Make species the first dimension
data_list$y <- aperm(data_list$y, c(3, 1, 2))
dimnames(data_list$y) <- list(species = dimnames(data_list$y)[[1]],
site = dimnames(data_list$y)[[2]],
week = dimnames(data_list$y)[[3]])
str(data_list)## List of 4
## $ y : num [1:8, 1:32, 1:36] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ species: chr [1:8] "Canis_latrans" "Procyon_lotor" "Dasypus_novemcinctus" "Lynx_rufus" ...
## .. ..$ site : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## .. ..$ week : chr [1:36] "2024-29" "2024-30" "2024-31" "2024-32" ...
## $ occ.covs:'data.frame': 32 obs. of 12 variables:
## ..$ Cogon_Patch_Size : num [1:32] -0.237 -0.446 -0.337 -0.308 0.039 ...
## ..$ VegetationDiversity : num [1:32] -0.273 0.455 1.619 -0.273 2.929 ...
## ..$ PostTreatmentDensities: num [1:32] 0.432 -0.729 0.432 2.169 1.13 ...
## ..$ Auth : num [1:32] -2.28 -2.28 -1.04 -1.04 -1.04 ...
## ..$ UTM_Zone : num [1:32] -0.473 -0.473 -0.473 -0.473 -0.473 ...
## ..$ CWD_presence : num [1:32] 1.86 -0.521 -0.521 -0.521 -0.521 ...
## ..$ Veg_shannon_index : num [1:32] 0.6829 0.0427 0.7279 -0.5991 1.1371 ...
## ..$ Avg_Cogongrass_Cover : num [1:32] -0.154 -0.708 0.308 2.045 1.121 ...
## ..$ total_shrub_cover : num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
## ..$ avg_veg_height : num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
## ..$ Tree_Density : num [1:32] -0.3629 -0.3564 -0.5111 3.5896 0.0958 ...
## ..$ Avg_Canopy_Cover : num [1:32] 0.1362 -0.0252 -0.9132 0.782 -1.9627 ...
## $ det.covs:List of 3
## ..$ shrub_cover: Named num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
## .. ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## ..$ veg_height : Named num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
## .. ..- attr(*, "names")= chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## ..$ week : num [1:32, 1:36] -0.613 -0.613 -0.613 -0.613 -0.613 ...
## $ coords : num [1:32, 1:2] 513432 512936 494295 491358 494388 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
## .. ..$ : chr [1:2] "X" "Y"Site Covariates
# Define detection formulas
det.null <- ~ 1
det.full <- ~ shrub_cover + veg_height + week
det.cover <- ~ shrub_cover + veg_height
det.week <- ~ week
det.week.quad <- ~ week + I(week^2)
det.full.quad <- ~ shrub_cover + veg_height + week + I(week^2)
# Define occupancy formulas
occ.null <- ~ 1
occ.full <- ~ #Cogon_Patch_Size +
Veg_shannon_index + total_shrub_cover +
Avg_Cogongrass_Cover +
#Tree_Density +
Avg_Canopy_Cover +
avg_veg_height + CWD_presence
occ.full.quad <- ~ #Cogon_Patch_Size +
Veg_shannon_index + total_shrub_cover +
Avg_Cogongrass_Cover + I(Avg_Cogongrass_Cover^2) +
#Tree_Density +
Avg_Canopy_Cover +
avg_veg_height + CWD_presence
occ.cover <- ~ Avg_Cogongrass_Cover + total_shrub_cover +
avg_veg_height + CWD_presence
occ.canopy <- ~ #Tree_Density +
Avg_Canopy_Cover + Avg_Cogongrass_Cover
occ.canopy.quad <- ~ #Tree_Density +
Avg_Canopy_Cover + Avg_Cogongrass_Cover + I(Avg_Cogongrass_Cover^2)
occ.move <- ~ #Cogon_Patch_Size +
Avg_Cogongrass_Cover + total_shrub_cover
occ.move.quad <- ~ #Cogon_Patch_Size +
Avg_Cogongrass_Cover + I(Avg_Cogongrass_Cover^2) +
total_shrub_cover
occ.forage <- ~ Veg_shannon_index + Avg_Cogongrass_Cover
occ.forage.quad <- ~ Veg_shannon_index + Avg_Cogongrass_Cover + I(Avg_Cogongrass_Cover^2)
occ.cogon <- ~ Avg_Cogongrass_Cover
occ.cogon.quad <- ~ Avg_Cogongrass_Cover + I(Avg_Cogongrass_Cover^2)Priors
priors <- list(
# Community-level effects
beta.comm.normal = list(mean = 0, var = 2.73),
alpha.comm.normal = list(mean = 0, var = 2.73),
# Species-level variance (hierarchical shrinkage)
tau.sq.beta.ig = list(2, 1),
tau.sq.alpha.ig = list(2, 1),
# Spatial variance
sigma.sq.ig = list(3, 2),
# Spatial decay
phi.unif = list(0.05, 3) # 0.05 (5km) local & 3 (300 km) large-scale
)Run the occupancy model
Model Comparison
waicOcc(ms_full_full, by.sp = FALSE)## elpd pD WAIC
## -975.4203 101.7356 2154.3117waicOcc(ms_null_null, by.sp = FALSE)## elpd pD WAIC
## -1040.0062 52.8479 2185.7083waicOcc(ms_cover_forage, by.sp = FALSE)## elpd pD WAIC
## -1000.82875 88.80023 2179.25797waicOcc(ms_cover_move, by.sp = FALSE)## elpd pD WAIC
## -999.42552 89.37232 2177.59568waicOcc(ms_cover_canopy, by.sp = FALSE)## elpd pD WAIC
## -989.0445 83.8199 2145.7288waicOcc(ms_cover_canopyQ, by.sp = FALSE)## elpd pD WAIC
## -985.21130 86.76208 2143.94677waicOcc(ms_weekQ_cogonQ, by.sp = FALSE)## elpd pD WAIC
## -1019.40873 74.94932 2188.71610waicOcc(ms_fullQ_moveQ, by.sp = FALSE)## elpd pD WAIC
## -981.3170 108.6137 2179.8613waicOcc(ms_fullQ_cover, by.sp = FALSE)## elpd pD WAIC
## -985.3002 114.8591 2200.3185waicOcc(ms_fullQ_forage, by.sp = FALSE)## elpd pD WAIC
## -989.0844 105.8171 2189.8031waicOcc(ms_fullQ_fullQ, by.sp = FALSE)## elpd pD WAIC
## -963.5634 111.9125 2150.9517Model Validation
This test explains how well the model fits that data at the community and species level. I believe 0.5 is the target p-value, though how far from this number is considered adequate, I do not know yet. I believe this is a good place to check when thinking about which species we include in the model.
ppc_ft <- ppcOcc(
ms_cover_canopyQ,
fit.stat = "freeman-tukey",
group = 1
)## Currently on species 1 out of 8## Currently on species 2 out of 8## Currently on species 3 out of 8## Currently on species 4 out of 8## Currently on species 5 out of 8## Currently on species 6 out of 8## Currently on species 7 out of 8## Currently on species 8 out of 8ppc_chisq <- ppcOcc(
ms_cover_canopyQ,
fit.stat = "chi-square",
group = 1
)## Currently on species 1 out of 8## Currently on species 2 out of 8## Currently on species 3 out of 8## Currently on species 4 out of 8## Currently on species 5 out of 8## Currently on species 6 out of 8## Currently on species 7 out of 8## Currently on species 8 out of 8summary(ppc_ft)##
## Call:
## ppcOcc(object = ms_cover_canopyQ, fit.stat = "freeman-tukey",
## group = 1)
##
## Samples per Chain: 10000
## Burn-in: 3000
## Thinning Rate: 2
## Number of Chains: 4
## Total Posterior Samples: 14000
##
## ----------------------------------------
## Community Level
## ----------------------------------------
## Bayesian p-value: 0.3655
##
## ----------------------------------------
## Species Level
## ----------------------------------------
## Canis_latrans Bayesian p-value: 0.5658
## Procyon_lotor Bayesian p-value: 0.0901
## Dasypus_novemcinctus Bayesian p-value: 1e-04
## Lynx_rufus Bayesian p-value: 0.5267
## Didelphis_virginiana Bayesian p-value: 0.3401
## Sylvilagus_floridanus Bayesian p-value: 0.4422
## Meleagris_gallopavo Bayesian p-value: 0.6181
## Sciurus_carolinensis Bayesian p-value: 0.3407
## Fit statistic: freeman-tukeysummary(ppc_chisq)##
## Call:
## ppcOcc(object = ms_cover_canopyQ, fit.stat = "chi-square", group = 1)
##
## Samples per Chain: 10000
## Burn-in: 3000
## Thinning Rate: 2
## Number of Chains: 4
## Total Posterior Samples: 14000
##
## ----------------------------------------
## Community Level
## ----------------------------------------
## Bayesian p-value: 0.3044
##
## ----------------------------------------
## Species Level
## ----------------------------------------
## Canis_latrans Bayesian p-value: 0.4676
## Procyon_lotor Bayesian p-value: 0.0028
## Dasypus_novemcinctus Bayesian p-value: 0
## Lynx_rufus Bayesian p-value: 0.4517
## Didelphis_virginiana Bayesian p-value: 0.2636
## Sylvilagus_floridanus Bayesian p-value: 0.2392
## Meleagris_gallopavo Bayesian p-value: 0.6877
## Sciurus_carolinensis Bayesian p-value: 0.3226
## Fit statistic: chi-squareProability of Direction- Occupancy
pd <- function(x) {
max(mean(x > 0), mean(x < 0))
}
pd_table <- apply(
ms_cover_canopyQ$beta.samples,
2,
pd
)
head(pd_table)## (Intercept)-Canis_latrans (Intercept)-Procyon_lotor
## 0.7610714 0.5065000
## (Intercept)-Dasypus_novemcinctus (Intercept)-Lynx_rufus
## 0.9718571 0.8350714
## (Intercept)-Didelphis_virginiana (Intercept)-Sylvilagus_floridanus
## 0.9902857 0.9268571beta_means <- colMeans(ms_cover_canopyQ$beta.samples)
beta_samps <- ms_cover_canopyQ$beta.samples
pd_vals <- apply(beta_samps, 2, pd)
beta_means <- colMeans(beta_samps)
results_species <- data.frame(
parameter = colnames(beta_samps),
mean = beta_means,
pd = pd_vals,
pd_95 = pd_vals >= 0.95
)
head(results_species)## parameter
## (Intercept)-Canis_latrans (Intercept)-Canis_latrans
## (Intercept)-Procyon_lotor (Intercept)-Procyon_lotor
## (Intercept)-Dasypus_novemcinctus (Intercept)-Dasypus_novemcinctus
## (Intercept)-Lynx_rufus (Intercept)-Lynx_rufus
## (Intercept)-Didelphis_virginiana (Intercept)-Didelphis_virginiana
## (Intercept)-Sylvilagus_floridanus (Intercept)-Sylvilagus_floridanus
## mean pd pd_95
## (Intercept)-Canis_latrans -0.414657186 0.7610714 FALSE
## (Intercept)-Procyon_lotor 0.003107734 0.5065000 FALSE
## (Intercept)-Dasypus_novemcinctus -1.096782013 0.9718571 TRUE
## (Intercept)-Lynx_rufus -0.704572513 0.8350714 FALSE
## (Intercept)-Didelphis_virginiana -1.530519435 0.9902857 TRUE
## (Intercept)-Sylvilagus_floridanus -0.958482335 0.9268571 FALSEbeta_comm <- ms_cover_canopyQ$beta.comm.samples
results_comm <- data.frame(
parameter = colnames(beta_comm),
mean = colMeans(beta_comm),
pd = apply(beta_comm, 2, pd)
)
results_comm$pd_95 <- results_comm$pd >= 0.95Probability of Direction- Detection
pd <- function(x) {
max(mean(x > 0), mean(x < 0))
}
alpha_samps <- ms_cover_canopyQ$alpha.samples
pd_vals_alpha <- apply(alpha_samps, 2, pd)
alpha_means <- colMeans(alpha_samps)
results_species_det <- data.frame(
parameter = colnames(alpha_samps),
mean = alpha_means,
pd = pd_vals_alpha,
pd_95 = pd_vals_alpha >= 0.95
)
head(results_species_det)## parameter mean
## (Intercept)-Canis_latrans (Intercept)-Canis_latrans -2.769030
## (Intercept)-Procyon_lotor (Intercept)-Procyon_lotor -2.307025
## (Intercept)-Dasypus_novemcinctus (Intercept)-Dasypus_novemcinctus -1.807891
## (Intercept)-Lynx_rufus (Intercept)-Lynx_rufus -3.581237
## (Intercept)-Didelphis_virginiana (Intercept)-Didelphis_virginiana -2.696502
## (Intercept)-Sylvilagus_floridanus (Intercept)-Sylvilagus_floridanus -3.133887
## pd pd_95
## (Intercept)-Canis_latrans 1 TRUE
## (Intercept)-Procyon_lotor 1 TRUE
## (Intercept)-Dasypus_novemcinctus 1 TRUE
## (Intercept)-Lynx_rufus 1 TRUE
## (Intercept)-Didelphis_virginiana 1 TRUE
## (Intercept)-Sylvilagus_floridanus 1 TRUEalpha_comm <- ms_cover_canopyQ$alpha.comm.samples
results_comm_det <- data.frame(
parameter = colnames(alpha_comm),
mean = colMeans(alpha_comm),
pd = apply(alpha_comm, 2, pd)
)
results_comm_det$pd_95 <- results_comm_det$pd >= 0.95
results_comm_det## parameter mean pd pd_95
## (Intercept) (Intercept) -2.78647331 1.0000000 TRUE
## shrub_cover shrub_cover 0.33689343 0.8709286 FALSE
## veg_height veg_height 0.07844567 0.6431429 FALSEPlot Species Occupancy
Batched Parameter Estimates
summary(ms_cover_canopyQ) # Summary of parameter estimates##
## Call:
## spMsPGOcc(occ.formula = occ.canopy.quad, det.formula = det.cover,
## data = data_list, priors = priors, cov.model = "exponential",
## NNGP = FALSE, n.batch = 1000, batch.length = 10, n.report = 100,
## n.burn = 3000, n.thin = 2, n.chains = 4)
##
## Samples per Chain: 10000
## Burn-in: 3000
## Thinning Rate: 2
## Number of Chains: 4
## Total Posterior Samples: 14000
## Run Time (min): 1.686
##
## ----------------------------------------
## Community Level
## ----------------------------------------
## Occurrence Means (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) -0.8348 0.4799 -1.7475 -0.8410 0.1304 1.0061 3039
## Avg_Canopy_Cover 1.3330 0.4323 0.5176 1.3191 2.2303 1.0033 3860
## Avg_Cogongrass_Cover -0.1475 0.4965 -1.1410 -0.1455 0.8286 1.0025 2580
## I(Avg_Cogongrass_Cover^2) 0.7980 0.4510 -0.0356 0.7813 1.7387 1.0025 1921
##
## Occurrence Variances (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) 0.8281 0.6418 0.2072 0.6550 2.4488 1.0094 4820
## Avg_Canopy_Cover 0.7977 0.6629 0.1958 0.6068 2.5021 1.0048 3796
## Avg_Cogongrass_Cover 0.6373 0.5183 0.1695 0.4973 1.9179 1.0073 5874
## I(Avg_Cogongrass_Cover^2) 0.6622 0.5877 0.1722 0.4955 2.1950 1.0072 2640
##
## Detection Means (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) -2.7865 0.298 -3.3724 -2.7880 -2.1852 1.0022 7055
## shrub_cover 0.3369 0.307 -0.2659 0.3350 0.9452 1.0021 6501
## veg_height 0.0784 0.233 -0.3852 0.0807 0.5371 1.0035 7592
##
## Detection Variances (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) 0.6163 0.3740 0.2192 0.5242 1.5639 1.0018 5522
## shrub_cover 0.6572 0.3914 0.2288 0.5589 1.6848 1.0033 6081
## veg_height 0.3673 0.2037 0.1398 0.3164 0.8927 1.0012 9850
##
## ----------------------------------------
## Species Level
## ----------------------------------------
## Occurrence (logit scale):
## Mean SD 2.5% 50%
## (Intercept)-Canis_latrans -0.4147 0.6191 -1.5843 -0.4304
## (Intercept)-Procyon_lotor 0.0031 0.6118 -1.1705 -0.0105
## (Intercept)-Dasypus_novemcinctus -1.0968 0.5913 -2.2857 -1.0902
## (Intercept)-Lynx_rufus -0.7046 0.7806 -2.1275 -0.7452
## (Intercept)-Didelphis_virginiana -1.5305 0.6819 -2.9644 -1.5023
## (Intercept)-Sylvilagus_floridanus -0.9585 0.6727 -2.2790 -0.9622
## (Intercept)-Meleagris_gallopavo -0.3392 0.7815 -1.7523 -0.3822
## (Intercept)-Sciurus_carolinensis -1.8740 0.7461 -3.4536 -1.8320
## Avg_Canopy_Cover-Canis_latrans 0.2814 0.5143 -0.7192 0.2821
## Avg_Canopy_Cover-Procyon_lotor 1.1639 0.5412 0.1637 1.1463
## Avg_Canopy_Cover-Dasypus_novemcinctus 1.3725 0.5649 0.3728 1.3366
## Avg_Canopy_Cover-Lynx_rufus 0.9873 0.7473 -0.3879 0.9577
## Avg_Canopy_Cover-Didelphis_virginiana 1.7583 0.6905 0.6061 1.6822
## Avg_Canopy_Cover-Sylvilagus_floridanus 2.1482 0.8632 0.7739 2.0281
## Avg_Canopy_Cover-Meleagris_gallopavo 1.5196 0.7210 0.2554 1.4661
## Avg_Canopy_Cover-Sciurus_carolinensis 1.8315 0.7405 0.6274 1.7485
## Avg_Cogongrass_Cover-Canis_latrans 0.0858 0.6652 -1.1803 0.0703
## Avg_Cogongrass_Cover-Procyon_lotor -0.2188 0.6375 -1.4617 -0.2134
## Avg_Cogongrass_Cover-Dasypus_novemcinctus 0.2066 0.6644 -1.0615 0.1940
## Avg_Cogongrass_Cover-Lynx_rufus -0.2705 0.7411 -1.7787 -0.2655
## Avg_Cogongrass_Cover-Didelphis_virginiana 0.3449 0.7216 -1.0083 0.3197
## Avg_Cogongrass_Cover-Sylvilagus_floridanus -0.7493 0.7480 -2.3211 -0.7166
## Avg_Cogongrass_Cover-Meleagris_gallopavo -0.5971 0.8491 -2.3696 -0.5583
## Avg_Cogongrass_Cover-Sciurus_carolinensis -0.0124 0.7171 -1.4114 -0.0226
## I(Avg_Cogongrass_Cover^2)-Canis_latrans 1.4361 0.8259 0.0982 1.3415
## I(Avg_Cogongrass_Cover^2)-Procyon_lotor 1.2205 0.7422 0.1007 1.1124
## I(Avg_Cogongrass_Cover^2)-Dasypus_novemcinctus 0.6207 0.5034 -0.2861 0.5903
## I(Avg_Cogongrass_Cover^2)-Lynx_rufus 1.2928 0.6594 0.2039 1.2205
## I(Avg_Cogongrass_Cover^2)-Didelphis_virginiana 0.2780 0.5499 -0.7033 0.2462
## I(Avg_Cogongrass_Cover^2)-Sylvilagus_floridanus 0.5932 0.5654 -0.3990 0.5590
## I(Avg_Cogongrass_Cover^2)-Meleagris_gallopavo 0.2881 0.8905 -1.2980 0.2035
## I(Avg_Cogongrass_Cover^2)-Sciurus_carolinensis 0.8451 0.5029 -0.0580 0.8162
## 97.5% Rhat ESS
## (Intercept)-Canis_latrans 0.8604 1.0011 4274
## (Intercept)-Procyon_lotor 1.2408 1.0040 4559
## (Intercept)-Dasypus_novemcinctus 0.0347 1.0030 4975
## (Intercept)-Lynx_rufus 0.9836 1.0098 2347
## (Intercept)-Didelphis_virginiana -0.2532 1.0004 4379
## (Intercept)-Sylvilagus_floridanus 0.3732 1.0023 4082
## (Intercept)-Meleagris_gallopavo 1.3053 1.0043 2714
## (Intercept)-Sciurus_carolinensis -0.5160 1.0022 3742
## Avg_Canopy_Cover-Canis_latrans 1.2975 1.0030 6472
## Avg_Canopy_Cover-Procyon_lotor 2.2726 1.0024 6728
## Avg_Canopy_Cover-Dasypus_novemcinctus 2.5971 1.0015 5573
## Avg_Canopy_Cover-Lynx_rufus 2.5498 1.0019 3117
## Avg_Canopy_Cover-Didelphis_virginiana 3.3557 1.0012 4119
## Avg_Canopy_Cover-Sylvilagus_floridanus 4.1489 1.0056 2727
## Avg_Canopy_Cover-Meleagris_gallopavo 3.1071 1.0018 4073
## Avg_Canopy_Cover-Sciurus_carolinensis 3.5361 1.0030 3163
## Avg_Cogongrass_Cover-Canis_latrans 1.4620 1.0016 4581
## Avg_Cogongrass_Cover-Procyon_lotor 1.0409 1.0028 5015
## Avg_Cogongrass_Cover-Dasypus_novemcinctus 1.5903 1.0023 4811
## Avg_Cogongrass_Cover-Lynx_rufus 1.1531 1.0022 3692
## Avg_Cogongrass_Cover-Didelphis_virginiana 1.8334 1.0009 4467
## Avg_Cogongrass_Cover-Sylvilagus_floridanus 0.6224 1.0006 3841
## Avg_Cogongrass_Cover-Meleagris_gallopavo 0.9665 1.0074 2150
## Avg_Cogongrass_Cover-Sciurus_carolinensis 1.4287 1.0003 3917
## I(Avg_Cogongrass_Cover^2)-Canis_latrans 3.2904 1.0002 2142
## I(Avg_Cogongrass_Cover^2)-Procyon_lotor 3.0318 1.0031 2059
## I(Avg_Cogongrass_Cover^2)-Dasypus_novemcinctus 1.6963 1.0011 4401
## I(Avg_Cogongrass_Cover^2)-Lynx_rufus 2.8077 1.0002 2882
## I(Avg_Cogongrass_Cover^2)-Didelphis_virginiana 1.5243 1.0013 2847
## I(Avg_Cogongrass_Cover^2)-Sylvilagus_floridanus 1.8269 1.0079 2685
## I(Avg_Cogongrass_Cover^2)-Meleagris_gallopavo 2.2669 1.0098 910
## I(Avg_Cogongrass_Cover^2)-Sciurus_carolinensis 1.9109 1.0005 4141
##
## Detection (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat
## (Intercept)-Canis_latrans -2.7690 0.1813 -3.1365 -2.7629 -2.4264 1.0003
## (Intercept)-Procyon_lotor -2.3070 0.1437 -2.6018 -2.3029 -2.0359 1.0004
## (Intercept)-Dasypus_novemcinctus -1.8079 0.1627 -2.1365 -1.8054 -1.4975 1.0006
## (Intercept)-Lynx_rufus -3.5812 0.3265 -4.2425 -3.5728 -2.9638 1.0105
## (Intercept)-Didelphis_virginiana -2.6965 0.2846 -3.2790 -2.6902 -2.1545 1.0006
## (Intercept)-Sylvilagus_floridanus -3.1339 0.2566 -3.6593 -3.1299 -2.6515 1.0038
## (Intercept)-Meleagris_gallopavo -3.8233 0.4392 -4.7133 -3.8176 -2.9684 1.0117
## (Intercept)-Sciurus_carolinensis -2.7728 0.3128 -3.4260 -2.7627 -2.1920 1.0027
## shrub_cover-Canis_latrans -0.2635 0.2175 -0.6870 -0.2645 0.1624 1.0086
## shrub_cover-Procyon_lotor 0.2494 0.1665 -0.0841 0.2546 0.5678 1.0027
## shrub_cover-Dasypus_novemcinctus 0.9859 0.3046 0.3979 0.9871 1.5818 1.0030
## shrub_cover-Lynx_rufus -0.2497 0.3532 -0.9395 -0.2508 0.4480 1.0043
## shrub_cover-Didelphis_virginiana 1.1357 0.3675 0.4446 1.1283 1.8729 1.0019
## shrub_cover-Sylvilagus_floridanus 0.4663 0.4052 -0.3344 0.4664 1.2535 1.0074
## shrub_cover-Meleagris_gallopavo -0.6221 0.3873 -1.3923 -0.6198 0.1395 1.0062
## shrub_cover-Sciurus_carolinensis 1.0483 0.3871 0.2984 1.0493 1.8141 1.0032
## veg_height-Canis_latrans -0.6377 0.1827 -0.9995 -0.6359 -0.2896 1.0037
## veg_height-Procyon_lotor 0.3551 0.1233 0.1157 0.3545 0.5988 1.0009
## veg_height-Dasypus_novemcinctus 0.2756 0.1408 0.0052 0.2757 0.5552 1.0005
## veg_height-Lynx_rufus 0.0822 0.2598 -0.4454 0.0882 0.5721 1.0041
## veg_height-Didelphis_virginiana 0.5004 0.2593 0.0129 0.4959 1.0251 1.0025
## veg_height-Sylvilagus_floridanus 0.1510 0.2588 -0.3486 0.1462 0.6567 1.0035
## veg_height-Meleagris_gallopavo -0.2008 0.4313 -1.0323 -0.2063 0.6600 1.0103
## veg_height-Sciurus_carolinensis 0.1368 0.2356 -0.3149 0.1334 0.6111 1.0000
## ESS
## (Intercept)-Canis_latrans 3500
## (Intercept)-Procyon_lotor 5148
## (Intercept)-Dasypus_novemcinctus 5569
## (Intercept)-Lynx_rufus 1604
## (Intercept)-Didelphis_virginiana 3051
## (Intercept)-Sylvilagus_floridanus 3229
## (Intercept)-Meleagris_gallopavo 1057
## (Intercept)-Sciurus_carolinensis 2927
## shrub_cover-Canis_latrans 4289
## shrub_cover-Procyon_lotor 6136
## shrub_cover-Dasypus_novemcinctus 4248
## shrub_cover-Lynx_rufus 2316
## shrub_cover-Didelphis_virginiana 2651
## shrub_cover-Sylvilagus_floridanus 2213
## shrub_cover-Meleagris_gallopavo 1507
## shrub_cover-Sciurus_carolinensis 3106
## veg_height-Canis_latrans 3789
## veg_height-Procyon_lotor 7420
## veg_height-Dasypus_novemcinctus 7933
## veg_height-Lynx_rufus 2985
## veg_height-Didelphis_virginiana 3820
## veg_height-Sylvilagus_floridanus 3325
## veg_height-Meleagris_gallopavo 1218
## veg_height-Sciurus_carolinensis 4085
##
## Spatial Covariance:
## Mean SD 2.5% 50% 97.5% Rhat ESS
## sigma.sq-Canis_latrans 1.1340 1.2681 0.2817 0.8067 4.0904 1.0247 1537
## sigma.sq-Procyon_lotor 1.1685 1.1554 0.2911 0.8265 4.4125 1.0205 1349
## sigma.sq-Dasypus_novemcinctus 1.2158 1.3396 0.2916 0.8337 4.6283 1.0193 1162
## sigma.sq-Lynx_rufus 1.0879 1.1693 0.2819 0.7992 3.5863 1.0766 1287
## sigma.sq-Didelphis_virginiana 0.9417 0.7785 0.2661 0.7244 2.9379 1.0056 1920
## sigma.sq-Sylvilagus_floridanus 0.9321 0.8017 0.2627 0.7221 2.9169 1.0094 1601
## sigma.sq-Meleagris_gallopavo 1.0452 1.0204 0.2824 0.7670 3.5113 1.0056 1448
## sigma.sq-Sciurus_carolinensis 0.8424 0.5909 0.2646 0.6811 2.3655 1.0059 2614
## phi-Canis_latrans 1.5288 0.8524 0.1366 1.5219 2.9314 1.0007 6701
## phi-Procyon_lotor 1.5234 0.8514 0.1290 1.5287 2.9213 1.0001 6115
## phi-Dasypus_novemcinctus 1.5389 0.8457 0.1267 1.5446 2.9244 1.0012 6129
## phi-Lynx_rufus 1.5134 0.8405 0.1194 1.5291 2.9115 1.0009 6393
## phi-Didelphis_virginiana 1.5305 0.8481 0.1297 1.5360 2.9234 1.0020 6449
## phi-Sylvilagus_floridanus 1.5186 0.8562 0.1303 1.5082 2.9280 1.0008 5909
## phi-Meleagris_gallopavo 1.5194 0.8497 0.1219 1.4987 2.9198 1.0018 6741
## phi-Sciurus_carolinensis 1.5414 0.8596 0.1167 1.5489 2.9355 1.0025 5971names(ms_cover_canopyQ)## [1] "rhat" "beta.comm.samples" "alpha.comm.samples"
## [4] "tau.sq.beta.samples" "tau.sq.alpha.samples" "beta.samples"
## [7] "alpha.samples" "theta.samples" "z.samples"
## [10] "psi.samples" "like.samples" "w.samples"
## [13] "X.p" "X.p.re" "X.re"
## [16] "ESS" "X" "y"
## [19] "call" "n.samples" "x.names"
## [22] "sp.names" "x.p.names" "theta.names"
## [25] "type" "coords" "cov.model.indx"
## [28] "n.post" "n.thin" "n.burn"
## [31] "n.chains" "pRE" "psiRE"
## [34] "run.time"str(ms_cover_canopyQ$beta.samples)## 'mcmc' num [1:14000, 1:32] -0.741 -0.113 -0.209 -1.137 0.156 ...
## - attr(*, "mcpar")= num [1:3] 1 14000 1
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:32] "(Intercept)-Canis_latrans" "(Intercept)-Procyon_lotor" "(Intercept)-Dasypus_novemcinctus" "(Intercept)-Lynx_rufus" ...mean(ms_cover_canopyQ$beta.samples[, 5] > 0)## [1] 0.009714286MCMCplot(ms_cover_canopyQ$beta.samples, ref_ovl = TRUE, ci = c(50, 95)) # Occupancy
MCMCplot(ms_cover_canopyQ$alpha.samples, ref_ovl = TRUE, ci = c(50, 95)) # Detection
## Occupancy
# Total number of parameters
n_params <- ncol(ms_cover_canopyQ$beta.samples)
# Number of parameters to plot at a time
chunk_size <- 10
# Split and plot
for (i in seq(1, n_params, by = chunk_size)) {
end <- min(i + chunk_size - 1, n_params)
param_names <- colnames(ms_cover_canopyQ$beta.samples)[i:end]
# Set filename
file_name <- paste0("MCMCplot_Occupancy_Params_", i, "_to_", end, ".png")
# Save plot to PNG
png(filename = file_name, width = 1200, height = 800, res = 150)
MCMCplot(ms_cover_canopyQ$beta.samples[, param_names],
ref_ovl = TRUE,
ci = c(50, 95),
main = paste0("Occupancy Parameters: ", i, " to ", end))
dev.off()
}
## Detection
# Number of parameters
n_params <- ncol(ms_cover_canopyQ$alpha.samples)
# Number of parameters to plot at a time
chunk_size <- 10
# Split and plot
for (i in seq(1, n_params, by = chunk_size)) {
end <- min(i + chunk_size - 1, n_params)
param_names <- colnames(ms_cover_canopyQ$alpha.samples)[i:end]
# Set filename
file_name <- paste0("MCMCplot_Detection_Params_", i, "_to_", end, ".png")
# Save plot to PNG
png(filename = file_name, width = 1200, height = 800, res = 150)
MCMCplot(ms_cover_canopyQ$alpha.samples[, param_names],
ref_ovl = TRUE,
ci = c(50, 95),
main = paste0("Detection Parameters: ", i, " to ", end))
dev.off()
}Plot Species level relationship with cogongrass cover
# Collecting scaled atrributes
cogon.mean <- attr(covariates_matrix, "scaled:center")["Avg_Cogongrass_Cover"]
cogon.sd <- attr(covariates_matrix, "scaled:scale")["Avg_Cogongrass_Cover"]
cogon.mean <- mean(CameraLoc$Avg_Cogongrass_Cover)
cogon.sd <- sd(CameraLoc$Avg_Cogongrass_Cover)
# Unscaled prediction values
cogon.pred.vals <- seq(0, 100, length.out = 100)
# Scale prediction values
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2
canopy.mean <- attr(covariates_matrix, "scaled:center")["Avg_Canopy_Cover"]
canopy.sd <- attr(covariates_matrix, "scaled:scale")["Avg_Canopy_Cover"]
pred.df <- as.matrix(data.frame(
intercept = 1,
Avg_Cogongrass_Cover = cogon.pred.scale,
I.Avg_Cogongrass_Cover.2. = cogon.pred.scale2,
Avg_Canopy_Cover = 0
))
# Mean coordinates
mean_coords <- matrix(
colMeans(ms_cover_canopyQ$coords),
nrow = nrow(pred.df),
ncol = ncol(ms_cover_canopyQ$coords),
byrow = TRUE
)
out.pred.with.spatial <- predict(ms_cover_canopyQ,
X.0 = pred.df,
coords.0 = mean_coords,
ignore.RE = FALSE) # Include spatial effects## ----------------------------------------
## Prediction description
## ----------------------------------------
## Spatial Multispecies Occupancy model with Polya-Gamma latent
## variable fit with 32 observations.
##
## Number of covariates 4 (including intercept if specified).
##
## Using the exponential spatial correlation model.
##
## Number of MCMC samples 14000.
##
## Predicting at 100 non-sampled locations.
##
##
## Source compiled with OpenMP support and model fit using 1 threads.
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 0 of 8, 0.00%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 1 of 8, 12.50%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 2 of 8, 25.00%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 3 of 8, 37.50%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 4 of 8, 50.00%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 5 of 8, 62.50%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 6 of 8, 75.00%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 7 of 8, 87.50%
## Species: 8 of 8, 100.00%# Extract the species-specific betas for Coyote
coyote_index <- which(ms_cover_canopyQ$sp.names == "Canis_latrans")
# Get column indices for the parameters we need
intercept_col <- paste0("(Intercept)-", "Canis_latrans")
cogon_col <- paste0("Avg_Cogongrass_Cover-", "Canis_latrans")
cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", "Canis_latrans")
canopy_col <- paste0("Avg_Canopy_Cover-", "Canis_latrans")
# Extract posterior samples for each coefficient
beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]
# Create prediction gradient (0 to 100%)
cogon.pred.vals <- seq(0, 100, length.out = 1000)
# Scale using same transformation
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2
n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)
logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
for(i in 1:n_samples) {
logit_psi[i, ] <- beta_intercept[i] +
beta_cogon[i] * cogon.pred.scale +
beta_cogon2[i] * cogon.pred.scale2 +
beta_canopy[i] * 0 # Setting to mean (0 on scaled)
}
# Convert to probability scale
psi_samples <- plogis(logit_psi)
# Calculate quantiles across MCMC samples for each prediction point
psi.coyote.quants <- apply(psi_samples, 2, quantile,
probs = c(0.025, 0.5, 0.975))
psi.coyote.dat <- data.frame(
Avg_Cogongrass_Cover = cogon.pred.vals,
psi.low = psi.coyote.quants[1, ],
psi.med = psi.coyote.quants[2, ],
psi.high = psi.coyote.quants[3, ]
)
ggplot(psi.coyote.dat, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
geom_line(size = 1.2) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Occupancy Probability (Canis latrans)") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
summary(psi.coyote.dat)## Avg_Cogongrass_Cover psi.low psi.med psi.high
## Min. : 0 Min. :0.1669 Min. :0.3938 Min. :0.6836
## 1st Qu.: 25 1st Qu.:0.2599 1st Qu.:0.5426 1st Qu.:0.8479
## Median : 50 Median :0.6118 Median :0.9885 Median :1.0000
## Mean : 50 Mean :0.5808 Mean :0.8140 Mean :0.9261
## 3rd Qu.: 75 3rd Qu.:0.8985 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :100 Max. :0.9798 Max. :1.0000 Max. :1.0000Community Occupancy
# Community-level (average across species)
comm_intercept <- ms_cover_canopyQ$beta.comm.samples[, "(Intercept)"]
comm_cogon <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Cogongrass_Cover"]
comm_cogon2 <- ms_cover_canopyQ$beta.comm.samples[, "I(Avg_Cogongrass_Cover^2)"]
comm_canopy <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Canopy_Cover"]
# Same prediction approach
logit_psi_comm <- matrix(NA, nrow = nrow(ms_cover_canopyQ$beta.comm.samples),
ncol = n_pred)
for(i in 1:nrow(ms_cover_canopyQ$beta.comm.samples)) {
logit_psi_comm[i, ] <- comm_intercept[i] +
comm_cogon[i] * cogon.pred.scale +
comm_cogon2[i] * cogon.pred.scale2 +
comm_canopy[i] * 0
}
psi_comm_samples <- plogis(logit_psi_comm)
psi.comm.quants <- apply(psi_comm_samples, 2, quantile,
probs = c(0.025, 0.5, 0.975))
psi.comm.dat <- data.frame(
Avg_Cogongrass_Cover = cogon.pred.vals,
psi.low = psi.comm.quants[1, ],
psi.med = psi.comm.quants[2, ],
psi.high = psi.comm.quants[3, ]
)
ggplot(psi.comm.dat, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
geom_line(size = 1.2) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Community Occupancy Probability") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100))
Loop through all species in the community
# Prediction gradient
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2
n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)
# List to store results
species_predictions <- list()
species_plots <- list()
# Loop through all species
for(sp in ms_cover_canopyQ$sp.names) {
# Get column names for this species
intercept_col <- paste0("(Intercept)-", sp)
cogon_col <- paste0("Avg_Cogongrass_Cover-", sp)
cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", sp)
canopy_col <- paste0("Avg_Canopy_Cover-", sp)
# Extract posterior samples
beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]
# Calculate linear predictor
logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
for(i in 1:n_samples) {
logit_psi[i, ] <- beta_intercept[i] +
beta_cogon[i] * cogon.pred.scale +
beta_cogon2[i] * cogon.pred.scale2 +
beta_canopy[i] * 0
}
# Convert to probability
psi_samples <- plogis(logit_psi)
# Calculate quantiles
psi.quants <- apply(psi_samples, 2, quantile, probs = c(0.025, 0.5, 0.975))
# Store results
species_predictions[[sp]] <- data.frame(
species = sp,
Avg_Cogongrass_Cover = cogon.pred.vals,
psi.low = psi.quants[1, ],
psi.med = psi.quants[2, ],
psi.high = psi.quants[3, ]
)
# Create plot
species_plots[[sp]] <- ggplot(species_predictions[[sp]],
aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
geom_line(size = 1.2) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Occupancy Probability",
title = gsub("_", " ", sp)) + # Replace underscores with spaces
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100))
}
# View individual plots
species_plots[["Canis_latrans"]]
species_plots[["Procyon_lotor"]]
# Combine all plots into a grid
all_species_plot <- wrap_plots(species_plots, ncol = 4)
all_species_plot
# Save combined plot
ggsave("all_species_cogongrass_response.png",
all_species_plot,
width = 12, height = 16, dpi = 300)
# Combine all predictions into one dataframe for faceted plot
all_predictions <- do.call(rbind, species_predictions)
# Create faceted plot
faceted_plot <- ggplot(all_predictions, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
geom_line(size = 1.2) +
facet_wrap(~ species, ncol = 2, labeller = labeller(species = function(x) gsub("_", " ", x))) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Occupancy Probability") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100)) +
theme(strip.text = element_text(face = "italic")) # Italicize species names
faceted_plot
ggsave("all_species_faceted.png",
faceted_plot,
width = 10, height = 12, dpi = 300)Species Lookup
SpeciesLookup <- tibble::tribble(
~Name, ~Scientific_Name, ~Common_Name,
"Canis_latrans", "Canis latrans", "Coyote",
"Dasypus_novemcinctus", "Dasypus novemcinctus", "Nine-banded Armadillo",
"Didelphis_virginiana", "Didelphis virginiana", "Virginia Opossum",
"Lynx_rufus", "Lynx rufus", "Bobcat",
"Meleagris_gallopavo", "Meleagris gallopavo", "Wild Turkey",
"Procyon_lotor", "Procyon lotor", "Raccoon",
"Sciurus_carolinensis", "Sciurus carolinensis", "Eastern Gray Squirrel",
"Sylvilagus_floridanus", "Sylvilagus floridanus", "Eastern Cottontail Rabbit",
)Varying cogongrass and canopy cover
# Prediction gradient for cogongrass
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2
# Define canopy cover scenarios with factor levels
canopy_scenarios <- data.frame(
scenario = factor(c("Low", "Average", "High"),
levels = c("Low", "Average", "High")),
canopy_scaled = c(-1, 0, 1),
color = c("#D55E00", "#009E73", "#0072B2")
)
n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)
# Lists for results
species_predictions <- list()
species_plots <- list()
# Loop through all species
for(sp in ms_cover_canopyQ$sp.names) {
# Get column names for this species
intercept_col <- paste0("(Intercept)-", sp)
cogon_col <- paste0("Avg_Cogongrass_Cover-", sp)
cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", sp)
canopy_col <- paste0("Avg_Canopy_Cover-", sp)
# Extract posterior samples
beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]
# Store predictions for each canopy scenario
scenario_predictions <- list()
for(j in 1:nrow(canopy_scenarios)) {
canopy_val <- canopy_scenarios$canopy_scaled[j]
# Calculate linear predictor
logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
for(i in 1:n_samples) {
logit_psi[i, ] <- beta_intercept[i] +
beta_cogon[i] * cogon.pred.scale +
beta_cogon2[i] * cogon.pred.scale2 +
beta_canopy[i] * canopy_val
}
# Convert to probability
psi_samples <- plogis(logit_psi)
# Calculate quantiles
psi.quants <- apply(psi_samples, 2, quantile, probs = c(0.025, 0.5, 0.975))
# Store results
scenario_predictions[[j]] <- data.frame(
species = sp,
scenario = canopy_scenarios$scenario[j],
Avg_Cogongrass_Cover = cogon.pred.vals,
psi.low = psi.quants[1, ],
psi.med = psi.quants[2, ],
psi.high = psi.quants[3, ]
)
}
# Combine all scenarios for this species
species_predictions[[sp]] <- do.call(rbind, scenario_predictions)
all_predictions <- all_predictions %>%
left_join(SpeciesLookup, by = c("species" = "Name"))
# Create plot for this species
species_plots[[sp]] <- ggplot(species_predictions[[sp]],
aes(x = Avg_Cogongrass_Cover, y = psi.med,
color = scenario, fill = scenario)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
geom_line(size = 1.2) +
scale_color_manual(values = canopy_scenarios$color) +
scale_fill_manual(values = canopy_scenarios$color) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Occupancy Probability",
title = gsub("_", " ", sp),
color = "Canopy Cover",
fill = "Canopy Cover") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100)) +
theme(legend.position = "bottom")
}
# View individual plots
species_plots[["Canis_latrans"]]
species_plots[["Procyon_lotor"]]
# Combine all plots into a grid
all_species_plot <- wrap_plots(species_plots, ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
all_species_plot
ggsave("all_species_cogongrass_canopy.png",
all_species_plot,
width = 14, height = 16, dpi = 300)
# Combine all predictions into one dataframe
all_predictions <- do.call(rbind, species_predictions)
all_predictions <- all_predictions %>%
left_join(
SpeciesLookup %>% dplyr::select(Name, Common_Name),
by = c("species" = "Name")
)
all_predictions$scenario <- factor(all_predictions$scenario,
levels = c("Low", "Average", "High"))
# Create faceted plot
faceted_plot <- ggplot(all_predictions,
aes(x = Avg_Cogongrass_Cover, y = psi.med,
color = scenario, fill = scenario)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
geom_line(size = 1.2) +
facet_wrap(~ Common_Name, ncol = 2) +
scale_color_manual(values = canopy_scenarios$color) +
scale_fill_manual(values = canopy_scenarios$color) +
theme_bw() +
labs(x = "Average Cogongrass Cover (%)",
y = "Occupancy Probability",
color = "Canopy Cover",
fill = "Canopy Cover") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100)) +
theme(strip.text = element_text(face = "bold"),
legend.position = "bottom")
faceted_plot
ggsave("all_species_faceted_canopy.png",
faceted_plot,
width = 12, height = 14, dpi = 600)Community Occupancy Plot- Cogongrass and Canopy Cover
# Prediction gradient for cogongrass
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2
# Define canopy cover scenarios
canopy_scenarios <- data.frame(
scenario = c("Low", "Average", "High"),
canopy_scaled = c(-1, 0, 1),
color = c("#D55E00", "#009E73", "#0072B2")
)
n_pred <- length(cogon.pred.vals)
# Get community-level parameters
comm_intercept <- ms_cover_canopyQ$beta.comm.samples[, "(Intercept)"]
comm_cogon <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Cogongrass_Cover"]
comm_cogon2 <- ms_cover_canopyQ$beta.comm.samples[, "I(Avg_Cogongrass_Cover^2)"]
comm_canopy <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Canopy_Cover"]
# Store predictions for each canopy scenario
community_predictions <- list()
for(j in 1:nrow(canopy_scenarios)) {
canopy_val <- canopy_scenarios$canopy_scaled[j]
# Calculate linear predictor
logit_psi_comm <- matrix(NA,
nrow = nrow(ms_cover_canopyQ$beta.comm.samples),
ncol = n_pred)
for(i in 1:nrow(ms_cover_canopyQ$beta.comm.samples)) {
logit_psi_comm[i, ] <- comm_intercept[i] +
comm_cogon[i] * cogon.pred.scale +
comm_cogon2[i] * cogon.pred.scale2 +
comm_canopy[i] * canopy_val
}
# Convert to probability
psi_comm_samples <- plogis(logit_psi_comm)
# Calculate quantiles
psi.comm.quants <- apply(psi_comm_samples, 2, quantile,
probs = c(0.025, 0.5, 0.975))
# Store results
community_predictions[[j]] <- data.frame(
scenario = canopy_scenarios$scenario[j],
Avg_Cogongrass_Cover = cogon.pred.vals,
psi.low = psi.comm.quants[1, ],
psi.med = psi.comm.quants[2, ],
psi.high = psi.comm.quants[3, ]
)
}
# Combine all scenarios
psi.comm.dat <- do.call(rbind, community_predictions)
# Reorder the scenario factor levels
psi.comm.dat$scenario <- factor(psi.comm.dat$scenario,
levels = c("Low", "Average", "High"))
community_plot <- ggplot(psi.comm.dat,
aes(x = Avg_Cogongrass_Cover, y = psi.med,
color = scenario, fill = scenario)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
geom_line(size = 1.5) +
scale_color_manual(values = canopy_scenarios$color) +
scale_fill_manual(values = canopy_scenarios$color) +
theme_bw(base_size = 14) +
labs(x = "Average Cogongrass Cover (%)",
y = "Community Occupancy Probability",
color = "Canopy Cover",
fill = "Canopy Cover") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(limits = c(0, 100)) +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5))
community_plot
ggsave("community_cogongrass_canopy.png",
community_plot,
width = 10, height = 7, dpi = 300)
write.csv(psi.comm.dat, "community_predictions_by_canopy.csv", row.names = FALSE)