plot_BAM_QPAD <-
function(spp) {
if (!exists(".BAMCOEFS"))
stop("Use 'load_BAM_QPAD()' to load estimates")
SPP <- getBAMspecieslist()
if (!(spp %in% SPP))
stop("species info not available")
jd <- seq(0.35, 0.55, 0.01) # TSSR
ts <- seq(-0.25, 0.5, 0.01) # JDAY
ls <- seq(0, 0.25, len=length(jd)) # DSLS
xp1 <- expand.grid(JDAY=jd,
TSSR=ts)
xp1$JDAY2 <- xp1$JDAY^2
xp1$TSSR2 <- xp1$TSSR^2
xp1$Jday <- xp1$JDAY * 365
xp1$Tssr <- xp1$TSSR * 24
xp2 <- expand.grid(DSLS=ls,
TSSR=ts)
xp2$DSLS2 <- xp2$DSLS^2
xp2$TSSR2 <- xp2$TSSR^2
xp2$Dsls <- xp2$DSLS * 365
xp2$Tssr <- xp2$TSSR * 24
Xp1 <- model.matrix(~., xp1)
Xp2 <- model.matrix(~., xp2)
if (getBAMversion() < 3) {
lc <- as.factor(1:5)
tr <- seq(0, 1, 0.01)
xq <- expand.grid(LCC=as.factor(lc),
TREE=tr)
} else {
lc <- factor(c("DecidMixed", "Conif", "Open", "Wet"),
c("DecidMixed", "Conif", "Open", "Wet"))
tr <- seq(0, 1, 0.1)
xq <- expand.grid(LCC4=lc, TREE=tr)
xq$LCC2 <- as.character(xq$LCC4)
xq$LCC2[xq$LCC2 %in% c("DecidMixed", "Conif")] <- "Forest"
xq$LCC2[xq$LCC2 %in% c("Open", "Wet")] <- "OpenWet"
xq$LCC2 <- factor(xq$LCC2, c("Forest", "OpenWet"))
}
Xq0 <- model.matrix(~., xq)
cfall <- exp(t(sapply(SPP, function(spp)
unlist(coefBAMspecies(spp, 0, 0)))))
t <- seq(0, 10, 0.1)
r <- seq(0, 4, 0.05)
pp <- sapply(SPP, function(spp) sra_fun(t, cfall[spp,1]))
qq <- sapply(SPP, function(spp) edr_fun(r, cfall[spp,2]))
## model weights
wp <- selectmodelBAMspecies(spp)$sra$wBIC
wq <- selectmodelBAMspecies(spp)$edr$wBIC
names(wp) <- rownames(selectmodelBAMspecies(spp)$sra)
names(wq) <- rownames(selectmodelBAMspecies(spp)$edr)
nsra <- selectmodelBAMspecies(spp)$sra$nobs[1]
nedr <- selectmodelBAMspecies(spp)$edr$nobs[1]
## covariate effects
mi <- bestmodelBAMspecies(spp, type="BIC")
cfi <- coefBAMspecies(spp, mi$sra, mi$edr)
vci <- vcovBAMspecies(spp, mi$sra, mi$edr)
Xp <- if (getBAMversion() > 2 && mi$sra %in% c("9","10","11","12","13","14"))
Xp2 else Xp1
if (getBAMversion() < 3)
colnames(Xp)[1] <- "INTERCEPT"
Xp <- Xp[,names(cfi$sra),drop=FALSE]
lphi1 <- drop(Xp %*% cfi$sra)
pmat <- matrix(exp(lphi1), length(jd), length(ts))
pmax <- sra_fun(10, max(exp(lphi1)))
pmat <- sra_fun(3, pmat)
pmax <- 1
if (getBAMversion() < 3)
colnames(Xq0)[1] <- "INTERCEPT"
Xq <- Xq0[,names(cfi$edr),drop=FALSE]
ltau1 <- drop(Xq %*% cfi$edr)
qmat <- matrix(exp(ltau1), length(lc), length(tr))
qmax <- edr_fun(0.5, max(exp(ltau1)))
qmat <- edr_fun(1, qmat)
qmax <- 1
op <- par(las=1, mfrow=c(3,2))
barplot(wp, space=0, col=grey(1-wp), border="grey", ylim=c(0,1),
main=paste0(spp, " (n=", nsra, ") v", getBAMversion()),
ylab="Model weight", xlab="Time-removal model ID")
barplot(wq, space=0, col=grey(1-wq), border="grey", ylim=c(0,1),
main=paste0(spp, " (n=", nedr, ") v", getBAMversion()),
ylab="Model weight", xlab="Distance sampling model ID")
plot(t, pp[,spp], type="n", ylim=c(0,1),
xlab="Point count duration (min)",
ylab="Probability of singing")
matlines(t, pp, col="grey", lwd=1, lty=1)
lines(t, pp[,spp], col=1, lwd=2)
plot(r*100, qq[,spp], type="n", ylim=c(0,1),
xlab="Point count radius (m)",
ylab="Probability of detection")
matlines(r*100, qq, col="grey", lwd=1, lty=1)
lines(r*100, qq[,spp], col=1, lwd=2)
abline(v=cfall[spp,2]*100, lty=2)
rug(cfall[,2]*100, side=1, col="grey")
box()
xval <- if (mi$sra %in% c("9","10","11","12","13","14"))
ls*365 else jd*365
image(xval, ts*24, pmat,
col = rev(grey(seq(0, pmax, len=12))),
xlab=ifelse(mi$sra %in% c("9","10","11","12","13","14"),
"Days since local springs", "Julian days"),
ylab="Hours since sunrise",
main=paste("Best model:", mi$sra))
box()
if (getBAMversion() < 3)
image(1:nlevels(xq$LCC), tr*100, qmat,
col = rev(grey(seq(0, qmax, len=12))), axes=FALSE,
xlab="Land cover types", ylab="Percent tree cover",
main=paste("Best model:", mi$edr))
if (getBAMversion() > 2)
image(1:nlevels(xq$LCC4), tr*100, qmat,
col = rev(grey(seq(0, qmax, len=12))), axes=FALSE,
xlab="Land cover types", ylab="Percent tree cover",
main=paste("Best model:", mi$edr))
if (getBAMversion() < 3)
axis(1, 1:5, c("DConif","DDecid","SConif","SDecid","Open"))
if (getBAMversion() > 2)
axis(1, 1:nlevels(xq$LCC4), levels(xq$LCC4))
axis(2)
box()
par(op)
invisible(NULL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.