diff options
author | Ken Kellner <kenkellner@users.noreply.github.com> | 2023-03-07 13:22:29 -0500 |
---|---|---|
committer | GitHub <noreply@github.com> | 2023-03-07 13:22:29 -0500 |
commit | 8c17748fdea81f61be92a87fdb0a13963bc67e8e (patch) | |
tree | be6907e2ad3702db5ac173de76ffc8177a730216 | |
parent | 8bf6bf43ca45808a5839e091534c6b4a4b494741 (diff) | |
parent | 304318512db5c6ad256324f2e506d4c90adba77e (diff) |
Merge pull request #247 from kenkellner/plotMarginal
Add plotEffects function
-rw-r--r-- | DESCRIPTION | 5 | ||||
-rw-r--r-- | NAMESPACE | 4 | ||||
-rw-r--r-- | R/plotEffects.R | 109 | ||||
-rw-r--r-- | man/plotEffects.Rd | 100 | ||||
-rw-r--r-- | tests/testthat/test_gdistremoval.R | 6 | ||||
-rw-r--r-- | tests/testthat/test_occuTTD.R | 6 | ||||
-rw-r--r-- | tests/testthat/test_plotEffects.R | 48 |
7 files changed, 274 insertions, 4 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 5454b14..5a04732 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.2.5.9010 -Date: 2022-12-29 +Version: 1.2.5.9011 +Date: 2023-03-07 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( @@ -46,6 +46,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R' 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R' 'nmixTTD.R' 'gdistremoval.R' + 'plotEffects.R' 'mixedModelTools.R' 'power.R' 'simulate.R' @@ -7,7 +7,7 @@ importFrom(stats, confint, fitted, coef, vcov, predict, update, profile, pnorm, qchisq, qnorm, quantile, rbinom, reshape, rmultinom, rnbinom, rpois, runif, sd, uniroot, update.formula, sigma) -importFrom(graphics, plot, hist, abline) +importFrom(graphics, plot, hist, abline, axis, lines, points, polygon, segments) importFrom(utils, head, read.csv) importFrom(grDevices, devAskNewPage, dev.interactive, palette.colors) importFrom(MASS, mvrnorm) @@ -51,7 +51,7 @@ exportMethods(backTransform, coef, confint, coordinates, fitted, getData, "siteCovs<-", summary, update, vcov, yearlySiteCovs, "yearlySiteCovs<-", "[", smoothed, projected, nonparboot, logLik, LRT, ranef, bup, crossVal, posteriorSamples, sigma, randomTerms, - optimizePenalty, unmarkedPowerList) + optimizePenalty, unmarkedPowerList, plotEffectsData, plotEffects) S3method("print", "unmarkedPostSamples") diff --git a/R/plotEffects.R b/R/plotEffects.R new file mode 100644 index 0000000..c77f5aa --- /dev/null +++ b/R/plotEffects.R @@ -0,0 +1,109 @@ +#type_to_covs <- function(umf, type){ +# if(type %in% c("state","psi","lam","lambda","sigma","dist")){ +# return(methods::slot(umf, "siteCovs")) +# } else if(type %in% c("det","rem","fp","b")){ +# return(methods::slot(umf, "obsCovs")) +# } else if(type %in% c("phi","transition","col","ext","gamma","omega","iota")){ +# return(methods::slot(umf, "yearlySiteCovs")) +# } +# return(NULL) +#} + +get_base_newdata <- function(object, type){ + covs <- get_orig_data(object, type) + out <- lapply(covs, function(x){ + if(is.numeric(x)){ + return(median(x, na.rm=TRUE)) + } else if(is.factor(x)){ + return(factor(levels(x)[1], levels=levels(x))) + } else { + stop("Unknown column type") + } + }) + as.data.frame(out) +} + +get_cov_seq <- function(covariate, object, type){ + cov_values <- get_orig_data(object, type)[[covariate]] + if(is.numeric(cov_values)){ + rng <- range(cov_values, na.rm=TRUE) + return(seq(rng[1], rng[2], length.out=100)) + } else if(is.factor(cov_values)){ + return(factor(levels(cov_values), levels=levels(cov_values))) + } else { + stop("Unknown covariate type") + } +} + +setGeneric("plotEffectsData", function(object, ...) standardGeneric("plotEffectsData")) + +setMethod("plotEffectsData", "unmarkedFit", + function(object, type, covariate, level=0.95, ...){ + + #umf <- umf_to_factor(object@data) + nd <- get_base_newdata(object, type) + if(! covariate %in% names(nd)){ + stop("Covariate not in this submodel", call.=FALSE) + } + values <- get_cov_seq(covariate, object, type) + nd <- nd[rep(1, length(values)),,drop=FALSE] + nd[[covariate]] <- values + + pr <- predict(object, type=type, newdata=nd, level=level, ...) + pr$covariate <- covariate + pr$covariateValue <- values + pr +}) + +setGeneric("plotEffects", function(object, ...) standardGeneric("plotEffects")) + +setMethod("plotEffects", "unmarkedFit", + function(object, type, covariate, level=0.95, ...){ + + # Get data for plot + plot_data <- plotEffectsData(object, type, covariate, level, ...) + + # Is covariate a factor? + is_factor <- is.factor(plot_data$covariateValue) + + # If not a factor (i.e., numeric) + if(!is_factor){ + # Setup basic plot structure + plot(plot_data$covariateValue, plot_data$Predicted, type='l', + xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)), + ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T))) + #Draw error ribbon + polygon(x=c(plot_data$covariateValue, rev(plot_data$covariateValue)), + y=c(plot_data$upper, rev(plot_data$lower)), + border=NA, col='gray90') + #Draw line on top of ribbon + lines(plot_data$covariateValue, plot_data$Predicted) + + # If covariate is a factor + } else { + # Get number of unique factor levels (groups) + ngroup <- nrow(plot_data) + # Calculate width of error bar tips based on # of factor levels + bw <- ngroup/10 + + # Add basic plot of points + plot(1:ngroup, plot_data$Predicted, + xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)), + xlim=c(1-bw/2,ngroup+bw/2), + xaxt='n', + ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T))) + + # Draw new x-axis with factor labels + axis(1, at=1:ngroup, labels=plot_data$covariateValue) + + # Draw error bars + for (i in 1:ngroup){ + segments(i, plot_data$lower[i], i, plot_data$upper[i]) + segments(i-bw/2, plot_data$lower[i], i+bw/2, plot_data$lower[i]) + segments(i-bw/2, plot_data$upper[i], i+bw/2, plot_data$upper[i]) + } + #Draw points on top of bars + points(1:ngroup, plot_data$Predicted, pch=19) + } + +}) diff --git a/man/plotEffects.Rd b/man/plotEffects.Rd new file mode 100644 index 0000000..d176d43 --- /dev/null +++ b/man/plotEffects.Rd @@ -0,0 +1,100 @@ +\name{plotEffects} +\alias{plotEffects} +\alias{plotEffects-methods} +\alias{plotEffects,unmarkedFit-method} +\alias{plotEffectsData} +\alias{plotEffectsData-methods} +\alias{plotEffectsData,unmarkedFit-method} + +\title{Plot marginal effects of covariates in unmarked models} + +\description{This function generates a plot visualizing the effects + of a single covariate on a parameter (e.g. occupancy, abundance) in an unmarked + model. If the covariate is numeric, the result is a line plot with an error + ribbon where the x-axis is the range of the covariate and the y-axis is the + predicted parameter value. If the covariate is an R factor (i.e., categorical), + the x-axis instead contains each unique value of the covariate. + + All covariates in the model besides the one being plotted are held either at + their median value (if they are numeric) or at their reference level (if they + are factors). + + Some types of unmarked models may require additional arguments, which are passed + to the matching \code{predict} method. For example, \code{unmarkedFitOccuMulti} + models require the \code{species} argument to be included in the function + call in order to work properly. + + If you want to customize a plot, the easiest approach is to get data + formatted for plotting using \code{plotEffectsData}, and use that. If + you want to see and/or modify the code used by \code{plotEffects} to generate + the default plots, run \code{getMethod("plotEffects", "unmarkedFit")} in + the R console. +} + +\usage{ +\S4method{plotEffects}{unmarkedFit}(object, type, covariate, level=0.95, ...) +\S4method{plotEffectsData}{unmarkedFit}(object, type, covariate, level=0.95, ...) +} + +\arguments{ + \item{object}{A fitted model inheriting class \code{unmarkedFit}} + \item{type}{Submodel in which the covariate of interest can be found, for + example \code{"state"} or \code{"det"}. This will depend on the fitted model} + \item{covariate}{The name of the covariate to be plotted, as a character string} + \item{level}{Confidence level for the error ribbons or bars} + \item{...}{Other arguments passed to the \code{predict} function, required + for some \code{unmarkedFit} types such as \code{unmarkedFitOccuMulti}} +} + +\value{A plot (\code{plotEffects} or a data frame (\code{plotEffectsData}) + containing values to be used in a plot. +} + +\author{Ken Kellner \email{contact@kenkellner.com}} + +\examples{ + +\dontrun{ + +# Simulate data and build an unmarked frame +set.seed(123) +dat_occ <- data.frame(x1=rnorm(500)) +dat_p <- data.frame(x2=rnorm(500*5)) + +y <- matrix(NA, 500, 5) +z <- rep(NA, 500) + +b <- c(0.4, -0.5, 0.3, 0.5) + +re_fac <- factor(sample(letters[1:5], 500, replace=T)) +dat_occ$group <- re_fac +re <- rnorm(5, 0, 1.2) +re_idx <- as.numeric(re_fac) + +idx <- 1 +for (i in 1:500){ + z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]])) + for (j in 1:5){ + y[i,j] <- z[i]*rbinom(1,1, + plogis(b[3] + b[4]*dat_p$x2[idx])) + idx <- idx + 1 + } +} + +umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p) + +# Fit model +(fm <- occu(~x2 ~x1 + group, umf)) + +# Plot marginal effects of various covariates +plotEffects(fm, "state", "x1") +plotEffects(fm, "state", "group") +plotEffects(fm, "det", "x2") + +# Get raw data used for a plot +plotEffectsData(fm, "state", "group") + +# See code used by plotEffects so you can edit it yourself and customize the plot +methods::getMethod("plotEffects", "unmarkedFit") +} +} diff --git a/tests/testthat/test_gdistremoval.R b/tests/testthat/test_gdistremoval.R index 8f453a5..3e205fb 100644 --- a/tests/testthat/test_gdistremoval.R +++ b/tests/testthat/test_gdistremoval.R @@ -343,6 +343,12 @@ test_that("gdistremoval can fit models",{ expect_is(fl, "unmarkedFitList") ms <- modSel(fl) expect_is(ms, "unmarkedModSel") + + #plotEffects + pdf(NULL) + plotEffects(fit, 'lambda', "sc1") + plotEffects(fit, 'rem', 'oc1') + dev.off() }) test_that("gdistremoval predict method works",{ diff --git a/tests/testthat/test_occuTTD.R b/tests/testthat/test_occuTTD.R index f152cde..0fb7d81 100644 --- a/tests/testthat/test_occuTTD.R +++ b/tests/testthat/test_occuTTD.R @@ -382,6 +382,12 @@ test_that("occuTTD can fit a dynamic model",{ expect_equivalent(dim(s[[1]]), c(100,4)) r <- residuals(fit) expect_equivalent(dim(r), c(100,4)) + pdf(NULL) + plotEffects(fit, 'psi', 'elev') + plotEffects(fit, 'col', 'forest') + plotEffects(fit, 'det', 'wind') + plotEffects(fit, 'det', 'obs') + dev.off() #Check ranef r <- ranef(fit) diff --git a/tests/testthat/test_plotEffects.R b/tests/testthat/test_plotEffects.R new file mode 100644 index 0000000..3aea3ec --- /dev/null +++ b/tests/testthat/test_plotEffects.R @@ -0,0 +1,48 @@ +context("plotEffects") + +skip_on_cran() + +set.seed(123) +dat_occ <- data.frame(x1=rnorm(500)) +dat_p <- data.frame(x2=rnorm(500*5)) + +y <- matrix(NA, 500, 5) +z <- rep(NA, 500) + +b <- c(0.4, -0.5, 0.3, 0.5) + +re_fac <- factor(sample(letters[1:5], 500, replace=T)) +dat_occ$group <- re_fac +re <- rnorm(5, 0, 1.2) +re_idx <- as.numeric(re_fac) + +idx <- 1 +for (i in 1:500){ + z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]])) + for (j in 1:5){ + y[i,j] <- z[i]*rbinom(1,1, + plogis(b[3] + b[4]*dat_p$x2[idx])) + idx <- idx + 1 + } +} + +umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p) + +fm <- occu(~x2 ~x1 + group, umf) + +test_that("plotEffects works", { + + pdf(NULL) + plotEffects(fm, "state", "x1") + plotEffects(fm, "state", "group") + plotEffects(fm, "det", "x2") + dev.off() + + expect_error(plotEffects(fm, "state", "x2")) + + dat <- plotEffectsData(fm, "state", "group") + + expect_true(inherits(dat, "data.frame")) + expect_equal(nrow(dat), 5) + +}) |