From 8e5780c74c80b3e5c5eb4fedb8a95f7cbd77f2cf Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 24 Mar 2020 20:48:30 -0400 Subject: [PATCH 1/3] Work on posteriorSamples method --- R/ranef.R | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/R/ranef.R b/R/ranef.R index 6e22b86c..9f4bae6f 100644 --- a/R/ranef.R +++ b/R/ranef.R @@ -1132,12 +1132,53 @@ setMethod("plot", c("unmarkedRanef", "missing"), function(x, y, ...) } }) +setClass("posteriorSamples", + representation(numSites="integer", + numPrimary="integer", + nsims="integer", + samples="integer") + ) + +setMethod("posteriorSamples", "unmarkedRanef", function(object, nsims) +{ + N <- dim(object@post)[1] + K <- dim(object@post)[2] + T <- dim(object@post)[3] + out <- array(NA, c(N, T, nsims)) + for (n in 1:N){ + for (t in 1:T){ + out[n, t, ] <- sample(0:(K-1), nsims, replace=TRUE, + prob=object@post[n,,t]) + } + } + new("posteriorSamples", numSites=N, numPrimary=T, nsims=nsims, + numPrimary=drop(out)) +}) +setMethod("posteriorSamples", "unmarkedFit", function(object, nsims) +{ + ran <- ranef(object) + posteriorSamples(ran, nsims) +}) +setMethod("show", "posteriorSamples", function(object) +{ + + tdim <- character(0) + if(object@numPrimary>1){ + tdim <- paste0("x ", object@numPrimary, " primary periods") + } + cat("Posterior samples from unmarked model") + cat(paste("\n"object@numSites, "sites", tdim, "x", object@nsims, "sims")) + cat("\nShowing first 5 sites. To see n sites, use print(object, n)") + + + +} From 1c8891693f239832ed27f80745ba03260dec8469 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 27 Mar 2020 15:55:46 -0400 Subject: [PATCH 2/3] Add predict method and tests --- DESCRIPTION | 8 ++-- NAMESPACE | 7 +++- R/posteriorSamples.R | 63 ++++++++++++++++++++++++++++ R/ranef.R | 57 ++++--------------------- inst/unitTests/runit.ranef.R | 81 ++++++++++++++++++++++++++++++++++++ man/posteriorSamples.Rd | 60 ++++++++++++++++++++++++++ man/predict-methods.Rd | 26 ++++++++++++ man/ranef-methods.Rd | 24 +++++++++++ 8 files changed, 272 insertions(+), 54 deletions(-) create mode 100644 R/posteriorSamples.R create mode 100644 man/posteriorSamples.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 53699ed5..ef74e899 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 0.13-2.9002 -Date: 2020-03-20 +Version: 0.13-2.9003 +Date: 2020-03-27 Type: Package Title: Models for Data from Unmarked Animals Author: Ian Fiske, Richard Chandler, David Miller, Andy Royle, Marc Kery, Jeff Hostetler, Rebecca Hutchinson, Adam Smith, Ken Kellner @@ -16,8 +16,8 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R' 'multinomPois.R' 'occu.R' 'occuRN.R' 'occuMulti.R' 'pcount.R' 'gmultmix.R' 'pcountOpen.R' 'gdistsamp.R' 'unmarkedFitList.R' 'unmarkedLinComb.R' 'ranef.R' 'boot.R' 'occuFP.R' 'gpcount.R' 'occuPEN.R' 'pcount.spHDS.R' - 'occuMS.R' 'occuTTD.R' 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' - 'distsampOpen.R' 'multmixOpen.R' + 'occuMS.R' 'occuTTD.R' 'distsampOpen.R' 'multmixOpen.R' + 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R' LinkingTo: Rcpp, RcppArmadillo SystemRequirements: GNU make URL: http://groups.google.com/group/unmarked, diff --git a/NAMESPACE b/NAMESPACE index 56a691a3..d57cd99f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,7 +36,8 @@ exportClasses(unmarkedFit, unmarkedFitOccu, unmarkedFitOccuFP, unmarkedFitDS, unmarkedFrameGMM, unmarkedFrameGDS, unmarkedFramePCO, unmarkedFrameGPC, unmarkedEstimate, unmarkedFitList, unmarkedModSel, unmarkedRanef, unmarkedFrameOccuMulti, - unmarkedFrameOccuMS, unmarkedCrossVal) + unmarkedFrameOccuMS, unmarkedCrossVal, + unmarkedPostSamples) # Methods exportMethods(backTransform, coef, confint, coordinates, fitted, getData, @@ -46,7 +47,9 @@ exportMethods(backTransform, coef, confint, coordinates, fitted, getData, projection, residuals, sampleSize, SE, show, simulate, siteCovs, "siteCovs<-", summary, update, vcov, yearlySiteCovs, "yearlySiteCovs<-", "[", smoothed, projected, nonparboot, logLik, - LRT, ranef, bup, crossVal) + LRT, ranef, bup, crossVal, posteriorSamples) + +S3method("print", "unmarkedPostSamples") # Constructors export(unmarkedEstimate, fitList, mapInfo, unmarkedFrame, diff --git a/R/posteriorSamples.R b/R/posteriorSamples.R new file mode 100644 index 00000000..3b84ef07 --- /dev/null +++ b/R/posteriorSamples.R @@ -0,0 +1,63 @@ +setGeneric("posteriorSamples", function(object, nsims, ...){ + standardGeneric("posteriorSamples") + }) + +setClass("unmarkedPostSamples", + representation(numSites="numeric", + numPrimary="numeric", + nsims="numeric", + samples="array") + ) + +setMethod("posteriorSamples", "unmarkedRanef", function(object, nsims=100, ...) +{ + + N <- dim(object@post)[1] + K <- dim(object@post)[2] + T <- dim(object@post)[3] + + out <- array(NA, c(N, T, nsims)) + + for (n in 1:N){ + for (t in 1:T){ + out[n, t, ] <- sample(0:(K-1), nsims, replace=TRUE, + prob=object@post[n,,t]) + } + } + new("unmarkedPostSamples", numSites=N, numPrimary=T, nsims=nsims, + samples=out) + +}) + +setMethod("posteriorSamples", "unmarkedFit", function(object, nsims=100, ...) +{ + ran <- ranef(object) + posteriorSamples(ran, nsims) +}) + +setMethod("show", "unmarkedPostSamples", function(object) +{ + + #tdim <- character(0) + #if(object@numPrimary>1){ + tdim <- paste0("x ", object@numPrimary, " primary periods") + #} + + cat("Posterior samples from unmarked model\n") + cat(paste(object@numSites, "sites", tdim, "x", object@nsims, "sims\n")) + cat(paste0("Showing first 5 sites and first 3 simulations\n", + "To see all samples, use print()\n")) + + print(object@samples[1:5,,1:3]) + +}) + +print.unmarkedPostSamples <- function(x, ...){ + print(x@samples) +} + +setMethod("[", c("unmarkedPostSamples","ANY","ANY","ANY"), + function(x, i, j, k) +{ + x@samples[i,j,k] +}) diff --git a/R/ranef.R b/R/ranef.R index 9f4bae6f..81512af6 100644 --- a/R/ranef.R +++ b/R/ranef.R @@ -1132,55 +1132,16 @@ setMethod("plot", c("unmarkedRanef", "missing"), function(x, y, ...) } }) -setClass("posteriorSamples", - representation(numSites="integer", - numPrimary="integer", - nsims="integer", - samples="integer") - ) - -setMethod("posteriorSamples", "unmarkedRanef", function(object, nsims) +setMethod("predict", "unmarkedRanef", function(object, func, nsims=100, ...) { - N <- dim(object@post)[1] - K <- dim(object@post)[2] - T <- dim(object@post)[3] + ps <- posteriorSamples(object, nsims)@samples + param <- apply(ps, 3, func) - out <- array(NA, c(N, T, nsims)) - - for (n in 1:N){ - for (t in 1:T){ - out[n, t, ] <- sample(0:(K-1), nsims, replace=TRUE, - prob=object@post[n,,t]) - } - } - new("posteriorSamples", numSites=N, numPrimary=T, nsims=nsims, - numPrimary=drop(out)) - -}) - -setMethod("posteriorSamples", "unmarkedFit", function(object, nsims) -{ - ran <- ranef(object) - posteriorSamples(ran, nsims) -}) - -setMethod("show", "posteriorSamples", function(object) -{ - - tdim <- character(0) - if(object@numPrimary>1){ - tdim <- paste0("x ", object@numPrimary, " primary periods") - } - cat("Posterior samples from unmarked model") - cat(paste("\n"object@numSites, "sites", tdim, "x", object@nsims, "sims")) - cat("\nShowing first 5 sites. To see n sites, use print(object, n)") + pr <- rowMeans(param, na.rm=TRUE) + se <- sqrt(apply(param, 1, function(x) stats::var(x, na.rm=TRUE))) + lower <- apply(param, 1, function(x) stats::quantile(x, 0.025, na.rm=TRUE)) + upper <- apply(param, 1, function(x) stats::quantile(x, 0.975, na.rm=TRUE)) - - -} - - - - - + data.frame(Predicted=pr, SE=se, lower=lower, upper=upper) +}) diff --git a/inst/unitTests/runit.ranef.R b/inst/unitTests/runit.ranef.R index 3142f25c..09b0b6c9 100644 --- a/inst/unitTests/runit.ranef.R +++ b/inst/unitTests/runit.ranef.R @@ -536,3 +536,84 @@ test.ranef.occuMulti <- function(){ ar <- as(re, "array") checkEqualsNumeric(colSums(ar),c(3.94470,6.055303),tol=1e-4) } + + + +#--------------------- predict -------------------------------- + +test.ranef.predict <- function(){ + + #Single-season model + set.seed(4564) + R <- 10 + J <- 5 + N <- rpois(R, 3) + y <- matrix(NA, R, J) + y[] <- rbinom(R*J, N, 0.5) + y[1,] <- NA + y[2,1] <- NA + K <- 15 + + umf <- unmarkedFramePCount(y=y) + fm <- pcount(~1 ~1, umf, K=K) + + re <- ranef(fm) + + ps <- posteriorSamples(re, nsim=10) + checkTrue(inherits(ps, "unmarkedPostSamples")) + #One is dropped bc of NA + checkEqualsNumeric(dim(ps@samples), c(9,1,10)) + + myfunc <- function(x){ + c(mean(x[1:4]), mean(x[5:9])) + } + + pr <- predict(re, fun=myfunc, nsim=10) + checkEqualsNumeric(nrow(pr), 2) + checkEqualsNumeric(as.numeric(pr[1,]), c(5.275,0.7115,4.1125,6.1937), + tol=1e-4) + + #Dynamic model + set.seed(7) + M <- 10 + J <- 3 + T <- 5 + lambda <- 5 + gamma <- 0.4 + omega <- 0.6 + p <- 0.5 + N <- matrix(NA, M, T) + y <- array(NA, c(M, J, T)) + S <- G <- matrix(NA, M, T-1) + N[,1] <- rpois(M, lambda) + y[,,1] <- rbinom(M*J, N[,1], p) + for(t in 1:(T-1)) { + S[,t] <- rbinom(M, N[,t], omega) + G[,t] <- rpois(M, gamma) + N[,t+1] <- S[,t] + G[,t] + y[,,t+1] <- rbinom(M*J, N[,t+1], p) + } + + # Prepare data + umf <- unmarkedFramePCO(y = matrix(y, M), numPrimary=T) + summary(umf) + + # Fit model and backtransform + m1 <- pcountOpen(~1, ~1, ~1, ~1, umf, K=20) + re1 <- ranef(m1) + + ps <- posteriorSamples(re1, nsim=10) + checkEqualsNumeric(dim(ps@samples), c(10,5,10)) + checkEqualsNumeric(ps@samples[1,,1],c(7,4,3,1,1)) + + myfunc <- function(x){ + apply(x, 2, function(x) c(mean(x[1:4]), mean(x[5:9]))) + } + + pr <- predict(re1, fun=myfunc, nsim=10) + checkEqualsNumeric(nrow(pr), 2*T) + checkEqualsNumeric(as.numeric(pr[1,]), c(3.525,0.2189,3.2500,3.7500), + tol=1e-4) + + +} diff --git a/man/posteriorSamples.Rd b/man/posteriorSamples.Rd new file mode 100644 index 00000000..0bd0ac74 --- /dev/null +++ b/man/posteriorSamples.Rd @@ -0,0 +1,60 @@ +\name{posteriorSamples} +\alias{posteriorSamples} +\alias{posteriorSamples-methods} +\alias{posteriorSamples,unmarkedRanef-method} +\alias{posteriorSamples,unmarkedFit-method} +\alias{unmarkedPostSamples-class} +\alias{show,unmarkedPostSamples-method} +\alias{[,unmarkedPostSamples,ANY,ANY,ANY-method} + +\title{Draw samples from the posterior predictive distribution} + +\description{ + Draw samples from the empirical Bayes posterior predictive distribution + derived from unmarked models or ranef objects +} + +\usage{ +\S4method{posteriorSamples}{unmarkedRanef}(object, nsims=100, ...) +\S4method{posteriorSamples}{unmarkedFit}(object, nsims=100, ...) +} + +\arguments{ + \item{object}{An object inheriting class \code{unmarkedRanef} or + \code{unmarkedFit}} + \item{nsims}{Number of draws to make from the posterior predictive distribution} + \item{...}{Other arguments} +} + +\value{\code{unmarkedPostSamples} object containing the draws from the + posterior predictive distribution. The draws are in the \code{@samples} slot. +} + +\author{Ken Kellner \email{contact@kenkellner.com}} + +\seealso{ + \code{\link{ranef}}, + \code{\link{predict}} +} + +\examples{ + +# Simulate data under N-mixture model +set.seed(4564) +R <- 20 +J <- 5 +N <- rpois(R, 10) +y <- matrix(NA, R, J) +y[] <- rbinom(R*J, N, 0.5) + +# Fit model +umf <- unmarkedFramePCount(y=y) +fm <- pcount(~1 ~1, umf, K=50) + +# Estimates of conditional abundance distribution at each site +(re <- ranef(fm)) + +#Draw from the posterior predictive distribution +(ppd <- posteriorSamples(re, nsims=100)) + +} diff --git a/man/predict-methods.Rd b/man/predict-methods.Rd index 8ee42a73..e757b192 100644 --- a/man/predict-methods.Rd +++ b/man/predict-methods.Rd @@ -1,5 +1,6 @@ \name{predict-methods} \docType{methods} +\alias{predict} \alias{predict-methods} \alias{predict,ANY-method} \alias{predict,unmarkedFit-method} @@ -14,6 +15,7 @@ \alias{predict,unmarkedFitPCO-method} \alias{predict,unmarkedFitDSO-method} \alias{predict,unmarkedFitList-method} +\alias{predict,unmarkedRanef-method} \title{ Methods for Function predict in Package `unmarked' } \description{ These methods return predicted values from fitted model objects. @@ -32,6 +34,30 @@ These methods return predicted values from fitted model objects. } \item{\code{signature(object = "unmarkedFitList")}}{ "type" depends upon the fitted models +} +\item{\code{signature(object = "unmarkedRanef")}}{ +Use this method to generate the empirical Bayes posterior predictive distribution +for functions of the random variables (latent abundance or occurrence). + +In addition to the output object from \code{ranef}, you must also supply a +custom function to argument \code{func}. The function must take as input a matrix +with dimensions M x T, where M is the number of sites and T is the number of +primary periods (T=1 for single-season models). The output of this function should +be a vector or matrix containing the derived parameters of interest. + +The output of predict is the mean value of each derived parameter, in +column-major order, along with the standard error and 95\% confidence interval. +See \code{\link{ranef}} for an example. + +You may also manually set the number of draws from the posterior predictive +distribution with argument \code{nsims}; the default is 100. + +Alternatively, you can use the \code{\link{posteriorSamples}} function on the +\code{ranef} output object to obtain the full posterior predictive distribution. +This is useful if you are having trouble designing your custom function or if +you want to obtain multiple different derived parameters from the same posterior +predictive distribution. + } }} \keyword{methods} diff --git a/man/ranef-methods.Rd b/man/ranef-methods.Rd index e4a32af9..2d2548bb 100644 --- a/man/ranef-methods.Rd +++ b/man/ranef-methods.Rd @@ -149,6 +149,30 @@ post.df <- as(re, "data.frame") head(post.df) post.arr <- as(re, "array") +#Generate posterior predictive distribution for a function +#of random variables using predict() + +#First, create a function that operates on a MxT matrix where +#M = nsites and T = n primary periods (1 in this case) +#Our function will generate mean abundance for sites 1-10 and sites 11-20 +myfunc <- function(x){ #x is 20 x 1 + + #Mean of first 10 sites + group1 <- mean(x[1:10,1]) + #Mean of sites 11-20 + group2 <- mean(x[11:20,1]) + + return(c(group1, group2)) + +} + +#Get mean, SE, and 95% confidence interval for the two derived parameters +#Using 100 draws from posterior predictive distribution +(predict(re, func=myfunc, nsims=100)) + +#Alternatively, you can return the posterior predictive distribution +#and run operations on it separately +(ppd <- posteriorSamples(re, nsims=100)) } \keyword{methods} From 64bc5a9ab57e458586b5fdfeb4049d9f3601c239 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 29 Mar 2020 17:44:53 -0400 Subject: [PATCH 3/3] Return all samples instead of summarizing --- R/ranef.R | 27 +++++++++++++++++++++------ inst/unitTests/runit.ranef.R | 15 +++++++-------- man/predict-methods.Rd | 12 ++++++++---- man/ranef-methods.Rd | 28 ++++++++++++++++++---------- 4 files changed, 54 insertions(+), 28 deletions(-) diff --git a/R/ranef.R b/R/ranef.R index 81512af6..e7963f42 100644 --- a/R/ranef.R +++ b/R/ranef.R @@ -1135,13 +1135,28 @@ setMethod("plot", c("unmarkedRanef", "missing"), function(x, y, ...) setMethod("predict", "unmarkedRanef", function(object, func, nsims=100, ...) { - ps <- posteriorSamples(object, nsims)@samples + ps <- posteriorSamples(object, nsims=nsims)@samples + s1 <- func(ps[,,1]) + nm <- names(s1) + row_nm <- rownames(s1) + col_nm <- colnames(s1) + + if(is.vector(s1)){ + out_dim <- c(length(s1), nsims) + } else{ + out_dim <- c(dim(s1), nsims) + } + param <- apply(ps, 3, func) - pr <- rowMeans(param, na.rm=TRUE) - se <- sqrt(apply(param, 1, function(x) stats::var(x, na.rm=TRUE))) - lower <- apply(param, 1, function(x) stats::quantile(x, 0.025, na.rm=TRUE)) - upper <- apply(param, 1, function(x) stats::quantile(x, 0.975, na.rm=TRUE)) + out <- array(param, out_dim) - data.frame(Predicted=pr, SE=se, lower=lower, upper=upper) + if(is.vector(s1)){ + rownames(out) <- nm + } else { + rownames(out) <- row_nm + colnames(out) <- col_nm + } + + drop(out) }) diff --git a/inst/unitTests/runit.ranef.R b/inst/unitTests/runit.ranef.R index 09b0b6c9..65f32d8b 100644 --- a/inst/unitTests/runit.ranef.R +++ b/inst/unitTests/runit.ranef.R @@ -565,14 +565,14 @@ test.ranef.predict <- function(){ checkEqualsNumeric(dim(ps@samples), c(9,1,10)) myfunc <- function(x){ - c(mean(x[1:4]), mean(x[5:9])) + c(gr1=mean(x[1:4]), gr2=mean(x[5:9])) } pr <- predict(re, fun=myfunc, nsim=10) - checkEqualsNumeric(nrow(pr), 2) - checkEqualsNumeric(as.numeric(pr[1,]), c(5.275,0.7115,4.1125,6.1937), - tol=1e-4) - + checkEqualsNumeric(dim(pr), c(2,10)) + checkEquals(rownames(pr), c("gr1","gr2")) + checkEqualsNumeric(as.numeric(pr[1,1:3]), c(6.0,4.0,5.75)) + #Dynamic model set.seed(7) M <- 10 @@ -611,9 +611,8 @@ test.ranef.predict <- function(){ } pr <- predict(re1, fun=myfunc, nsim=10) - checkEqualsNumeric(nrow(pr), 2*T) - checkEqualsNumeric(as.numeric(pr[1,]), c(3.525,0.2189,3.2500,3.7500), - tol=1e-4) + checkEqualsNumeric(dim(pr), c(2,5,10)) + checkEqualsNumeric(pr[1,1:3,1], c(3.5,2.5,1.5)) } diff --git a/man/predict-methods.Rd b/man/predict-methods.Rd index e757b192..fcc3b21a 100644 --- a/man/predict-methods.Rd +++ b/man/predict-methods.Rd @@ -45,13 +45,17 @@ with dimensions M x T, where M is the number of sites and T is the number of primary periods (T=1 for single-season models). The output of this function should be a vector or matrix containing the derived parameters of interest. -The output of predict is the mean value of each derived parameter, in -column-major order, along with the standard error and 95\% confidence interval. -See \code{\link{ranef}} for an example. - You may also manually set the number of draws from the posterior predictive distribution with argument \code{nsims}; the default is 100. +The output of \code{predict} will be a vector or array with one more dimension +than the output of the function supplied \code{func}, corresponding to the number +of draws requested \code{nsims}. For example, if \code{func} +outputs a scalar, the output of \code{predict} will be a vector with length +equal to \code{nsims}. If \code{func} outputs a 3x2 matrix, the output of +\code{predict} will be an array with dimensions 3x2x\code{nsims}. +See \code{\link{ranef}} for an example. + Alternatively, you can use the \code{\link{posteriorSamples}} function on the \code{ranef} output object to obtain the full posterior predictive distribution. This is useful if you are having trouble designing your custom function or if diff --git a/man/ranef-methods.Rd b/man/ranef-methods.Rd index 2d2548bb..1cb62184 100644 --- a/man/ranef-methods.Rd +++ b/man/ranef-methods.Rd @@ -152,23 +152,31 @@ post.arr <- as(re, "array") #Generate posterior predictive distribution for a function #of random variables using predict() -#First, create a function that operates on a MxT matrix where -#M = nsites and T = n primary periods (1 in this case) +#First, create a function that operates on a vector of +#length M (if you fit a single-season model) or a matrix of +#dimensions MxT (if a dynamic model), where +#M = nsites and T = n primary periods #Our function will generate mean abundance for sites 1-10 and sites 11-20 -myfunc <- function(x){ #x is 20 x 1 +myfunc <- function(x){ #x will be length 20 since M=20 #Mean of first 10 sites - group1 <- mean(x[1:10,1]) + group1 <- mean(x[1:10]) #Mean of sites 11-20 - group2 <- mean(x[11:20,1]) - - return(c(group1, group2)) + group2 <- mean(x[11:20]) + + #Naming elements of the output is optional but helpful + return(c(group1=group1, group2=group2)) } -#Get mean, SE, and 95% confidence interval for the two derived parameters -#Using 100 draws from posterior predictive distribution -(predict(re, func=myfunc, nsims=100)) +#Get 100 samples of the values calculated in your function +(pr <- predict(re, func=myfunc, nsims=100)) + +#Summarize posterior +data.frame(mean=rowMeans(pr), + se=apply(pr, 1, stats::sd), + lower=apply(pr, 1, stats::quantile, 0.025), + upper=apply(pr, 1, stats::quantile, 0.975)) #Alternatively, you can return the posterior predictive distribution #and run operations on it separately