In pbs-assess/gfranges: Groundfish Ranges
knitr::opts_chunk$set(
fig.width = 11, fig.height = 8.5,
echo = FALSE, warning = FALSE, message = FALSE
)
options(scipen = 999)
options(digits = 4)
Set species
# # Species run so far...
# species <- "Arrowtooth Flounder"
species <- "Pacific Cod"
# species <- "Sablefish"
# species <- "Silvergray Rockfish"
# species <- "Lingcod"
# species <- "North Pacific Spiny Dogfish" # note: using all data for maturity thresholds
# species <- "Quillback Rockfish"
# species <- "Pacific Ocean Perch"
# species <- "Yelloweye Rockfish"
Get species-specific optimal values from mountain plots for most abundant year
spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species)))
params <- readRDS(paste0("data/optimal-values-by-species.rds"))
params <- params[rev(order(as.Date(params$date))), ]
params <- params[!duplicated(params$species), ]
optimal_do <- round(params[params$species == species, ]$optimal_do, digits = 2)
optimal_temp <- round(params[params$species == species, ]$optimal_temp, digits = 2)
optimal_values <- c(optimal_do, optimal_temp)
Set vector type
climate <- "temperature"
# climate <- "do"
# climate <- "do and temperature"
Make climate vectors for biannual changes between surveys
trend <- FALSE
start_time <- 2012
end_time <- 2015
temp_threshold <- c(1)
do_threshold <- c(0.1)
OR
Make climate vectors for average changes and lower threshold values
# trend <- TRUE
# temp_trend_threshold <- c(0.25)
# do_trend_threshold <- c(0.2)
Set Region
# region <- "West Coast Vancouver Island"
region <- "both odd year surveys"
# region <- "West Coast Haida Gwaii"
Run all subsequent code...
library(dplyr)
library(ggplot2)
library(gfplot)
# library(gfdata)
library(sdmTMB)
#library(itsadug)
library(gfranges)
if (region == "both odd year surveys") {
survey <- c("SYN QCS", "SYN HS")
model_ssid <- c(1, 3)
ssid_string <- paste0(model_ssid, collapse = "n")
trend_indices_do <- c(1, 1, 1, 2, 2)
trend_indices_temp <- c(1, 1, 1, 2, 2)
years <- NULL
}
if (region == "West Coast Vancouver Island") {
survey <- c("SYN WCVI")
model_ssid <- c(4)
ssid_string <- paste0(model_ssid, collapse = "n")
trend_indices_do <- c(1, 1, 1, 2, 2)
trend_indices_temp <- c(1, 1, 1, 2, 2, 2)
years <- NULL
}
if (region == "West Coast Haida Gwaii") {
survey <- c("SYN WCHG")
model_ssid <- c(16)
ssid_string <- paste0(model_ssid, collapse = "n")
trend_indices_do <- c(1, 1, 1, 2, 2)
trend_indices_temp <- c(1, 1, 1, 2, 2)
years <- NULL
}
# Set default cell sizes and scale
input_cell_size <- 2
scale_fac <- 1
delta_t_step <- 2
skip_time <- NULL
if (climate == "temperature") {
variable_names <- c("temp")
min_thresholds <- c(Inf)
#delta_threshold <- temp_threshold
optimal_values <- optimal_temp
if (trend) {
max_thresholds <- temp_trend_threshold
start_time <- 2008 # change if can add back in prior 2008 data
end_time <- NULL
delta_t_total <- 6 # change if can add back in prior 2008 data
indices <- trend_indices_temp
} else {
max_thresholds <- temp_threshold
indices <- c(1, 2)
delta_t_total <- 2
}
}
if (climate == "do") {
variable_names <- c("do_est")
max_thresholds <- c(Inf)
#delta_threshold <- -1 * (do_threshold)
optimal_values <- optimal_do
if (trend) {
min_thresholds <- do_trend_threshold
start_time <- 2008
skip_time <- 2016
end_time <- NULL
delta_t_total <- 6
indices <- trend_indices_do
} else {
min_thresholds <- do_threshold
indices <- c(1, 2)
delta_t_total <- 2
}
}
if (climate == "do and temperature") {
variable_names <- c("do_est", "temp")
if (trend) {
min_thresholds <- c(do_trend_threshold, Inf)
max_thresholds <- c(Inf, temp_trend_threshold)
start_time <- 2008
skip_time <- 2016
end_time <- NULL
delta_t_total <- 6
indices <- trend_indices_do
} else {
min_thresholds <- c(do_threshold, Inf)
max_thresholds <- c(Inf, temp_threshold)
indices <- c(1, 2)
delta_t_total <- 2
}
}
min_string <- paste0(min_thresholds, collapse = "-")
max_string <- paste0(max_thresholds, collapse = "-")
optimal_string <- paste0(optimal_values, collapse = "-")
climate_string <- gsub(" ", "-", gsub("\\/", "-", tolower(climate)))
Calculate climate vectors
rm(vocc)
try({
vocc <- readRDS(paste0(
"data/vocc-", ssid_string, "-", climate_string, "-",
min_string, "to", max_string, "-", start_time, "-", delta_t_total, ".rds"
))
})
if (!exists("vocc")) {
pred_temp_do <- readRDS(here::here("analysis/tmb-sensor-explore/models/predicted-DO.rds"))
if (length(variable_names) > 1) {
climate_predictions <- replicate(length(variable_names), pred_temp_do, simplify = FALSE)
} else {
climate_predictions <- pred_temp_do
}
starttime <- Sys.time()
vocc <- make_vector_data(climate_predictions,
variable_names = variable_names,
ssid = model_ssid,
start_time = start_time,
end_time = end_time,
skip_time = skip_time,
input_cell_size = input_cell_size,
scale_fac = scale_fac,
delta_t_total = delta_t_total,
delta_t_step = delta_t_step,
indices = indices,
min_thresholds = min_thresholds,
max_thresholds = max_thresholds,
round_fact = 10
)
endtime <- Sys.time()
time_vocc <- round(starttime - endtime)
vocc$var_1_min <- min_thresholds[1]
vocc$var_2_min <- min_thresholds[2]
vocc$var_1_max <- max_thresholds[1]
vocc$var_2_max <- max_thresholds[2]
saveRDS(vocc, file = paste0(
"data/vocc-", ssid_string, "-", climate_string, "-",
min_string, "to", max_string, "-", start_time, "-", delta_t_total, ".rds"
))
}
# glimpse(vocc)
Plot all possible vectors
if (climate == "temperature") {
vocc$diff <- vocc$units_per_decade / 10 * delta_t_total
vocc_plot <- plot_vocc(vocc,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle(climate)
}
if (climate == "do") {
vocc$diff <- vocc$units_per_decade / 10 * delta_t_total
vocc_plot <- plot_vocc(vocc,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
high_fill = "Steel Blue 4",
low_fill = "Red 3",
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle(climate)
}
if (climate == "do and temperature") {
vocc$temp_diff <- vocc$temp.units_per_decade / 10 * delta_t_total
vocc_plot_temp <- plot_vocc(vocc,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "temp_diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle("temperature")
vocc$do_diff <- vocc$do_est.units_per_decade / 10 * delta_t_total
vocc_plot_do <- plot_vocc(vocc,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "do_diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
high_fill = "Steel Blue 4",
low_fill = "Red 3",
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle("do")
vocc_plot <- gridExtra::grid.arrange(
vocc_plot_temp, vocc_plot_do,
ncol = 2,
top = grid::textGrob(paste("do of ", optimal_values[1], "less", min_thresholds[1], "or temperature of ", optimal_values[2], "plus", max_thresholds[2], " over ", start_time, "-", start_time + delta_t_total, ""))
)
}
vocc_plot
Filter vectors to match optimal range for chosen species
vectors <- trim_vector_data(vocc, variable_names, optimal_values, min_thresholds, max_thresholds,
cell_size = input_cell_size, dist_intercept = input_cell_size/2,
max_dist = 120
)
# if (length(vectors[, 2]) == 0) {
# other_coefs <- readRDS(paste0("data/", climate_string, "-vocc-by-spp-sdmTMB-coefs.RDS"))
# coefs <- other_coefs[1, ]
# coefs$var <- "NA"
# coefs$beta <- "NA"
# coefs$se <- "NA"
# coefs$model <- "NA"
# coefs$species <- spp
# coefs$lifestage <- "NA"
# coefs$climate <- climate
# coefs$min_thresholds <- min_string
# coefs$max_thresholds <- max_string
# coefs$optimum <- optimal_string
# coefs$start_time <- start_time
# coefs$timespan <- delta_t_total
# coefs$region <- region
# coefs$ssid_string <- ssid_string
#
# coefs <- rbind(other_coefs, coefs) %>% distinct()
# saveRDS(coefs, file = paste0("data/", climate_string, "-vocc-by-spp-sdmTMB-coefs.RDS"))
#
# stop("No vectors. Conditions saved to 'data/vocc-byspp-sdmTMB-coefs'.")
# }
# glimpse(vocc)
# glimpse(vectors)
Plot remaining climate vectors for regions where change exceeds thresholds
if (climate == "temperature") {
vectors$diff <- vectors$units_per_decade / 10 * delta_t_total
vocc_plot <- plot_vocc(vectors,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle(climate)
}
if (climate == "do") {
vectors$diff <- vectors$units_per_decade / 10 * delta_t_total
vocc_plot <- plot_vocc(vectors,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
high_fill = "Steel Blue 4",
low_fill = "Red 3",
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle(climate)
}
if (climate == "do and temperature") {
vectors$temp_diff <- vectors$temp.units_per_decade / 10 * delta_t_total
vocc_plot_temp <- plot_vocc(vectors,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "temp_diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle("temperature")
vectors$do_diff <- vectors$do_est.units_per_decade / 10 * delta_t_total
vocc_plot_do <- plot_vocc(vectors,
vec_aes = "distance",
min_vec_plotted = input_cell_size,
NA_label = ".",
fill_col = "do_diff",
fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""),
raster_alpha = 1,
white_zero = TRUE,
high_fill = "Steel Blue 4",
low_fill = "Red 3",
vec_alpha = 0.25,
axis_lables = FALSE,
transform_col = no_trans
) + ggtitle("do")
vocc_plot <- gridExtra::grid.arrange(
vocc_plot_temp, vocc_plot_do,
ncol = 2,
top = grid::textGrob(paste("DO of ", optimal_values[1], "less", min_thresholds[1], "or temperature of ", optimal_values[2], "plus", max_thresholds[2], " over ", start_time, "-", start_time + delta_t_total, ""))
)
}
vocc_plot
Add species biomass data
Get all biomass predictions
priors <- FALSE
# priors <- TRUE
covariates <- "+muddy+any_rock"
covs <- gsub("\\+", "-", covariates)
rm(adult_predictions)
rm(imm_predictions)
rm(predictions)
try({
adult_predictions <- readRDS(paste0(
"data/", spp,
"/all-predictions-", spp, covs, "-mat-biomass-prior-", priors, ".rds"
))
})
try({
imm_predictions <- readRDS(paste0(
"data/", spp,
"/all-predictions-", spp, covs, "-imm-biomass-prior-", priors, ".rds"
))
})
try({
predictions <- readRDS(paste0(
"data/", spp,
"/all-predictions-", spp, covs, "-total-biomass-prior-", priors, ".rds"
))
})
if (!exists("adult_predictions")) adult_predictions <- NULL
if (!exists("imm_predictions")) imm_predictions <- NULL
if (!exists("predictions")) predictions <- NULL
biomass_list <- list(
adult = adult_predictions,
imm = imm_predictions,
total = predictions
)
biomass_list <- biomass_list[!sapply(biomass_list, is.null)]
# glimpse(biomass_list)
Extract biomass values to match spatial and temporal range of climate vectors
alldata <- lapply(biomass_list, function(i) {
biomass_change <- make_trend_data(i,
ssid = model_ssid,
start_time = start_time,
skip_time = skip_time,
end_time = end_time,
input_cell_size = input_cell_size,
scale_fac = scale_fac,
delta_t_total = delta_t_total,
delta_t_step = delta_t_step,
indices = indices,
variable_names = "est_exp"
)
biomass_covs <- i %>% select(X, Y, depth, trawled, muddy, any_rock, omega_s) %>% distinct()
biomass_change <- biomass_change %>%
rename(start_density = est_exp_s, end_density = est_exp_e) %>%
dplyr::select(icell, x, y, start_density, end_density)
# dplyr::select(-N, -time_step, - units_per_decade, -slope)
biomass <- left_join(biomass_change, biomass_covs, by = c("x" = "X", "y"="Y"))
# keep just one row for each cell that exceeded the end climate threshold
source_cells <- vectors %>% distinct(icell, .keep_all = TRUE) %>%
select(-target_X, -target_Y, -target_values, -tid)
source_biomass <- left_join(source_cells, biomass, by = c("icell", "x", "y")) %>%
mutate(cell_type = "source") %>% filter(!is.na(start_density))
# identify all possible control cells, assume distance = 0 given
control_cells <- anti_join(vocc, vectors, by = c("icell", "x", "y")) %>%
distinct(icell, .keep_all = TRUE) %>%
mutate(distance = 0) %>%
select(-target_X, -target_Y, -target_values, -tid)
control_biomass <- left_join(control_cells, biomass, by = c("icell", "x", "y")) %>%
mutate(cell_type = "control") %>% filter(!is.na(start_density))
alldata <- rbind(source_biomass, control_biomass)
alldata <- mutate(alldata, treat = ifelse(cell_type=="source", 1, 0)) %>%
select(-var_2_min, -var_2_max) %>%
rename(vect_dist = distance)
alldata$treat <- as.factor(alldata$treat)
alldata
})
lapply(alldata, function(i){ggplot(i, aes(x=cell_type, y=log(start_density))) + geom_boxplot()})
Select matched control cells
matched <- lapply(alldata, function(alldata) {
## trim non-overlapping values to reduce size of dataset
control_min <- quantile(alldata$start_density[alldata$treat=="0"], 0.01)
source_min <- quantile(alldata$start_density[alldata$treat=="1"], 0.01)
# control_max <- quantile(alldata$start_density[alldata$treat=="0"], 1)
# source_max <- quantile(alldata$start_density[alldata$treat=="1"], 1)
group_min <- max(control_min, source_min)
# group_max <- min(control_max, source_max)
alldata <- filter(alldata, start_density >= group_min) #%>% filter(start_density <= group_max)
# need to be standardized? seems not
alldata$depth2 <- alldata$depth^2
alldata$start_density.st <- round(arm::rescale(log(alldata$start_density))*5)/5
alldata$start_density_st <- round(arm::rescale(log(alldata$start_density)), 1)
alldata$depth.st <- arm::rescale(log(alldata$depth))
alldata$muddy.st <- arm::rescale(alldata$muddy)
alldata$any_rock.st <- arm::rescale(alldata$any_rock)
#alldata$omega_s.st <- arm::rescale(alldata$omega_s)
## Using optimization ... doesn't work
# if(
# class(
# try(
# fit_formula <- paste0("treat ~ start_density_st + depth.st + muddy.st + any_rock.st")
# # fit_formula <- paste0("treat ~ start_density_st")
# # fit_formula <- paste0("treat ~ start_density_st + depth.st")
# # # for (j in (variable_names)) {
# # # fit_formula <- paste0("", fit_formula, " + ", j, "_s")
# # # }
# # options("optmatch_max_problem_size" = Inf)
# matched <- MatchIt::matchit(as.formula(fit_formula), method = "optimal", data = alldata)
# , silent = TRUE)
# ) == "try-error"){
# Using exact match and rounding
fit_formula <- paste0("treat ~ depth.st + muddy.st + any_rock.st + trawled + omega_s + x + y") #
# for (j in (variable_names)) {
# fit_formula <- paste0("", fit_formula, " + ", j, "_s")
# }
# Choose data order
# alldata <- alldata[(order(alldata$y)), ]
# if (order == "increasing") alldata <- alldata[rev(order(alldata$start_density)), ]
# if (order == "decreasing") alldata <- alldata[rev(order(alldata$start_density)), ]
matched <- MatchIt::matchit(as.formula(fit_formula), exact = c("start_density.st"),
method = "nearest", m.order = "random",
distance = "probit",
discard="control",
data = alldata)
#}
matched
})
Print diagnostics for matching
lapply(matched, function(i){
summary(i)
plot(i)
})
Prep matched data for analysis
vocc_by_sp <- lapply(matched, function(matched) {
mdata <- MatchIt::match.data(matched, distance = "pscore")
sources <- mdata[row.names(matched$match.matrix), ]
sources$origobs <- row.names(sources)
sources$matchobs <- matched$match.matrix[,1]
controls <- mdata[matched$match.matrix, ]
controls$origobs <- row.names(controls)
controls$matchobs <- matched$match.matrix[,1]
# if selecting 2 matches
# controls2 <- mdata[matched$match.matrix, ]
# controls2$origobs <- row.names(controls)
# controls2$matchobs <- matched$match.matrix[,2]
mdata <- rbind(sources, controls) %>% tidyr::drop_na() %>% filter(matchobs!="-1")
matched_time0 <- mdata %>% mutate(after = 0, density = start_density)
#%>% dplyr::select(-start_density, -end_density)
matched_time1 <- mdata %>% mutate(after = 1, density = end_density)
#%>% dplyr::select(-start_density, -end_density)
for (j in seq_along(variable_names)) {
climate <- paste0(variable_names[j])
matched_time0 <- matched_time0 %>%
mutate((!!climate) := (UQ(rlang::sym(paste0(variable_names[j], "_s")))))
matched_time1 <- matched_time1 %>%
mutate((!!climate) := (UQ(rlang::sym(paste0(variable_names[j], "_e")))))
}
vocc_by_sp <- rbind(matched_time0, matched_time1)
# calculate standardized climate variables
for (j in seq_along(variable_names)) {
mean.var <- paste0("mean.", variable_names[j])
sd.var <- paste0("sd.", variable_names[j])
std.var <- paste0("std_", variable_names[j])
std.var2 <- paste0("std_", variable_names[j], "2")
vocc_by_sp <- mutate(
vocc_by_sp,
(!!mean.var) := mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE),
(!!sd.var) := sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE),
(!!std.var) := ((UQ(rlang::sym(paste0(variable_names[j]))) -
mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE)) /
(sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE))*2),
(!!std.var2) := (( UQ(rlang::sym(paste0(variable_names[j]))) -
mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE) ) /
(sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE)*2))^2
)
}
# rename spatial coordinates for inclusion in models and plots without naming conflicts
data <- vocc_by_sp %>%
rename(X = x, Y = y)
#%>% rename(distance = vect_dist) %>% select(-start_density.st, -depth.st, -muddy.st, -any_rock.st, -omega_s.st)
# other variable prep
data$log_density <- log(data$density)
data$vect_dist2 <- data$vect_dist^2
data <- mutate(data, time = ifelse(after==1, "end", "begin"))
data$matchobs <- as.factor(data$matchobs)
data$origobs <- as.factor(data$origobs)
data
})
Adults
data <- vocc_by_sp [["adult"]]
# View(data)
Plot matches
ggplot(data,
#ggplot(filter(data, as.numeric(matchobs) <= 550),
aes(y=Y, x=X, shape=as.factor(treat), colour=as.factor(matchobs))) +
geom_point(alpha=0.8, size = 1.5) +
scale_shape_manual(values= c(0, 15)) +
scale_colour_viridis_d(option="C", guide = FALSE) +
theme_bw()
lapply(variable_names, function(j) {
ggplot(data, aes(y=log_density, x=UQ(rlang::sym(paste(j))), colour=time)) +
geom_point(alpha=0.35) +
scale_colour_viridis_d(option="C") +
theme_bw()
})
ggplot(data, aes(y=log_density, x=depth, colour=cell_type)) +
geom_point(alpha=0.35) +
scale_colour_viridis_d(option="C") +
facet_wrap(~time) +
theme_bw()
ggplot(data, aes(log_density, fill = cell_type)) +
geom_density(alpha=0.5) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)+
theme_bw()
ggplot(data, aes(log_density, fill = cell_type)) +
geom_histogram(alpha=0.45) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)+
theme_bw()
ggplot(data, aes(y=log_density, x=cell_type, fill=cell_type)) +
geom_boxplot(alpha=0.5) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)+
theme_bw()
Basic BACI
for (j in (variable_names)) {
baci_formula <- paste0("log_density~after + cell_type + after:cell_type + n_targets + vect_dist + vect_dist2 + std_", j, "+ std_", j, "2 + (1|icell) + (1|matchobs)")
}
m <- lmerTest::lmer(as.formula(baci_formula), data=data)
summary(m)
sjPlot::plot_model(m, type = "std2", order.terms=c(1,2,8,3,4,5,6,7), title = paste(species, "adult"))
GAMM
for (j in (variable_names)) {
gamm_formula <- paste0("log_density ~ time * cell_type + s(vect_dist) + s(std_", j, ") + s(X,Y)")
}
gam4 <- gamm4::gamm4(as.formula(gamm_formula), random = ~(1|matchobs) + (1|origobs), data=data) # + s(X,Y, bs="re")
summary(gam4$gam) ## summary of gam
plot(gam4$gam,pages=1)
plot_parametric(gam4$gam, pred=list(time=c("begin","end"), cell_type=c("control", "source")), parametricOnly = TRUE, main="", xlab ="log_density"
#, groups = cell_type, gcolor = my_cols, color = my_cols["time"]
)
Immatures
idata <- vocc_by_sp[["imm"]]
#View(data)
Plot matches
ggplot(data,
#ggplot(filter(data, as.numeric(matchobs) <= 550),
aes(y=Y, x=X, shape=as.factor(treat), colour=matchobs)) +
geom_point(alpha=0.4, size = 2) +
scale_shape_manual(values= c(0, 15)) +
scale_colour_viridis_d(option="C", guide = FALSE) +
theme_bw()
lapply(variable_names, function(j) {
ggplot(idata, aes(y=log_density, x=UQ(rlang::sym(paste(j))), colour=time)) +
geom_point(alpha=0.35) +
scale_colour_viridis_d(option="C")
})
ggplot(idata, aes(log_density, fill = cell_type)) +
geom_density(alpha=0.5) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)
ggplot(idata, aes(log_density, fill = cell_type)) +
geom_histogram(alpha=0.5) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)
ggplot(idata, aes(y=log_density, x=cell_type, fill=cell_type)) +
geom_boxplot(alpha=0.5) +
scale_fill_viridis_d(option="C") +
facet_wrap(~time)
Basic BACI
mi <- lmerTest::lmer(as.formula(baci_formula), data=idata)
summary(mi)
sjPlot::plot_model(mi, type = "std2", order.terms=c(1,2,8,3,4,5,6,7), title = paste(species, "immature"))
GAMM
# for (j in (variable_names)) {
# gamm_formula <- paste0("log_density ~ time * cell_type + s(vect_dist) + s(std_", j, ") + s(X,Y)")
# }
igam4 <- gamm4::gamm4(as.formula(gamm_formula), random = ~(1|matchobs) + (1|origobs), data=idata) # + s(X,Y, bs="re")
summary(igam4$gam) ## summary of gam
plot(igam4$gam,pages=1)
plot_parametric(igam4$gam, pred=list(time=c("begin","end"), cell_type=c("control", "source")), parametricOnly = TRUE, main="", xlab ="log_density"
#, groups = cell_type, gcolor = my_cols, color = my_cols["time"]
)
Misc exploring gamm
# summary(igam4$mer) ## underlying mixed model
# anova(igam4$gam)
# itsadug::plot_parametric(gam4$gam, pred=list(cell_type=c("control", "source"), after=c("0","1")), parametricOnly = TRUE, main="")
# mi <- mgcv::gamm(log_density~ as.factor(after) + cell_type + after:cell_type + std_temp + std_temp2 + vect_dist + vect_dist2 + s(as.factor(idata$icell), bs="re") + s(as.factor(idata$matchobs), bs="re"), data=idata)
# summary(mi)
#
# mi <- mgcv::gamm(log_density~ as.factor(after) * cell_type + s(std_temp) + s(vect_dist) + s(as.factor(idata$matchobs), bs="re") , data=idata) # + s(X,Y, bs="re")
# summary(mi)
# mi
#
#
# idata$after <- as.factor(idata$after)
# idata$matchobs <- as.factor(idata$matchobs)
#
# m4b <- gamm4::gamm4(log_density~ after * cell_type + s(std_temp) + s(vect_dist) + #s(X) + s(Y) +
# s(X,Y), random = ~(1|matchobs), data=idata) # + s(X,Y, bs="re")
#
#
# plot(m4b$gam,pages=1)
#
# summary(m4b$gam) ## summary of gam
# summary(m4b$mer) ## underlying mixed model
# anova(m4b$gam)
#
# gamtabs(m1)
#
# #library(itsadug)
# plot_parametric(m4b$gam, pred=list(cell_type=c("control", "source"), after=c("0","1")))
#
# library(mgcv)
pbs-assess/gfranges documentation built on Dec. 13, 2021, 4:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( fig.width = 11, fig.height = 8.5, echo = FALSE, warning = FALSE, message = FALSE ) options(scipen = 999) options(digits = 4)
Set species
# # Species run so far... # species <- "Arrowtooth Flounder" species <- "Pacific Cod" # species <- "Sablefish" # species <- "Silvergray Rockfish" # species <- "Lingcod" # species <- "North Pacific Spiny Dogfish" # note: using all data for maturity thresholds # species <- "Quillback Rockfish" # species <- "Pacific Ocean Perch" # species <- "Yelloweye Rockfish"
Get species-specific optimal values from mountain plots for most abundant year
spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species))) params <- readRDS(paste0("data/optimal-values-by-species.rds")) params <- params[rev(order(as.Date(params$date))), ] params <- params[!duplicated(params$species), ] optimal_do <- round(params[params$species == species, ]$optimal_do, digits = 2) optimal_temp <- round(params[params$species == species, ]$optimal_temp, digits = 2) optimal_values <- c(optimal_do, optimal_temp)
Set vector type
climate <- "temperature" # climate <- "do" # climate <- "do and temperature"
Make climate vectors for biannual changes between surveys
trend <- FALSE start_time <- 2012 end_time <- 2015 temp_threshold <- c(1) do_threshold <- c(0.1)
OR
Make climate vectors for average changes and lower threshold values
# trend <- TRUE # temp_trend_threshold <- c(0.25) # do_trend_threshold <- c(0.2)
Set Region
# region <- "West Coast Vancouver Island" region <- "both odd year surveys" # region <- "West Coast Haida Gwaii"
Run all subsequent code...
library(dplyr) library(ggplot2) library(gfplot) # library(gfdata) library(sdmTMB) #library(itsadug) library(gfranges)
if (region == "both odd year surveys") { survey <- c("SYN QCS", "SYN HS") model_ssid <- c(1, 3) ssid_string <- paste0(model_ssid, collapse = "n") trend_indices_do <- c(1, 1, 1, 2, 2) trend_indices_temp <- c(1, 1, 1, 2, 2) years <- NULL } if (region == "West Coast Vancouver Island") { survey <- c("SYN WCVI") model_ssid <- c(4) ssid_string <- paste0(model_ssid, collapse = "n") trend_indices_do <- c(1, 1, 1, 2, 2) trend_indices_temp <- c(1, 1, 1, 2, 2, 2) years <- NULL } if (region == "West Coast Haida Gwaii") { survey <- c("SYN WCHG") model_ssid <- c(16) ssid_string <- paste0(model_ssid, collapse = "n") trend_indices_do <- c(1, 1, 1, 2, 2) trend_indices_temp <- c(1, 1, 1, 2, 2) years <- NULL }
# Set default cell sizes and scale input_cell_size <- 2 scale_fac <- 1 delta_t_step <- 2 skip_time <- NULL if (climate == "temperature") { variable_names <- c("temp") min_thresholds <- c(Inf) #delta_threshold <- temp_threshold optimal_values <- optimal_temp if (trend) { max_thresholds <- temp_trend_threshold start_time <- 2008 # change if can add back in prior 2008 data end_time <- NULL delta_t_total <- 6 # change if can add back in prior 2008 data indices <- trend_indices_temp } else { max_thresholds <- temp_threshold indices <- c(1, 2) delta_t_total <- 2 } } if (climate == "do") { variable_names <- c("do_est") max_thresholds <- c(Inf) #delta_threshold <- -1 * (do_threshold) optimal_values <- optimal_do if (trend) { min_thresholds <- do_trend_threshold start_time <- 2008 skip_time <- 2016 end_time <- NULL delta_t_total <- 6 indices <- trend_indices_do } else { min_thresholds <- do_threshold indices <- c(1, 2) delta_t_total <- 2 } } if (climate == "do and temperature") { variable_names <- c("do_est", "temp") if (trend) { min_thresholds <- c(do_trend_threshold, Inf) max_thresholds <- c(Inf, temp_trend_threshold) start_time <- 2008 skip_time <- 2016 end_time <- NULL delta_t_total <- 6 indices <- trend_indices_do } else { min_thresholds <- c(do_threshold, Inf) max_thresholds <- c(Inf, temp_threshold) indices <- c(1, 2) delta_t_total <- 2 } } min_string <- paste0(min_thresholds, collapse = "-") max_string <- paste0(max_thresholds, collapse = "-") optimal_string <- paste0(optimal_values, collapse = "-") climate_string <- gsub(" ", "-", gsub("\\/", "-", tolower(climate)))
Calculate climate vectors
rm(vocc) try({ vocc <- readRDS(paste0( "data/vocc-", ssid_string, "-", climate_string, "-", min_string, "to", max_string, "-", start_time, "-", delta_t_total, ".rds" )) }) if (!exists("vocc")) { pred_temp_do <- readRDS(here::here("analysis/tmb-sensor-explore/models/predicted-DO.rds")) if (length(variable_names) > 1) { climate_predictions <- replicate(length(variable_names), pred_temp_do, simplify = FALSE) } else { climate_predictions <- pred_temp_do } starttime <- Sys.time() vocc <- make_vector_data(climate_predictions, variable_names = variable_names, ssid = model_ssid, start_time = start_time, end_time = end_time, skip_time = skip_time, input_cell_size = input_cell_size, scale_fac = scale_fac, delta_t_total = delta_t_total, delta_t_step = delta_t_step, indices = indices, min_thresholds = min_thresholds, max_thresholds = max_thresholds, round_fact = 10 ) endtime <- Sys.time() time_vocc <- round(starttime - endtime) vocc$var_1_min <- min_thresholds[1] vocc$var_2_min <- min_thresholds[2] vocc$var_1_max <- max_thresholds[1] vocc$var_2_max <- max_thresholds[2] saveRDS(vocc, file = paste0( "data/vocc-", ssid_string, "-", climate_string, "-", min_string, "to", max_string, "-", start_time, "-", delta_t_total, ".rds" )) } # glimpse(vocc)
Plot all possible vectors
if (climate == "temperature") { vocc$diff <- vocc$units_per_decade / 10 * delta_t_total vocc_plot <- plot_vocc(vocc, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle(climate) } if (climate == "do") { vocc$diff <- vocc$units_per_decade / 10 * delta_t_total vocc_plot <- plot_vocc(vocc, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, high_fill = "Steel Blue 4", low_fill = "Red 3", vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle(climate) } if (climate == "do and temperature") { vocc$temp_diff <- vocc$temp.units_per_decade / 10 * delta_t_total vocc_plot_temp <- plot_vocc(vocc, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "temp_diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle("temperature") vocc$do_diff <- vocc$do_est.units_per_decade / 10 * delta_t_total vocc_plot_do <- plot_vocc(vocc, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "do_diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, high_fill = "Steel Blue 4", low_fill = "Red 3", vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle("do") vocc_plot <- gridExtra::grid.arrange( vocc_plot_temp, vocc_plot_do, ncol = 2, top = grid::textGrob(paste("do of ", optimal_values[1], "less", min_thresholds[1], "or temperature of ", optimal_values[2], "plus", max_thresholds[2], " over ", start_time, "-", start_time + delta_t_total, "")) ) } vocc_plot
Filter vectors to match optimal range for chosen species
vectors <- trim_vector_data(vocc, variable_names, optimal_values, min_thresholds, max_thresholds, cell_size = input_cell_size, dist_intercept = input_cell_size/2, max_dist = 120 ) # if (length(vectors[, 2]) == 0) { # other_coefs <- readRDS(paste0("data/", climate_string, "-vocc-by-spp-sdmTMB-coefs.RDS")) # coefs <- other_coefs[1, ] # coefs$var <- "NA" # coefs$beta <- "NA" # coefs$se <- "NA" # coefs$model <- "NA" # coefs$species <- spp # coefs$lifestage <- "NA" # coefs$climate <- climate # coefs$min_thresholds <- min_string # coefs$max_thresholds <- max_string # coefs$optimum <- optimal_string # coefs$start_time <- start_time # coefs$timespan <- delta_t_total # coefs$region <- region # coefs$ssid_string <- ssid_string # # coefs <- rbind(other_coefs, coefs) %>% distinct() # saveRDS(coefs, file = paste0("data/", climate_string, "-vocc-by-spp-sdmTMB-coefs.RDS")) # # stop("No vectors. Conditions saved to 'data/vocc-byspp-sdmTMB-coefs'.") # } # glimpse(vocc) # glimpse(vectors)
Plot remaining climate vectors for regions where change exceeds thresholds
if (climate == "temperature") { vectors$diff <- vectors$units_per_decade / 10 * delta_t_total vocc_plot <- plot_vocc(vectors, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle(climate) } if (climate == "do") { vectors$diff <- vectors$units_per_decade / 10 * delta_t_total vocc_plot <- plot_vocc(vectors, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, high_fill = "Steel Blue 4", low_fill = "Red 3", vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle(climate) } if (climate == "do and temperature") { vectors$temp_diff <- vectors$temp.units_per_decade / 10 * delta_t_total vocc_plot_temp <- plot_vocc(vectors, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "temp_diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle("temperature") vectors$do_diff <- vectors$do_est.units_per_decade / 10 * delta_t_total vocc_plot_do <- plot_vocc(vectors, vec_aes = "distance", min_vec_plotted = input_cell_size, NA_label = ".", fill_col = "do_diff", fill_label = paste("Change from\n", start_time, "to", start_time + delta_t_total, ""), raster_alpha = 1, white_zero = TRUE, high_fill = "Steel Blue 4", low_fill = "Red 3", vec_alpha = 0.25, axis_lables = FALSE, transform_col = no_trans ) + ggtitle("do") vocc_plot <- gridExtra::grid.arrange( vocc_plot_temp, vocc_plot_do, ncol = 2, top = grid::textGrob(paste("DO of ", optimal_values[1], "less", min_thresholds[1], "or temperature of ", optimal_values[2], "plus", max_thresholds[2], " over ", start_time, "-", start_time + delta_t_total, "")) ) } vocc_plot
Add species biomass data
Get all biomass predictions
priors <- FALSE # priors <- TRUE covariates <- "+muddy+any_rock" covs <- gsub("\\+", "-", covariates) rm(adult_predictions) rm(imm_predictions) rm(predictions) try({ adult_predictions <- readRDS(paste0( "data/", spp, "/all-predictions-", spp, covs, "-mat-biomass-prior-", priors, ".rds" )) }) try({ imm_predictions <- readRDS(paste0( "data/", spp, "/all-predictions-", spp, covs, "-imm-biomass-prior-", priors, ".rds" )) }) try({ predictions <- readRDS(paste0( "data/", spp, "/all-predictions-", spp, covs, "-total-biomass-prior-", priors, ".rds" )) }) if (!exists("adult_predictions")) adult_predictions <- NULL if (!exists("imm_predictions")) imm_predictions <- NULL if (!exists("predictions")) predictions <- NULL biomass_list <- list( adult = adult_predictions, imm = imm_predictions, total = predictions ) biomass_list <- biomass_list[!sapply(biomass_list, is.null)] # glimpse(biomass_list)
Extract biomass values to match spatial and temporal range of climate vectors
alldata <- lapply(biomass_list, function(i) { biomass_change <- make_trend_data(i, ssid = model_ssid, start_time = start_time, skip_time = skip_time, end_time = end_time, input_cell_size = input_cell_size, scale_fac = scale_fac, delta_t_total = delta_t_total, delta_t_step = delta_t_step, indices = indices, variable_names = "est_exp" ) biomass_covs <- i %>% select(X, Y, depth, trawled, muddy, any_rock, omega_s) %>% distinct() biomass_change <- biomass_change %>% rename(start_density = est_exp_s, end_density = est_exp_e) %>% dplyr::select(icell, x, y, start_density, end_density) # dplyr::select(-N, -time_step, - units_per_decade, -slope) biomass <- left_join(biomass_change, biomass_covs, by = c("x" = "X", "y"="Y")) # keep just one row for each cell that exceeded the end climate threshold source_cells <- vectors %>% distinct(icell, .keep_all = TRUE) %>% select(-target_X, -target_Y, -target_values, -tid) source_biomass <- left_join(source_cells, biomass, by = c("icell", "x", "y")) %>% mutate(cell_type = "source") %>% filter(!is.na(start_density)) # identify all possible control cells, assume distance = 0 given control_cells <- anti_join(vocc, vectors, by = c("icell", "x", "y")) %>% distinct(icell, .keep_all = TRUE) %>% mutate(distance = 0) %>% select(-target_X, -target_Y, -target_values, -tid) control_biomass <- left_join(control_cells, biomass, by = c("icell", "x", "y")) %>% mutate(cell_type = "control") %>% filter(!is.na(start_density)) alldata <- rbind(source_biomass, control_biomass) alldata <- mutate(alldata, treat = ifelse(cell_type=="source", 1, 0)) %>% select(-var_2_min, -var_2_max) %>% rename(vect_dist = distance) alldata$treat <- as.factor(alldata$treat) alldata })
lapply(alldata, function(i){ggplot(i, aes(x=cell_type, y=log(start_density))) + geom_boxplot()})
Select matched control cells
matched <- lapply(alldata, function(alldata) { ## trim non-overlapping values to reduce size of dataset control_min <- quantile(alldata$start_density[alldata$treat=="0"], 0.01) source_min <- quantile(alldata$start_density[alldata$treat=="1"], 0.01) # control_max <- quantile(alldata$start_density[alldata$treat=="0"], 1) # source_max <- quantile(alldata$start_density[alldata$treat=="1"], 1) group_min <- max(control_min, source_min) # group_max <- min(control_max, source_max) alldata <- filter(alldata, start_density >= group_min) #%>% filter(start_density <= group_max) # need to be standardized? seems not alldata$depth2 <- alldata$depth^2 alldata$start_density.st <- round(arm::rescale(log(alldata$start_density))*5)/5 alldata$start_density_st <- round(arm::rescale(log(alldata$start_density)), 1) alldata$depth.st <- arm::rescale(log(alldata$depth)) alldata$muddy.st <- arm::rescale(alldata$muddy) alldata$any_rock.st <- arm::rescale(alldata$any_rock) #alldata$omega_s.st <- arm::rescale(alldata$omega_s) ## Using optimization ... doesn't work # if( # class( # try( # fit_formula <- paste0("treat ~ start_density_st + depth.st + muddy.st + any_rock.st") # # fit_formula <- paste0("treat ~ start_density_st") # # fit_formula <- paste0("treat ~ start_density_st + depth.st") # # # for (j in (variable_names)) { # # # fit_formula <- paste0("", fit_formula, " + ", j, "_s") # # # } # # options("optmatch_max_problem_size" = Inf) # matched <- MatchIt::matchit(as.formula(fit_formula), method = "optimal", data = alldata) # , silent = TRUE) # ) == "try-error"){ # Using exact match and rounding fit_formula <- paste0("treat ~ depth.st + muddy.st + any_rock.st + trawled + omega_s + x + y") # # for (j in (variable_names)) { # fit_formula <- paste0("", fit_formula, " + ", j, "_s") # } # Choose data order # alldata <- alldata[(order(alldata$y)), ] # if (order == "increasing") alldata <- alldata[rev(order(alldata$start_density)), ] # if (order == "decreasing") alldata <- alldata[rev(order(alldata$start_density)), ] matched <- MatchIt::matchit(as.formula(fit_formula), exact = c("start_density.st"), method = "nearest", m.order = "random", distance = "probit", discard="control", data = alldata) #} matched })
Print diagnostics for matching
lapply(matched, function(i){ summary(i) plot(i) })
Prep matched data for analysis
vocc_by_sp <- lapply(matched, function(matched) { mdata <- MatchIt::match.data(matched, distance = "pscore") sources <- mdata[row.names(matched$match.matrix), ] sources$origobs <- row.names(sources) sources$matchobs <- matched$match.matrix[,1] controls <- mdata[matched$match.matrix, ] controls$origobs <- row.names(controls) controls$matchobs <- matched$match.matrix[,1] # if selecting 2 matches # controls2 <- mdata[matched$match.matrix, ] # controls2$origobs <- row.names(controls) # controls2$matchobs <- matched$match.matrix[,2] mdata <- rbind(sources, controls) %>% tidyr::drop_na() %>% filter(matchobs!="-1") matched_time0 <- mdata %>% mutate(after = 0, density = start_density) #%>% dplyr::select(-start_density, -end_density) matched_time1 <- mdata %>% mutate(after = 1, density = end_density) #%>% dplyr::select(-start_density, -end_density) for (j in seq_along(variable_names)) { climate <- paste0(variable_names[j]) matched_time0 <- matched_time0 %>% mutate((!!climate) := (UQ(rlang::sym(paste0(variable_names[j], "_s"))))) matched_time1 <- matched_time1 %>% mutate((!!climate) := (UQ(rlang::sym(paste0(variable_names[j], "_e"))))) } vocc_by_sp <- rbind(matched_time0, matched_time1) # calculate standardized climate variables for (j in seq_along(variable_names)) { mean.var <- paste0("mean.", variable_names[j]) sd.var <- paste0("sd.", variable_names[j]) std.var <- paste0("std_", variable_names[j]) std.var2 <- paste0("std_", variable_names[j], "2") vocc_by_sp <- mutate( vocc_by_sp, (!!mean.var) := mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE), (!!sd.var) := sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE), (!!std.var) := ((UQ(rlang::sym(paste0(variable_names[j]))) - mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE)) / (sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE))*2), (!!std.var2) := (( UQ(rlang::sym(paste0(variable_names[j]))) - mean(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE) ) / (sd(UQ(rlang::sym(paste0(variable_names[j]))), na.rm = TRUE)*2))^2 ) } # rename spatial coordinates for inclusion in models and plots without naming conflicts data <- vocc_by_sp %>% rename(X = x, Y = y) #%>% rename(distance = vect_dist) %>% select(-start_density.st, -depth.st, -muddy.st, -any_rock.st, -omega_s.st) # other variable prep data$log_density <- log(data$density) data$vect_dist2 <- data$vect_dist^2 data <- mutate(data, time = ifelse(after==1, "end", "begin")) data$matchobs <- as.factor(data$matchobs) data$origobs <- as.factor(data$origobs) data })
Adults
data <- vocc_by_sp [["adult"]] # View(data)
Plot matches
ggplot(data, #ggplot(filter(data, as.numeric(matchobs) <= 550), aes(y=Y, x=X, shape=as.factor(treat), colour=as.factor(matchobs))) + geom_point(alpha=0.8, size = 1.5) + scale_shape_manual(values= c(0, 15)) + scale_colour_viridis_d(option="C", guide = FALSE) + theme_bw() lapply(variable_names, function(j) { ggplot(data, aes(y=log_density, x=UQ(rlang::sym(paste(j))), colour=time)) + geom_point(alpha=0.35) + scale_colour_viridis_d(option="C") + theme_bw() }) ggplot(data, aes(y=log_density, x=depth, colour=cell_type)) + geom_point(alpha=0.35) + scale_colour_viridis_d(option="C") + facet_wrap(~time) + theme_bw() ggplot(data, aes(log_density, fill = cell_type)) + geom_density(alpha=0.5) + scale_fill_viridis_d(option="C") + facet_wrap(~time)+ theme_bw() ggplot(data, aes(log_density, fill = cell_type)) + geom_histogram(alpha=0.45) + scale_fill_viridis_d(option="C") + facet_wrap(~time)+ theme_bw() ggplot(data, aes(y=log_density, x=cell_type, fill=cell_type)) + geom_boxplot(alpha=0.5) + scale_fill_viridis_d(option="C") + facet_wrap(~time)+ theme_bw()
Basic BACI
for (j in (variable_names)) { baci_formula <- paste0("log_density~after + cell_type + after:cell_type + n_targets + vect_dist + vect_dist2 + std_", j, "+ std_", j, "2 + (1|icell) + (1|matchobs)") } m <- lmerTest::lmer(as.formula(baci_formula), data=data) summary(m)
sjPlot::plot_model(m, type = "std2", order.terms=c(1,2,8,3,4,5,6,7), title = paste(species, "adult"))
GAMM
for (j in (variable_names)) { gamm_formula <- paste0("log_density ~ time * cell_type + s(vect_dist) + s(std_", j, ") + s(X,Y)") } gam4 <- gamm4::gamm4(as.formula(gamm_formula), random = ~(1|matchobs) + (1|origobs), data=data) # + s(X,Y, bs="re") summary(gam4$gam) ## summary of gam plot(gam4$gam,pages=1)
plot_parametric(gam4$gam, pred=list(time=c("begin","end"), cell_type=c("control", "source")), parametricOnly = TRUE, main="", xlab ="log_density" #, groups = cell_type, gcolor = my_cols, color = my_cols["time"] )
Immatures
idata <- vocc_by_sp[["imm"]] #View(data)
Plot matches
ggplot(data, #ggplot(filter(data, as.numeric(matchobs) <= 550), aes(y=Y, x=X, shape=as.factor(treat), colour=matchobs)) + geom_point(alpha=0.4, size = 2) + scale_shape_manual(values= c(0, 15)) + scale_colour_viridis_d(option="C", guide = FALSE) + theme_bw() lapply(variable_names, function(j) { ggplot(idata, aes(y=log_density, x=UQ(rlang::sym(paste(j))), colour=time)) + geom_point(alpha=0.35) + scale_colour_viridis_d(option="C") }) ggplot(idata, aes(log_density, fill = cell_type)) + geom_density(alpha=0.5) + scale_fill_viridis_d(option="C") + facet_wrap(~time) ggplot(idata, aes(log_density, fill = cell_type)) + geom_histogram(alpha=0.5) + scale_fill_viridis_d(option="C") + facet_wrap(~time) ggplot(idata, aes(y=log_density, x=cell_type, fill=cell_type)) + geom_boxplot(alpha=0.5) + scale_fill_viridis_d(option="C") + facet_wrap(~time)
Basic BACI
mi <- lmerTest::lmer(as.formula(baci_formula), data=idata) summary(mi)
sjPlot::plot_model(mi, type = "std2", order.terms=c(1,2,8,3,4,5,6,7), title = paste(species, "immature"))
GAMM
# for (j in (variable_names)) { # gamm_formula <- paste0("log_density ~ time * cell_type + s(vect_dist) + s(std_", j, ") + s(X,Y)") # } igam4 <- gamm4::gamm4(as.formula(gamm_formula), random = ~(1|matchobs) + (1|origobs), data=idata) # + s(X,Y, bs="re") summary(igam4$gam) ## summary of gam plot(igam4$gam,pages=1)
plot_parametric(igam4$gam, pred=list(time=c("begin","end"), cell_type=c("control", "source")), parametricOnly = TRUE, main="", xlab ="log_density" #, groups = cell_type, gcolor = my_cols, color = my_cols["time"] )
Misc exploring gamm
# summary(igam4$mer) ## underlying mixed model # anova(igam4$gam) # itsadug::plot_parametric(gam4$gam, pred=list(cell_type=c("control", "source"), after=c("0","1")), parametricOnly = TRUE, main="") # mi <- mgcv::gamm(log_density~ as.factor(after) + cell_type + after:cell_type + std_temp + std_temp2 + vect_dist + vect_dist2 + s(as.factor(idata$icell), bs="re") + s(as.factor(idata$matchobs), bs="re"), data=idata) # summary(mi) # # mi <- mgcv::gamm(log_density~ as.factor(after) * cell_type + s(std_temp) + s(vect_dist) + s(as.factor(idata$matchobs), bs="re") , data=idata) # + s(X,Y, bs="re") # summary(mi) # mi # # # idata$after <- as.factor(idata$after) # idata$matchobs <- as.factor(idata$matchobs) # # m4b <- gamm4::gamm4(log_density~ after * cell_type + s(std_temp) + s(vect_dist) + #s(X) + s(Y) + # s(X,Y), random = ~(1|matchobs), data=idata) # + s(X,Y, bs="re") # # # plot(m4b$gam,pages=1) # # summary(m4b$gam) ## summary of gam # summary(m4b$mer) ## underlying mixed model # anova(m4b$gam) # # gamtabs(m1) # # #library(itsadug) # plot_parametric(m4b$gam, pred=list(cell_type=c("control", "source"), after=c("0","1"))) # # library(mgcv)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.