Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
fig.align = "center",
fig.width = 6, fig.height = 6,
cache = FALSE,
collapse = TRUE,
comment = "#>",
highlight = TRUE
)
## ----Willow-tit, echo=FALSE, out.width="\\textwidth", fig.cap="(ref:cap-Willow-tit)"----
knitr::include_graphics("figures/Poecile_montanus.jpg")
## ----libraries----------------------------------------------------------------
# Load libraries
library(sp)
library(raster)
library(hSDM)
## ----Kery2010-data------------------------------------------------------------
# Load Kéry et al. 2010 data
data(data.Kery2010, package="hSDM")
head(data.Kery2010)
## ----normalizing--------------------------------------------------------------
# Normalized variables
elev.mean <- mean(data.Kery2010$elevation)
elev.sd <- sd(data.Kery2010$elevation)
juldate.mean <- mean(c(data.Kery2010$juldate1,
data.Kery2010$juldate2,
data.Kery2010$juldate3),na.rm=TRUE)
juldate.sd <- sd(c(data.Kery2010$juldate1,
data.Kery2010$juldate2,
data.Kery2010$juldate3),na.rm=TRUE)
data.Kery2010$elevation <- (data.Kery2010$elevation-elev.mean)/elev.sd
data.Kery2010$juldate1 <- (data.Kery2010$juldate1-juldate.mean)/juldate.sd
data.Kery2010$juldate2 <- (data.Kery2010$juldate2-juldate.mean)/juldate.sd
data.Kery2010$juldate3 <- (data.Kery2010$juldate3-juldate.mean)/juldate.sd
## ----sites, out.width="\\textwidth", fig.cap="(ref:cap-sites)", fig.width=9, fig.height=6----
# Landscape and observation sites
sites.sp <- SpatialPointsDataFrame(coords=data.Kery2010[c("coordx","coordy")],
data=data.Kery2010[,-c(1,2)])
xmin <- min(data.Kery2010$coordx)
xmax <- max(data.Kery2010$coordx)
ymin <- min(data.Kery2010$coordy)
ymax <- max(data.Kery2010$coordy)
res_spatial_cell <- 10
ncol <- ceiling((xmax-xmin)/res_spatial_cell)
nrow <- ceiling((ymax-ymin)/res_spatial_cell)
Xmax <- xmin+ncol*res_spatial_cell
Ymax <- ymin+nrow*res_spatial_cell
ext <- extent(c(xmin, Xmax,
ymin, Ymax))
landscape <- raster(ncols=ncol,nrows=nrow,ext)
values(landscape) <- runif(ncell(landscape),0,1)
landscape.po <- rasterToPolygons(landscape)
par(mar=c(0.1,0.1,0.1,0.1))
plot(landscape.po)
points(sites.sp, pch=16, cex=1, col="black")
## ----Neighborhood-------------------------------------------------------------
# Neighborhood
# Rasters must be projected to correctly compute the neighborhood
crs(landscape) <- '+proj=utm +zone=1'
# Cell for each site
cells <- extract(landscape,sites.sp,cell=TRUE)[,1]
# Neighborhood matrix
ncells <- ncell(landscape)
neighbors.mat <- adjacent(landscape, cells=c(1:ncells), directions=8,
pairs=TRUE, sorted=TRUE)
# Number of neighbors by cell
n.neighbors <- as.data.frame(table(as.factor(neighbors.mat[,1])))[,2]
# Adjacent cells
adj <- neighbors.mat[,2]
## ----Arranging-data-----------------------------------------------------------
# Arranging data
# data.obs
nsite <- length(data.Kery2010$coordx)
count <- c(data.Kery2010$count1,data.Kery2010$count2,data.Kery2010$count3)
juldate <- c(data.Kery2010$juldate1,data.Kery2010$juldate2,
data.Kery2010$juldate3)
site <- rep(1:nsite,3)
data.obs <- data.frame(count,juldate,site)
data.obs <- data.obs[!is.na(data.obs$juldate),]
# data.suit
data.suit <- data.Kery2010[c("coordx","coordy","elevation")]
data.suit$cells <- cells
data.suit <- data.suit[-139,] # Removing site 139 with no juldate
## ----Kery2010-pois-mod, cache=TRUE--------------------------------------------
# hSDM.poisson
data.pois <- data.obs
data.pois$elevation <- data.suit$elevation[as.numeric(as.factor(data.obs$site))]
mod.Kery2010.pois <- hSDM.poisson(counts=data.pois$count,
suitability=~elevation+I(elevation^2),
data=data.pois, beta.start=0)
## ----Kery2010-pois-mcmc, results="markup"-------------------------------------
# Outputs
summary(mod.Kery2010.pois$mcmc)
## ----Kery2010-pois-pred-------------------------------------------------------
# Predictions
npred <- 100
nsamp <- dim(mod.Kery2010.pois$mcmc)[1]
# Abundance-elevation
elev.seq <- seq(500,3000,length.out=npred)
elev.seq.n <- (elev.seq-elev.mean)/elev.sd
beta <- as.matrix(mod.Kery2010.pois$mcmc[,1:3])
tbeta <- t(beta)
X <- matrix(c(rep(1,npred),elev.seq.n,elev.seq.n^2),ncol=3)
N <- matrix(NA,nrow=nsamp,ncol=npred)
for (i in 1:npred) {
N[,i] <- exp(X[i,] %*% tbeta)
}
N.est.pois <- apply(N,2,mean)
N.q1.pois <- apply(N,2,quantile,0.025)
N.q2.pois <- apply(N,2,quantile,0.975)
## ----Kery2010-Nmix-mod, cache=TRUE--------------------------------------------
# hSDM.Nmixture
mod.Kery2010.Nmix <- hSDM.Nmixture(# Observations
counts=data.obs$count,
observability=~juldate+I(juldate^2),
site=data.obs$site,
data.observability=data.obs,
# Habitat
suitability=~elevation+I(elevation^2),
data.suitability=data.suit,
# Predictions
suitability.pred=NULL,
# Chains
burnin=10000, mcmc=5000, thin=5,
# Starting values
beta.start=0,
gamma.start=0,
# Priors
mubeta=0, Vbeta=1.0E6,
mugamma=0, Vgamma=1.0E6,
# Various
seed=1234, verbose=1,
save.p=0, save.N=0)
## ----Kery2010-Nmix-mcmc, results="markup"-------------------------------------
# Outputs
summary(mod.Kery2010.Nmix$mcmc)
## ----Kery2010-Nmix-pred-------------------------------------------------------
# Predictions
nsamp <- dim(mod.Kery2010.Nmix$mcmc)[1]
# Abundance-elevation
beta <- as.matrix(mod.Kery2010.Nmix$mcmc[,1:3])
tbeta <- t(beta)
N <- matrix(NA,nrow=nsamp,ncol=npred)
for (i in 1:npred) {
N[,i] <- exp(X[i,] %*% tbeta)
}
N.est.Nmix <- apply(N,2,mean)
N.q1.Nmix <- apply(N,2,quantile,0.025)
N.q2.Nmix <- apply(N,2,quantile,0.975)
# Detection-Julian date
juldate.seq <- seq(100,200,length.out=npred)
juldate.seq.n <- (juldate.seq-juldate.mean)/juldate.sd
gamma <- as.matrix(mod.Kery2010.Nmix$mcmc[,4:6])
tgamma <- t(gamma)
W <- matrix(c(rep(1,npred),juldate.seq.n,juldate.seq.n^2),ncol=3)
delta <- matrix(NA,nrow=nsamp,ncol=npred)
for (i in 1:npred) {
delta[,i] <- inv.logit(X[i,] %*% tgamma)
}
delta.est.Nmix <- apply(delta,2,mean)
delta.q1.Nmix <- apply(delta,2,quantile,0.025)
delta.q2.Nmix <- apply(delta,2,quantile,0.975)
## ----Kery2010-Nmix-iCAR-mod, cache=TRUE---------------------------------------
# hSDM.Nmixture.iCAR
mod.Kery2010.Nmix.iCAR <- hSDM.Nmixture.iCAR(# Observations
counts=data.obs$count,
observability=~juldate+I(juldate^2),
site=data.obs$site,
data.observability=data.obs,
# Habitat
suitability=~elevation+I(elevation^2),
data.suitability=data.suit,
# Spatial structure
spatial.entity=data.suit$cells,
n.neighbors=n.neighbors, neighbors=adj,
# Chains
burnin=20000, mcmc=10000, thin=10,
# Starting values
beta.start=0,
gamma.start=0,
Vrho.start=1,
# Priors
mubeta=0, Vbeta=1.0E6,
mugamma=0, Vgamma=1.0E6,
priorVrho="1/Gamma",
shape=1, rate=1,
# Various
seed=1234, verbose=1,
save.rho=0, save.p=0, save.N=0)
## ----Kery2010-Nmix-iCAR-mcmc,results="markup"---------------------------------
summary(mod.Kery2010.Nmix.iCAR$mcmc)
## ----Kery2010-Nmix-iCAR-spatial-effects, out.width="\\textwidth", fig.width=9, fig.height=6, fig.cap="(ref:cap-spatial)"----
# Spatial random effects
rho.pred <- mod.Kery2010.Nmix.iCAR$rho.pred
r.rho.pred <- rasterFromXYZ(cbind(coordinates(landscape),rho.pred))
plot(r.rho.pred)
# Mean abundance by site
ma <- apply(sites.sp@data[,3:5],1,mean,na.rm=TRUE)
points(sites.sp,pch=16,cex=0.5)
points(sites.sp,pch=1,cex=ma/2)
## ----Kery2010-Nmix-iCAR-pred--------------------------------------------------
# Predictions
nsamp <- dim(mod.Kery2010.Nmix.iCAR$mcmc)[1]
# Abundance-elevation
beta <- as.matrix(mod.Kery2010.Nmix.iCAR$mcmc[,1:3])
tbeta <- t(beta)
N <- matrix(NA,nrow=nsamp,ncol=npred)
# Simplified way of obtaining samples for rho
rho.samp <- sample(rho.pred,nsamp,replace=TRUE)
for (i in 1:npred) {
N[,i] <- exp(X[i,] %*% tbeta + rho.samp)
}
N.est.Nmix.iCAR <- apply(N,2,mean)
N.q1.Nmix.iCAR <- apply(N,2,quantile,0.025)
N.q2.Nmix.iCAR <- apply(N,2,quantile,0.975)
# Detection-Julian date
gamma <- as.matrix(mod.Kery2010.Nmix.iCAR$mcmc[,4:6])
tgamma <- t(gamma)
delta <- matrix(NA,nrow=nsamp,ncol=npred)
for (i in 1:npred) {
delta[,i] <- inv.logit(X[i,] %*% tgamma)
}
delta.est.Nmix.iCAR <- apply(delta,2,mean)
delta.q1.Nmix.iCAR <- apply(delta,2,quantile,0.025)
delta.q2.Nmix.iCAR <- apply(delta,2,quantile,0.975)
## ----Kery2010-comp-abundance, fig.cap="(ref:cap-abundance)"-------------------
# Expected abundance - Elevation
par(mar=c(4,4,1,1),cex=1.4,tcl=+0.5)
plot(elev.seq,N.est.pois,type="l",
xlim=c(500,3000),
ylim=c(0,10),
lwd=2,
xlab="Elevation (m a.s.l.)",
ylab="Expected abundance",
axes=FALSE)
#lines(elev.seq,N.q1.pois,lty=3,lwd=1)
#lines(elev.seq,N.q2.pois,lty=3,lwd=1)
axis(1,at=seq(500,3000,by=500),labels=seq(500,3000,by=500))
axis(2,at=seq(0,10,by=2),labels=seq(0,10,by=2))
# Nmix
lines(elev.seq,N.est.Nmix,lwd=2,col="red")
#lines(elev.seq,N.q1.Nmix,lty=3,lwd=1,col="red")
#lines(elev.seq,N.q2.Nmix,lty=3,lwd=1,col="red")
# Nmix.iCAR
lines(elev.seq,N.est.Nmix.iCAR,lwd=2,col="dark green")
#lines(elev.seq,N.q1.Nmix.iCAR,lty=3,lwd=1,col="dark green")
#lines(elev.seq,N.q2.Nmix.iCAR,lty=3,lwd=1,col="dark green")
## ----Kery2010-comp-detection, fig.cap="(ref:cap-detection)"-------------------
# Detection probability - Julian date
par(mar=c(4,4,1,1),cex=1.4,tcl=+0.5)
plot(juldate.seq,delta.est.Nmix,type="l",
xlim=c(100,200),
ylim=c(0,1),
lwd=2,
col="red",
xlab="Julian date",
ylab="Detection probability",
axes=FALSE)
lines(juldate.seq,delta.q1.Nmix,lty=3,lwd=1,col="red")
lines(juldate.seq,delta.q2.Nmix,lty=3,lwd=1,col="red")
axis(1,at=seq(100,200,by=20),labels=seq(100,200,by=20))
axis(2,at=seq(0,1,by=0.2),labels=seq(0,1,by=0.2))
# Nmix.iCAR
lines(juldate.seq,delta.est.Nmix.iCAR,lwd=2,col="dark green")
lines(juldate.seq,delta.q1.Nmix.iCAR,lty=3,lwd=1,col="dark green")
lines(juldate.seq,delta.q2.Nmix.iCAR,lty=3,lwd=1,col="dark green")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.