# 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))
# Best Unbiased Predictors
bup(re, stat="mean") # Posterior mean
bup(re, stat="mode") # Posterior mode
confint(re, level=0.9) # 90% CI
# Plots
plot(re, subset=site %in% c(1:10), layout=c(5, 2), xlim=c(-1,20))
# Compare estimates to truth
sum(N)
sum(bup(re))
# Extract all values in convenient formats
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))
Run the code above in your browser using DataLab