Skip to content

Commit 690a7b2

Browse files
author
marlonecobos
committed
changes to make models with and without replicates work with all other functions.
1 parent da83398 commit 690a7b2

File tree

4 files changed

+120
-78
lines changed

4 files changed

+120
-78
lines changed

‎R/fit_selected.R‎

Lines changed: 44 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -124,7 +124,7 @@ fit_selected <- function(calibration_results,
124124
m_ids <- calibration_results$selected_models$ID
125125
algorithm <- calibration_results$algorithm
126126

127-
# Fitting models over multiple replicates
127+
# Fitting models if they have multiple replicates
128128
if (n_replicates > 1) {
129129
if (verbose) {
130130
message("Fitting replicates...")
@@ -224,6 +224,7 @@ fit_selected <- function(calibration_results,
224224
if (verbose) {
225225
message("\nFitting full models...")
226226
}
227+
227228
# Full models grid setup
228229
n_models <- length(m_ids)
229230
dfgrid <- expand.grid(models = m_ids, replicates = 1)
@@ -278,6 +279,11 @@ fit_selected <- function(calibration_results,
278279
}
279280
}
280281

282+
# Stop the cluster for full models
283+
if (parallel) {
284+
parallel::stopCluster(cl)
285+
}
286+
281287
# Assign names to full models
282288
names(full_models) <- paste0("Model_", m_ids)
283289

@@ -286,10 +292,6 @@ fit_selected <- function(calibration_results,
286292
best_models[[i]]$Full_model <- full_models[[i]]
287293
}
288294

289-
# Stop the cluster for full models
290-
if (parallel) {
291-
parallel::stopCluster(cl)
292-
}
293295

294296
# Compute thresholds for predictions
295297
occ <- calibration_results$calibration_data[
@@ -298,7 +300,7 @@ fit_selected <- function(calibration_results,
298300
# Predictions and consensus for occurrences
299301
p_occ <- lapply(names(best_models), function(x) {
300302
m_x <- best_models[[x]]
301-
if (any(grepl("Rep", names(m_x)))) {
303+
if (n_replicates > 1) {
302304
m_x$Full_model <- NULL
303305
}
304306

@@ -308,33 +310,59 @@ fit_selected <- function(calibration_results,
308310
type = type))
309311
} else if (algorithm == "glm") {
310312
p_r <- sapply(m_x, function(i) suppressWarnings(
311-
as.numeric(predict_glm_mx(model = i,
312-
newdata = occ, type = type))
313+
as.numeric(predict_glm_mx(model = i, newdata = occ, type = type))
313314
))
314315
}
315316

316-
p_mean <- apply(p_r, 1, mean, na.rm = TRUE)
317-
p_median <- apply(p_r, 1, median, na.rm = TRUE)
317+
if (n_replicates > 1) {
318+
p_mean <- apply(p_r, 1, mean, na.rm = TRUE)
319+
p_median <- apply(p_r, 1, median, na.rm = TRUE)
318320

319-
list(mean = p_mean, median = p_median)
321+
list(mean = p_mean, median = p_median, rep = p_r)
322+
} else {
323+
list(Full_model = p_r[, 1])
324+
}
320325
})
321326

322327
names(p_occ) <- names(best_models)
323328

324329
# Calculate consensus across models
325-
mean_consensus <- apply(sapply(p_occ, function(x) x$mean), 1,
326-
mean, na.rm = TRUE)
327-
median_consensus <- apply(sapply(p_occ, function(x) x$median), 1,
328-
median, na.rm = TRUE)
330+
if (length(p_occ) == 1) {
331+
if (n_replicates > 1) {
332+
mean_consensus <- p_occ[[1]]$mean
333+
median_consensus <- p_occ[[1]]$median
334+
} else {
335+
mean_consensus <- p_occ[[1]]$Full_model
336+
median_consensus <- p_occ[[1]]$Full_model
337+
}
338+
339+
} else {
340+
if (n_replicates > 1) {
341+
mean_consensus <- apply(sapply(p_occ, function(x) x$mean), 1,
342+
mean, na.rm = TRUE)
343+
median_consensus <- apply(do.call(cbind, lapply(p_occ, `[[`, "rep")), 1,
344+
median, na.rm = TRUE)
345+
} else {
346+
mean_consensus <- apply(sapply(p_occ, function(x) x$Full_model), 1,
347+
mean, na.rm = TRUE)
348+
median_consensus <- apply(sapply(p_occ, function(x) x$Full_model), 1,
349+
median, na.rm = TRUE)
350+
}
351+
}
352+
353+
p_occ <- lapply(p_occ, function(x) x[names(x) != "rep"])
354+
329355
consensus <- list(mean = mean_consensus, median = median_consensus)
356+
330357
p_occ <- c(p_occ, list(consensus = consensus))
331358

332359
# Calculate thresholds
333360
p_thr <- lapply(p_occ, function(model) {
334361
lapply(model, calc_thr,
335362
thr = calibration_results$summary$omission_rate_thr / 100)
336363
})
337-
#Append type of predictions
364+
365+
# Append type of predictions
338366
p_thr$type <- type
339367

340368
#Prepare final data

‎R/helpers_project_selected_glmnetx.R‎

Lines changed: 10 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -68,9 +68,11 @@ var_models_rep_by_gcm <- function(path) {
6868
mean_replicates <- terra::rast(lapply(1:n_replicates, function(n) {
6969
rep_n <- terra::mean(rast(lapply(r_x, function(x) x[[n]])))
7070
}))
71-
var_rep_x <- terra::app(mean_replicates, "var")} else {
72-
r_x <- terra::rast(model_files)
73-
var_rep_x <- terra::app(r_x, "var")}
71+
var_rep_x <- terra::app(mean_replicates, "var")
72+
} else {
73+
r_x <- terra::rast(model_files)
74+
var_rep_x <- terra::app(r_x, "var")
75+
}
7476
names(var_rep_x) <- basename(path)
7577
return(var_rep_x)
7678
}
@@ -80,7 +82,8 @@ var_models_model_by_gcm <- function(path, consensus) {
8082
full.names = TRUE))
8183
r_x <- r_x[[sapply(r_x, function(r) names(r) == consensus)]]
8284
if (terra::nlyr(r_x) > 1) {
83-
var_x <- terra::app(r_x, "var") } else {
85+
var_x <- terra::app(r_x, "var")
86+
} else {
8487
var_x <- r_x * 0
8588
}
8689
return(var_x)
@@ -124,7 +127,7 @@ check_pred_scenarios <- function(projection_data, out_dir) {
124127
#Present
125128
if ("Present" %in% sc) {
126129
#Create folder
127-
present_dir <- file.path(out_dir, "Present/")
130+
present_dir <- file.path(out_dir, "Present")
128131
present_sc <- names(projection_data[["Present"]])
129132
suppressWarnings({
130133
d_present <- data.frame(
@@ -140,7 +143,7 @@ check_pred_scenarios <- function(projection_data, out_dir) {
140143
#Past
141144
if ("Past" %in% sc) {
142145
#Create folder
143-
past_dir <- file.path(out_dir, "Past/")
146+
past_dir <- file.path(out_dir, "Past")
144147
#Get grid of projections
145148
df_past <- do.call(
146149
rbind,
@@ -172,7 +175,7 @@ check_pred_scenarios <- function(projection_data, out_dir) {
172175
####Project to Future scenarios####
173176
if ("Future" %in% sc) {
174177
#Create folder
175-
future_dir <- file.path(out_dir, "Future/")
178+
future_dir <- file.path(out_dir, "Future")
176179

177180
#Create grid of time-ssp-gcm
178181
df_future <- do.call(

‎R/projection_variability.R‎

Lines changed: 49 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -6,22 +6,23 @@
66
#' replicates, model parameterizations, and general circulation models (GCMs).
77
#'
88
#' @usage
9-
#' projection_variability(model_projections, by_replicate = TRUE, by_gcm = TRUE,
10-
#' by_model = TRUE, consensus = "median",
11-
#' write_files = FALSE, output_dir = NULL,
12-
#' return_rasters = TRUE, progress_bar = FALSE,
13-
#' verbose = TRUE, overwrite = FALSE)
9+
#' projection_variability(model_projections, from_replicates = TRUE,
10+
#' from_parameters = TRUE, from_gcms = TRUE,
11+
#' consensus = "median", write_files = FALSE,
12+
#' output_dir = NULL, return_rasters = TRUE,
13+
#' progress_bar = FALSE, verbose = TRUE,
14+
#' overwrite = FALSE)
1415
#'
1516
#' @param model_projections a `model_projections` object generated by the
1617
#' \code{\link{project_selected}}() function. This object contains the file
1718
#' paths to the raster projection results and the thresholds used for binarizing
1819
#' the predictions.
19-
#' @param by_replicate (logical) whether to compute the variance originating
20+
#' @param from_replicates (logical) whether to compute the variance originating
2021
#' from replicates.
21-
#' @param by_gcm (logical) whether to compute the variance originating from
22-
#' general circulation models (GCMs)
23-
#' @param by_model (logical) whether to compute the variance originating from
22+
#' @param from_parameters (logical) whether to compute the variance originating from
2423
#' model parameterizations.
24+
#' @param from_gcms (logical) whether to compute the variance originating from
25+
#' general circulation models (GCMs)
2526
#' @param consensus (character) (character) the consensus measure to use for
2627
#' calculating changes. Available options are 'mean', 'median', 'range', and
2728
#' 'stdev' (standard deviation). Default is 'median'.
@@ -108,19 +109,19 @@
108109
#' out_dir = out_dir)
109110
#'
110111
#' # Step 5: Compute variance from distinct sources
111-
#' v <- projection_variability(model_projections = p, by_replicate = FALSE)
112+
#' v <- projection_variability(model_projections = p, from_replicates = FALSE)
112113
#'
113-
#' #terra::plot(v$Present$by_replicate) # Variance from replicates, present projection
114-
#' terra::plot(v$Present$by_model) # From models
115-
#' #terra::plot(v$`Future_2041-2060_ssp126`$by_replicate) # From replicates in future projection
116-
#' terra::plot(v$`Future_2041-2060_ssp126`$by_model) # From models
117-
#' terra::plot(v$`Future_2041-2060_ssp126`$by_gcm) # From GCMs
114+
#' #terra::plot(v$Present$from_replicates) # Variance from replicates, present projection
115+
#' terra::plot(v$Present$from_parameters) # From models with distinct parameters
116+
#' #terra::plot(v$`Future_2041-2060_ssp126`$from_replicates) # From replicates in future projection
117+
#' terra::plot(v$`Future_2041-2060_ssp126`$from_parameters) # From models
118+
#' terra::plot(v$`Future_2041-2060_ssp585`$from_gcms) # From GCMs
118119

119120

120121
projection_variability <- function(model_projections,
121-
by_replicate = TRUE,
122-
by_gcm = TRUE,
123-
by_model = TRUE,
122+
from_replicates = TRUE,
123+
from_parameters = TRUE,
124+
from_gcms = TRUE,
124125
consensus = "median", # MAKE IT WORK WHEN CONSENSUS IS FULL MODEL
125126
write_files = FALSE,
126127
output_dir = NULL,
@@ -140,13 +141,13 @@ projection_variability <- function(model_projections,
140141
}
141142

142143
if (length(consensus) > 1) {
143-
stop("Argument 'consensus' must be a unique value.",
144-
"\nAvailable options are: 'median', 'range', 'mean' or 'stdev'.")
144+
stop("Argument 'consensus' must be a single value.",
145+
"\nOptions are: 'median' or 'mean'.")
145146
}
146-
consensus_out <- setdiff(consensus, c("median", "range", "mean", "stdev"))
147+
consensus_out <- setdiff(consensus, c("median", "mean"))
147148
if (length(consensus_out) > 0) {
148149
stop("Invalid 'consensus' provided.",
149-
"\nAvailable options are: 'median', 'range', 'mean' or 'stdev'.")
150+
"\nOptions are: 'median' or 'mean'.")
150151
}
151152

152153
if (write_files & is.null(output_dir)) {
@@ -167,7 +168,10 @@ projection_variability <- function(model_projections,
167168
if (write_files) {
168169
out_dir <- file.path(output_dir, "variance")
169170
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
170-
} else {out_dir <- NULL}
171+
} else {
172+
out_dir <- NULL
173+
}
174+
171175
#### Get data ####
172176
d <- model_projections[["paths"]]
173177

@@ -183,9 +187,6 @@ projection_variability <- function(model_projections,
183187

184188
####Iteration over combinations####
185189
res <- lapply(1:nrow(uc), function(z) {
186-
187-
#To test
188-
#z = 1
189190
time <- uc$Time[z]
190191
period <- uc$Period[z]
191192
scenario <- uc$Scenario[z]
@@ -199,7 +200,7 @@ projection_variability <- function(model_projections,
199200
paths <- d_p$output_path
200201

201202
#### By replicate ####
202-
if (by_replicate) {
203+
if (from_replicates) {
203204
if (verbose) {
204205
message("\nCalculating variability from distinct replicates: scenario ",
205206
z, " of ", nrow(uc))
@@ -208,42 +209,51 @@ projection_variability <- function(model_projections,
208209

209210
#### By replicates ####
210211
# Get variance of replicates in each gcm, than get the average across gcms
211-
var_rep_by_gcm <- terra::rast(lapply(paths, var_models_rep_by_gcm))
212-
var_rep <- terra::mean(var_rep_by_gcm)
213-
} else {#End of by_replicate
212+
var_rep <- terra::rast(lapply(paths, var_models_rep_by_gcm))
213+
var_rep <- terra::mean(var_rep)
214+
names(var_rep) <- "from_replicates"
215+
} else {#End of from_replicates
214216
var_rep <- NULL
215217
}
216218

217219
#### By Model ####
218-
if (by_model) {
220+
if (from_parameters) {
219221
if (verbose) {
220222
message("Calculating variability from distinct models: scenario ",
221223
z, " of ", nrow(uc))
222224
}
223225

224226
# Get variance of models in each gcm, than get the average
225-
var_model_by_gcm <- terra::rast(lapply(paths, var_models_model_by_gcm, consensus))
226-
var_model <- terra::mean(var_model_by_gcm)
227-
names(var_model) <- "by_model"
227+
if (names(model_projections$thresholds[[1]])[1] == "Full_model") {
228+
var_model <- terra::rast(lapply(paths, var_models_model_by_gcm,
229+
"Full_model"))
230+
} else {
231+
var_model <- terra::rast(lapply(paths, var_models_model_by_gcm,
232+
consensus))
233+
}
234+
235+
var_model <- terra::mean(var_model)
236+
names(var_model) <- "from_parameters"
228237
} else { #End of by model
229238
var_model <- NULL
230239
}
231240

232241

233242
####By GCM####
234-
if (by_gcm & period != "Present") {
243+
if (from_gcms & period != "Present") {
235244
if (verbose) {
236245
message("Calculating variability from distinct GCMs: scenario ",
237246
z, " of ", nrow(uc))
238247
}
239248

240249
var_gcm <- var_models_across_gcm(paths = paths, consensus = consensus)
241-
names(var_gcm) <- "by_gcm"
250+
names(var_gcm) <- "from_gcms"
242251
} else {
243-
var_gcm <- NULL}#End of by_gcm
252+
var_gcm <- NULL
253+
}#End of from_gcms
244254

245-
all_var <- terra::rast(c("by_replicate" = var_rep, "by_model" = var_model,
246-
"by_gcm" = var_gcm))
255+
all_var <- terra::rast(c("from_replicates" = var_rep, "from_parameters" = var_model,
256+
"from_gcms" = var_gcm))
247257

248258

249259
#Write results

0 commit comments

Comments
 (0)