Multi-species occupancy models

ecology tutorial

occupancy

hierarchical models

Bayesian statistics

MCMC

community ecology

Fit a community (multi-species) occupancy model in base R with a hand-coded Gibbs sampler, and watch partial pooling stabilise sparse-species estimates.

Author

Tidy Ecology

Published

2026-05-22

Modified

2026-09-28

Corrected 27 September 2026: the figure description said independent estimates scatter widely at low detection counts and pooled ones cluster at the community mean; it now says the two mostly agree and pooling mainly moves the fits pinned near a boundary. On 28 September the figure title followed: it said rare-species estimates are pulled toward the community, but the two species seen at only two sites are pulled past the community mean to about 0.1.

When you survey a whole assemblage, most species are rare. A handful turn up almost everywhere; the long tail is seen at a few sites, on a few visits. Fit an occupancy model to each species on its own and the common ones behave, but the sparse ones fall apart: the likelihood has almost nothing to work with.

A community occupancy model treats the species as draws from a shared distribution. Each species keeps its own occupancy and detection probability, but those probabilities are tied together by community-level means and spreads. The rare species then borrow strength from the assemblage instead of standing alone. This tutorial builds the model from scratch in base R, samples it with a short Metropolis-within-Gibbs routine, and compares the pooled estimates against fitting every species separately.

A community with a long rare tail

We simulate 32 species across 90 sites, each visited 4 times. Species-level occupancy and detection are drawn on the logit scale from community distributions, so a few species are genuinely widespread and many are scarce.

library(ggplot2); library(dplyr); library(tidyr)
logit <- qlogis; expit <- plogis

# brand palette + theme
te_ink<-"#16241d"; te_body<-"#2c3a31"; te_forest<-"#275139"; te_faint<-"#5d6b61"
te_sage<-"#93a87f"; te_paper<-"#f5f4ee"; te_line<-"#dad9ca"; te_brick<-"#b5534e"
theme_te <- function(base_size=12){
  theme_minimal(base_size=base_size) +
    theme(plot.background=element_rect(fill=te_paper,colour=NA),
          panel.background=element_rect(fill="#f5f4ee",colour=NA),
          panel.grid.minor=element_blank(),
          panel.grid.major=element_line(colour=te_line,linewidth=0.3),
          axis.title=element_text(colour=te_body), axis.text=element_text(colour=te_faint),
          plot.title=element_text(colour=te_ink,face="bold"),
          plot.subtitle=element_text(colour=te_faint,size=rel(0.9)),
          legend.text=element_text(colour=te_faint,size=rel(0.8)))
}

set.seed(143)
S <- 32L; R <- 90L; K <- 4L                 # species, sites, visits
mu_lpsi <- logit(0.40); sd_lpsi <- 1.30     # community occupancy (logit scale)
mu_lp   <- logit(0.27); sd_lp   <- 0.95     # community detection (logit scale)

lpsi_t <- rnorm(S, mu_lpsi, sd_lpsi)        # true species occupancy logits
lp_t   <- rnorm(S, mu_lp,   sd_lp)          # true species detection logits
psi_t  <- expit(lpsi_t); p_t <- expit(lp_t)

z_t <- matrix(rbinom(S*R, 1, psi_t), nrow=S, ncol=R)          # occupancy state
y   <- matrix(rbinom(S*R, K, ifelse(z_t==1, p_t, 0)), S, R)   # detection counts
det_sites <- rowSums(y > 0)                                    # sites each species was seen at

Detection frequencies run from 2 to 80 sites, and 7 species were detected at fewer than ten sites. Every species was seen at least once here; the separate problem of species that are never detected is taken up in estimating species richness with data augmentation.

Fitting every species on its own

The single-season occupancy likelihood collapses over the unknown occupancy state at each site: a site with a detection must be occupied, and a site with no detection is either unoccupied or occupied but missed. We maximise that likelihood species by species.

nll_ss <- function(par, yi){
  psi <- expit(par[1]); p <- expit(par[2]); det <- yi > 0
  ll <- sum(log(psi) + dbinom(yi[det], K, p, log=TRUE)) +
        sum((!det) * log(psi*(1-p)^K + (1-psi)))
  -ll
}
ss_psi <- ss_p <- rep(NA_real_, S); boundary <- logical(S)
for (s in 1:S){
  op <- optim(c(0,0), nll_ss, yi=y[s,], method="BFGS")
  ss_psi[s] <- expit(op$par[1]); ss_p[s] <- expit(op$par[2])
  boundary[s] <- ss_psi[s]>0.995 || ss_psi[s]<0.005 || ss_p[s]>0.995 || ss_p[s]<0.005
}
sum(boundary)   # species whose independent estimate is pinned at a boundary

[1] 5

5 of the 32 species land on a boundary of the parameter space. The rarest, detected at just 2 sites, comes back with occupancy 0.98 and detection 0.0057.

The community model

Give species s an occupancy logit lpsi[s] and a detection logit lp[s], and let those come from community distributions with means mu_psi, mu_p and standard deviations sd_psi, sd_p. We sample the whole set with random-walk Metropolis updates for the species effects and the hyperparameters, drawing the latent occupancy state at the undetected sites as a Gibbs step, so every species update sees a complete occupancy matrix.

ni <- 18000L; nb <- 6000L; nt <- 3L
lpsi <- logit(pmin(pmax((det_sites+0.5)/(R+1), 0.02), 0.98)); lp <- rep(logit(0.3), S)
mps <- 0; sps <- 1; mp <- logit(0.3); sp <- 1               # hyperparameters
t_lpsi<-0.30; t_lp<-0.26; t_mu<-0.14; t_lsd<-0.14           # proposal scales
keep <- seq(nb+1, ni, by=nt); nk <- length(keep)
PSI <- matrix(NA,nk,S); P <- matrix(NA,nk,S); HYP <- matrix(NA,nk,4); kk <- 0L

for (it in 1:ni){
  psi <- expit(lpsi); p <- expit(lp)
  # sum of occupied sites, drawing the latent state at undetected sites (Gibbs)
  pz <- psi*(1-p)^K / (psi*(1-p)^K + (1-psi))
  z  <- ifelse(y==0, matrix(rbinom(S*R,1,pz),S,R), 1L); sumz <- rowSums(z)
  # species occupancy logits (Metropolis)
  prop <- lpsi + rnorm(S,0,t_lpsi); psip <- expit(prop)
  a <- sumz*log(psip)+(R-sumz)*log(1-psip)+dnorm(prop,mps,sps,log=TRUE) -
       (sumz*log(psi)+(R-sumz)*log(1-psi)+dnorm(lpsi,mps,sps,log=TRUE))
  au <- log(runif(S)) < a; lpsi[au] <- prop[au]
  # species detection logits (Metropolis), over occupied sites
  p <- expit(lp); occ <- z==1; Sy <- rowSums(y*occ); SKy <- rowSums((K-y)*occ)
  prop <- lp + rnorm(S,0,t_lp); pp <- expit(prop)
  a <- Sy*log(pp)+SKy*log(1-pp)+dnorm(prop,mp,sp,log=TRUE) -
       (Sy*log(p)+SKy*log(1-p)+dnorm(lp,mp,sp,log=TRUE))
  au <- log(runif(S)) < a; lp[au] <- prop[au]
  # community means and spreads (Metropolis; half-normal prior on the spreads)
  pr <- mps+rnorm(1,0,t_mu)
  if (log(runif(1)) < sum(dnorm(lpsi,pr,sps,log=TRUE))-sum(dnorm(lpsi,mps,sps,log=TRUE))) mps <- pr
  pr <- exp(log(sps)+rnorm(1,0,t_lsd))
  if (log(runif(1)) < sum(dnorm(lpsi,mps,pr,log=TRUE))+dnorm(pr,0,2.5,log=TRUE)+log(pr)-
                      sum(dnorm(lpsi,mps,sps,log=TRUE))-dnorm(sps,0,2.5,log=TRUE)-log(sps)) sps <- pr
  pr <- mp+rnorm(1,0,t_mu)
  if (log(runif(1)) < sum(dnorm(lp,pr,sp,log=TRUE))-sum(dnorm(lp,mp,sp,log=TRUE))) mp <- pr
  pr <- exp(log(sp)+rnorm(1,0,t_lsd))
  if (log(runif(1)) < sum(dnorm(lp,mp,pr,log=TRUE))+dnorm(pr,0,2.5,log=TRUE)+log(pr)-
                      sum(dnorm(lp,mp,sp,log=TRUE))-dnorm(sp,0,2.5,log=TRUE)-log(sp)) sp <- pr
  if (it>nb && ((it-nb) %% nt == 0)){ kk<-kk+1L
    PSI[kk,]<-expit(lpsi); P[kk,]<-expit(lp); HYP[kk,]<-c(mps,sps,mp,sp) }
}
comm_psi <- colMeans(PSI); hyp <- colMeans(HYP)
c(community_occupancy = round(expit(hyp[1]),3), community_detection = round(expit(hyp[3]),3))

community_occupancy community_detection 
              0.514               0.280

The community mean occupancy is 0.514, close to the assemblage average of 0.494 that generated the data. The among-species spread is the weak point: the posterior for the occupancy standard deviation runs from 1.45 to 2.73 on the logit scale, against the 1.30 that generated the data and the 1.51 spread of this particular draw. Variance components are the hardest part of a hierarchical fit to recover from imperfect detections, so treat the community spread as the least trustworthy number in the output rather than as a result.

Pooling in action

The effect is largest where the independent fit breaks down. For most species the two estimates nearly agree; the long segments belong to the fits pinned near a boundary, above all the two species seen at only two sites, and the community model pulls those well away from the boundary.

comm_mean <- expit(hyp[1])
df  <- data.frame(det=det_sites, ind=ss_psi, comm=comm_psi)
seg <- data.frame(det=df$det, y1=df$ind, y2=df$comm)
ggplot() +
  geom_hline(yintercept=comm_mean, linetype="dashed", colour=te_faint) +
  geom_segment(data=seg, aes(det, y1, xend=det, yend=y2), colour=te_sage, linewidth=0.4) +
  geom_point(data=df, aes(det, ind, shape="Independent (per species)"), colour=te_brick, size=2) +
  geom_point(data=df, aes(det, comm, shape="Community (pooled)"), colour=te_forest, size=2) +
  scale_shape_manual(NULL, values=c("Independent (per species)"=1,"Community (pooled)"=16)) +
  scale_x_continuous(transform="sqrt", breaks=c(2,5,10,20,40,80)) +
  labs(x="sites where the species was detected (sqrt scale)", y="occupancy estimate (psi)",
       title="Pooling mostly moves the estimates pinned at a boundary",
       subtitle="dashed line: community mean; segments join the two estimates for each species") +
  theme_te() + theme(legend.position = "inside", legend.position.inside = c(0.75,0.15),
                     legend.background=element_rect(fill="#ffffff",colour=te_line))

Occupancy estimate against the number of sites where each species was detected, on a sqrt-scaled x axis, with a hollow independent estimate and a filled pooled estimate for each species joined by a segment, and a dashed community mean near 0.51. For most species the two estimates nearly coincide and rise with detections, from below 0.1 at three sites to about 0.9 at eighty. The long segments are the fits pinned near a boundary: the two species seen at only two sites, whose independent estimate of 0.98 drops to about 0.1 when pooled, and five species between 22 and 66 sites, whose independent estimate of 1.00 drops to between about 0.78 and 0.93. — Figure 1: Occupancy estimates against detection frequency, fitted per species and pooled across the community.

Averaged over all species, the mean absolute error in occupancy is 0.125 for the independent fits against 0.067 for the community model. The estimates are less scattered, with a spread of 0.311 against 0.341.

The model also reports the distribution of species itself, with the per-species estimates sitting underneath it.

xx <- seq(-6, 5, length=400)
dens <- data.frame(psi=expit(xx)); lp_ax <- logit(dens$psi)
dens$dp <- dnorm(lp_ax, hyp[1], hyp[2]) / (dens$psi*(1-dens$psi))    # change of variable to psi scale
sp <- data.frame(psi=comm_psi)
ggplot() +
  geom_area(data=subset(dens, psi>1e-3 & psi<0.999), aes(psi, dp), fill=te_sage, alpha=0.35) +
  geom_line(data=subset(dens, psi>1e-3 & psi<0.999), aes(psi, dp), colour=te_forest, linewidth=0.7) +
  geom_rug(data=sp, aes(psi), colour=te_brick, length=unit(0.05,"npc")) +
  geom_point(data=sp, aes(psi, y=0), colour=te_brick, size=1.6, alpha=0.8) +
  labs(x="species occupancy (psi)", y="community density",
       title="The estimated distribution of species occupancy",
       subtitle="curve: fitted community (logit-normal); ticks: per-species posterior means") +
  theme_te()

A shaded density curve from the fitted community distribution, with species occupancy on the x axis and community density on the y axis, and rug ticks for the individual species along the bottom. — Figure 2: The fitted community distribution of occupancy, with per-species posterior means beneath it.

What the pooling does and does not do

Partial pooling is a compromise, not a free lunch. It lets each species borrow from the community distribution in proportion to how little that species tells us on its own, which here mainly moved the fits pinned near a boundary rather than drawing every sparse species to the community mean, and it shrinks a genuinely unusual rare species toward the crowd. The assumption that species are exchangeable draws from one distribution is doing real work; if the assemblage is really two guilds with very different occupancy, a single community distribution will misrepresent both, and that is worth checking (see checking a community occupancy model).

The among-species variance is not recovered well, and that is the usual pattern with imperfect detection. From here you can add site covariates with species-specific responses and let a trait explain the variation among species, which is the subject of community covariates and species traits.

References

Dorazio & Royle 2005 Journal of the American Statistical Association 100(470):389-398 (10.1198/016214505000000015)
Dorazio, Royle, Soderstrom & Glimskar 2006 Ecology 87(4):842-854 (10.1890/0012-9658(2006)87[842:ESRAAB]2.0.CO;2)
Kery & Royle 2008 Journal of Applied Ecology 45(2):589-598 (10.1111/j.1365-2664.2007.01441.x)
Zipkin, DeWan & Royle 2009 Journal of Applied Ecology 46(4):815-822 (10.1111/j.1365-2664.2009.01664.x)
Gelman, Meng & Stern 1996 Statistica Sinica 6(4):733-807
Kery & Royle 2016. Applied Hierarchical Modeling in Ecology, Volume 1. Academic Press. ISBN 978-0-12-801378-6
Royle & Dorazio 2008. Hierarchical Modeling and Inference in Ecology. Academic Press. ISBN 978-0-12-374097-7

--- title: "Multi-species occupancy models" description: "Fit a community (multi-species) occupancy model in base R with a hand-coded Gibbs sampler, and watch partial pooling stabilise sparse-species estimates." date: "2026-05-22 14:00" categories: [ecology tutorial, R, occupancy, hierarchical models, Bayesian statistics, MCMC, community ecology] image: thumbnail.png image-alt: "Occupancy estimates for many species, with per-species values pinned near one pulled well down by the community model." date-modified: "2026-09-28" --- *Corrected 27 September 2026: the figure description said independent estimates scatter widely at low detection counts and pooled ones cluster at the community mean; it now says the two mostly agree and pooling mainly moves the fits pinned near a boundary. On 28 September the figure title followed: it said rare-species estimates are pulled toward the community, but the two species seen at only two sites are pulled past the community mean to about 0.1.* When you survey a whole assemblage, most species are rare. A handful turn up almost everywhere; the long tail is seen at a few sites, on a few visits. Fit an occupancy model to each species on its own and the common ones behave, but the sparse ones fall apart: the likelihood has almost nothing to work with. A community occupancy model treats the species as draws from a shared distribution. Each species keeps its own occupancy and detection probability, but those probabilities are tied together by community-level means and spreads. The rare species then borrow strength from the assemblage instead of standing alone. This tutorial builds the model from scratch in base R, samples it with a short Metropolis-within-Gibbs routine, and compares the pooled estimates against fitting every species separately. ## A community with a long rare tail We simulate `r 32` species across `r 90` sites, each visited `r 4` times. Species-level occupancy and detection are drawn on the logit scale from community distributions, so a few species are genuinely widespread and many are scarce. ```{r} #| label: setup #| message: false #| warning: false library(ggplot2); library(dplyr); library(tidyr) logit <- qlogis; expit <- plogis # brand palette + theme te_ink<-"#16241d"; te_body<-"#2c3a31"; te_forest<-"#275139"; te_faint<-"#5d6b61" te_sage<-"#93a87f"; te_paper<-"#f5f4ee"; te_line<-"#dad9ca"; te_brick<-"#b5534e" theme_te <- function(base_size=12){ theme_minimal(base_size=base_size) + theme(plot.background=element_rect(fill=te_paper,colour=NA), panel.background=element_rect(fill="#f5f4ee",colour=NA), panel.grid.minor=element_blank(), panel.grid.major=element_line(colour=te_line,linewidth=0.3), axis.title=element_text(colour=te_body), axis.text=element_text(colour=te_faint), plot.title=element_text(colour=te_ink,face="bold"), plot.subtitle=element_text(colour=te_faint,size=rel(0.9)), legend.text=element_text(colour=te_faint,size=rel(0.8))) } ``` ```{r} #| label: simulate set.seed(143) S <- 32L; R <- 90L; K <- 4L # species, sites, visits mu_lpsi <- logit(0.40); sd_lpsi <- 1.30 # community occupancy (logit scale) mu_lp <- logit(0.27); sd_lp <- 0.95 # community detection (logit scale) lpsi_t <- rnorm(S, mu_lpsi, sd_lpsi) # true species occupancy logits lp_t <- rnorm(S, mu_lp, sd_lp) # true species detection logits psi_t <- expit(lpsi_t); p_t <- expit(lp_t) z_t <- matrix(rbinom(S*R, 1, psi_t), nrow=S, ncol=R) # occupancy state y <- matrix(rbinom(S*R, K, ifelse(z_t==1, p_t, 0)), S, R) # detection counts det_sites <- rowSums(y > 0) # sites each species was seen at ``` Detection frequencies run from `r min(rowSums(y>0))` to `r max(rowSums(y>0))` sites, and `r sum(rowSums(y>0) < 10)` species were detected at fewer than ten sites. Every species was seen at least once here; the separate problem of species that are never detected is taken up in [estimating species richness with data augmentation](../estimating-species-richness-augmentation/). ## Fitting every species on its own The single-season occupancy likelihood collapses over the unknown occupancy state at each site: a site with a detection must be occupied, and a site with no detection is either unoccupied or occupied but missed. We maximise that likelihood species by species. ```{r} #| label: independent-mle nll_ss <- function(par, yi){ psi <- expit(par[1]); p <- expit(par[2]); det <- yi > 0 ll <- sum(log(psi) + dbinom(yi[det], K, p, log=TRUE)) + sum((!det) * log(psi*(1-p)^K + (1-psi))) -ll } ss_psi <- ss_p <- rep(NA_real_, S); boundary <- logical(S) for (s in 1:S){ op <- optim(c(0,0), nll_ss, yi=y[s,], method="BFGS") ss_psi[s] <- expit(op$par[1]); ss_p[s] <- expit(op$par[2]) boundary[s] <- ss_psi[s]>0.995 || ss_psi[s]<0.005 || ss_p[s]>0.995 || ss_p[s]<0.005 } sum(boundary) # species whose independent estimate is pinned at a boundary ``` `r sum(boundary)` of the `r 32` species land on a boundary of the parameter space. The rarest, detected at just `r sort(det_sites)[1]` sites, comes back with occupancy `r round(ss_psi[order(det_sites)][1],2)` and detection `r signif(ss_p[order(det_sites)][1],2)`. ## The community model Give species `s` an occupancy logit `lpsi[s]` and a detection logit `lp[s]`, and let those come from community distributions with means `mu_psi`, `mu_p` and standard deviations `sd_psi`, `sd_p`. We sample the whole set with random-walk Metropolis updates for the species effects and the hyperparameters, drawing the latent occupancy state at the undetected sites as a Gibbs step, so every species update sees a complete occupancy matrix. ```{r} #| label: community-mcmc ni <- 18000L; nb <- 6000L; nt <- 3L lpsi <- logit(pmin(pmax((det_sites+0.5)/(R+1), 0.02), 0.98)); lp <- rep(logit(0.3), S) mps <- 0; sps <- 1; mp <- logit(0.3); sp <- 1 # hyperparameters t_lpsi<-0.30; t_lp<-0.26; t_mu<-0.14; t_lsd<-0.14 # proposal scales keep <- seq(nb+1, ni, by=nt); nk <- length(keep) PSI <- matrix(NA,nk,S); P <- matrix(NA,nk,S); HYP <- matrix(NA,nk,4); kk <- 0L for (it in 1:ni){ psi <- expit(lpsi); p <- expit(lp) # sum of occupied sites, drawing the latent state at undetected sites (Gibbs) pz <- psi*(1-p)^K / (psi*(1-p)^K + (1-psi)) z <- ifelse(y==0, matrix(rbinom(S*R,1,pz),S,R), 1L); sumz <- rowSums(z) # species occupancy logits (Metropolis) prop <- lpsi + rnorm(S,0,t_lpsi); psip <- expit(prop) a <- sumz*log(psip)+(R-sumz)*log(1-psip)+dnorm(prop,mps,sps,log=TRUE) - (sumz*log(psi)+(R-sumz)*log(1-psi)+dnorm(lpsi,mps,sps,log=TRUE)) au <- log(runif(S)) < a; lpsi[au] <- prop[au] # species detection logits (Metropolis), over occupied sites p <- expit(lp); occ <- z==1; Sy <- rowSums(y*occ); SKy <- rowSums((K-y)*occ) prop <- lp + rnorm(S,0,t_lp); pp <- expit(prop) a <- Sy*log(pp)+SKy*log(1-pp)+dnorm(prop,mp,sp,log=TRUE) - (Sy*log(p)+SKy*log(1-p)+dnorm(lp,mp,sp,log=TRUE)) au <- log(runif(S)) < a; lp[au] <- prop[au] # community means and spreads (Metropolis; half-normal prior on the spreads) pr <- mps+rnorm(1,0,t_mu) if (log(runif(1)) < sum(dnorm(lpsi,pr,sps,log=TRUE))-sum(dnorm(lpsi,mps,sps,log=TRUE))) mps <- pr pr <- exp(log(sps)+rnorm(1,0,t_lsd)) if (log(runif(1)) < sum(dnorm(lpsi,mps,pr,log=TRUE))+dnorm(pr,0,2.5,log=TRUE)+log(pr)- sum(dnorm(lpsi,mps,sps,log=TRUE))-dnorm(sps,0,2.5,log=TRUE)-log(sps)) sps <- pr pr <- mp+rnorm(1,0,t_mu) if (log(runif(1)) < sum(dnorm(lp,pr,sp,log=TRUE))-sum(dnorm(lp,mp,sp,log=TRUE))) mp <- pr pr <- exp(log(sp)+rnorm(1,0,t_lsd)) if (log(runif(1)) < sum(dnorm(lp,mp,pr,log=TRUE))+dnorm(pr,0,2.5,log=TRUE)+log(pr)- sum(dnorm(lp,mp,sp,log=TRUE))-dnorm(sp,0,2.5,log=TRUE)-log(sp)) sp <- pr if (it>nb && ((it-nb) %% nt == 0)){ kk<-kk+1L PSI[kk,]<-expit(lpsi); P[kk,]<-expit(lp); HYP[kk,]<-c(mps,sps,mp,sp) } } comm_psi <- colMeans(PSI); hyp <- colMeans(HYP) c(community_occupancy = round(expit(hyp[1]),3), community_detection = round(expit(hyp[3]),3)) ``` The community mean occupancy is `r round(expit(hyp[1]),3)`, close to the assemblage average of `r round(mean(psi_t),3)` that generated the data. The among-species spread is the weak point: the posterior for the occupancy standard deviation runs from `r round(quantile(HYP[,2],0.025),2)` to `r round(quantile(HYP[,2],0.975),2)` on the logit scale, against the `r sprintf("%.2f", sd_lpsi)` that generated the data and the `r round(sd(lpsi_t),2)` spread of this particular draw. Variance components are the hardest part of a hierarchical fit to recover from imperfect detections, so treat the community spread as the least trustworthy number in the output rather than as a result. ## Pooling in action The effect is largest where the independent fit breaks down. For most species the two estimates nearly agree; the long segments belong to the fits pinned near a boundary, above all the two species seen at only two sites, and the community model pulls those well away from the boundary. ```{r} #| label: fig-shrinkage #| fig-cap: "Occupancy estimates against detection frequency, fitted per species and pooled across the community." #| fig-alt: "Occupancy estimate against the number of sites where each species was detected, on a sqrt-scaled x axis, with a hollow independent estimate and a filled pooled estimate for each species joined by a segment, and a dashed community mean near 0.51. For most species the two estimates nearly coincide and rise with detections, from below 0.1 at three sites to about 0.9 at eighty. The long segments are the fits pinned near a boundary: the two species seen at only two sites, whose independent estimate of 0.98 drops to about 0.1 when pooled, and five species between 22 and 66 sites, whose independent estimate of 1.00 drops to between about 0.78 and 0.93." comm_mean <- expit(hyp[1]) df <- data.frame(det=det_sites, ind=ss_psi, comm=comm_psi) seg <- data.frame(det=df$det, y1=df$ind, y2=df$comm) ggplot() + geom_hline(yintercept=comm_mean, linetype="dashed", colour=te_faint) + geom_segment(data=seg, aes(det, y1, xend=det, yend=y2), colour=te_sage, linewidth=0.4) + geom_point(data=df, aes(det, ind, shape="Independent (per species)"), colour=te_brick, size=2) + geom_point(data=df, aes(det, comm, shape="Community (pooled)"), colour=te_forest, size=2) + scale_shape_manual(NULL, values=c("Independent (per species)"=1,"Community (pooled)"=16)) + scale_x_continuous(transform="sqrt", breaks=c(2,5,10,20,40,80)) + labs(x="sites where the species was detected (sqrt scale)", y="occupancy estimate (psi)", title="Pooling mostly moves the estimates pinned at a boundary", subtitle="dashed line: community mean; segments join the two estimates for each species") + theme_te() + theme(legend.position = "inside", legend.position.inside = c(0.75,0.15), legend.background=element_rect(fill="#ffffff",colour=te_line)) ``` Averaged over all species, the mean absolute error in occupancy is `r round(mean(abs(ss_psi-psi_t)),3)` for the independent fits against `r round(mean(abs(comm_psi-psi_t)),3)` for the community model. The estimates are less scattered, with a spread of `r round(sd(comm_psi),3)` against `r round(sd(ss_psi),3)`. The model also reports the distribution of species itself, with the per-species estimates sitting underneath it. ```{r} #| label: fig-distribution #| fig-cap: "The fitted community distribution of occupancy, with per-species posterior means beneath it." #| fig-alt: "A shaded density curve from the fitted community distribution, with species occupancy on the x axis and community density on the y axis, and rug ticks for the individual species along the bottom." xx <- seq(-6, 5, length=400) dens <- data.frame(psi=expit(xx)); lp_ax <- logit(dens$psi) dens$dp <- dnorm(lp_ax, hyp[1], hyp[2]) / (dens$psi*(1-dens$psi)) # change of variable to psi scale sp <- data.frame(psi=comm_psi) ggplot() + geom_area(data=subset(dens, psi>1e-3 & psi<0.999), aes(psi, dp), fill=te_sage, alpha=0.35) + geom_line(data=subset(dens, psi>1e-3 & psi<0.999), aes(psi, dp), colour=te_forest, linewidth=0.7) + geom_rug(data=sp, aes(psi), colour=te_brick, length=unit(0.05,"npc")) + geom_point(data=sp, aes(psi, y=0), colour=te_brick, size=1.6, alpha=0.8) + labs(x="species occupancy (psi)", y="community density", title="The estimated distribution of species occupancy", subtitle="curve: fitted community (logit-normal); ticks: per-species posterior means") + theme_te() ``` ## What the pooling does and does not do Partial pooling is a compromise, not a free lunch. It lets each species borrow from the community distribution in proportion to how little that species tells us on its own, which here mainly moved the fits pinned near a boundary rather than drawing every sparse species to the community mean, and it shrinks a genuinely unusual rare species toward the crowd. The assumption that species are exchangeable draws from one distribution is doing real work; if the assemblage is really two guilds with very different occupancy, a single community distribution will misrepresent both, and that is worth checking (see [checking a community occupancy model](../checking-a-community-occupancy-model/)). The among-species variance is not recovered well, and that is the usual pattern with imperfect detection. From here you can add site covariates with species-specific responses and let a trait explain the variation among species, which is the subject of [community covariates and species traits](../community-covariates-and-traits/). ## References - Dorazio & Royle 2005 Journal of the American Statistical Association 100(470):389-398 (10.1198/016214505000000015) - Dorazio, Royle, Soderstrom & Glimskar 2006 Ecology 87(4):842-854 (10.1890/0012-9658(2006)87[842:ESRAAB]2.0.CO;2) - Kery & Royle 2008 Journal of Applied Ecology 45(2):589-598 (10.1111/j.1365-2664.2007.01441.x) - Zipkin, DeWan & Royle 2009 Journal of Applied Ecology 46(4):815-822 (10.1111/j.1365-2664.2009.01664.x) - Gelman, Meng & Stern 1996 Statistica Sinica 6(4):733-807 - Kery & Royle 2016. Applied Hierarchical Modeling in Ecology, Volume 1. Academic Press. ISBN 978-0-12-801378-6 - Royle & Dorazio 2008. Hierarchical Modeling and Inference in Ecology. Academic Press. ISBN 978-0-12-374097-7 ## Related tutorials - [Estimating species richness with data augmentation](../estimating-species-richness-augmentation/) - [Community covariates and species traits](../community-covariates-and-traits/) - [Checking a community occupancy model](../checking-a-community-occupancy-model/) - [Single-season occupancy model](../single-season-occupancy-model/)