# If necessary install packages using the three lines of code below, these packages will be used throughout the Occupancy analysis.
#when saving items as a csv or needing to install packages get rid of the # symbol as this prevents packages being reinstalled every time this code is attempted to run.
#install.packages("lubridate", lib = "C:/path_to_your_directory")
#install.packages("readr", repos = "https://cran.r-project.org")
#install.packages(c("readr", "lubridate", "dplyr", "knitr", "kableExtra", "raster", "nngeo", "elevatr", "sf", "tidyr", "ggplot2"))
library(readr)
library(lubridate)
library(dplyr)
library(knitr)
library(kableExtra)
deployment <- read_csv("deployments.csv",
col_types = cols(start_date = col_character()))
#The format of the data set doesn't conform with what is required to run the anlysis so converting it to a POSIXct format to ensure a standardised value for time.
# Convert start_date from DD/MM/YYYY HH:MM to POSIXct
deployment$start_date <- dmy_hm(deployment$start_date) # telling R it is currently in the day, month, year - hours, minutes format.
deployment$end_date <- dmy_hm(deployment$end_date)
# Format start_date to the desired format YYYY-MM-DD HH:MM:SS
deployment$start_date <- format(deployment$start_date, "%Y-%m-%d %H:%M:%S")
deployment$end_date <- format(deployment$end_date, "%Y-%m-%d %H:%M:%S")
camtraps <- read.csv("images.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
# Convert start_date from DD/MM/YYYY HH:MM to POSIXct
camtraps$timestamp <- dmy_hm(camtraps$timestamp) # telling R it is currently in the day, month, year - hours, minutes format.
# change the time stamp so it is compatible with the deployments data set.
camtraps$timestamp <- format(camtraps$timestamp, "%Y-%m-%d %H:%M:%S")1 Step 1 Upload the data and load necessary packages:
2 Step 2 Add necesarry variables:
3 Step 2:1 add the variable elevation for each deployment
# due to rugged landscape elevation may be a variable that alters the probability of detection. the elevatr package which is utilized below uses meters as the default unit of measurement.
# install and load the package "elevatr"
#install.packages("elevatr")
library(elevatr)
# for this analysis to work the dataset must have longitude and latitude which in this case it does!
coordinates <- data.frame(x = deployment$longitude, y = deployment$latitude)
# This line retrieves the elevation data, the code src = "aws" retrieves the data which requires an active internet connection.
elevation_data <- get_elev_point(coordinates, prj = "EPSG:4326", src = "aws")
# Add the elevation (in meters) as a new column to the deployment dataset
deployment$elevation <- elevation_data$elevation4 Step 2:2 add the distance to rivers/streams, water bodies, roads, urban areas and the NDVI for each deployment:
#for this step to work the coordinate system used needs to be changed from ESPG:4326 which is the standard for gps units however it measures distance in degrees rather than metres which is no use for when it comes to analyzing the data.
#so we first need to convert the coordinate system used which will aid in the analysis of the
#install.packages(c("raster", "terra", "sf", "dplyr", "nngeo"))
#package that handles raster data
library(raster)
#load the necessary packages for calculating distances and handling spatial data.
library(nngeo)
#install.packages("sf")
library(sf)
# Load the NDVI Tiff as a raster
ndvi_raster <- raster("NDVI.tif")
# Convert deployment data to an sf object
deployments_sf <- st_as_sf(deployment, coords = c("longitude", "latitude"), crs = 4326)
ndvi_values <- raster::extract(ndvi_raster, deployments_sf)
print(head(ndvi_values))
# Add NDVI values to the SpatialPointsDataFrame
deployments_sf$NDVI <- ndvi_values
# as the maximum value for the NDVI before it has been scaled properly is 255 dividing it by that will provide us with the actual values of NDVI which has to fall between 0 - 1.
deployments_sf$NDVI <- deployments_sf$NDVI / 255
# using the below code we can inspect the data to make sure it actually worked properly
View(deployments_sf)
# Replace NA values with 0 in the NDVI column
deployments_sf$NDVI[is.na(deployments_sf$NDVI)] <- 0
# re project the data to the coordinate system used most commonly in California which is ESPG:32610
deployments_utm <- st_transform(deployments_sf, crs = 32610)
# load the shapefiles into R
roads <- st_read("roads.shp")
rivers <- st_read("rivers_streams.shp")
urban <- st_read("urban_areas.shp")
waterbodies <- st_read("large_waterways.shp")
# Reproject these layers to ensure they are compatible with the deployment data.
roads <- st_transform(roads, crs = 32610)
rivers <- st_transform(rivers, crs = 32610)
urban <- st_transform(urban, crs = 32610)
waterbodies <- st_transform(waterbodies, crs = 32610)
#after reprojection we can now calculate the distance of each deployment to the nearest feature for each of the shapefiles previously loaded in.
# Find the nearest feature st_nn calculates what the nearest neighbor is. for each deployment k = 1 specifies that it is only the nearest feature we want, this prevents R from calculating the proximity of every feature and just the nearest.
#this step may take a little long
#nearest river
nearest_rivers <- st_nn(deployments_utm, rivers, k = 1, returnDist = TRUE)
deployments_sf$distance_to_rivers <- unlist(nearest_rivers$dist)
# Nearest road
nearest_roads <- st_nn(deployments_utm, roads, k = 1, returnDist = TRUE)
deployments_sf$distance_to_roads <- unlist(nearest_roads$dist)
# Nearest urban area
nearest_urban <- st_nn(deployments_utm, urban, k = 1, returnDist = TRUE)
deployments_sf$distance_to_urban <- unlist(nearest_urban$dist)
# Nearest waterbody
nearest_waterbody <- st_nn(deployments_utm, waterbodies, k = 1, returnDist = TRUE)
deployments_sf$distance_to_waterbody <- unlist(nearest_waterbody$dist)Although filling the NDVI column NA’s with 0’s will likely alter the results slightly the models where not working to their best efficiency and getting rid of over 52 deployments because there where NA’s. in addition to this the models output showed that there was a large confidence interval which is why i came back and altered the this code
#After adding generating all of the variables that require the deployments to be a shapefile we need to convert it back to a CSV so we can save it as a new dataset.
#if you inspect the dataset you will see theres a column called geometry this is no longer required so we would want to get rid of this column to simplify the dataframe.
#this line tells R to get rid of the geometry column converting it to a dataframe and leaving the variables that where calculated in the dataframe.
deployments_df <- deployments_sf %>%
mutate(longitude = st_coordinates(.)[, 1],
latitude = st_coordinates(.)[, 2]) %>%
st_drop_geometry()
#this saves the dataframe as a csv
write.csv(deployments_df, "deployment_with_distances.csv", row.names = FALSE)##step 2:3 integrate the RIA into the dataframe
Integrate the RIA data previously created to add an extra variable. For this occupancy analysis i will be using the file previously produced called trapping_events_with_threshold.csv. the reason the RIA with the threshold will be used it to prevent multiple captures of the same individual altering the results.
# Load necessary libraries
library(dplyr)
library(tidyr)
# Read in the data for ria
#deployments_df <- read.csv("deployment_with_distances.csv")
trapping_events_df <- read.csv("trapping_events_with_threshold.csv")
# Pivot the trapping events data to wide format with species as columns
ria_wide <- trapping_events_df %>%
dplyr::select(deployment_id, full_species, RIA) %>%
tidyr::pivot_wider(names_from = full_species, values_from = RIA, values_fill = 0)
# Filter only the desired species
ria_filtered <- ria_wide %>%
dplyr::select(deployment_id, `Tamiasciurus douglasii`, `Ursus americanus`, `Odocoileus hemionus`)
# Join the filtered RIA data with the deployments DataFrame
final_df <- deployments_df %>%
left_join(ria_filtered, by = "deployment_id")
# Write the resulting DataFrame to a new CSV file
write.csv(final_df, "deployment+variables.csv", row.names = FALSE)5 Step 3 create a dataframe with just covariates
# Create a new dataset with the columns that will be used for the occupancy analysis
deployment_subset <- final_df[, c("deployment_id","longitude","latitude","start_date","end_date","elevation","distance_to_roads","distance_to_rivers","distance_to_urban","distance_to_waterbody","Tamiasciurus douglasii","Ursus americanus","Odocoileus hemionus","NDVI")]
# Get the first 10 rows to ensure the filtering worked
first_10_rows <- head(deployment_subset, 10)
# using the kable function to display the result in a table
kable(first_10_rows, format = "html", caption = "Deployment Data with 3 covariates") %>%
kable_styling(full_width = FALSE, position = "left")
# save the new dataframe as a csv h
write.csv(deployment_subset,
"deployment+covariates.csv",
row.names = FALSE)
deployment_data <- read.csv("deployment+covariates.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)6 Step 4 inputting the effort as a covariate:
library(dplyr)
library(lubridate)
# convert the start and end date from POSIXct format to date format
deployment_data <- deployment_data %>%
mutate(start_date = as.Date(start_date, format = "%Y-%m-%d"),
end_date = as.Date(end_date, format = "%Y-%m-%d"))
#convert camtraps to a POSIXct format for later use
camtraps <- camtraps %>%
mutate(timestamp = as.POSIXct(timestamp, format = "%Y-%m-%d %H:%M:%S"))
# Get the earliest start date from the deployments dataset
overall_start_date <- as.Date(min(deployment_data$start_date))
# Get the latest end date from the deployments dataset
overall_end_date <- as.Date(max(deployment_data$end_date))
# Define a period length of 10 days
period_length <- 10
# Calculate the total number of periods required to cover the entire range
total_periods <- ceiling(as.numeric(difftime(overall_end_date, overall_start_date, units = "days")) / period_length)
# create unique column names for each column called SO1, SO2, ect
so_columns <- paste0("SO", 1:total_periods)
# create a dataframe to fit the calculated data sets above
effort_df <- data.frame(deployment_id = character(),
matrix(NA, nrow = 0, ncol = total_periods),
start_date = as.Date(character()),
end_date = as.Date(character()),
stringsAsFactors = FALSE)
# Set the column names for the periods
colnames(effort_df)[2:(total_periods+1)] <- so_columns
# tell r to go through every deployment to calculate the effort which for each period which was previously defined
for (index in 1:nrow(deployment_data)) {
deployment_id <- deployment_data$deployment_id[index]
start_date <- deployment_data$start_date[index]
end_date <- deployment_data$end_date[index]
total_days <- as.numeric(difftime(end_date, start_date, units = "days")) + 1
# Initialize an array to hold the days for each period
period_days <- rep(NA, total_periods) # Adjust the number of periods dynamically
# god through each period and fill it out
for (period in 0:(total_periods-1)) {
period_start <- overall_start_date + days(period * period_length)
period_end <- period_start + days(period_length - 1)
# Check if the deployment overlaps with the current period to ensure r has calculated the period correctly
if (period_end >= start_date && period_start <= end_date) {
active_days <- max(0, as.numeric(difftime(min(end_date, period_end), max(start_date, period_start), units = "days")) + 1)
period_days[period + 1] <- as.integer(active_days)
}
}
# Combine results into the effort_df
effort_df <- rbind(effort_df, data.frame(deployment_id, t(period_days), start_date, end_date, stringsAsFactors = FALSE))
}
# Rename columns to match the SO format previously made
colnames(effort_df)[2:(total_periods+1)] <- so_columns
# remove any duplicate columns that may have been produced
effort_df <- effort_df %>%
dplyr::select(deployment_id, all_of(so_columns), start_date, end_date)
# Organize the dataset by the earliest date to the latest
effort_df <- effort_df %>%
arrange(start_date)
# Save the effort_df as a CSV file
write.csv(effort_df,
"effort_data.csv",
row.names = FALSE)
# Get the first 5 rows
first_5_rows <- head(effort_df, 10)
# Display the first five rows in a table to sample the data
kable(first_5_rows, format = "html", caption = "Effort data over 105 days, split by 10 days ocassions") %>%
kable_styling(full_width = FALSE, position = "left")7 Step 5.1 creating a detection matrix:
#load the necessary packages
library(lubridate)
library(dplyr)
#convert the deployment data from a date format back to the POSIXct format
deployment_data <- deployment_data %>%
mutate(start_date = as.Date(start_date, format = "%Y-%m-%d"),
end_date = as.Date(end_date, format = "%Y-%m-%d"))
camtraps <- camtraps %>%
mutate(timestamp = as.POSIXct(timestamp, format = "%Y-%m-%d %H:%M:%S"))
#get the start and end data similarly to when calculating the effort
overall_start_date <- as.Date(min(deployment_data$start_date))
overall_end_date <- as.Date(max(deployment_data$end_date))
#define the period length of 10 days
period_length <- 10
#calculated the number of 5 day periods that would fall within the start and end date
total_periods <- ceiling(as.numeric(difftime(overall_end_date, overall_start_date, units = "days")) / period_length)
#like with the effort data create unique title for each column
so_columns <- paste0("SO", 1:total_periods)
#select the desired columns from the deployment data creating a vector named squirrel detection. repeat this step for each of the desired species. Black bear and Mule deer.
squirrel_detection <- deployment_data %>%
dplyr::select(deployment_id, longitude, latitude, elevation, distance_to_roads, distance_to_rivers, distance_to_urban, distance_to_waterbody, Tamiasciurus.douglasii, Ursus.americanus, Odocoileus.hemionus, NDVI)
bear_detection <- deployment_data %>%
dplyr::select(deployment_id, longitude, latitude, elevation, distance_to_roads, distance_to_rivers, distance_to_urban, distance_to_waterbody, Tamiasciurus.douglasii, Ursus.americanus, Odocoileus.hemionus, NDVI)
deer_detection <- deployment_data %>%
dplyr::select(deployment_id, longitude, latitude, elevation, distance_to_roads, distance_to_rivers, distance_to_urban, distance_to_waterbody, Tamiasciurus.douglasii, Ursus.americanus, Odocoileus.hemionus, NDVI)8 Step 5.2 create the detection matrix for Douglas’s squirrel.
# Prepare the SO columns for the detection data and fill them with NAs.
for (col in so_columns) {
squirrel_detection[[col]] <- NA
}
#similarly to the effort data we now need to tell r to go through each column and row and see wether or not a squirrel was detected.
for (i in 1:nrow(squirrel_detection)) {
current_deployment_id <- squirrel_detection$deployment_id[i] # Current deployment id
# Get start and end date for the current deployment
start_date <- deployment_data$start_date[deployment_data$deployment_id == current_deployment_id][1]
end_date <- deployment_data$end_date[deployment_data$deployment_id == current_deployment_id][1]
# go through each period from S01 to the last column to generate the detection matrix
for (period in 1:total_periods) {
# Calculate the period for the time between the start and end date
period_start <- overall_start_date + days((period - 1) * period_length)
period_end <- period_start + days(period_length - 1)
# Check if the deployment overlaps with the current period that has been set to ensure it was created correctly
if (max(as.Date(period_start), as.Date(start_date)) <= min(as.Date(period_end), as.Date(end_date))) {
# go through the data sets looking specifically for the current deployment and species (Douglas squirrel)
squirrel_count <- camtraps %>%
filter(deployment_id == current_deployment_id,
common_name == "Douglas's Squirrel",
as.POSIXct(timestamp) >= as.POSIXct(period_start) &
as.POSIXct(timestamp) <= as.POSIXct(period_end))
# Count the occurrences that douglass squirrel was caputured
squirrel_count_num <- nrow(squirrel_count)
# Assign the count to detection dataframe
squirrel_detection[[paste0("SO", period)]][i] <- ifelse(squirrel_count_num > 0, 1, 0)
} else {
squirrel_detection[[paste0("SO", period)]][i] <- NA
}
}
}write.csv(squirrel_detection,
"squirrel_detection.csv",
row.names = FALSE)
# Get the first 10 rows
first_10_rows_squirrel <- head(squirrel_detection, 10)
# Create a nice table for visualization
kable(first_10_rows_squirrel, format = "html", caption = "Detection matrix for squirrel over 105 days, split by 10-day occasions") %>%
kable_styling(full_width = FALSE, position = "left")9 Step 5.3 create the detection matrix for the Black bear
# Prepare the SO columns for the detection data and fill them with NAs.
for (col in so_columns) {
bear_detection[[col]] <- NA
}
#similarly to the effort data we now need to tell r to go through each column and row and see whether or not a Bear was detected.
for (i in 1:nrow(bear_detection)) {
current_deployment_id <- bear_detection$deployment_id[i] # Current deployment id
# Get start and end date for the current deployment
start_date <- deployment_data$start_date[deployment_data$deployment_id == current_deployment_id][1]
end_date <- deployment_data$end_date[deployment_data$deployment_id == current_deployment_id][1]
# go through each period from S01 to the last column to generate the detection matrix
for (period in 1:total_periods) {
# Calculate the period for the time between the start and end date
period_start <- overall_start_date + days((period - 1) * period_length)
period_end <- period_start + days(period_length - 1)
# Check if the deployment overlaps with the current period that has been set to ensure it was created correctly
if (max(as.Date(period_start), as.Date(start_date)) <= min(as.Date(period_end), as.Date(end_date))) {
# go through the data sets looking specifically for the current deployment and species (Black Bear)
bear_count <- camtraps %>%
filter(deployment_id == current_deployment_id,
common_name == "American Black Bear",
as.POSIXct(timestamp) >= as.POSIXct(period_start) &
as.POSIXct(timestamp) <= as.POSIXct(period_end))
# Count the occurrences that Black Bear was caputured
bear_count_num <- nrow(bear_count)
# Assign the count to detection dataframe
bear_detection[[paste0("SO", period)]][i] <- ifelse(bear_count_num > 0, 1, 0)
} else {
bear_detection[[paste0("SO", period)]][i] <- NA
}
}
}write.csv(bear_detection,
"bear_detection.csv",
row.names = FALSE)
# Get the first 10 rows
first_10_rows_bear <- head(bear_detection, 10)
# Create a nice table for visualization
kable(first_10_rows_bear, format = "html", caption = "Detection matrix for Black Bear over 105 days, split by 10-day occasions") %>%
kable_styling(full_width = FALSE, position = "left")10 Step 5.4 create the detection matrix the Mule deer
# Prepare the SO columns for the detection data and fill them with NAs.
for (col in so_columns) {
deer_detection[[col]] <- NA
}
#similarly to the effort data we now need to tell r to go through each column and row and see whether or not a Deer was detected.
for (i in 1:nrow(deer_detection)) {
current_deployment_id <- deer_detection$deployment_id[i] # Current deployment id
# Get start and end date for the current deployment
start_date <- deployment_data$start_date[deployment_data$deployment_id == current_deployment_id][1]
end_date <- deployment_data$end_date[deployment_data$deployment_id == current_deployment_id][1]
# go through each period from S01 to the last column to generate the detection matrix
for (period in 1:total_periods) {
# Calculate the period for the time between the start and end date
period_start <- overall_start_date + days((period - 1) * period_length)
period_end <- period_start + days(period_length - 1)
# Check if the deployment overlaps with the current period that has been set to ensure it was created correctly
if (max(as.Date(period_start), as.Date(start_date)) <= min(as.Date(period_end), as.Date(end_date))) {
# go through the data sets looking specifically for the current deployment and species (Mule Deer)
deer_count <- camtraps %>%
filter(deployment_id == current_deployment_id,
common_name == "Mule Deer",
as.POSIXct(timestamp) >= as.POSIXct(period_start) &
as.POSIXct(timestamp) <= as.POSIXct(period_end))
# Count the occurrences that Mule deer was caputured
deer_count_num <- nrow(deer_count)
# Assign the count to detection dataframe
deer_detection[[paste0("SO", period)]][i] <- ifelse(deer_count_num > 0, 1, 0)
} else {
deer_detection[[paste0("SO", period)]][i] <- NA
}
}
}write.csv(deer_detection,
"deer_detection.csv",
row.names = FALSE)
# Get the first 10 rows
first_10_rows_deer <- head(deer_detection, 10)
# Create a nice table for visualization
kable(first_10_rows_deer, format = "html", caption = "Detection matrix for Mule Deer over 105 days, split by 10-day occasions") %>%
kable_styling(full_width = FALSE, position = "left")11 Step 6 conducting the occupancy analysis
12 Step 6.1 download necesarry packages
# create a list of the required packages that will be used for the occupancy analysis.
list.of.packages <- c(
"unmarked",
"tidyr",
"dplyr",
"vegan",
"AICcmodavg",
"readr",
"ggplot2",
"gridExtra",
"kableExtra",
"knitr")
# check if they are in the library if not install them and load them up.
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages,repos = "http://cran.us.r-project.org")
lapply(list.of.packages, require, character.only = TRUE)13 Step 6.2.1 read in the detection history and covariates for the Squirrel.
to do this for each of the three species after every new vector created i will be adding an s at the end to represent squirrel, b at the end to represent the bears and d at the end to represent deers to ensure nothing is being overwritten when running different parts of the code.
# load in the squirrel detection data and the effort dataframe
squirrel_detection <- read.csv("squirrel_detection.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
effort_df <- read.csv("effort_data.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
# Match the effort dataframe to squirrel the squirrel detection data based on the deployment_id, to make sure detections, covariates and effort to ensure they align and match.
effort_aligneds <- effort_df[match(squirrel_detection$deployment_id, effort_df$deployment_id), ]
# Create data frame with detection histories columns 13 -> 93 (SO1-SO82); remember that we only make use of 32 sessions, 8 per season
y.crosss <- squirrel_detection[,13:92]
# Create data frame with site covariates, which are found in columns 4 to 12
siteCovss <- squirrel_detection[, 4:12]
# Standardize covariates
siteCovss$elevation <- scale(siteCovss$elevation)
siteCovss$distance_to_roads <- scale(siteCovss$distance_to_roads)
siteCovss$distance_to_rivers <- scale(siteCovss$distance_to_rivers)
siteCovss$distance_to_urban <- scale(siteCovss$distance_to_urban)
siteCovss$distance_to_waterbody <- scale(siteCovss$distance_to_waterbody)
# Standardize RIA data as this will be used in the analysis
siteCovss$Tamiasciurus.douglasii <- scale(siteCovss$Tamiasciurus.douglasii)
siteCovss$Ursus.americanus <- scale(siteCovss$Ursus.americanus)
siteCovss$Odocoileus.hemionus <- scale(siteCovss$Odocoileus.hemionus)
# Create a list for yearly site covariates named summer, autumn, winter, and spring
#seasons <- list(season = rep(c("summer", "autumn", "winter", "spring")))
#seasons <- matrix(seasons, nrow(squirrel_detection),4,byrow=TRUE)
# Define seasons as a matrix with each column representing a season
season_matrixs <- matrix(c("summer", "autumn", "winter", "spring"),
nrow = nrow(y.crosss), # as many rows as deployments
ncol = 4, # same number than seasons
byrow = TRUE)
# Convert to a list for yearlySiteCovs
seasons_lists <- list(season = season_matrixs)
# Read in & create data frame with effort (columns 2-81) ; covering the 80 sessions
efforts <- effort_aligneds[, 2:81]
# Check dimensions (should match)
dim(y.crosss)[1] 300 80dim(efforts)[1] 300 80dim(season_matrixs)[1] 300 414 Step 6.2.2 read in the detection history and covariates for the Bear.
# load in the bear detection data and the effort dataframe
bear_detection <- read.csv("bear_detection.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
effort_df <- read.csv("effort_data.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
# As the effort is the same for every species as there was not unique deployments for each species we can reuse the squirrels effort dataset for the bears
# Match the effort dataframe to squirrel the squirrel detection data based on the deployment_id, to make sure detections, covariates and effort to ensure they align and match.
effort_alignedb <- effort_df[match(bear_detection$deployment_id, effort_df$deployment_id), ]
# Create data frame with detection histories columns 13 -> 92 (SO1-SO82); remember that we only make use of 82 sessions, 20 per season
y.crossb <- bear_detection[,13:92]
# Create data frame with site covariates, which are found in columns 4 to 12
siteCovsb <- bear_detection[, 4:12]
# Standardize covariates
siteCovsb$elevation <- scale(siteCovsb$elevation)
siteCovsb$distance_to_roads <- scale(siteCovsb$distance_to_roads)
siteCovsb$distance_to_rivers <- scale(siteCovsb$distance_to_rivers)
siteCovsb$distance_to_urban <- scale(siteCovsb$distance_to_urban)
siteCovsb$distance_to_waterbody <- scale(siteCovsb$distance_to_waterbody)
# Standardize RIA data as this will be used in the analysis
siteCovsb$Tamiasciurus.douglasii <- scale(siteCovsb$Tamiasciurus.douglasii)
siteCovsb$Ursus.americanus <- scale(siteCovsb$Ursus.americanus)
siteCovsb$Odocoileus.hemionus <- scale(siteCovsb$Odocoileus.hemionus)
# Create a list for yearly site covariates named summer, autumn, winter, and spring
#seasons <- list(season = rep(c("summer", "autumn", "winter", "spring")))
#seasons <- matrix(seasons, nrow(squirrel_detection),4,byrow=TRUE)
# Define seasons as a matrix with each column representing a season
season_matrixb <- matrix(c("summer", "autumn", "winter", "spring"),
nrow = nrow(y.crossb), # as many rows as deployments
ncol = 4, # same number than seasons
byrow = TRUE)
# Convert to a list for yearlySiteCovs
seasons_listb <- list(season = season_matrixb)
# Read in & create data frame with effort (columns 2-81) ; covering the 80 sessions. Using the same effort data frame as squirrels.
effortb <- effort_alignedb[, 2:81]
# Check dimensions (should match)
dim(y.crossb)[1] 300 80dim(effortb)[1] 300 80dim(season_matrixb)[1] 300 415 Step 6.2.3 read in the detection history and covariates for the Deer.
# load in the deer detection data and the effort dataframe
deer_detection <- read.csv("deer_detection.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
effort_df <- read.csv("effort_data.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
# Match the effort dataframe to the deer detection data based on the deployment_id, to make sure detections, covariates and effort to ensure they align and match.
effort_alignedd <- effort_df[match(deer_detection$deployment_id, effort_df$deployment_id), ]
# Create data frame with detection histories columns 13 -> 92 (SO1-SO82); remember that we only make use of 82 sessions, 20 per season
y.crossd <- deer_detection[,13:92]
# Create data frame with site covariates, which are found in columns 4 to 12
siteCovsd <- deer_detection[, 4:12]
# Standardize covariates
siteCovsd$elevation <- scale(siteCovsd$elevation)
siteCovsd$distance_to_roads <- scale(siteCovsd$distance_to_roads)
siteCovsd$distance_to_rivers <- scale(siteCovsd$distance_to_rivers)
siteCovsd$distance_to_urban <- scale(siteCovsd$distance_to_urban)
siteCovsd$distance_to_waterbody <- scale(siteCovsd$distance_to_waterbody)
# Standardize RIA data as this will be used in the analysis
siteCovsd$Tamiasciurus.douglasii <- scale(siteCovsd$Tamiasciurus.douglasii)
siteCovsd$Ursus.americanus <- scale(siteCovsd$Ursus.americanus)
siteCovsd$Odocoileus.hemionus <- scale(siteCovsd$Odocoileus.hemionus)
# Create a list for yearly site covariates named summer, autumn, winter, and spring
#seasons <- list(season = rep(c("summer", "autumn", "winter", "spring")))
#seasons <- matrix(seasons, nrow(squirrel_detection),4,byrow=TRUE)
# Define seasons as a matrix with each column representing a season
season_matrixd <- matrix(c("summer", "autumn", "winter", "spring"),
nrow = nrow(y.crossd), # as many rows as deployments
ncol = 4, # same number than seasons
byrow = TRUE)
# Convert to a list for yearlySiteCovs
seasons_listd <- list(season = season_matrixd)
# Read in & create data frame with effort (columns 2-83) ; covering the 80 sessions. Using the same effort data frame as squirrels.
effortd <- effort_alignedd[, 2:81]
# Check dimensions (should match)
dim(y.crossd)[1] 300 80dim(effortd)[1] 300 80dim(season_matrixd)[1] 300 416 Step 6.3.1 create an unmarked dataframe which contains all of the dataframes for the squirrels
#install.packages("unmarked")
library(unmarked)Warning: package 'unmarked' was built under R version 4.4.2# create single unmarked data frame with detection histories, site covs and sample covs #
# Create the unmarkedMultFrame object for occupancy modeling
squirrel <- unmarkedMultFrame(
y = y.crosss, # Dataframe for the detection matrix
siteCovs = siteCovss,
yearlySiteCovs = seasons_lists,
obsCovs = list(efforts = efforts),
numPrimary = 4 # Number of primary survey periods (4 seasons)
)Warning: yearlySiteCovs contains characters. Converting them to factors.# Get a breakdown of the number of sites, detections, non-detections, and information on each covariate
summary(squirrel)unmarkedFrame Object
300 sites
Maximum number of observations per site: 80
Mean number of observations per site: 4.45
Number of primary survey periods: 4
Number of secondary survey periods: 20
Sites with at least one detection: 213
Tabulation of y observations:
0 1 <NA>
633 703 22664
Site-level covariates:
elevation.V1 distance_to_roads.V1 distance_to_rivers.V1
Min. :-1.5721628 Min. :-1.2017612 Min. :-1.7262304
1st Qu.:-0.8287921 1st Qu.:-0.9040452 1st Qu.:-0.7611969
Median :-0.1010265 Median :-0.2487094 Median :-0.1029578
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.0000000
3rd Qu.: 0.7430680 3rd Qu.: 0.6231422 3rd Qu.: 0.5935525
Max. : 1.9262191 Max. : 2.5650240 Max. : 2.8381711
distance_to_urban.V1 distance_to_waterbody.V1 Tamiasciurus.douglasii.V1
Min. :-2.0993534 Min. :-2.0665421 Min. :-0.255594
1st Qu.:-0.6239979 1st Qu.:-0.7740383 1st Qu.:-0.253983
Median : 0.0150121 Median :-0.0485365 Median :-0.225527
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
3rd Qu.: 0.8093207 3rd Qu.: 0.8341673 3rd Qu.:-0.095593
Max. : 2.1511238 Max. : 2.2903886 Max. :14.548294
Ursus.americanus.V1 Odocoileus.hemionus.V1 NDVI
Min. :-0.441652 Min. :-0.498041 Min. :0.0000
1st Qu.:-0.441652 1st Qu.:-0.483937 1st Qu.:0.4471
Median :-0.402400 Median :-0.321742 Median :0.6627
Mean : 0.000000 Mean : 0.000000 Mean :0.5744
3rd Qu.:-0.063850 3rd Qu.: 0.061415 3rd Qu.:0.8235
Max. : 6.201775 Max. : 9.299492 Max. :0.9804
Observation-level covariates:
efforts
Min. : 1.00
1st Qu.: 8.00
Median :10.00
Mean : 8.59
3rd Qu.:10.00
Max. :10.00
NA's :22664
Yearly-site-level covariates:
season
autumn:300
spring:300
summer:300
winter:300 17 Step 6.3.2 Create an unmarked dataframe which contains all of the dataframes for the Bears
# create single unmarked data frame with detection histories, site covs and sample covs
# Create the unmarkedMultFrame object for occupancy modeling
bear <- unmarkedMultFrame(
y = y.crossb, # Dataframe for the detection matrix
siteCovs = siteCovsb,
yearlySiteCovs = seasons_listb,
obsCovs = list(effortb = effortb),
numPrimary = 4 # Number of primary survey periods (4 seasons)
)Warning: yearlySiteCovs contains characters. Converting them to factors.# Get a breakdown of the number of sites, detections, non-detections, and information on each covariate
summary(bear)unmarkedFrame Object
300 sites
Maximum number of observations per site: 80
Mean number of observations per site: 4.45
Number of primary survey periods: 4
Number of secondary survey periods: 20
Sites with at least one detection: 161
Tabulation of y observations:
0 1 <NA>
1000 336 22664
Site-level covariates:
elevation.V1 distance_to_roads.V1 distance_to_rivers.V1
Min. :-1.5721628 Min. :-1.2017612 Min. :-1.7262304
1st Qu.:-0.8287921 1st Qu.:-0.9040452 1st Qu.:-0.7611969
Median :-0.1010265 Median :-0.2487094 Median :-0.1029578
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.0000000
3rd Qu.: 0.7430680 3rd Qu.: 0.6231422 3rd Qu.: 0.5935525
Max. : 1.9262191 Max. : 2.5650240 Max. : 2.8381711
distance_to_urban.V1 distance_to_waterbody.V1 Tamiasciurus.douglasii.V1
Min. :-2.0993534 Min. :-2.0665421 Min. :-0.255594
1st Qu.:-0.6239979 1st Qu.:-0.7740383 1st Qu.:-0.253983
Median : 0.0150121 Median :-0.0485365 Median :-0.225527
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
3rd Qu.: 0.8093207 3rd Qu.: 0.8341673 3rd Qu.:-0.095593
Max. : 2.1511238 Max. : 2.2903886 Max. :14.548294
Ursus.americanus.V1 Odocoileus.hemionus.V1 NDVI
Min. :-0.441652 Min. :-0.498041 Min. :0.0000
1st Qu.:-0.441652 1st Qu.:-0.483937 1st Qu.:0.4471
Median :-0.402400 Median :-0.321742 Median :0.6627
Mean : 0.000000 Mean : 0.000000 Mean :0.5744
3rd Qu.:-0.063850 3rd Qu.: 0.061415 3rd Qu.:0.8235
Max. : 6.201775 Max. : 9.299492 Max. :0.9804
Observation-level covariates:
effortb
Min. : 1.00
1st Qu.: 8.00
Median :10.00
Mean : 8.59
3rd Qu.:10.00
Max. :10.00
NA's :22664
Yearly-site-level covariates:
season
autumn:300
spring:300
summer:300
winter:300 # Create the unmarkedMultFrame object for occupancy modeling
deer <- unmarkedMultFrame(
y = y.crossd, # Dataframe for the detection matrix
siteCovs = siteCovsd,
yearlySiteCovs = seasons_listd,
obsCovs = list(effortd = effortd),
numPrimary = 4 # Number of primary survey periods (4 seasons)
)Warning: yearlySiteCovs contains characters. Converting them to factors.# Get a breakdown of the number of sites, detections, non-detections, and information on each covariate
summary(deer)unmarkedFrame Object
300 sites
Maximum number of observations per site: 80
Mean number of observations per site: 4.45
Number of primary survey periods: 4
Number of secondary survey periods: 20
Sites with at least one detection: 215
Tabulation of y observations:
0 1 <NA>
794 542 22664
Site-level covariates:
elevation.V1 distance_to_roads.V1 distance_to_rivers.V1
Min. :-1.5721628 Min. :-1.2017612 Min. :-1.7262304
1st Qu.:-0.8287921 1st Qu.:-0.9040452 1st Qu.:-0.7611969
Median :-0.1010265 Median :-0.2487094 Median :-0.1029578
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.0000000
3rd Qu.: 0.7430680 3rd Qu.: 0.6231422 3rd Qu.: 0.5935525
Max. : 1.9262191 Max. : 2.5650240 Max. : 2.8381711
distance_to_urban.V1 distance_to_waterbody.V1 Tamiasciurus.douglasii.V1
Min. :-2.0993534 Min. :-2.0665421 Min. :-0.255594
1st Qu.:-0.6239979 1st Qu.:-0.7740383 1st Qu.:-0.253983
Median : 0.0150121 Median :-0.0485365 Median :-0.225527
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
3rd Qu.: 0.8093207 3rd Qu.: 0.8341673 3rd Qu.:-0.095593
Max. : 2.1511238 Max. : 2.2903886 Max. :14.548294
Ursus.americanus.V1 Odocoileus.hemionus.V1 NDVI
Min. :-0.441652 Min. :-0.498041 Min. :0.0000
1st Qu.:-0.441652 1st Qu.:-0.483937 1st Qu.:0.4471
Median :-0.402400 Median :-0.321742 Median :0.6627
Mean : 0.000000 Mean : 0.000000 Mean :0.5744
3rd Qu.:-0.063850 3rd Qu.: 0.061415 3rd Qu.:0.8235
Max. : 6.201775 Max. : 9.299492 Max. :0.9804
Observation-level covariates:
effortd
Min. : 1.00
1st Qu.: 8.00
Median :10.00
Mean : 8.59
3rd Qu.:10.00
Max. :10.00
NA's :22664
Yearly-site-level covariates:
season
autumn:300
spring:300
summer:300
winter:300 18 Step 6.4 plot the detection
# plot of detection/non-detection from this we can see the main three deployment times and we can see that the squirrel seemed to have been detected the most, then the deer and then the bear which coincides with what the RIA previously informed us.
plot(squirrel)
plot(bear)
plot(deer)19 Step 6.5.1 test for the correlation between the covariates
as the covariates and the values they have are the same for every deployment for all three species there is no point in running multiple correlation tests so for the code below only the siteCovss will be used which was initially created for the Squirrel.
#To test the correlation between all of the variables running this code will produce correlations between each variable.
cor(siteCovss[, c("elevation", "distance_to_roads", "distance_to_rivers", "distance_to_urban", "distance_to_waterbody", "Tamiasciurus.douglasii", "Ursus.americanus", "Odocoileus.hemionus", "NDVI")])
#after running this code it is clear there is no real correlation that exceeds the threshold of 0.7 or -0.7 which means theres no significant correlation between any two covariates.
#However some of the covariates are different in their ecological role and may influence habitat preference so testing the covariates with the higher correlations that both might influnce the data in the same way is useful to see if they need to be taken out of the analysis
# Test correlation between distance to roads and distance to urban areas
cor.test(siteCovss$distance_to_roads, siteCovss$distance_to_urban)
# Test correlation between elevation and distance to roads
cor.test(siteCovss$elevation, siteCovss$distance_to_roads)
# Test correlation between elevation and distance to waterbody
cor.test(siteCovss$elevation, siteCovss$distance_to_waterbody)
# Test correlation between distance to urban and distance to waterbody
cor.test(siteCovss$distance_to_urban, siteCovss$distance_to_waterbody)
# Test correlation between distance to roads and distance to rivers
cor.test(siteCovss$distance_to_roads, siteCovss$distance_to_rivers)20 Step 6.6.1 create a range of models to fit the squirrel data
# Fit the model to the squirrel data the squirrel model named m0
m0 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~1, # detection constant
data = squirrel)
# Extract parameter estimates
occupancy_est <- backTransform(m0, type = "psi") # Occupancy probability
colonization_est <- backTransform(m0, type = "col") # Colonisation probability
extinction_est <- backTransform(m0, type = "ext") # Extinction probability
detection_est <- backTransform(m0, type = "det") # Detection probability
# Print the estimates
cat("Occupancy probability (psi):", occupancy_est@estimate, "\n")Occupancy probability (psi): 0.7818954 cat("Colonization probability (gamma):", colonization_est@estimate, "\n")Colonization probability (gamma): 0.5840802 cat("Extinction probability (epsilon):", extinction_est@estimate, "\n")Extinction probability (epsilon): 0.2326869 cat("Detection probability (p):", detection_est@estimate, "\n")Detection probability (p): 0.7157974 21 Step 6.6.2 create a range of models to fit the bear data
# Fit the model to the deer data the squirrel model named m1
m1 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~1, # detection constant
data = bear)
# Extract parameter estimates
occupancy_est <- backTransform(m1, type = "psi") # Occupancy probability
colonization_est <- backTransform(m1, type = "col") # Colonisation probability
extinction_est <- backTransform(m1, type = "ext") # Extinction probability
detection_est <- backTransform(m1, type = "det") # Detection probability
# Print the estimates
cat("Occupancy probability (psi):", occupancy_est@estimate, "\n")Occupancy probability (psi): 0.6660325 cat("Colonization probability (gamma):", colonization_est@estimate, "\n")Colonization probability (gamma): 0.5771848 cat("Extinction probability (epsilon):", extinction_est@estimate, "\n")Extinction probability (epsilon): 0.3542619 cat("Detection probability (p):", detection_est@estimate, "\n")Detection probability (p): 0.3817947 22 Step 6.6.3 create a range of models to fit the deer data
# Fit the model to the deer data the deer model named m2
m2 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~1, # detection constant
data = deer)
# Extract parameter estimates
occupancy_est <- backTransform(m2, type = "psi") # Occupancy probability
colonization_est <- backTransform(m2, type = "col") # Colonisation probability
extinction_est <- backTransform(m2, type = "ext") # Extinction probability
detection_est <- backTransform(m2, type = "det") # Detection probability
# Print the estimates
cat("Occupancy probability (psi):", occupancy_est@estimate, "\n")Occupancy probability (psi): 0.8254695 cat("Colonization probability (gamma):", colonization_est@estimate, "\n")Colonization probability (gamma): 0.5717693 cat("Extinction probability (epsilon):", extinction_est@estimate, "\n")Extinction probability (epsilon): 0.1900285 cat("Detection probability (p):", detection_est@estimate, "\n")Detection probability (p): 0.5123236 23 Step 6.7.1 fit the models by detection probability for the squirrels
# model with detection dependent on effort
ms1 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~efforts, # probability of detection is affected by effort
data = squirrel)
ms2 <- colext(~1, ~1, ~1, ~elevation, data = squirrel)
ms3 <- colext(~1, ~1, ~1, ~efforts+elevation, data = squirrel)
# compare the three models made above
dlist<-fitList(Null = m0,ms1=ms1, ms2=ms2, ms3=ms3)
selmod<-modSel(dlist,nullmod="Null")
selmod
# as the model produced a lower AIC value than the base model the combonation of elevation and effort makes this model the best fitting one, with an AIC weight of nearly one and an AIC value of 626.2 it fits better than the other models that where simulated.
#keeping these 2 variables in the model for further evaluation will be important seeing which further covariates contribute to the model to wittle out the covariates that do not contribute as much.
ms4 <- colext(~1, ~1, ~1, ~efforts+elevation+distance_to_rivers, data = squirrel)
ms5 <- colext(~1, ~1, ~1, ~efforts+elevation+distance_to_roads, data = squirrel)
ms6 <- colext(~1, ~1, ~1, ~efforts+elevation+Ursus.americanus, data = squirrel)
ms7 <- colext(~1, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii, data = squirrel)
ms8 <- colext(~1, ~1, ~1, ~efforts+elevation+Odocoileus.hemionus, data = squirrel)
ms9 <- colext(~1, ~1, ~1, ~efforts+elevation+distance_to_urban, data = squirrel)
ms10 <- colext(~1, ~1, ~1, ~efforts+elevation+distance_to_waterbody, data = squirrel)
# compare the three models made above
dlist<-fitList(Null = m0,ms1=ms1, ms2=ms2, ms3=ms3, ms4=ms4, ms5=ms5, ms6=ms6, ms7=ms7, ms8=ms8, ms9=ms9, ms10=ms10)
selmod<-modSel(dlist,nullmod="Null")
selmod
# create a model of the top 3 covariates from the previous comparison which where the RIA for the douglas squirrel and black bear and the distance to rivers.
#although distance to roads is relevent to our overall question in regards to forest fires it doesnt support the models being almost near the bottom for AIC.
ms11 <- colext(~1, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers, data = squirrel)
dlist<-fitList ( m0=m0, ms7=ms7, ms6=ms6, ms4=ms4, ms11=ms11)
selmod<-modSel(dlist,nullmod= "m0")
selmod
coef(ms11)
# from this we can see that effort, elevation, the ria for douglas squirrel and black bear and distance to rivers correlate the detection.
# testing for variables that best fit the occupuncy will now be undertaken.
ms1 <- colext(~elevation, # occupancy no longer constant
~1, # colonisation constant
~1, # extinction constant
~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, # probability of detection is affected by effort, tamiasciurus.douglisii, ursus.americanus, distance to rivers, elevation
data = squirrel)
ms2 <- colext(~distance_to_waterbody, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms3 <- colext(~distance_to_waterbody+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms4 <- colext(~distance_to_rivers, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms5 <- colext(~distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms6 <- colext(~distance_to_roads+distance_to_waterbody, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms7 <- colext(~distance_to_roads+distance_to_rivers, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms8 <- colext(~distance_to_roads+distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms9 <- colext(~distance_to_roads+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms10 <- colext(~distance_to_waterbody+distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms11 <- colext(~distance_to_waterbody+distance_to_rivers, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms12 <- colext(~distance_to_rivers+distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms13 <- colext(~distance_to_rivers+distance_to_roads, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms14 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_rivers, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms15 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_roads, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms16 <- colext(~distance_to_waterbody+distance_to_urban+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms17 <- colext(~distance_to_rivers+distance_to_urban+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms18 <- colext(~distance_to_rivers+distance_to_roads+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms19 <- colext(~distance_to_rivers+distance_to_waterbody+elevation, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms20 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms21 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms22 <- colext(~distance_to_urban+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms23 <- colext(~distance_to_urban+distance_to_rivers+elevation+distance_to_roads, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
ms24 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers+elevation, data = squirrel)
dlist<-fitList ( m0=m0, ms1=ms1, ms2=ms2, ms3=ms3, ms4=ms4, ms5=ms5, ms6=ms6, ms7=ms7, ms8=ms8, ms9=ms9, ms10=ms10, ms11=ms11, ms12=ms12, ms13=ms13, ms14=ms14, ms15=ms15, ms16=ms16, ms17=ms17, ms18=ms18, ms19=ms19, ms20=ms20, ms21=ms21, ms22=ms22, ms23=ms23, ms24=ms24)
selmod<-modSel(dlist,nullmod= "m0")
selmod
coef(ms7)
confint(ms7, type="psi")
# even though the ms7 has a large confidence interval meaning and doesn't fall into the 95% it has the highest AIC out of all the other variables meaning they are the best fit for predicting the occupancy. the variables in ms7 are distance to rivers and distance to roads.
ms1 <- colext(~distance_to_rivers+distance_to_roads, # occupancy no longer constant and is represented by the distance to rivers and roads
~season, # colonisation constant
~1, # extinction constant
~efforts+elevation+Ursus.americanus+distance_to_rivers,# probability of detection is affected by effort, tamiasciurus.douglisii, ursus.americanus, distance to rivers, elevation
data = squirrel)
ms2 <- colext(~distance_to_roads+distance_to_rivers, ~Tamiasciurus.douglasii, ~1, ~efforts+elevation+Ursus.americanus+distance_to_rivers, data = squirrel)
ms3 <- colext(~distance_to_roads+distance_to_rivers, ~Tamiasciurus.douglasii+season, ~1, ~efforts+elevation+Ursus.americanus+distance_to_rivers, data = squirrel)
ms4 <- colext(~distance_to_roads+distance_to_rivers, ~Tamiasciurus.douglasii+season, ~season, ~efforts+elevation+Ursus.americanus+distance_to_rivers, data = squirrel)
ms5 <- colext(~distance_to_roads+distance_to_rivers+NDVI, ~Tamiasciurus.douglasii+season, ~season, ~efforts+elevation+Ursus.americanus+distance_to_rivers, data = squirrel)
ms6 <- colext(~distance_to_roads+distance_to_rivers, ~Tamiasciurus.douglasii, ~Tamiasciurus.douglasii, ~efforts+elevation+Ursus.americanus+distance_to_rivers, data = squirrel)
dlist<-fitList ( m0=m0, ms1=ms1, ms2=ms2, ms3=ms3, ms4=ms4, ms5=ms5, ms6=ms6)
selmod<-modSel(dlist,nullmod= "m0")
selmod
coef(ms6)
confint(ms6, type="col")
confint(ms6, type="ext")
#ms6 being the best fit makes sense as the RIA would help predict the colonisation and extinction of the animals in this area therefore because it fits best and makes most ecological sense this will be used in the final model
ms1 <- colext(~distance_to_roads+distance_to_rivers+NDVI, ~Tamiasciurus.douglasii, ~Tamiasciurus.douglasii, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers, data = squirrel)
dlist<-fitList ( m0=m0, ms1=ms1, ms2=ms2, ms3=ms3)
selmod<-modSel(dlist,nullmod= "m0")
selmod
confint(ms1, type="psi") # for occupancy (distance_to_road+distance_to_rivers)
confint(ms1, type="col") # for colonisation (Tamiascurus.douglasii)
confint(ms1, type="ext") # for extinction (Tamiascurus.douglasii)
confint(ms1, type="det") # for detection (efforts+elevation+Tamiascurus.douglasii+Ursus.americanus+distance_to_rivers)ms1 <- colext(~distance_to_roads+distance_to_rivers+NDVI, ~Tamiasciurus.douglasii, ~Tamiasciurus.douglasii, ~efforts+elevation+Tamiasciurus.douglasii+Ursus.americanus+distance_to_rivers, data = squirrel)
confint(ms1, type="psi") # for occupancy (distance_to_road+distance_to_rivers) 0.025 0.975
psi(Int) 0.1996443 3.61986016
psi(distance_to_roads) -0.5357932 0.99159641
psi(distance_to_rivers) -1.2273510 0.02803148
psi(NDVI) -2.7990149 2.06152152confint(ms1, type="col") # for colonisation (Tamiascurus.douglasii) 0.025 0.975
col(Int) -48.21704 254.5764
col(Tamiasciurus.douglasii) -184.88899 1031.7958confint(ms1, type="ext") # for extinction (Tamiascurus.douglasii) 0.025 0.975
ext(Int) -364.9714 18.42883
ext(Tamiasciurus.douglasii) -1440.9300 70.32044confint(ms1, type="det") # for detection (efforts+elevation+Tamiascurus.douglasii+Ursus.americanus+distance_to_rivers) 0.025 0.975
p(Int) -2.878036187 -1.1740297
p(efforts) 0.243619896 0.4420617
p(elevation) 0.254377147 0.8711799
p(Tamiasciurus.douglasii) 1.212099636 3.1653867
p(Ursus.americanus) 0.064691084 1.2030422
p(distance_to_rivers) 0.004795157 0.538724524 Step 6.7.2 fit the models by detection probability for the Bears
# model with detection dependent on effort
mb1 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~effortb, # probability of detection is affected by effort
data = bear)
mb2 <- colext(~1, ~1, ~1, ~elevation, data = bear)
mb3 <- colext(~1, ~1, ~1, ~effortb + elevation, data = bear)
# compare the three models made above
dlist<-fitList(Null = m1,mb1=mb1, mb2=mb2, mb3=mb3)
selmod<-modSel(dlist,nullmod="Null")
selmod
# as the model produced a lower AIC value than the base model the combonation of elevation and effort makes this model the best fitting one, with an AIC weight of nearly one and an AIC value of 626.2 it fits better than the other models that where simulated.
#keeping these 2 variables in the model for further evaluation will be important seeing which further covariates contribute to the model to wittle out the covariates that do not contribute as much.
mb4 <- colext(~1, ~1, ~1, ~effortb+elevation+distance_to_rivers, data = bear)
mb5 <- colext(~1, ~1, ~1, ~effortb+elevation+distance_to_roads, data = bear)
mb6 <- colext(~1, ~1, ~1, ~effortb+elevation+Ursus.americanus, data = bear)
mb7 <- colext(~1, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii, data = bear)
mb8 <- colext(~1, ~1, ~1, ~effortb+elevation+Odocoileus.hemionus, data = bear)
mb9 <- colext(~1, ~1, ~1, ~effortb+elevation+distance_to_urban, data = bear)
mb10 <- colext(~1, ~1, ~1, ~effortb+elevation+distance_to_waterbody, data = bear)
# compare the three models made above
dlist<-fitList(Null = m1,mb1=mb1, mb2=mb2, mb3=mb3, mb4=mb4, mb5=mb5, mb6=mb6, mb7=mb7, mb8=mb8, mb9=mb9, mb10=mb10)
selmod<-modSel(dlist,nullmod="Null")
selmod
# create a model of the top 3 covariates from the previous comparison which where the RIA for the douglas squirrel and black bear and the distance to rivers.
#although distance to roads is relevent to our overall question in regards to forest fires it doesnt support the models being almost near the bottom for AIC.
mb11 <- colext(~1, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban, data = bear)
# the reason i have colated all of these variables as each of them on their own have a higher AIC than model 3 which was the highest out of the origincal 3 models that was produced. we can now compare how all of the variables influence the AIC in comparison to when they are on their own.
mb12 <- colext(~1, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
dlist<-fitList(Null = mb7, mb6=mb6, mb9=mb9, mb8=mb8, mb11=mb11, mb12=mb12)
selmod<-modSel(dlist,nullmod="Null")
selmod
# The top 3 models are mb6, mb11,and mb12 normally the lowest AIC would be the choice however after inspecting the variables which make up these three mb12 is most relevant to the research question which is exploring wildfires effect on the occupancy of the Black bear so although mb6 and mb12 very slightly fits better the presence of the extra variables can provide useful insight into their effect on the detection and occupancy of the bears.
coef(mb12)
mb1 <- colext(~1, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb2 <- colext(~distance_to_waterbody, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb3 <- colext(~distance_to_waterbody+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb4 <- colext(~distance_to_rivers, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb5 <- colext(~distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb6 <- colext(~distance_to_roads+distance_to_waterbody, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb7 <- colext(~distance_to_roads+distance_to_rivers, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb8 <- colext(~distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb9 <- colext(~distance_to_roads+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb10 <- colext(~distance_to_waterbody+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb11 <- colext(~distance_to_waterbody+distance_to_rivers, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb12 <- colext(~distance_to_rivers+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb13 <- colext(~distance_to_rivers+distance_to_roads, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb14 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_rivers, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb15 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_roads, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb16 <- colext(~distance_to_waterbody+distance_to_urban+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb17 <- colext(~distance_to_rivers+distance_to_urban+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb18 <- colext(~distance_to_rivers+distance_to_roads+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb19 <- colext(~distance_to_rivers+distance_to_waterbody+elevation, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb20 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb21 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb22 <- colext(~distance_to_urban+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb23 <- colext(~distance_to_urban+distance_to_rivers+elevation+distance_to_roads, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb24 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb3=mb3, mb4=mb4, mb5=mb5, mb6=mb6, mb7=mb7, mb8=mb8, mb9=mb9, mb10=mb10, mb11=mb11, mb12=mb12, mb13=mb13, mb14=mb14, mb15=mb15, mb16=mb16, mb17=mb17, mb18=mb18, mb19=mb19, mb20=mb20, mb21=mb21, mb22=mb22, mb23=mb23, mb24=mb24)
selmod<-modSel(dlist,nullmod= "m1")
selmod
coef(mb24)
#due to the fact the same covariates are fitting best when present in both occupancy and and detection to avoid multicolinearity plotting the models that fit best depending on the AIC will determine where the variables are will be used either in detection or occupancy.
mb1 <- colext(~distance_to_rivers+distance_to_waterbody+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
mb2 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
#these two models will compare which model best fits the data determined by wether elevation is a predictor of occupuncy or detection.
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb24=mb24)
selmod<-modSel(dlist,nullmod= "m1")
selmod
#although the model fits better when elevation is helping predict occupancy and detection, mb2 which is elevation predicting occupancy has a higher AIC out of the two therefore making it a better fit for predicting occupancy.
#repeating the above steps for distance to rivers and distance to urban as they both appear in predicting occupancy and detection.
mb1 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban, data = bear)
mb2 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb24=mb24)
selmod<-modSel(dlist,nullmod= "m1")
selmod
#distance to rivers fits detection better
mb1 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb2 <- colext(~distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_urban+distance_to_rivers, data = bear)
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb24=mb24)
selmod<-modSel(dlist,nullmod= "m1")
selmod
#distance to urban fits occupancy better
# next is to see what variables best predict extinction and colonisation. only a couple variables really make sense in these columns and they are the RIA of the Black bear and season
#first colonisation.
mb1 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season+Ursus.americanus, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb2 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~Ursus.americanus, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb3 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season, ~1, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb3=mb3)
selmod<-modSel(dlist,nullmod= "m1")
selmod
#adding the season to extinction makes sense especially for bears as during the winter months if bears do not 'bulk' up enough it can result in death for them.
mb4 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season+Ursus.americanus, ~season, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb5 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season+Ursus.americanus, ~Ursus.americanus, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb6 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~Ursus.americanus, ~season+Ursus.americanus, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb7 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season+Ursus.americanus, ~season+Ursus.americanus, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
mb8 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban+NDVI, ~season+Ursus.americanus, ~Ursus.americanus, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
dlist<-fitList ( m1=m1, mb1=mb1, mb2=mb2, mb3=mb3, mb4=mb4, mb5=mb5, mb6=mb6, mb7=mb7, mb8=mb8)
selmod<-modSel(dlist,nullmod= "m1")
selmod
#as i will not be using all of the variables in the final anlysis to simplify the model i will remove the variables that do not fit in the final anlysis
mb1 <- colext(~distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~season+Ursus.americanus, ~Ursus.americanus, ~effortb+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = bear)
confint(mb1, type="psi") # for occupancy (distance_to_waterbody+distance_to_roads+distance_to_urban+elevation)
confint(mb1, type="col") # for colonisation (Ursus.americanus RIA for Black bears)
confint(mb1, type="ext") # for extinction (season+Ursus.americanus)
confint(mb1, type="det") # for detection (efforts+Tamiascurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers)
#simplifying the model is helpful for the predictions to work accurately at the end the confint command helps identify which variables are to weak to keep in the model This can be seen by looking the the large confidence intervals that are shown for example in the detection part of the confint
mb1 <- colext(~elevation+distance_to_roads+distance_to_urban, ~Ursus.americanus, ~Ursus.americanus, ~effortb+Ursus.americanus+Odocoileus.hemionus, data = bear)mb1 <- colext(~elevation+distance_to_roads+distance_to_urban, ~Ursus.americanus, ~Ursus.americanus, ~effortb+Ursus.americanus+Odocoileus.hemionus, data = bear)
confint(mb1, type="psi") # for occupancy (distance_to_roads+distance_to_urban+elevation) 0.025 0.975
psi(Int) 0.2933946 3.9854337
psi(elevation) -4.0169708 -0.9033902
psi(distance_to_roads) -2.2164368 0.5436605
psi(distance_to_urban) -1.9841822 0.4699215confint(mb1, type="col") # for colonisation (Ursus.americanus RIA for Black bears) 0.025 0.975
col(Int) -109.9871 307.3407
col(Ursus.americanus) -265.0873 743.3713confint(mb1, type="ext") # for extinction (season+Ursus.americanus) 0.025 0.975
ext(Int) -309.8429 89.91842
ext(Ursus.americanus) -716.7929 204.51803confint(mb1, type="det") # for detection (efforts+Ursus.americanus+Odocoileus.hemionus) 0.025 0.975
p(Int) -3.13895850 -1.2439602
p(effortb) 0.06367968 0.2672774
p(Ursus.americanus) 0.43880913 0.9979814
p(Odocoileus.hemionus) -0.12956649 0.374989625 Step 6.7.3 fit the models by detection probability for the Deers
# model with detection dependent on effort
md1 <- colext(~1, # occupancy constant
~1, # colonisation constant
~1, # extinction constant
~effortd, # probability of detection is affected by effort
data = deer)
md2 <- colext(~1, ~1, ~1, ~elevation, data = deer)
md3 <- colext(~1, ~1, ~1, ~effortd + elevation, data = deer)
# compare the three models made above
dlist<-fitList(Null = m2,md1=md1, md2=md2, md3=md3)
selmod<-modSel(dlist,nullmod="Null")
selmod
md4 <- colext(~1, ~1, ~1, ~effortd+elevation+distance_to_rivers, data = deer)
md5 <- colext(~1, ~1, ~1, ~effortd+elevation+distance_to_roads, data = deer)
md6 <- colext(~1, ~1, ~1, ~effortd+elevation+Ursus.americanus, data = deer)
md7 <- colext(~1, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii, data = deer)
md8 <- colext(~1, ~1, ~1, ~effortd+elevation+Odocoileus.hemionus, data = deer)
md9 <- colext(~1, ~1, ~1, ~effortd+elevation+distance_to_urban, data = deer)
md10 <- colext(~1, ~1, ~1, ~effortd+elevation+distance_to_waterbody, data = deer)
# compare the three models made above
dlist<-fitList(Null = m2,md1=md1, md2=md2, md3=md3, md4=md4, md5=md5, md6=md6, md7=md7, md8=md8, md9=md9, md10=md10)
selmod<-modSel(dlist,nullmod="Null")
selmod
# create a model of the top covariates from the previous comparison which where the RIA for the douglas squirrel, mule deer and black bear and the distance to rivers, distance_to_waterbody.
md11 <- colext(~1, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers+distance_to_waterbody, data = deer)
# the reason i have colated all of these variables as each of them on their own have a higher AIC than model 3 which was the highest out of the origincal 3 models that was produced. we can now compare how all of the variables influence the AIC in comparison to when they are on their own.
dlist<-fitList(Null = md7, md6=md6, md10=md10, md8=md8, md4=md4, md11=md11)
selmod<-modSel(dlist,nullmod="Null")
selmod
# model 11 which is the culmination of all of the covariates which gave the highest AIC and compared to when they are on their own they offer a higher AIC when they are all together however they look teh culumitave weight of 0.04 which is something to bare in mind when running the test and looking at the results.
# as the coefficients for some of the variables in md11 are rather low and contain 0s this tells us that we can further refine teh test by dropping variables that are similar.
#one variable is waterbodys if we drop the waterbodies variable and test the AIC to see if it had a large impact or not.
md12 <- colext(~1, ~1, ~1, ~effortd + elevation + Tamiasciurus.douglasii + Ursus.americanus + Odocoileus.hemionus + distance_to_rivers, data = deer)
md13 <- colext(~1, ~1, ~1, ~effortd + elevation + Ursus.americanus + Odocoileus.hemionus + distance_to_rivers, data = deer)
md14 <- colext(~1, ~1, ~1, ~effortd + Tamiasciurus.douglasii + Ursus.americanus + Odocoileus.hemionus + distance_to_rivers, data = deer)
dlist_refined <- fitList(Null = md11, md11 = md11, md12 = md12, md13=md13, md14=md14)
selmod_refined <- modSel(dlist_refined, nullmod = "Null")
selmod_refined
# after removing waterbodies we get a slightly better AIC therefore we can remove it from the overall model to avoid multicolinearity which is where multiple variables follow the same trend which can cause issues with the occupuncy model.
coef(md12)
# model 12 is the best fit overall having dropped the distance to waterbodies variable it had a lower AIC and will avoid any mulitcolinearity with distance to rivers however if the rivers prove
md1 <- colext(~elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md2 <- colext(~distance_to_waterbody, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md3 <- colext(~distance_to_waterbody+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md4 <- colext(~distance_to_rivers, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md5 <- colext(~distance_to_urban, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md6 <- colext(~distance_to_roads+distance_to_waterbody, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md7 <- colext(~distance_to_roads+distance_to_rivers, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md8 <- colext(~distance_to_roads+distance_to_urban, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md9 <- colext(~distance_to_roads+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md10 <- colext(~distance_to_waterbody+distance_to_urban, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md11 <- colext(~distance_to_waterbody+distance_to_rivers, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md12 <- colext(~distance_to_rivers+distance_to_urban, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md13 <- colext(~distance_to_rivers+distance_to_roads, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md14 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_rivers, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md15 <- colext(~distance_to_waterbody+distance_to_urban+distance_to_roads, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md16 <- colext(~distance_to_waterbody+distance_to_urban+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md17 <- colext(~distance_to_rivers+distance_to_urban+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md18 <- colext(~distance_to_rivers+distance_to_roads+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md19 <- colext(~distance_to_rivers+distance_to_waterbody+elevation, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md20 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_urban, ~1, ~1, ~effortd + elevation + Tamiasciurus.douglasii + Ursus.americanus + Odocoileus.hemionus + distance_to_rivers, data = deer)
md21 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md22 <- colext(~distance_to_urban+distance_to_waterbody+elevation+distance_to_roads, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md23 <- colext(~distance_to_urban+distance_to_rivers+elevation+distance_to_roads, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md24 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+distance_to_urban, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md25 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~1, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
dlist<-fitList ( m2=m2, md1=md1, md2=md2, md3=md3, md4=md4, md5=md5, md6=md6, md7=md7, md8=md8, md9=md9, md10=md10, md11=md11, md12=md12, md13=md13, md14=md14, md15=md15, md16=md16, md17=md17, md18=md18, md19=md19, md20=md20, md21=md21, md22=md22, md23=md23, md24=md24, md25=md25)
selmod<-modSel(dlist,nullmod= "m2")
selmod
#although model 21 fit best, model 25 contains the NDVI values which are values that represent the presence, density and health of the vegetation using near infrared light and red light from satellite imagery to produce number between 0 and 1. higher numbers are better and show, denser and healthier flora coverage and low numbers show the opposite.
#this is a useful and commonly used tool which is why it will be included in the final model.
#next is to see what variables best fit the role of predicting extinction and colonisation.
md1 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~season, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md2 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md3 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus+season, ~1, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
dlist<-fitList ( m2=m2, md1=md1, md2=md2, md3=md3)
selmod<-modSel(dlist,nullmod= "m2")
selmod
#after running the above code we can see that model 2 ran best and has the highest AIC leaving the RIA for the mule deer to be left to predict colonisation.
md1 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus, ~season, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md2 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus, ~Odocoileus.hemionus, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
md3 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus, ~season+Odocoileus.hemionus, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
dlist<-fitList ( m2=m2, md1=md1, md2=md2, md3=md3)
selmod<-modSel(dlist,nullmod= "m2")
selmod
#again model 2 had the lowest AIC which has the RIA for mule deer account for the mule deers extinction.
md1 <- colext(~distance_to_rivers+distance_to_waterbody+elevation+distance_to_roads+NDVI, ~Odocoileus.hemionus, ~Odocoileus.hemionus, ~effortd+elevation+Tamiasciurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
confint(md1, type="psi") # for occupancy (distance_to_rivers+distance_to_waterbody+distance_to_roads+elevation+NDVI)
confint(md1, type="col") # for colonisation (Odocoileus.hemionus RIA for Mule Deer)
confint(md1, type="ext") # for extinction (Odocoileus.hemionus)
confint(md1, type="det") # for detection (effortd+elevation+Tamiascurus.douglasii+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers)
#simplifying the model is helpful for the predictions to work accuratly at the end the confint command helps identify which variables are to weak to keep in the model. This can be seen by looking the the large confidence intervals that are shown for example in the detection part of the confint, Tamiasciurus Douglassi which is the RIA of the douglas squirrel, this has a large range which overlaps the value 0 indicating there is no clear relationship between them. this is the case as well for the NDVI and Distance to Roads in the occupuncy section.
md1 <- colext(~distance_to_rivers, ~Odocoileus.hemionus, ~season, ~effortd+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)md1 <- colext(~distance_to_rivers, ~Odocoileus.hemionus, ~season, ~effortd+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers, data = deer)
confint(md1, type="psi") # for occupancy (distance_to_rivers) 0.025 0.975
psi(Int) -61.93975 81.90976
psi(distance_to_rivers) -30.05366 28.83469confint(md1, type="col") # for colonisation (Odocoileus.hemionus RIA for Mule Deer) 0.025 0.975
col(Int) -555.2439 535.8151
col(Odocoileus.hemionus) -186.0853 191.3249confint(md1, type="ext") # for extinction (Odocoileus.hemionus) 0.025 0.975
ext(Int) -787.5971 763.8319
ext(seasonsummer) -445.3415 447.5415
ext(seasonwinter) -268.1292 268.1292confint(md1, type="det") # for detection (effortd+Ursus.americanus+Odocoileus.hemionus+distance_to_rivers) 0.025 0.975
p(Int) -2.92104690 -1.32619811
p(effortd) 0.10747377 0.28201560
p(Ursus.americanus) 0.03495282 0.46611223
p(Odocoileus.hemionus) 1.04089233 1.73353617
p(distance_to_rivers) -0.47456665 -0.0769379526 Step 6.8.1 Plot the previously created occupuncy model for the Douglas squirrel
Although the variables that where used to predict occupancy revealed that they didn’t impact it significantly
#first using distance to roads and rivers as a predictor for occupancy
# Load necessary libraries
library(ggplot2)Warning: package 'ggplot2' was built under R version 4.4.2library(dplyr)Warning: package 'dplyr' was built under R version 4.4.2
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union# Generate a range of values for distance to roads and distance to rivers
min_distance_to_roads <- min(siteCovss$distance_to_roads, na.rm = TRUE)
max_distance_to_roads <- max(siteCovss$distance_to_roads, na.rm = TRUE)
distance_to_roads_vals <- seq(min_distance_to_roads, max_distance_to_roads, length.out = 100)
min_distance_to_rivers <- min(siteCovss$distance_to_rivers, na.rm = TRUE)
max_distance_to_rivers <- max(siteCovss$distance_to_rivers, na.rm = TRUE)
distance_to_rivers_vals <- seq(min_distance_to_rivers, max_distance_to_rivers, length.out = 100)
# Calculate the mean values for other covariates
mean_effort <- mean(rowMeans(efforts, na.rm = TRUE), na.rm = TRUE)
mean_NDVI <- mean(siteCovss$NDVI, na.rm = TRUE)
mean_distance_to_roads <- mean(siteCovss$distance_to_roads, na.rm = TRUE)
mean_distance_to_rivers <- mean(siteCovss$distance_to_rivers, na.rm = TRUE)
# Create data frames for predictions
roads_data <- data.frame(
distance_to_roads = distance_to_roads_vals,
distance_to_rivers = mean_distance_to_rivers, # Fix rivers for this plot
NDVI = mean_NDVI,
season = "default_season",
effort = mean_effort
)
rivers_data <- data.frame(
distance_to_roads = mean_distance_to_roads, # Fix roads for this plot
distance_to_rivers = distance_to_rivers_vals,
NDVI = mean_NDVI,
season = "default_season",
effort = mean_effort
)
# Make predictions for both variables
roads_preds <- predict(ms1, newdata = roads_data, type = "psi")
rivers_preds <- predict(ms1, newdata = rivers_data, type = "psi")
# Add predictions and confidence intervals to data frames
roads_data$Occupancy <- roads_preds$Predicted
roads_data$CI_Lower <- roads_preds$Predicted - 1.96 * roads_preds$SE
roads_data$CI_Upper <- roads_preds$Predicted + 1.96 * roads_preds$SE
rivers_data$Occupancy <- rivers_preds$Predicted
rivers_data$CI_Lower <- rivers_preds$Predicted - 1.96 * rivers_preds$SE
rivers_data$CI_Upper <- rivers_preds$Predicted + 1.96 * rivers_preds$SE
# Plot the effect of distance to roads
roads_plot <- ggplot(roads_data, aes(x = distance_to_roads, y = Occupancy)) +
geom_line(color = "blue") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "blue") +
labs(
title = "Effect of Distance to Roads on Occupancy Probability",
x = "Distance to Roads",
y = "Occupancy Probability"
) +
theme_minimal()
# Plot the effect of distance to rivers
rivers_plot <- ggplot(rivers_data, aes(x = distance_to_rivers, y = Occupancy)) +
geom_line(color = "green") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "green") +
labs(
title = "Effect of Distance to Rivers on Occupancy Probability",
x = "Distance to Rivers",
y = "Occupancy Probability"
) +
theme_minimal()
# Display the plots
print(roads_plot)
print(rivers_plot)
after plotting the occupancy probability we can move onto the probability of detection using the model called ms1 the reason a prediction cannot be made for colonisation and extinction is because the variables produced both do not fit into those roles and do not produce statistically valid results.
# Calculate mean values for other detection covariates
mean_elevation <- mean(siteCovss$elevation, na.rm = TRUE)
mean_distance_to_rivers <- mean(siteCovss$distance_to_rivers, na.rm = TRUE)
mean_Tamiasciurus <- mean(siteCovss$Tamiasciurus.douglasii, na.rm = TRUE)
mean_Ursus <- mean(siteCovss$Ursus.americanus, na.rm = TRUE)
# Generate effort values for prediction
effort_vals <- seq(min(efforts, na.rm = TRUE), max(efforts, na.rm = TRUE), length.out = 100)
# Create a data frame for predictions
effort_data <- data.frame(
efforts = effort_vals,
elevation = mean_elevation,
distance_to_rivers = mean_distance_to_rivers,
Tamiasciurus.douglasii = mean_Tamiasciurus,
Ursus.americanus = mean_Ursus
)
# Make predictions for detection probability with confidence intervals
detection_preds <- predict(ms1, newdata = effort_data, type = "det", se.fit = TRUE)
# Add predictions and confidence intervals to data frame
effort_data$Detection_Probability <- detection_preds$fit
# Add predictions and confidence intervals to effort_data
effort_data$Detection_Probability <- detection_preds$Predicted
effort_data$CI_Lower <- detection_preds$lower
effort_data$CI_Upper <- detection_preds$upper
# Plot the effect of effort on detection probability
detection_plot <- ggplot(effort_data, aes(x = efforts, y = Detection_Probability)) +
geom_line(color = "blue") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "blue") +
labs(
title = "Effect of Effort on Detection Probability",
x = "Effort",
y = "Detection Probability"
) +
theme_minimal()
# Display the plot
print(detection_plot)
#using the @ sign instead of the $ sign traditionally used as the squirrel dataset is an unmarked dataframe meaning the data is organised into slots rather than columns. for some reason it didnt work the same for these variables as it did for the others.
site_data <- squirrel@siteCovs
obs_data <- squirrel@obsCovs
# Generate a sequence of elevation values
elevation_vals <- seq(min(site_data$elevation), max(site_data$elevation), length.out = 100)
ria_squirrel <- seq(min(site_data$Tamiasciurus.douglasii), max(site_data$Tamiasciurus.douglasii), length.out = 100)
ria_bear <- seq(min(site_data$Ursus.americanus), max(site_data$Ursus.americanus), length.out = 100)
# Create a data frame for predictions
elevation_data <- data.frame(
efforts = mean(obs_data$efforts, na.rm = TRUE), # Mean of efforts
elevation = elevation_vals, # Sequence of elevation values
Tamiasciurus.douglasii = mean(site_data$Tamiasciurus.douglasii, na.rm = TRUE), # Mean of species presence
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE), # Mean of species presence
distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE) # Mean of distance to rivers
)
# Create a data frame for predictions
doug_ria <- data.frame(
efforts = mean(obs_data$efforts, na.rm = TRUE), # Mean of efforts
elevation = mean(site_data$elevation, na.rm = TRUE), # Mean of elevation
Tamiasciurus.douglasii = ria_squirrel, # sequence of ria for douglas squirrel
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE), # Mean of species presence
distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE) # Mean of distance to rivers
)
# Create a data frame for predictions
bear_ria <- data.frame(
efforts = mean(obs_data$efforts, na.rm = TRUE), # Mean of efforts
elevation = mean(site_data$elevation, na.rm = TRUE), # Mean of elevation
Tamiasciurus.douglasii = mean(site_data$Tamiasciurus.douglasii, na.rm = TRUE), # Mean of the RIA for douglas squirrel
Ursus.americanus = ria_bear, # data sequence of the black bear ria
distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE) # Mean of distance to rivers
)
elevation_preds <- predict(ms1, newdata = elevation_data, type = "det", se.fit = TRUE)
ria_doug_preds <- predict(ms1, newdata = doug_ria, type = "det", se.fit = TRUE)
ria_bear_preds <- predict(ms1, newdata = bear_ria, type = "det", se.fit = TRUE)
# Extract predictions and confidence intervals for elevation
elevation_data$Detection_Probability <- elevation_preds$Predicted
elevation_data$CI_Lower <- elevation_preds$Predicted - 1.96 * elevation_preds$SE
elevation_data$CI_Upper <- elevation_preds$Predicted + 1.96 * elevation_preds$SE
#for the squirrel ria
doug_ria$Detection_Probability <- ria_doug_preds$Predicted
doug_ria$CI_Lower <- ria_doug_preds$Predicted - 1.96 * ria_doug_preds$SE
doug_ria$CI_Upper <- ria_doug_preds$Predicted + 1.96 * ria_doug_preds$SE
#and finally for the bears ria
bear_ria$Detection_Probability <- ria_bear_preds$Predicted
bear_ria$CI_Lower <- ria_bear_preds$Predicted - 1.96 * ria_bear_preds$SE
bear_ria$CI_Upper <- ria_bear_preds$Predicted + 1.96 * ria_bear_preds$SE
library(ggplot2)
ggplot(elevation_data, aes(x = elevation, y = Detection_Probability)) +
geom_line(color = "green") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "green") +
labs(title = "Effect of Elevation on Detection Probability",
x = "Elevation",
y = "Detection Probability") +
theme_minimal()
ggplot(doug_ria, aes(x = Tamiasciurus.douglasii, y = Detection_Probability)) +
geom_line(color = "red") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "red") +
labs(title = "Effect of Douglas's Squirrel abundance on Detection Probability",
x = "Douglas's Squirrel RIA",
y = "Detection Probability") +
theme_minimal()
ggplot(bear_ria, aes(x = Ursus.americanus, y = Detection_Probability)) +
geom_line(color = "purple") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "purple") +
labs(title = "Effect of American Black Bears abundance on Detection Probability",
x = "American Black Bear RIA",
y = "Detection Probability") +
theme_minimal()
27 Step 6.8.2 Plot the Occupancy model for the American Black Bear
site_data <- bear@siteCovs # Extract site-level covariates
obs_data <- bear@obsCovs # Extract observational-level covariate
# Generate a sequence of elevation values
elevation_vals <- seq(min(site_data$elevation), max(site_data$elevation), length.out = 100)
distance_to_roads_vals <- seq(min(site_data$distance_to_roads), max(site_data$distance_to_roads), length.out = 100)
distance_to_urban_vals <- seq(min(site_data$distance_to_urban), max(site_data$distance_to_urban), length.out = 100)
# Create a data frame for predictions whilst calculating the mean
elevation_data <- data.frame(
distance_to_waterbody = mean(site_data$distance_to_waterbody, na.rm = TRUE),
elevation = elevation_vals,
distance_to_roads = mean(site_data$distance_to_roads, na.rm = TRUE),
distance_to_urban = mean(site_data$distance_to_urban, na.rm = TRUE))
# Create a data frame for predictions
roads_data <- data.frame(
distance_to_urban = mean(site_data$distance_to_urban, na.rm = TRUE),
elevation = mean(site_data$elevation, na.rm = TRUE),
distance_to_roads = distance_to_roads_vals,
distance_to_waterbody = mean(site_data$distance_to_waterbody, na.rm = TRUE))
# Create a data frame for predictions
urban_data <- data.frame(
distance_to_urban = distance_to_urban_vals,
elevation = mean(site_data$elevation, na.rm = TRUE),
distance_to_roads = mean(site_data$distance_to_roads, na.rm = TRUE),
distance_to_waterbody = mean(site_data$distance_to_waterbody, na.rm = TRUE))
#introducing model mb1 to the newly created dataframes for the predictions
elevation_preds <- predict(mb1, newdata = elevation_data, type = "psi", se.fit = TRUE)
roads_preds <- predict(mb1, newdata = roads_data, type = "psi", se.fit = TRUE)
urban_preds <- predict(mb1, newdata = urban_data, type = "psi", se.fit = TRUE)
# Extract predictions and confidence intervals for elevation
elevation_data$Detection_Probability <- elevation_preds$Predicted
elevation_data$CI_Lower <- elevation_preds$Predicted - 1.96 * elevation_preds$SE
elevation_data$CI_Upper <- elevation_preds$Predicted + 1.96 * elevation_preds$SE
#for the Distance to roads
roads_data$Detection_Probability <- roads_preds$Predicted
roads_data$CI_Lower <- roads_preds$Predicted - 1.96 * roads_preds$SE
roads_data$CI_Upper <- roads_preds$Predicted + 1.96 * roads_preds$SE
#and finally distance to urban
urban_data$Detection_Probability <- urban_preds$Predicted
urban_data$CI_Lower <- urban_preds$Predicted - 1.96 * urban_preds$SE
urban_data$CI_Upper <- urban_preds$Predicted + 1.96 * urban_preds$SE
library(ggplot2)
ggplot(elevation_data, aes(x = elevation, y = Detection_Probability)) +
geom_line(color = "blue") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "blue") +
labs(title = "Effect of Elevation on Occupancy Probability",
x = "Elevation",
y = "Detection Probability") +
theme_minimal()
ggplot(roads_data, aes(x = distance_to_roads, y = Detection_Probability)) +
geom_line(color = "green") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "green") +
labs(title = "Effect of Distance to roads on Occupancy Probability",
x = "Distance to roads",
y = "Detection Probability") +
theme_minimal()
ggplot(urban_data, aes(x = distance_to_urban, y = Detection_Probability)) +
geom_line(color = "red") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "red") +
labs(title = "Effect of Distance to Urban areas on Occupancy Probability",
x = "Distance to urban",
y = "Detection Probability") +
theme_minimal()
site_data <- bear@siteCovs # Extract site-level covariates
obs_data <- bear@obsCovs # Extract observational-level covariate
# Generate a sequence of elevation values
bear_vals <- seq(min(site_data$Ursus.americanus), max(site_data$Ursus.americanus), length.out = 100)
deer_vals <- seq(min(site_data$Odocoileus.hemionus), max(site_data$Odocoileus.hemionus), length.out = 100)
effort_vals <- seq(min(obs_data$effortb, na.rm = TRUE), max(obs_data$effortb, na.rm = TRUE), length.out = 100)
# Create a data frame for predictions whilst calculating the mean
bear_data <- data.frame(
Ursus.americanus = bear_vals,
effortb = mean(obs_data$effortb, na.rm = TRUE),
Odocoileus.hemionus = mean(site_data$Odocoileus.hemionus, na.rm = TRUE))
# Create a data frame for predictions
deer_data <- data.frame(
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE),
effortb = mean(obs_data$effortb, na.rm = TRUE),
Odocoileus.hemionus = deer_vals)
# Create a data frame for predictions
effort_data <- data.frame(
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE),
effortb = effort_vals,
Odocoileus.hemionus = mean(site_data$Odocoileus.hemionus, na.rm = TRUE))
#introducing model mb1 to the newly created dataframes for the predictions
bear_preds <- predict(mb1, newdata = bear_data, type = "det", se.fit = TRUE)
deer_preds <- predict(mb1, newdata = deer_data, type = "det", se.fit = TRUE)
effort_preds <- predict(mb1, newdata = effort_data, type = "det", se.fit = TRUE)
# Extract predictions and confidence intervals for elevation
bear_data$Detection_Probability <- bear_preds$Predicted
bear_data$CI_Lower <- bear_preds$Predicted - 1.96 * bear_preds$SE
bear_data$CI_Upper <- bear_preds$Predicted + 1.96 * bear_preds$SE
#for the Distance to roads
deer_data$Detection_Probability <- deer_preds$Predicted
deer_data$CI_Lower <- deer_preds$Predicted - 1.96 * deer_preds$SE
deer_data$CI_Upper <- deer_preds$Predicted + 1.96 * deer_preds$SE
#and finally distance to urban
effort_data$Detection_Probability <- effort_preds$Predicted
effort_data$CI_Lower <- effort_preds$Predicted - 1.96 * effort_preds$SE
effort_data$CI_Upper <- effort_preds$Predicted + 1.96 * effort_preds$SE
library(ggplot2)
ggplot(bear_data, aes(x = Ursus.americanus, y = Detection_Probability)) +
geom_line(color = "blue") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "blue") +
labs(title = "Effect of Black Bear RIA on Detection Probability",
x = "RIA of Black bear",
y = "Detection Probability") +
theme_minimal()
ggplot(deer_data, aes(x = Odocoileus.hemionus, y = Detection_Probability)) +
geom_line(color = "green") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "green") +
labs(title = "Effect of Mule Deer RIA on Detection Probability",
x = "RIA of Mule Deer",
y = "Detection Probability") +
theme_minimal()
ggplot(effort_data, aes(x = effortb, y = Detection_Probability)) +
geom_line(color = "red") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "red") +
labs(title = "Effect of effort on Detection Probability",
x = "Effort",
y = "Detection Probability") +
theme_minimal()
28 Step 6.8.3 Plot the Occupancy model for the Mule Deer
Due to large intervals and variables that coincide with one another significantly for the deer, the detection was the only thing that was able to be predicted however this can still inform us about how to create future conservation projects on these species.
site_data <- deer@siteCovs # Extract site-level covariates
obs_data <- deer@obsCovs # Extract observational-level covariate
# Generate a sequence of elevation values
bear_vals <- seq(min(site_data$Ursus.americanus), max(site_data$Ursus.americanus), length.out = 100)
deer_vals <- seq(min(site_data$Odocoileus.hemionus), max(site_data$Odocoileus.hemionus), length.out = 100)
effort_vals <- seq(min(obs_data$effortd, na.rm = TRUE), max(obs_data$effortd, na.rm = TRUE), length.out = 100)
rivers_vals <- seq(min(site_data$distance_to_rivers, na.rm = TRUE), max(site_data$distance_to_rivers, na.rm = TRUE), length.out = 100)
# Create a data frame for predictions whilst calculating the mean
bear_data <- data.frame(
Ursus.americanus = bear_vals,
effortd = mean(obs_data$effortd, na.rm = TRUE),
Odocoileus.hemionus = mean(site_data$Odocoileus.hemionus, na.rm = TRUE), distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE), distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE))
# Create a data frame for predictions
deer_data <- data.frame(
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE),
effortd = mean(obs_data$effortd, na.rm = TRUE),
Odocoileus.hemionus = deer_vals, distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE))
# Create a data frame for predictions
effort_data <- data.frame(
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE),
effortd = effort_vals,
Odocoileus.hemionus = mean(site_data$Odocoileus.hemionus, na.rm = TRUE), distance_to_rivers = mean(site_data$distance_to_rivers, na.rm = TRUE))
river_data <- data.frame(
Ursus.americanus = mean(site_data$Ursus.americanus, na.rm = TRUE),
effortd = mean(obs_data$effortd, na.rm = TRUE),
Odocoileus.hemionus = mean(site_data$Odocoileus.hemionus, na.rm = TRUE), distance_to_rivers = rivers_vals )
#introducing model mb1 to the newly created dataframes for the predictions
bear_preds <- predict(md1, newdata = bear_data, type = "det", se.fit = TRUE)
deer_preds <- predict(md1, newdata = deer_data, type = "det", se.fit = TRUE)
effort_preds <- predict(md1, newdata = effort_data, type = "det", se.fit = TRUE)
river_preds <- predict(md1, newdata = river_data, type = "det", se.fit = TRUE)
# Extract predictions and confidence intervals for deer RIA
deer_data$Detection_Probability <- deer_preds$Predicted
deer_data$CI_Lower <- deer_preds$Predicted - 1.96 * deer_preds$SE
deer_data$CI_Upper <- deer_preds$Predicted + 1.96 * deer_preds$SE
# do the same for effort
effort_data$Detection_Probability <- effort_preds$Predicted
effort_data$CI_Lower <- effort_preds$Predicted - 1.96 * effort_preds$SE
effort_data$CI_Upper <- effort_preds$Predicted + 1.96 * effort_preds$SE
# and lastly distance to rivers
river_data$Detection_Probability <- river_preds$Predicted
river_data$CI_Lower <- river_preds$Predicted - 1.96 * river_preds$SE
river_data$CI_Upper <- river_preds$Predicted + 1.96 * river_preds$SE
library(ggplot2)
ggplot(deer_data, aes(x = Odocoileus.hemionus, y = Detection_Probability)) +
geom_line(color = "green") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "green") +
labs(title = "Effect of Mule Deer RIA on Detection Probability",
x = "RIA of Mule Deer",
y = "Detection Probability") +
theme_minimal()
ggplot(effort_data, aes(x = effortd, y = Detection_Probability)) +
geom_line(color = "red") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "red") +
labs(title = "Effect of effort on Detection Probability",
x = "Effort",
y = "Detection Probability") +
theme_minimal()
ggplot(river_data, aes(x = distance_to_rivers, y = Detection_Probability)) +
geom_line(color = "blue") +
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), alpha = 0.2, fill = "blue") +
labs(title = "Effect of distance to rivers on Detection Probability",
x = "Distance to rivers",
y = "Detection Probability") +
theme_minimal()