Merge pull request #170 from kenkellner/predict_ranef

Adds predict method for unmarkedRanef class

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,

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,

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]
+})

R/ranef.R

@@ -1132,14 +1132,31 @@ setMethod("plot", c("unmarkedRanef", "missing"), function(x, y, ...)
     }
 })
 
+setMethod("predict", "unmarkedRanef", function(object, func, nsims=100, ...)
+{
 
+  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)
 
+  out <- array(param, out_dim)
+
+  if(is.vector(s1)){
+    rownames(out) <- nm
+  } else {
+    rownames(out) <- row_nm
+    colnames(out) <- col_nm
+  }
 
-
-
-
-
-
-
+  drop(out)
+})

inst/unitTests/runit.ranef.R

@@ -536,3 +536,83 @@ 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(gr1=mean(x[1:4]), gr2=mean(x[5:9]))
+  }
+
+  pr <- predict(re, fun=myfunc, nsim=10)
+  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
+  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(dim(pr), c(2,5,10))
+  checkEqualsNumeric(pr[1,1:3,1], c(3.5,2.5,1.5))
+
+
+}

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))
+
+}

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,34 @@ 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. 
+
+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 
+you want to obtain multiple different derived parameters from the same posterior
+predictive distribution.
+
 }
 }}
 \keyword{methods}

man/ranef-methods.Rd

@@ -149,6 +149,38 @@ 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 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 will be length 20 since M=20
+
+  #Mean of first 10 sites
+  group1 <- mean(x[1:10])
+  #Mean of sites 11-20
+  group2 <- mean(x[11:20])
+
+  #Naming elements of the output is optional but helpful
+  return(c(group1=group1, group2=group2))
+
+}
+
+#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
+(ppd <- posteriorSamples(re, nsims=100))
 
 }
 \keyword{methods}