#load packages library(amt) library(momentuHMM) library(tidyverse) library(lubridate) library(patchwork) library(purrr) library(dplyr) library(viridis) #### Dispersal #### # Cali Data # Subset the data for each Study Site Cali_dispersal_data <- readRDS("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/California/Dop6/masters/dispersal/Cali_dispersal_covar.rds") Nevada_dispersal_data <- readRDS("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/Nevada/Dop6/masters/dispersal/Nevada_dispersal_covar.rds") # Combine the three data frames into one combined_data <- rbind(Cali_dispersal_data, Nevada_dispersal_data) data2 <- na.omit(combined_data) unique(data2$id) data <- combined_data data2[data2$id == "M198",] names(data) summary(is.na(data2)) data2$study_site <-(as.factor(data2$study_site)) nd_amt1 <- data2 %>% nest(rsteps = -id) #nest data by individual # nested Cali nd_amt1_Cali <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "California",) })) nd_amt1_Cali # Remove any ids that dont have data for that behavior nd_amt1_Cali <- nd_amt1_Cali %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Cali # nested nevada nd_amt1_Nevada <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "Nevada",) })) nd_amt1_Nevada # Remove any ids that dont have data for that behavior nd_amt1_Nevada <- nd_amt1_Nevada %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Nevada names(data2) # Disperal state model issa1_Cali <- nd_amt1_Cali %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + I(sc_elev_end^2) + #Terrain Ruggedness Index(TRI) sc_tri_end + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) # Dispersal state model issa1_Nevada <- nd_amt1_Nevada %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + I(sc_elev_end^2) + #Terrain Ruggedness Index(TRI) sc_tri_end + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) ################# issa1_Cali_tidy <- issa1_Cali %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) issa1_Nevada_tidy <- issa1_Nevada %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) # Make dataframe issa1_Cali_df <- issa1_Cali_tidy %>% unnest(cols = tidy_res) issa1_Nevada_df <- issa1_Nevada_tidy %>% unnest(cols = tidy_res) # Organize beta coefficients and create Confidence Interval for # a one unit change in the covariate for each term. issa1_Cali_iw <- issa1_Cali_df %>% nest(data = c(-term)) issa1_Nevada_iw <- issa1_Nevada_df %>% nest(data = c(-term)) all.eff <- bind_rows( issa1_Cali_df %>% mutate(study_site = "California"), issa1_Nevada_df %>% mutate(study_site = "Nevada") ) %>% nest(data = c(-term)) all.eff$data[[1]] %>% count(study_site) %>% arrange(n) ### all.eff_iw <- all.eff %>% mutate(iw = lapply(data, function(x) { mod <- lm(estimate ~ study_site, data = x, weights = 1 / (std.error)^2) return(mod) })) %>% mutate(pred = lapply(iw, function(x) { pred <- predict(x, newdata = expand_grid(study_site = factor(c("California","Nevada"))), se.fit = TRUE) est <- expand_grid(study_site = factor(c("California","Nevada"))) %>% mutate(mean = pred$fit, lwr = pred$fit - 1.96 * pred$se.fit, upr = pred$fit + 1.96 * pred$se.fit) return(est) })) all.eff_iw$term %>% unique() %>% sort() all.eff_iw$iw all.eff_iw$pred iw_results_dis <- all.eff_iw %>% dplyr::select(term, pred) %>% unnest(cols = pred) #rename covariates iw_results_dis_renamed <- iw_results_dis %>% mutate(term = recode( term, "sc_log_Dist_Developed_end" = "Distance to Developed", "sc_log_Dist_Haycrop_end" = "Distance to Hay and Crop", "sc_log_fourwd_priv_rds_end" = "Distance to 4WD Roads", "sc_elev_end" = "Elevation", "I(sc_elev_end^2)" = "Quadratic Elevation", "sc_tri_end" = "Terrain Ruggedness Index", "sc_log_Dist_forest_end" = "Distance to Forest", "sc_log_Dist_shrub_end" = "Distance to Shrub", "sc_log_Dist_Water_end" = "Distance to Water", "log_sl_" = "Log Step Length", "cos_ta_" = "Cosine Turn Angle", "cos_ta_:sc_log_Dist_Developed_start" = "Distance to Developed:TA", "log_sl_:sc_log_Dist_Developed_start" = "Distance to Developed:SL", "cos_ta_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:TA", "log_sl_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:SL", "cos_ta_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:TA", "log_sl_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:SL", "log_sl_:cos_ta_" = "Log Step Length:Cos Turn Angle", )) identical_cols <- all( identical(iw_results_dis$mean, iw_results_dis_renamed$mean), identical(iw_results_dis$lwr, iw_results_dis_renamed$lwr), identical(iw_results_dis$study_site, iw_results_dis_renamed$study_site), identical(iw_results_dis$upr, iw_results_dis_renamed$upr) ) identical_cols # Add a new column indicating if zero falls between lower and upper bounds iw_results_dis_renamed$Significant <- iw_results_dis_renamed$lwr <= 0 & iw_results_dis_renamed$upr >= 0 # Replace TRUE with " " and FALSE with "Sig" in the zero_between_bounds column iw_results_dis_renamed$Significant <- ifelse(iw_results_dis_renamed$Significant, "---","Sig") iw_results_dis_renamed # Extract rows with Significant value of "Sig" significant_rows_dis <- iw_results_dis_renamed[iw_results_dis_renamed$Significant == "Sig", ] unique(significant_rows_dis$term) #write.csv(iw_results_dis_renamed, "C:/Users/jaran/Desktop/USU_Thesis_Project/Output/THEGOOOOODS/Dispersal_model_ouput.csv") #ggsave("Dispersal_sig_RSS.png", plot = piw_results.sig, width = 10, height = 4, dpi = 750) p.iw_results_all_dis<-iw_results_dis_renamed %>% #filter(!term %in% c("cos_ta", "log_sl","log_sl:cos_ta")) %>% ggplot(aes(x = mean, y = term, color = study_site)) + geom_point(position = position_dodge(width = 0.7), size = 2) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7)) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed")+ labs( x = "log-RSS", y = "Covariate") + theme_classic() + theme(strip.background = element_blank())+ scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) p.iw_results_all_dis #ggsave("Dispersal_all_RSS.png", plot = p.iw_results_all_dis, width = 10, height = 4, dpi = 750) ################################################################################################################## #### Exploratory #### # Cali Data # Subset the data for each Study Site Cali_Exploratory_data <- read.csv("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/California/Dop6/masters/Exploratory/Cali_Exploratory_covar.csv") Nevada_Exploratory_data <- read.csv("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/Nevada/Dop6/masters/Exploratory/Nevada_Exploratory_covar.csv") # Combine the three data frames into one combined_data <- rbind(Cali_Exploratory_data, Nevada_Exploratory_data) data2 <- na.omit(combined_data) unique(data2$id) data <- combined_data data2[data2$id == "M198",] names(data) summary(is.na(data2)) data2$study_site <-(as.factor(data2$study_site)) nd_amt1 <- data2 %>% nest(rsteps = -id) # nested Cali nd_amt1_Cali <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "California",) })) nd_amt1_Cali # Remove any ids that dont have data for that behavior nd_amt1_Cali <- nd_amt1_Cali %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Cali # nested Exploratory nd_amt1_Nevada <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "Nevada",) })) nd_amt1_Nevada # Remove any ids that dont have data for that behavior nd_amt1_Nevada <- nd_amt1_Nevada %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Nevada names(data2) # Exploratory state model issa1_Cali <- nd_amt1_Cali %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + #Terrain Ruggedness Index(TRI) sc_tri_end + I(sc_elev_end^2) + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) # Exploratory State model issa1_Nevada <- nd_amt1_Nevada %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + I(sc_elev_end^2) + #Terrain Ruggedness Index(TRI) sc_tri_end + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) ################# issa1_Cali_tidy <- issa1_Cali %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) issa1_Nevada_tidy <- issa1_Nevada %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) # Create Dataframe issa1_Cali_df <- issa1_Cali_tidy %>% unnest(cols = tidy_res) issa1_Nevada_df <- issa1_Nevada_tidy %>% unnest(cols = tidy_res) # Organize beta coefficients and create Confidence Interval for # a one unit change in the covariate for each term. issa1_Cali_iw <- issa1_Cali_df %>% nest(data = c(-term)) issa1_Nevada_iw <- issa1_Nevada_df %>% nest(data = c(-term)) all.eff <- bind_rows( issa1_Cali_df %>% mutate(study_site = "California"), issa1_Nevada_df %>% mutate(study_site = "Nevada") ) %>% nest(data = c(-term)) all.eff$data[[1]] %>% count(study_site) %>% arrange(n) ### all.eff_iw <- all.eff %>% mutate(iw = lapply(data, function(x) { mod <- lm(estimate ~ study_site, data = x, weights = 1 / (std.error)^2) return(mod) })) %>% mutate(pred = lapply(iw, function(x) { pred <- predict(x, newdata = expand_grid(study_site = factor(c("California","Nevada"))), se.fit = TRUE) est <- expand_grid(study_site = factor(c("California","Nevada"))) %>% mutate(mean = pred$fit, lwr = pred$fit - 1.96 * pred$se.fit, upr = pred$fit + 1.96 * pred$se.fit) return(est) })) all.eff_iw$term %>% unique() %>% sort() all.eff_iw$iw all.eff_iw$pred iw_results_exp <- all.eff_iw %>% dplyr::select(term, pred) %>% unnest(cols = pred) #rename covariates iw_results_exp_renamed <- iw_results_exp %>% mutate(term = recode( term, "sc_log_Dist_Developed_end" = "Distance to Developed", "sc_log_Dist_Haycrop_end" = "Distance to Hay and Crop", "sc_log_fourwd_priv_rds_end" = "Distance to 4WD Roads", "sc_elev_end" = "Elevation", "I(sc_elev_end^2)" = "Quadratic Elevation", "sc_tri_end" = "Terrain Ruggedness Index", "sc_log_Dist_forest_end" = "Distance to Forest", "sc_log_Dist_shrub_end" = "Distance to Shrub", "sc_log_Dist_Water_end" = "Distance to Water", "log_sl_" = "Log Step Length", "cos_ta_" = "Cosine Turn Angle", "cos_ta_:sc_log_Dist_Developed_start" = "Distance to Developed:TA", "log_sl_:sc_log_Dist_Developed_start" = "Distance to Developed:SL", "cos_ta_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:TA", "log_sl_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:SL", "cos_ta_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:TA", "log_sl_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:SL", "log_sl_:cos_ta_" = "Log Step Length:Cos Turn Angle", )) identical_cols <- all( identical(iw_results_exp$mean, iw_results_exp_renamed$mean), identical(iw_results_exp$lwr, iw_results_exp_renamed$lwr), identical(iw_results_exp$study_site, iw_results_exp_renamed$study_site), identical(iw_results_exp$upr, iw_results_exp_renamed$upr) ) identical_cols # Add a new column indicating if zero falls between lower and upper bounds iw_results_exp_renamed$Significant <- iw_results_exp_renamed$lwr <= 0 & iw_results_exp_renamed$upr >= 0 # Replace TRUE with " " and FALSE with "Sig" in the zero_between_bounds column iw_results_exp_renamed$Significant <- ifelse(iw_results_exp_renamed$Significant, "---","Sig") iw_results_exp_renamed #write.csv(iw_results_exp_renamed, "C:/Users/jaran/Desktop/USU_Thesis_Project/Output/THEGOOOOODS/Exploratory_model_ouput.csv") # Extract rows with Significant value of "Sig" significant_rows_exp <- iw_results_exp_renamed[iw_results_exp_renamed$Significant == "Sig", ] unique(significant_rows_exp$term) # significant covariates piw_results.sig<-iw_results_exp_renamed %>% filter(term %in% c(significant_rows_exp$term)) %>% ggplot(aes(x = mean, y = term, color = study_site)) + geom_point(position = position_dodge(width = 0.7), size = 2) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7)) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed") + labs( x = "log-RSS", y = "Covariate") + theme_classic() + theme(strip.background = element_blank()) + scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) piw_results.sig #ggsave("Exploratory_sig_RSS.png", plot = piw_results.sig, width = 10, height = 4, dpi = 750) p.iw_results_all_exp<-iw_results_exp_renamed %>% filter(!term %in% c("cos_ta", "log_sl","log_sl:cos_ta")) %>% ggplot(aes(x = mean, y = term, color = study_site)) + geom_point(position = position_dodge(width = 0.7), size = 2) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7)) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed")+ labs( x = "log-RSS", y = "Covariate") + theme_classic() + theme(strip.background = element_blank())+ scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) p.iw_results_all_exp #ggsave("Exploratory_all_RSS.png", plot = p.iw_results_all_exp, width = 10, height = 4, dpi = 750) ################################################################################################33 #### THR #### # Cali Data # Subset the data for each Study Site Cali_THR_data <- read.csv("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/California/Dop6/masters/THR/Cali_THR_covar.csv") Nevada_THR_data <- read.csv("C:/Users/jaran/Desktop/USU_Thesis_Project/Cleaned_Data/CleanedData/Nevada/Dop6/masters/THR/Nevada_THR_covar.csv") # Combine the three data frames into one combined_data <- rbind(Cali_THR_data, Nevada_THR_data) data2 <- na.omit(combined_data) unique(data2$id) data <- combined_data data2[data2$id == "M198",] names(data) summary(is.na(data2)) data2$study_site <-(as.factor(data2$study_site)) nd_amt1 <- data2 %>% nest(rsteps = -id) # nested Cali nd_amt1_Cali <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "California",) })) nd_amt1_Cali # Remove any ids that dont have data for that behavior nd_amt1_Cali <- nd_amt1_Cali %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Cali # nested THR nd_amt1_Nevada <- nd_amt1 %>% mutate(rsteps = lapply(rsteps, FUN = function(x) { d <- x %>% filter(study_site == "Nevada",) })) nd_amt1_Nevada # Remove any ids that dont have data for that behavior nd_amt1_Nevada <- nd_amt1_Nevada %>% filter(map_lgl(rsteps, ~any(!is.na(.)))) nd_amt1_Nevada names(data2) # THR State model issa1_Cali <- nd_amt1_Cali %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + I(sc_elev_end^2) + #Terrain Ruggedness Index(TRI) sc_tri_end + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) # THR State model issa1_Nevada <- nd_amt1_Nevada %>% mutate(issa = lapply(rsteps, FUN = function (x) { res <- NA res <- try(fit_issf(data = x, formula = case_ ~ ### Anthropogenic Influence #Dist to Developed sc_log_Dist_Developed_end + sc_log_Dist_Developed_start:cos_ta_ + sc_log_Dist_Developed_start:log_sl_ + #Dist to Hay and Crop land sc_log_Dist_Haycrop_end + sc_log_Dist_Haycrop_start:cos_ta_ + sc_log_Dist_Haycrop_start:log_sl_ + # Dist to fourwd rds sc_log_fourwd_priv_rds_end + sc_log_fourwd_priv_rds_start:cos_ta_ + sc_log_fourwd_priv_rds_start:log_sl_ + ### Topography # Elevation sc_elev_end + I(sc_elev_end^2) + #Terrain Ruggedness Index(TRI) sc_tri_end + ### Valuable Lion Habitat # Dist to forest sc_log_Dist_forest_end + # Dist to shrub sc_log_Dist_shrub_end + ### Water # Dist to Water sc_log_Dist_Water_end + log_sl_ + cos_ta_ + cos_ta_:log_sl_ + strata(sid), model = TRUE, method="efron")) })) ################# issa1_Cali_tidy <- issa1_Cali %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) issa1_Nevada_tidy <- issa1_Nevada %>% mutate(tidy_res = lapply(issa, function(x) { if(inherits(x, "fit_clogit")) { res <- broom::tidy(x$model) return(res) } else { return(NA) }})) %>% filter(!is.na(tidy_res)) #Create Dataframe issa1_Cali_df <- issa1_Cali_tidy %>% unnest(cols = tidy_res) issa1_Nevada_df <- issa1_Nevada_tidy %>% unnest(cols = tidy_res) # Organize beta coefficients and create Confidence Interval for # a one unit change in the covariate for each term. issa1_Cali_iw <- issa1_Cali_df %>% nest(data = c(-term)) issa1_Nevada_iw <- issa1_Nevada_df %>% nest(data = c(-term)) all.eff <- bind_rows( issa1_Cali_df %>% mutate(study_site = "California"), issa1_Nevada_df %>% mutate(study_site = "Nevada") ) %>% nest(data = c(-term)) all.eff$data[[1]] %>% count(study_site) %>% arrange(n) ### all.eff_iw <- all.eff %>% mutate(iw = lapply(data, function(x) { mod <- lm(estimate ~ study_site, data = x, weights = 1 / (std.error)^2) return(mod) })) %>% mutate(pred = lapply(iw, function(x) { pred <- predict(x, newdata = expand_grid(study_site = factor(c("California","Nevada"))), se.fit = TRUE) est <- expand_grid(study_site = factor(c("California","Nevada"))) %>% mutate(mean = pred$fit, lwr = pred$fit - 1.96 * pred$se.fit, upr = pred$fit + 1.96 * pred$se.fit) return(est) })) all.eff_iw$term %>% unique() %>% sort() all.eff_iw$iw all.eff_iw$pred iw_results_thr <- all.eff_iw %>% dplyr::select(term, pred) %>% unnest(cols = pred) # rename covariates iw_results_thr_renamed <- iw_results_thr %>% mutate(term = recode( term, "sc_log_Dist_Developed_end" = "Distance to Developed", "sc_log_Dist_Haycrop_end" = "Distance to Hay and Crop", "sc_log_fourwd_priv_rds_end" = "Distance to 4WD Roads", "sc_elev_end" = "Elevation", "I(sc_elev_end^2)" = "Quadratic Elevation", "sc_tri_end" = "Terrain Ruggedness Index", "sc_log_Dist_forest_end" = "Distance to Forest", "sc_log_Dist_shrub_end" = "Distance to Shrub", "sc_log_Dist_Water_end" = "Distance to Water", "log_sl_" = "Log Step Length", "cos_ta_" = "Cosine Turn Angle", "cos_ta_:sc_log_Dist_Developed_start" = "Distance to Developed:TA", "log_sl_:sc_log_Dist_Developed_start" = "Distance to Developed:SL", "cos_ta_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:TA", "log_sl_:sc_log_Dist_Haycrop_start" = "Distance to Hay and Crop:SL", "cos_ta_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:TA", "log_sl_:sc_log_fourwd_priv_rds_start" = "Distance to 4WD Roads:SL", "log_sl_:cos_ta_" = "Log Step Length:Cos Turn Angle", )) identical_cols <- all( identical(iw_results_thr$mean, iw_results_thr_renamed$mean), identical(iw_results_thr$lwr, iw_results_thr_renamed$lwr), identical(iw_results_thr$study_site, iw_results_thr_renamed$study_site), identical(iw_results_thr$upr, iw_results_thr_renamed$upr) ) identical_cols # Add a new column indicating if zero falls between lower and upper bounds iw_results_thr_renamed$Significant <- iw_results_thr_renamed$lwr <= 0 & iw_results_thr_renamed$upr >= 0 # Replace TRUE with " " and FALSE with "Sig" in the zero_between_bounds column iw_results_thr_renamed$Significant <- ifelse(iw_results_thr_renamed$Significant, "---","Sig") iw_results_thr_renamed #write.csv(iw_results_thr_renamed, "C:/Users/jaran/Desktop/USU_Thesis_Project/Output/THEGOOOOODS/THR_model_ouput.csv") # Extract rows with Significant value of "Sig" significant_rows_thr <- iw_results_thr_renamed[iw_results_thr_renamed$Significant == "Sig", ] unique(significant_rows_thr$term) # significant covariates piw_results.sig<-iw_results_thr_renamed %>% filter(term %in% c(significant_rows_thr$term)) %>% ggplot(aes(x = mean, y = term, color = study_site)) + geom_point(position = position_dodge(width = 0.7), size = 2) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7)) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed") + labs( x = "log-RSS", y = "Covariate") + theme_classic() + theme(strip.background = element_blank()) + scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) piw_results.sig #ggsave("THR_sig_RSS.png", plot = piw_results.sig, width = 10, height = 4, dpi = 750) p.iw_results_all_thr<-iw_results_thr_renamed %>% filter(!term %in% c("cos_ta", "log_sl","log_sl:cos_ta")) %>% ggplot(aes(x = mean, y = term, color = study_site)) + geom_point(position = position_dodge(width = 0.7), size = 2) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7)) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed")+ labs( x = "log-RSS", y = "Covariate") + theme_classic() + theme(strip.background = element_blank())+ scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) p.iw_results_all_thr #ggsave("THR_all_RSS.png", plot = p.iw_results_all_thr, width = 10, height = 4, dpi = 750) ### Make one single plot with all behaviors in them ############################################### #### p.iw_results_all_dis p.iw_results_all_exp p.iw_results_all_thr # Assuming p.iw_results_all_dis, p.iw_results_all_exp, and p.iw_results_all_thr are your data data_dis <- iw_results_dis_renamed data_exp <- iw_results_exp_renamed data_thr <- iw_results_thr_renamed # Load the gridExtra package library(gridExtra) library(ggplot2) library(cowplot) library(ggbreak) library(ggsignif) # Your existing code plot_dis <- iw_results_dis_renamed %>% filter(!term %in% c("Cosine Turn Angle", "Log Step Length", "Log Step Length:Cos Turn Angle")) %>% ggplot(aes(x = mean, y = term, color = study_site, alpha = Significant)) + geom_point(position = position_dodge(width = 0.7), size = 2.5) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7), size = .8) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed") + labs(x = "log-RSS", y = "") + theme_classic() + theme(strip.background = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank()) + scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) + guides(alpha = "none", color = "none") + scale_alpha_manual(values = c("---" = .3, "Sig" = 1)) # Adjust alpha values as needed plot_dis # ############################################# plot_exp <- iw_results_exp_renamed %>% filter(!term %in% c("Cosine Turn Angle", "Log Step Length","Log Step Length:Cos Turn Angle")) %>% ggplot(aes(x = mean, y = term, color = study_site, alpha = Significant)) + geom_point(position = position_dodge(width = 0.7), size = 2.5) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7), size = .8) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed") + labs(x = "", y = "") + theme_classic() + theme(strip.background = element_blank()) + scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue")) + guides(alpha = "none", color = "none") + lims(x = c(-1,1)) + # Set x-axis limits scale_alpha_manual(values = c("---" = .3, "Sig" = 1)) plot_thr <- iw_results_thr_renamed %>% filter(!term %in% c("Cosine Turn Angle", "Log Step Length", "Log Step Length:Cos Turn Angle")) %>% ggplot(aes(x = mean, y = term, color = study_site, alpha = Significant)) + geom_point(position = position_dodge(width = 0.7), size = 2.5) + geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.6, position = position_dodge(width = 0.7), size = .8) + facet_grid() + geom_vline(xintercept = 0, linetype = "dashed") + labs(x = "", y = "") + theme_classic() + theme(strip.background = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank()) + scale_color_manual(name = "Study Site", values = c("California" = "red", "Nevada" = "blue"), labels = c("California" = "Protected", "Nevada" = "Hunted")) + # Change the labels guides(alpha = "none") + scale_alpha_manual(values = c("---" = 0.3, "Sig" = 1)) # Combine plots into a 1x3 grid with centered titles all_behave <- grid.arrange( ggdraw() + draw_label("Exploratory", size = 12, hjust = -.8), ggdraw() + draw_label("Departure", size = 12, hjust = 0.6), ggdraw() + draw_label("Transient Home Range", size = 12, hjust = .6), plot_exp, plot_dis, plot_thr, ncol = 3, heights = c(0.1, 3.2) ) plot_exp <- plot_exp + theme(axis.text.x = element_text(size = 14), # Adjust x-axis font size axis.text.y = element_text(size = 14)) # Adjust y-axis font size plot_dis <- plot_dis + theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14)) # Adjust y-axis font size plot_thr <- plot_thr + theme(axis.text.x = element_text(size = 14), legend.text = element_text(size = 16), # Adjust legend text size legend.title = element_text(size = 16)) # Adjust legend title size) # Adjust y-axis font size # Combine plots into a 1x3 grid with centered titles all_behave <- grid.arrange( ggdraw() + draw_label("Exploratory", size = 14, hjust = -.925), # Center the title ggdraw() + draw_label("Departure", size = 14, hjust = 0.59), # Center the title ggdraw() + draw_label("Transient Home Range", size = 14, hjust = .65), # Center the title plot_exp, plot_dis, plot_thr, ncol = 3, heights = c(0.1, 3.2) ) ggsave("all_behave_opacity.png", plot = all_behave, width = 16, height = 8, dpi = 750)