Nothing
#' @title Plot cell-specific decay curve for carcass persistence
#'
#' @description Produce the figure panel for a specific cell (factor
#' level combination) including the specific fitted decay curves.
#'
#' @param model model of class cpm
#'
#' @param specificCell name of the specific cell to plot
#'
#' @param col color to use
#'
#' @param axis_x logical of whether or not to plot the x axis
#'
#' @param axis_y logical of whether or not to plot the y axis
#'
#' @export
#'
cpmCPCellPlot <- function(model, specificCell, col, axis_y = TRUE,
axis_x = TRUE){
# Juniper's original function
dist <- model$dist
CL <- model$CL
cellwise <- model$cell_ab
cellNames <- model$cells[ , "CellNames"]
whichCarcs <- which(model$carcCell == specificCell)
observations <- model$observations[whichCarcs, ]
ncarc <- nrow(observations)
whichSpecificCell <- which(cellNames == specificCell)
a <- cellwise[whichSpecificCell, "pda_median"]
b <- cellwise[whichSpecificCell, "pdb_median"]
abs <- rcp(n = 1000, model = model, type = "ppersist")
as <- abs[[whichSpecificCell]][ , "pda"]
bs <- abs[[whichSpecificCell]][ , "pdb"]
xVals <- model$observations[model$observations!= Inf]
max_x <- max(xVals, na.rm = TRUE)
max_x <- ceiling(max_x / 10) * 10
pred_x <- seq(0, max_x, length.out = max_x * 10)
pts <- pred_x[2:length(pred_x)]
pta0 <- rep(0, length(pts))
pta1 <- rep(0.000001, length(pts))
CP <- ppersist(pda = as, pdb = bs, dist = dist,
t_arrive0 = pta0, t_arrive1 = pta1, t_search = pts)
pred_y <- apply(CP, 1, median)
pred_yl <- apply(CP, 1, quantile, probs = (1 - CL) / 2)
pred_yu <- apply(CP, 1, quantile, probs = 1 - (1 - CL) / 2)
t1 <- observations[ , 1]
t2 <- observations[ , 2]
event <- rep(3, length(t1))
event[which(!is.finite(t2))] <- 0
event[which(t1 == t2)] <- 1
t1[which(t1 == 0)] <- 0.0001
survobj <- survival::Surv(t1, t2, event, "interval")
form <- formula("survobj ~ 1")
smod <- tryCatch(
survival::survfit(form, data = observations),
error = function(e) NULL,
warning = function(w) NULL
)
plot(smod, ylim = c(0, 1), xlim = c(0, max_x), # this gives the persistence probabilities
xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "o", lwd = c(2, 1, 1)
)
axis(1, las = 1, cex.axis = 0.9, at = seq(0, max_x, by = 10), labels = axis_x)
axis(2, las = 1, cex.axis = 0.9, at = seq(0, 1, 0.2), labels = axis_y)
points(pts, pred_y, type = "l", col = col, lwd = 3)
points(pts, pred_yl, type = "l", col = col, lwd = 2, lty = 3)
points(pts, pred_yu, type = "l", col = col, lwd = 2, lty = 3)
countText <- paste("N = ", ncarc, sep = "")
text(x = max_x, y = 0.95, countText, cex = 0.75, xpd = TRUE, adj = 1)
text(max_x, 1.02, specificCell, adj = 1, cex = 0.75, font = 2)
}
#' @title Plot results of a single CP model
#'
#' @description Plot a single \code{\link{cpm}} model
#'
#' @param x model of class cpm
#'
#' @param col color to use
#'
#' @param ... to be passed down
#'
#' @examples
#' data(wind_RP)
#' mod <- cpm(formula_l = l ~ Season, formula_s = s ~ Season,
#' data = wind_RP$CP, left = "LastPresent", right = "FirstAbsent")
#' plot(mod)
#'
#' @export
#'
plot.cpm <- function(x, col = "black", ...){
# loop through cells
# create header
model <- x
name_l <- format(model$formula_l)
name_s <- format(model$formula_s)
modelName <- paste(name_l, "; ", name_s, sep = "")
if (model$dist == "exponential"){
modelName <- name_l
}
par(mar = c(0, 0, 0, 0), fig = c(0, 1, 0.95, 1))
plot(1, 1, type = 'n', bty = 'n', xaxt = 'n', yaxt = 'n', xlab = "",
ylab = "", ylim = c(0, 1), xlim = c(0, 1))
points(c(0.01, 0.06), c(0.25, 0.25), type = 'l', lwd = 2, col = col)
text(x = 0.07, y = 0.3, "= Median", adj = 0, cex = 0.9)
points(c(0.2, 0.25), c(0.25, 0.25), type = 'l', lwd = 2, lty = 3, col = col)
text(x = 0.26, y = 0.3, "= Confidence Bounds", adj = 0)
labelsText <- paste(model$predictors, collapse = ".")
text_label <- paste("Labels: ", labelsText, sep = "")
text_model <- paste("Model: ", modelName, sep = "")
text(x = 0.58, y = 0.3, text_label, adj = 0, cex = 0.75)
text(x = 0.58, y = 0.7, text_model, adj = 0, cex = 0.75)
par(fig = c(0, 1, 0, 1), mar = c(1, 1, 1, 1), new = TRUE)
plot(1,1, type = 'n', bty = 'n', xaxt = 'n', yaxt = 'n', xlab = "",
ylab = ""
)
mtext(side = 1, "Time", line = -0.25, cex = 1.5)
mtext(side = 2, "Carcass Persistence", line = -0.25, cex = 1.5)
ncell <- model$ncell
cellNames <- model$cells[ , "CellNames"]
nmatrix_col <- min(3, ncell)
nmatrix_row <- ceiling(ncell / nmatrix_col)
figxspace <- 0.925 / nmatrix_col
figyspace <- 0.875 / nmatrix_row
x1 <- rep(figxspace * ((1:nmatrix_col) - 1), nmatrix_row) + 0.05
x2 <- rep(figxspace * ((1:nmatrix_col)), nmatrix_row) + 0.05
y1 <- rep(figyspace * ((nmatrix_row:1) - 1), each = nmatrix_col) + 0.04
y2 <- rep(figyspace * ((nmatrix_row:1)), each = nmatrix_col) + 0.04
bottomCells <- seq(ncell - (nmatrix_col - 1), ncell, 1)
leftCells <- which(1:ncell %% nmatrix_col == 1)
if (length(leftCells) == 0){
leftCells <- 1
}
for (celli in 1:ncell){
par(mar = c(2.5, 2, 0, 0))
par(fig = c(x1[celli], x2[celli], y1[celli], y2[celli]), new = T)
specificCell <- cellNames[celli]
axis_x <- FALSE
axis_y <- FALSE
if (celli %in% bottomCells){
axis_x <- TRUE
}
if (celli %in% leftCells){
axis_y <- TRUE
}
cpmCPCellPlot(model, specificCell, col, axis_y, axis_x)
}
par(.par_default)
}
#' @title Plot results of a set of CP models
#'
#' @description Produce a set of figures for a set of CP models, as fit by
#' \code{\link{cpmSet}}
#'
#' @param x pk model set of class pkmSet
#'
#' @param specificModel the name(s) or index number(s) of specific model(s) to
#' restrict the plot
#'
#' @param cols named vector of the colors to use for the distributions
#'
#' @param ... to be passed down
#'
#' @examples
#' data(wind_RP)
#' mod <- cpmSet(formula_l = l ~ Season, formula_s = s ~ Season,
#' data = wind_RP$CP, left = "LastPresent", right = "FirstAbsent")
#' \donttest{plot(mod)}
#'
#' @export
#'
plot.cpmSet <- function(x, specificModel = NULL, cols = NULL, ...){
modelSet <- tidyModelSetCP(x)
cells_set <- modelSetCells(modelSet)
ncell <- nrow(cells_set)
preds_set <- modelSetPredictors(modelSet)
if (!is.null(specificModel) &&
(!"character" %in% class(specificModel) || anyNA(specificModel))){
stop("specific model must be specified by name (character string) or number")
}
if (is.null(specificModel)){
specMods <- names(modelSet)
} else {
for (spi in specificModel){
if (is.null(modelSet[[spi]])){ # spi not in modelSet (maybe exponential?)
if (!grepl("exponential", spi))
stop("specific CP model not in model set")
tmp_spc <- strsplit(spi, ";")[[1]]
tmp_spc[3] <- " NULL"
if (is.null(modelSet[paste(tmp_spc, collapse = ";")]))
stop("specific CP model not in model set")
}
}
specMods <- specificModel
}
dist_set <- NULL # which distributions are included?
# the possibilities are limited to names in the .cols_CP vector (internal data)
for (di in names(.cols_CP))
if(any(grepl(di, names(modelSet)))) dist_set <- c(dist_set, di)
if (is.null(cols)){
cols_CP <- .cols_CP[dist_set]
} else if (length(dist_set) == 1 & length(cols) == 1){
cols_CP <- cols
names(cols_CP) <- dist_set
} else {
if (!all(dist_set %in% names(cols))){
stop("'cols' must be NULL or a named vector of colors to be used in ",
"plotting the distributions. Names are the names of the distributions ",
"to associate with each respective color in the vector.")
}
cols_CP <- .cols_CP[dist_set]
}
preds_set <- modelSetPredictors(modelSet)
n_col <- ifelse(length(preds_set) == 0, 1,
length(unique(cells_set[ , preds_set[1]])))
n_row <- nrow(cells_set)/n_col
for (spi in specMods){
if (spi != specMods[1]){ # more than one specific & more than one fig to draw
devAskNewPage(TRUE)
}
par(oma = c(5, 7.5, ifelse(nrow(cells_set) > 1, 9, 6), ifelse(n_row > 1, 4, 0.5)),
mar = c(0, 0, 0, 0))
par(mfrow = c(n_row, n_col))
CPFig(modelSet, spi, cells_set, preds_set, cols_CP)
}
devAskNewPage(FALSE)
par(.par_default)
}
# draw a CP figure graph
# NOTE: This is an internal utility function called by plot.cpmSet
#' @keywords internal
CPFig <- function(modelSet, specificModel, cells_set, preds_set, cols_CP){
n_col <- ifelse(length(preds_set) == 0, 1,
length(unique(cells_set[ , preds_set[1]])))
n_row <- nrow(cells_set)/n_col
distr_featured <- sub(";", "", strsplit(specificModel, " ")[[1]][2])
if (distr_featured == "exponential" && grepl("NULL", specificModel)){
modelSet <- modelSet[grepl("exponential", names(modelSet))]
cols_CP <- cols_CP["exponential"]
}
if (distr_featured == "exponential" && length(cols_CP) == 1){
dists <- "exponential"
forms <- strsplit(specificModel, "; ")[[1]]
exind <- grep(paste0(forms[2], ";"), names(modelSet), fixed = TRUE)
model_spc <- modelSet[[exind]]
names(modelSet)[[exind]] <- specificModel
} else {
# reorder distributions to plot featured distr last
dists <- c(names(cols_CP), distr_featured)
dists <- dists[-min(which(dists == distr_featured))]
if (("exponential" %in% dists) & length(dists) > 1){
# change "NULL" to "s ~ <whatever's needed>"
forms <- strsplit(specificModel, "; ")[[1]]
exind <- intersect(
grep("exponential", names(modelSet)),
grep(paste0(forms[2], ";"), names(modelSet), fixed = TRUE))
names(modelSet)[exind] <-
sub("NULL", strsplit(specificModel, "; ")[[1]][3], names(modelSet)[exind])
}
if (distr_featured == "exponential" & length(dists > 1)){
specificModel <- sub("exponential", dists[1], specificModel)
}
model_spc <- modelSet[[specificModel]]
}
# draw the cells
for (celli in 1:nrow(cells_set)){
specificCell <- cells_set[celli, "CellNames"]
if (nrow(cells_set) == 1){
observations <- model_spc$observations
} else {
carcCells <- apply(data.frame(model_spc$data[ , preds_set]),
MARGIN = 1, FUN = paste, collapse = ".")
observations <- model_spc$observations[carcCells == specificCell, ]
}
ncarc <- nrow(observations)
cind <- which(cells_set$CellNames == specificCell)
cell_spc <- matchCells(model_spc, modelSet)[cind]
if ("exponential" %in% dists)
cell_exp <- matchCells(modelSet[[exind]], modelSet)[cind]
xVals <- observations[observations != Inf]
max_x <- max(xVals, na.rm = TRUE)
max_x <- ceiling(max_x / 10) * 10
t1 <- observations[ , 1]
t2 <- observations[ , 2]
event <- rep(3, length(t1))
event[which(!is.finite(t2))] <- 0
event[which(t1 == t2)] <- 1
t1[which(t1 == 0)] <- 0.0001
smod <- tryCatch(
survival::survfit(survival::Surv(t1, t2, event, type = "interval") ~ 1),
error = function(e) NULL,
warning = function(w) NULL
)
plot(smod, ylim = c(0, 1), xlim = c(0, max_x), xlab = "", ylab = "",
axes = FALSE, lwd = c(2, 1, 1))
box()
if (celli <= n_col & length(preds_set) > 0)
mtext(side = 3, line = 0.3, text = cells_set[celli, preds_set[1]])
if (celli %% n_col == 0 && length(preds_set) == 2)
mtext(side = 4, line = 0.3, text = cells_set[celli, preds_set[2]])
if (celli %% n_col == 1 | n_col == 1){
axis(2, las = 1, at = seq(0, 1, by = 0.2))
axis(2, las = 1, at = seq(0, 1, by = 0.1), labels = FALSE, tcl = -0.3)
}
if (celli > n_col * (n_row - 1))
axis(1)
xvals <- seq(0.001, max_x, length = 500)
text(x = max_x, y = 0.95, paste0("N = ", ncarc), adj = 1)
for (distr in dists){
lwd <- 2 + 2 * (distr == distr_featured)
model_line <- gsub(model_spc$dist, distr, specificModel)
mod <- modelSet[[model_line]]
if (!is.null(mod)){
if (distr == "exponential"){
crow <- which(mod$cell_ab$cell == cell_exp)
lines(xvals, 1 - pexp(xvals, rate = 1/mod$cell_ab[crow, "pdb_median"]),
col = cols_CP[distr], lwd = lwd)
} else {
cind <- which(mod$cell_ab$cell == cell_spc)
pda <- mod$cell_ab[cind, "pda_median"]
pdb <- mod$cell_ab[cind, "pdb_median"]
if (distr == "weibull"){
lines(xvals, 1 - pweibull(xvals, shape = pda, scale = pdb),
lwd = lwd, col = cols_CP[distr])
} else if (distr == "loglogistic"){
lines(xvals, 1 - pllogis(xvals, pda = pda, pdb = pdb),
lwd = lwd, col = cols_CP[distr])
} else if (distr == "lognormal"){
lines(xvals, 1 - plnorm(xvals, meanlog = pdb, sdlog = sqrt(pda)),
lwd = lwd, col = cols_CP[distr])
}
}
}
}
}
mtext(side = 1, line = 3.3, text = "Time", outer = TRUE, cex = 1.4)
mtext(side = 2, line = 5.2, text = "Carcass Persistence", outer = TRUE, cex = 1.4)
mtext(side = 2, line = 3.3, text = "Pr(persist > t)", cex = 1.2, outer = TRUE,)
if (length(preds_set) > 0)
mtext(side = 3, line = 2.5, text = preds_set[1], cex = 1.2, outer = TRUE)
if (length(preds_set) == 2)
mtext(side = 4, line = 2.5, text = preds_set[2], cex = 1.2, outer = TRUE)
par(new = TRUE, mar = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mfrow = c(1, 1) )
plot(0, type = "n", axes = FALSE)
if (length(dists) > 1)
dists <- dists[c(length(dists), 1:(length(dists) - 1))]
lwd <- rep(2, length(dists))
lwd[1] <- 4
if (grepl("exponential", specificModel) & length(dists) == 1){
# exponential model only is plotted
legend("topleft",
xpd = TRUE,
legend = "exponential", lty = 1, lwd = lwd, ncol = 2, bty = "n",
cex = 0.85, col = cols_CP["exponential"],
title = paste("model:",
gsub("~ 1", "constant", paste(strsplit(specificModel, ";")[[1]][2],
collapse = ";"))),
title.adj = 0)
} else {
legend("topleft", xpd = TRUE, legend = dists, lty = 1, col = cols_CP[dists],
lwd = lwd, ncol = 2, bty = "n", cex = 0.85,
title = paste("model:",
gsub("~ 1", "~ constant",
paste(strsplit(specificModel, ";")[[1]][-1], collapse = ";"))),
title.adj = 0)
}
}
#' @title Plot results of a single CP model in a set
#'
#' @description Produce a figures for a specific CP model, as fit by
#' \code{\link{cpmSet}}
#'
#' @param modelSet cp model set of class \code{cpmSet}
#'
#' @param specificModel the name of the specific model for the plot
#'
#' @param cols named vector of the colors to use for the distributions
#'
#' @return a plot
#'
#' @export
#'
plotCPFigure <- function(modelSet, specificModel, cols = CPcols()){
plotCPHeader(modelSet, specificModel, cols)
plotCPCells(modelSet, specificModel, cols)
}
#' @title Plot the cellwise results of a single model in a set of CP models
#'
#' @description Produce a set of cellwise figures for a specific CP model, as
#' fit by \code{\link{cpmSet}}
#'
#' @param modelSet cp model set of class cpmSet
#'
#' @param specificModel the name of the specific model for the plot
#'
#' @param cols named vector of the colors to use for the distributions
#'
#' @return a plot
#'
#' @export
#'
plotCPCells <- function(modelSet, specificModel, cols){
## need cells_set, preds_set
preds_set <- modelSetPredictors(modelSet)
cells_set <- modelSetCells(modelSet)
par(fig = c(0, 1, 0, 0.925), new = TRUE, mar = c(0, 0, 0, 0),
oma = c(4, 4, 4, 0.5))
plot(1, 1, type = "n", bty = "n", xaxt = "n", yaxt = "n", xlab = "",
ylab = "")
mtext(side = 1, "Time", line = -0.25, cex = 1.5)
mtext(side = 2, "Carcass Persistence", line = -0.25, cex = 1.5)
n_col <- length(unique(cells_set[ , preds_set[1]]))
n_row <- nrow(cells_set)/n_col
H <- .res * (.header + n_row * .panel_H)
W <- min(1200, .res * .panel_W * max(2, n_col))
cell_W <- round(W/n_col)
box_W <- round(W/2)
# cells
if (nrow(cells_set) > 1) omd = c(.res/W, 1 - .res/W,
.res * .footer/H, 1 - .res * (.header + .box_H + 0.75 * .buffer)/H)
if (nrow(cells_set) == 1) omd = c(.res/W, 0.96,
.res * .footer/H, 1 - .res * (.header + .box_H + 0.2 * .buffer)/H)
par(omd = omd, mar = c(0, 0, 0, 0), xaxt = "n", yaxt = "n")
par(mfrow = c(n_row, n_col))
cells_set <- modelSetCells(modelSet)
ncell <- nrow(cells_set)
cellNames <- cells_set[ , "CellNames"]
nmatrix_col <- min(3, ncell)
nmatrix_row <- ceiling(ncell / nmatrix_col)
figxspace <- 0.925 / nmatrix_col
figyspace <- 0.875 / nmatrix_row
x1 <- rep(figxspace * ((1:nmatrix_col) - 1), nmatrix_row) + 0.04
x2 <- rep(figxspace * ((1:nmatrix_col)), nmatrix_row) + 0.04
y1 <- rep(figyspace * ((nmatrix_row:1) - 1), each = nmatrix_col) + 0.03
y2 <- rep(figyspace * ((nmatrix_row:1)), each = nmatrix_col) + 0.03
bottomCells <- seq(ncell - (nmatrix_col - 1), ncell, 1)
leftCells <- which(1:ncell %% nmatrix_col == 1)
if (length(leftCells) == 0){
leftCells <- 1
}
for (celli in 1:ncell){
par(mar = c(2.5, 2, 0, 0))
par(fig = c(x1[celli], x2[celli], y1[celli], y2[celli]), new = TRUE)
specificCell <- cellNames[celli]
axis_x <- FALSE
axis_y <- FALSE
if (celli %in% bottomCells){
axis_x <- TRUE
}
if (celli %in% leftCells){
axis_y <- TRUE
}
axes <- c("x" = axis_x, "y" = axis_y)
cpmSetSpecCPCellPlot(modelSet, specificModel, specificCell, cols, axes)
}
}
#' @title Plot cell-specific decay curve for carcass persistence
#'
#' @description Produce the figure panel for a specific cell (factor
#' level combination) including the specific fitted decay curves.
#'
#' @param modelSet modelSet of class cpmSet
#'
#' @param specificModel name of the specific submodel to plot
#'
#' @param specificCell name of the specific cell to plot
#'
#' @param cols named vector of the colors to use for the distributions
#'
#' @param axes named vector of logical values indicating whether or not to
#' plot the x axis and the y axis
#'
#' @return a specific cell plot panel
#'
#' @export
#'
cpmSetSpecCPCellPlot <- function(modelSet, specificModel, specificCell,
cols, axes){
model_spec <- modelSet[[specificModel]]
dist <- model_spec$dist
CL <- model_spec$CL
cellNames_set <- modelSetCells(modelSet)[ , "CellNames"]
preds_set <- modelSetPredictors(modelSet)
carcCells <- apply(data.frame(model_spec$data[ , preds_set]),
MARGIN = 1, FUN = paste, collapse = ".")
whichCarcs <- which(carcCells == specificCell)
if (specificCell == "all"){
whichCarcs <- 1:length(carcCells)
}
observations <- model_spec$observations[whichCarcs, ]
ncarc <- nrow(observations)
cellMatch_spec <- matchCells(model_spec, modelSet)
whichSpecificCell_spec <- cellMatch_spec[cellNames_set == specificCell]
xVals <- model_spec$observations[model_spec$observations!=Inf]
max_x <- max(xVals, na.rm = TRUE)
max_x <- ceiling(max_x / 10) * 10
t1 <- observations[ , 1]
t2 <- observations[ , 2]
event <- rep(3, length(t1))
event[which(!is.finite(t2))] <- 0
event[which(t1 == t2)] <- 1
t1[which(t1 == 0)] <- 0.0001
smod <- tryCatch(
survival::survfit(survival::Surv(t1, t2, event, type = "interval") ~ 1),
error = function(e) NULL,
warning = function(w) NULL
)
plot(smod, ylim = c(0, 1), xlim = c(0, max_x),
xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "o", lwd = c(2, 1, 1)
)
axis(1, las = 1, cex.axis = 0.9, at = seq(0, max_x, by = 10),
labels = axes["x"]
)
axis(2, las = 1, cex.axis = 0.9, at = seq(0, 1, 0.2), labels = axes["y"])
countText <- paste("N = ", ncarc, sep = "")
text(x = max_x, y = 0.95, countText, cex = 0.75, xpd = TRUE, adj = 1)
text(max_x, 1.02, specificCell, adj = 1, cex = 0.75, font = 2)
nmodelsInSet <- length(modelSet)
modelSetNames <- names(modelSet)
modelSetNames_components <- unlist(strsplit(modelSetNames, "; "))
modelSetNames_matrix <- matrix(modelSetNames_components, ncol = 3, byrow = TRUE)
colnames(modelSetNames_matrix) <- c("dist", "form_l", "form_s")
modelMatches <- vector("list", length = nmodelsInSet)
for (modelsInSeti in 1:nmodelsInSet){
mod_l <- gsub(" ", "", modelSetNames_matrix[modelsInSeti, 2])
mod_s <- gsub(" ", "", modelSetNames_matrix[modelsInSeti, 3])
modComp_l <- gsub(" ", "", modelSetNames_matrix[ , 2])
modComp_s <- gsub(" ", "", modelSetNames_matrix[ , 3])
matches <- which(modComp_l == mod_l & modComp_s == mod_s)
modelMatches[[modelsInSeti]] <- matches
}
whichSpecificModel <- which(names(modelSet) == specificModel)
modsToPlot <- modelMatches[[whichSpecificModel]]
modsToPlot <- modsToPlot[-which(modsToPlot == whichSpecificModel)]
modsToPlot <- c(modsToPlot, whichSpecificModel)
nmodsToPlot <- length(modsToPlot)
for (modsToPloti in 1:nmodsToPlot){
model_i <- modelSet[[modsToPlot[modsToPloti]]]
dist_i <- model_i$dist
cellwise_i <- model_i$cell_ab
cellNames_i <- cellwise_i[ , "cell"]
cellMatch_i <- matchCells(model_i, modelSet)
reducedCell_i <- cellMatch_i[cellNames_set == specificCell]
whichSpecificCell <- which(cellNames_i == reducedCell_i)
a <- cellwise_i[whichSpecificCell, "pda_median"]
b <- cellwise_i[whichSpecificCell, "pdb_median"]
pred_x <- seq(0, max_x, length.out = max_x * 10)
pts <- pred_x[2:length(pred_x)]
pta0 <- rep(0, length(pts))
pta1 <- rep(0.000001, length(pts))
pred_y <- t(ppersist(pda = a, pdb = b, dist = dist_i,
t_arrive0 = pta0, t_arrive1 = pta1, t_search = pts))
col_i <- cols[dist_i]
lwd <- 2
if (modsToPloti == nmodsToPlot){
lwd <- 5
}
points(pts, pred_y, type = "l", col = col_i, lwd = lwd)
}
}
#' @title The CP plot header
#'
#' @description Produce the header for a CP plot
#'
#' @param modelSet cp model set of class cpmSet
#'
#' @param specificModel the name of the specific model for the plot
#'
#' @param cols named vector of the colors to use for the distributions
#'
#' @return a plot
#'
#' @export
#'
plotCPHeader <- function(modelSet, specificModel, cols = CPcols()){
modelSetNames <- names(modelSet)
nmodelsInSet <- length(modelSetNames)
whichSpecificModel <- which(names(modelSet) == specificModel)
modelSetNames_components <- unlist(strsplit(modelSetNames, "; "))
modelSetNames_matrix <- matrix(modelSetNames_components, ncol = 3,
byrow = TRUE)
colnames(modelSetNames_matrix) <- c("dist", "form_l", "form_s")
distsIncluded <- gsub("dist: ", "", unique(modelSetNames_matrix[ , "dist"]))
par(fig = c(0, 1, 0.925, 1), mar = c(0, 0, 0, 0))
plot(1,1, type = "n", bty = "n", xaxt = "n", yaxt = "n", xlab = "",
ylab = "", ylim = c(0, 1), xlim = c(0, 1))
x1 <- c(0, 0.29, 0, 0.29)
y1 <- c(0.4, 0.4, 0.0, 0.0)
x2 <- c(0.04, 0.33, 0.04, 0.33)
y2 <- c(0.6, 0.6, 0.2, 0.2)
xt <- c(0.05, 0.34, 0.05, 0.34)
yt <- c(0.5, 0.5, 0.1, 0.1)
ndists <- length(distsIncluded)
for (disti in 1:ndists){
rect(x1[disti], y1[disti], x2[disti], y2[disti], border = NA,
col = cols[distsIncluded[disti]]
)
distName <- distsIncluded[disti]
distName <- paste(toupper(substring(distName, 1, 1)),
substring(distName, 2), sep = "", collapse = " ")
distText <- paste0("= ", gsub("Log", "Log-", distName))
text(x = xt[disti], y = yt[disti], adj = 0, cex = 0.8, distText)
}
labelsText <- paste(modelSetPredictors(modelSet), collapse = ".")
text_label <- paste("Labels: ", labelsText, sep = "")
forms <- modelSetNames_matrix[whichSpecificModel, c("form_l", "form_s")]
modelsText <- paste(forms, collapse = "; ")
modelsText <- gsub("~ 1", "~ constant", modelsText)
text_model <- paste("Model: ", modelsText, sep = "")
text(x = 0.59, y = 0.1, text_label, adj = 0, cex = 0.75)
text(x = 0.59, y = 0.5, text_model, adj = 0, cex = 0.75)
}
#' @title Error check a specific model selection for a CP plot
#'
#' @description Make sure it's available and good, update the name for usage
#'
#' @param modelSet cp model set of class cpmSet
#'
#' @param specificModel the name of the specific model for the plot
#'
#' @return updated name of the model to use
#'
#' @export
#'
checkSpecificModelCP <- function(modelSet, specificModel){
if (!is.null(specificModel) && anyNA(specificModel)){
stop(
"specificModel must be NULL or a vector of model names or positions.",
"\nNAs not allowed."
)
}
if (length(specificModel) > 0){
if (is.numeric(specificModel)){
if (anyNA(specificModel)){
warning("specificModel cannot be NA. NA models removed.")
specificModel <- specificModel[!is.na(specificModel)]
if (length(specificModel) == 0){
stop("No valid specificModel")
}
}
if (any(specificModel > length(modelSet))){
stop(paste0("there are only ", length(modelSet), " model choices."))
}
specificModel <- names(modelSet)[specificModel]
} else{
specificModel <- gsub("NULL", "s ~ 1", specificModel)
}
if (any(specificModel %in% names(modelSet)) == FALSE){
stop("Selected model not in set. To see options use names(modelSet).")
}
modNames <- specificModel
for (modi in modNames){
if (cpmFail(modelSet[[modi]])){
stop("specificModel ", modi, " is not a well-fit model")
}
}
} else{
specificModel <- names(cpmSetFailRemove(modelSet))
}
return(specificModel)
}
#' @title Expand a CP model set for plotting
#'
#' @description Expand the exponential models across the other distributions
#'
#' @param modelSet cp model set of class cpmSet
#'
#' @return updated model set
#'
#' @export
#'
expandModelSetCP <- function(modelSet){
modelSetNames <- names(modelSet)
nmodelsInSet <- length(modelSetNames)
modelSetNames_components <- unlist(strsplit(modelSetNames, "; "))
modelSetNames_matrix <- matrix(modelSetNames_components, ncol = 3,
byrow = TRUE
)
colnames(modelSetNames_matrix) <- c("dist", "form_l", "form_s")
whichExp <- which(modelSetNames_matrix[ , "dist"] == "dist: exponential")
modelSetNames_matrix[whichExp, "form_s"] <- "s ~ 1"
if (length(whichExp) > 0){
uniqueForm_s <- unique(modelSetNames_matrix[ , "form_s"])
names(modelSet)[whichExp] <- gsub("NULL", "s ~ 1",
names(modelSet)[whichExp]
)
if (length(uniqueForm_s) > 1){
uniqueForm_s <- uniqueForm_s[-which(uniqueForm_s == "s ~ 1")]
nuniqueForm_s <- length(uniqueForm_s)
for (uniqueForm_si in 1:nuniqueForm_s){
replacement <- uniqueForm_s[uniqueForm_si]
newEntries <- modelSetNames_matrix[whichExp, ]
newEntries[ , "form_s"] <- replacement
modelSetNames_matrix <- rbind(modelSetNames_matrix, newEntries)
newEntries <- modelSet[whichExp]
names(newEntries) <- gsub("s ~ 1", replacement, names(newEntries))
modelSet <- c(modelSet, newEntries)
}
}
}
return(modelSet)
}
#' @title Tidy a CP model set
#'
#' @description Remove bad fit models
#'
#' @param modelSet cp model set of class cpmSet
#'
#' @return a trimmed model set
#'
#' @export
#'
tidyModelSetCP <- function(modelSet){
modelSet <- cpmSetFailRemove(modelSet)
modelSet <- modelSet[order(sapply(modelSet, "[[", "AICc"))]
class(modelSet) <- c("cpmSet", "list")
return(modelSet)
}
#' @title Produce a named vector of standard CP plot colors
#'
#' @description Produce a named vector of standard CP plot colors
#'
#' @export
#'
CPcols <- function(){
c(exponential = rgb(0.80, 0.38, 0.56),
weibull = rgb(1.00, 0.76, 0.15),
loglogistic = rgb(0.00, 1.00, 1.00),
lognormal = rgb(0.00, 0.41, 0.55))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.