# Plot the density of predicted values for presences and absences

January 15, 2020
By

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The ‘predDensity’ function, included in the modEvA package since version 1.5 (currently available on R-Forge), produces a histogram and/or a kernel density plot of predicted values for observed presences and absences in a binomial GLM:

```predDensity <- function (model = NULL, obs = NULL, pred = NULL, separate = TRUE,
type = c("both"), legend.pos = "topright")
{
if (!is.null(model)) {
if (!("glm" %in% class(model)) || family(model)\$family !=
"binomial")
stop("'model' must be of class 'glm' and family 'binomial'.")
if (!is.null(obs))
message("Argument 'obs' ignored in favour of 'model'.")
if (!is.null(pred))
message("Argument 'pred' ignored in favour of 'model'.")
obs <- model\$y
pred <- model\$fitted.values
}
if (is.null(obs)) {
if (is.null(pred))
stop("You must provide either 'model' or 'pred'.")
separate <- FALSE
obs <- sample(c(0, 1), length(pred), replace = TRUE)
}
else {
if (length(obs) != length(pred))
stop("'obs' and 'pred' must have the same length.")
}
pred0 <- pred[obs == 0]
pred1 <- pred[obs == 1]
type <- match.arg(type, c("histogram", "density", "both"))
rslt <- vector("list")
if (type %in% c("density", "both")) {
if (!separate) {
dens <- density(pred)
xrange <- range(dens\$x, finite = TRUE)
yrange <- range(dens\$y, finite = TRUE)
rslt[["density"]] <- dens
}
else {
dens0 <- density(pred0)
dens1 <- density(pred1)
xrange <- range(dens0\$x, dens1\$x, finite = TRUE)
yrange <- range(dens0\$y, dens1\$y, finite = TRUE)
rslt[["density_obs1"]] <- dens1
rslt[["density_obs0"]] <- dens0
}
plot(x = xrange, y = yrange, xlab = "Predicted value",
ylab = "Density", type = "n")
}
if (type %in% c("histogram", "both")) {
hist0 <- hist(pred0, plot = FALSE)
hist1 <- hist(pred1, plot = FALSE)
if (type == "histogram") {
yrange <- range(hist0\$density, hist1\$density, finite = TRUE)
plot(x = c(0, 1), y = yrange, type = "n", xlab = "Predicted value",
ylab = "Density")
}
if (!separate) {
histogram <- hist(c(pred0, pred1), freq = FALSE,
col = "grey20", add = TRUE)
rslt[["histogram"]] <- histogram
}
else {
hist(pred1, freq = FALSE, col = "grey20", add = TRUE)
hist(pred0, freq = FALSE, col = "darkgrey", density = 40,
angle = 45, add = TRUE)
rslt[["histogram_obs1"]] <- hist1
rslt[["histogram_obs0"]] <- hist0
if (legend.pos != "n" && type == "histogram")
legend(legend.pos, legend = c("absences", "presences"),
fill = c("darkgrey", "grey20"), border = NA,
density = c(40, NA), bty = "n")
}
}
if (type %in% c("density", "both")) {
if (!separate) {
lines(dens, col = "black", lwd = 2)
}
else {
lines(dens1, col = "black", lwd = 2)
lines(dens0, col = "darkgrey", lty = 5, lwd = 2)
if (legend.pos != "n" && type == "density")
legend(legend.pos, legend = c("absences", "presences"),
col = c("darkgrey", "black"), lty = c(5, 1),
bty = "n")
if (legend.pos != "n" && type == "both")
legend(legend.pos, legend = c("absences", "presences"),
fill = c("darkgrey", "grey20"), border = NA,
lty = c(5, 1), col = c("darkgrey", "grey15"),
density = c(40, NA), bty = "n")
}
}
return(rslt)
}
```

Usage example:

```install.packages("modEvA", repos="http://R-Forge.R-project.org")
library(modEvA)
data(rotif.mods)
predDensity(model = rotif.mods\$models[[3]])
```

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.