Introduction to the multiocc package
library(interp)
library(MCMCpack)
#> Loading required package: coda
#> Loading required package: MASS
#>
#> Attaching package: 'MASS'
#> The following object is masked from 'package:interp':
#>
#> area
#> ##
#> ## Markov Chain Monte Carlo Package (MCMCpack)
#> ## Copyright (C) 2003-2026 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
#> ##
#> ## Support provided by the U.S. National Science Foundation
#> ## (Grants SES-0350646 and SES-0350613)
#> ##
library(tmvtnorm)
#> Loading required package: mvtnorm
#> Loading required package: Matrix
#> Loading required package: stats4
#> Loading required package: gmm
#> Loading required package: sandwich
library(truncnorm)
library(multiocc)
library(MASS)
library(corrplot)
#> corrplot 0.95 loaded
library(fields)
#> Loading required package: spam
#> Spam version 2.11-4 (2026-05-28) is loaded.
#> Type 'help( Spam)' or 'demo( spam)' for a short introduction
#> and overview of this package.
#> Help for individual functions is also obtained by adding the
#> suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
#>
#> Attaching package: 'spam'
#> The following object is masked from 'package:stats4':
#>
#> mle
#> The following object is masked from 'package:Matrix':
#>
#> det
#> The following objects are masked from 'package:mvtnorm':
#>
#> rmvnorm, rmvt
#> The following objects are masked from 'package:base':
#>
#> backsolve, forwardsolve
#> Loading required package: viridisLite
#> Loading required package: RColorBrewer
#>
#> Try help(fields) to get started.data(detection)
data(occupancy)
data(coords)
DataNames <- list("species"=colnames(detection)[4:9],
"detection"=c("duration"),"occupancy"=c("forest","elev"))
model.input <- multioccbuild(detection, occupancy, coords, DataNames, threshold = 15000)
#> Warning: Rows in detection with missing covariates have been removed for purposes of fitting the model, but the site/season combination is retained in occupancy and therefore predictions will be outputted.Perform some exploratory data analysis
par(mfrow=c(1,3))
hist(occupancy$forest, main="", xlab="Forest")
hist(occupancy$elev, main="", xlab="Elevation")
hist(detection$duration, main="", xlab="Duration")
par(mfrow=c(3,2), mar=c(3,3,3,1))
quilt.plot(coords[,2:3], occupancy$forest[1:267], main="Forest Cover", zlim=c(-1.5,3))
fit <- Tps(coords[,2:3], occupancy$forest[1:267])
out <- predictSurface(fit, df=100)
image.plot(out, main="Forest Cover (interpolated)", zlim=c(-1.5,2))
quilt.plot(coords[,2:3], occupancy$elev[1:267], main="Elevation", zlim=c(-1.5,3.5))
fit <- Tps(coords[,2:3], occupancy$elev[1:267])
out <- predictSurface(fit, df=100)
image.plot(out, main="Elevation (interpolated)", zlim=c(-1.5,2))
quilt.plot(coords[,2:3], detection$duration[1:267], main="Duration", zlim=c(-2.5,3))
fit <- Tps(coords[,2:3], detection$duration[1:267])
out <- predictSurface(fit, df=100)
image.plot(out, main="Duration (Survey 1)", zlim=c(-2.5,2.5))
A short run for demonstration purposes
## Shorter run for demonstration purposes.
## library(tmvtnorm)
mcmc.out <- GibbsSampler(M.iter=10, M.burn=1, M.thin=1, model.input, q=10, sv=FALSE)
#> | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%A longer run (not executed) for scientific results
Can summarize output
summary(mcmc.out$samples$alpha)
#>
#> Iterations = 1:9
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 9
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> Great.tit Int 0.648677 0.07970 0.026566 0.054599
#> Great.tit forest -0.115320 0.03093 0.010311 0.010311
#> Great.tit elev -0.133691 0.03411 0.011371 0.011371
#> Blue.tit Int 0.417075 0.08384 0.027948 0.058228
#> Blue.tit forest -0.085162 0.02538 0.008460 0.015758
#> Blue.tit elev -0.153027 0.04177 0.013922 0.028239
#> Coal.tit Int 0.868624 0.11564 0.038545 0.077637
#> Coal.tit forest -0.012489 0.01144 0.003812 0.003812
#> Coal.tit elev -0.110031 0.05658 0.018862 0.039379
#> Crested.tit Int 0.578194 0.11629 0.038763 0.072414
#> Crested.tit forest 0.004438 0.03424 0.011415 0.013387
#> Crested.tit elev -0.080836 0.03547 0.011824 0.026574
#> Marsh.tit Int 0.387941 0.07746 0.025821 0.054045
#> Marsh.tit forest -0.116182 0.01598 0.005327 0.005327
#> Marsh.tit elev -0.172140 0.02995 0.009985 0.011275
#> Willow.tit Int 0.103039 0.07472 0.024906 0.058038
#> Willow.tit forest 0.064760 0.02804 0.009347 0.018535
#> Willow.tit elev -0.015727 0.01781 0.005937 0.005937
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> Great.tit Int 0.50794 0.60128 0.677308 0.70894 0.728827
#> Great.tit forest -0.14045 -0.13542 -0.125805 -0.10704 -0.053995
#> Great.tit elev -0.16719 -0.15148 -0.150857 -0.12338 -0.067662
#> Blue.tit Int 0.27868 0.37188 0.416328 0.48987 0.518170
#> Blue.tit forest -0.12098 -0.09658 -0.087575 -0.07419 -0.042026
#> Blue.tit elev -0.20897 -0.17679 -0.171426 -0.12586 -0.086640
#> Coal.tit Int 0.66894 0.79805 0.917427 0.95421 0.986211
#> Coal.tit forest -0.02727 -0.01608 -0.013525 -0.01141 0.007578
#> Coal.tit elev -0.18168 -0.14651 -0.119468 -0.06109 -0.027047
#> Crested.tit Int 0.36833 0.52737 0.602434 0.66961 0.706158
#> Crested.tit forest -0.04033 -0.01954 -0.001832 0.02918 0.050734
#> Crested.tit elev -0.13753 -0.10536 -0.072367 -0.04625 -0.041101
#> Marsh.tit Int 0.26006 0.35013 0.388444 0.44819 0.477173
#> Marsh.tit forest -0.13552 -0.12761 -0.117085 -0.10815 -0.092252
#> Marsh.tit elev -0.22398 -0.18978 -0.161597 -0.15409 -0.143039
#> Willow.tit Int -0.02129 0.06448 0.133367 0.16257 0.169323
#> Willow.tit forest 0.02981 0.03292 0.071681 0.07914 0.101467
#> Willow.tit elev -0.04129 -0.02369 -0.018627 -0.01017 0.014129
summary(mcmc.out$samples$rho)
#>
#> Iterations = 1:9
#> Thinning interval = 1
#> Number of chains = 1
#> Sample size per chain = 9
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> Great.tit rho 0.8866 0.10108 0.03369 0.07059
#> Blue.tit rho 0.6517 0.22922 0.07641 0.17392
#> Coal.tit rho 0.2193 0.10948 0.03649 0.03649
#> Crested.tit rho 0.9169 0.07267 0.02422 0.02422
#> Marsh.tit rho 0.8144 0.11236 0.03745 0.07863
#> Willow.tit rho 0.8865 0.09861 0.03287 0.03287
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> Great.tit rho 0.73002 0.8530 0.8825 0.9704 0.9939
#> Blue.tit rho 0.29048 0.5919 0.6017 0.8467 0.9125
#> Coal.tit rho 0.08922 0.1075 0.2362 0.3191 0.3556
#> Crested.tit rho 0.76777 0.9142 0.9304 0.9547 0.9763
#> Marsh.tit rho 0.68038 0.7202 0.7799 0.9101 0.9592
#> Willow.tit rho 0.69707 0.8703 0.9138 0.9441 0.9842Visualize correlation matrix
par(mfrow=c(1,1), mar=c(3,3,3,1))
sigout <- mcmc.out$samples$sig
Sig <- matrix(colMeans(sigout),6,6)
SpeciesCor <- cov2cor(Sig)
rownames(SpeciesCor) <- DataNames$species
colnames(SpeciesCor) <- DataNames$species
corrplot::corrplot(SpeciesCor)
Make predictions from fitted model
y.agg1 <- aggregate(model.input$y[,1], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot1 <- 1*(y.agg1$x>0)
y.agg2 <- aggregate(model.input$y[,2], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot2 <- 1*(y.agg2$x>0)
y.agg3 <- aggregate(model.input$y[,3], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot3 <- 1*(y.agg3$x>0)
y.agg4 <- aggregate(model.input$y[,4], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot4 <- 1*(y.agg4$x>0)
y.agg5 <- aggregate(model.input$y[,5], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot5 <- 1*(y.agg5$x>0)
y.agg6 <- aggregate(model.input$y[,6], by=list(model.input$detection.info$siteID,
model.input$detection.info$season), FUN=sum, na.rm=TRUE)
y.plot6 <- 1*(y.agg6$x>0)
for (yr in c(1,4,7,10)){
print(yr)
range <- which(model.input$occupancy.info$season == yr)
psiout <- mcmc.out$samples$psi
#pout <- mcmc.out$p
dim(psiout)
psi1 <- apply(psiout[,0*2670+range],2,mean)
psi2 <- apply(psiout[,1*2670+range],2,mean)
psi3 <- apply(psiout[,2*2670+range],2,mean)
psi4 <- apply(psiout[,3*2670+range],2,mean)
psi5 <- apply(psiout[,4*2670+range],2,mean)
psi6 <- apply(psiout[,5*2670+range],2,mean)
par(mfrow=c(3,2), mar=c(1,3,3,1))
fit <- Tps(coords[1:267,2:3], psi1)
out <- predictSurface(fit, df=100)
image.plot(out, main="Great Tit", zlim=c(-0.01,1.01))
mtext(paste("Year",yr), side=3, line=-2, outer=TRUE)
y.plot1.in <- y.plot1[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot1.in==1),2:3])
fit <- Tps(coords[1:267,2:3], psi2)
out <- predictSurface(fit, df=100)
image.plot(out, main="Blue Tit", zlim=c(-0.01,1.01))
y.plot2.in <- y.plot2[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot2.in==1),2:3])
fit <- Tps(coords[1:267,2:3], psi3)
out <- predictSurface(fit, df=100)
image.plot(out, main="Coal Tit", zlim=c(-0.01,1.01))
y.plot3.in <- y.plot3[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot3.in==1),2:3])
fit <- Tps(coords[1:267,2:3], psi4)
out <- predictSurface(fit, df=100)
image.plot(out, main="Crested Tit", zlim=c(-0.01,1.01))
y.plot4.in <- y.plot4[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot4.in==1),2:3])
fit <- Tps(coords[1:267,2:3], psi5)
out <- predictSurface(fit, df=100)
image.plot(out, main="Marsh Tit", zlim=c(-0.01,1.01))
y.plot5.in <- y.plot5[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot5.in==1),2:3])
fit <- Tps(coords[1:267,2:3], psi6)
out <- predictSurface(fit, df=100)
image.plot(out, main="Willow Tit", zlim=c(-0.01,1.01))
y.plot6.in <- y.plot6[which(model.input$occupancy.info$season ==yr)]
points(coords[which(y.plot6.in==1),2:3])
}
#> [1] 1
#> [1] 4
#> [1] 7
#> [1] 10