Created
September 8, 2018 07:40
-
-
Save dill/1f601a4ad43a62cb327456f6bbbb49a4 to your computer and use it in GitHub Desktop.
Calculate the average abundance estimates and uncertainty over multiple time periods
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# average predictions | |
library(Distance) | |
library(dsm) | |
# download data from http://workshops.distancesampling.org/stand-intermed-2018/practicals/spermwhale.RData | |
load("spermwhale.RData") | |
# make sure we get the same results every time | |
set.seed(123) | |
# fit a simple model | |
df_hn <- ds(data=dist, truncation=6000, key="hn", adjustment=NULL) | |
dsm_nb_xy_ms <- dsm(count~#s(x,y, bs="ts") + | |
#s(Depth, bs="ts") + | |
s(DistToCAS, bs="ts") + | |
s(SST, bs="ts"),# + | |
#s(EKE, bs="ts"),# + | |
#s(NPP, bs="ts"), | |
df_hn, segs, obs, | |
family=nb()) | |
summary(dsm_nb_xy_ms) | |
plot(dsm_nb_xy_ms, page=1) | |
## extra data generation | |
# now create some extra prediction data to pretend that we | |
# have temporal variation in SST | |
# don't need to do this for real data of course | |
pred <- rbind(predgrid, predgrid) | |
pred$SST[(nrow(predgrid)+1):nrow(pred)] <- | |
pred$SST[(nrow(predgrid)+1):nrow(pred)] + | |
rnorm((nrow(predgrid))) | |
# check that did something | |
# plot(pred$SST[1:nrow(predgrid)], pred$SST[(nrow(predgrid)+1):nrow(pred)], | |
# xlab="old", ylab="new") | |
## end extra data generation | |
## here are some checks | |
# abundance in the "old" area | |
N_old <- predict(dsm_nb_xy_ms, pred[1:nrow(predgrid), ]) | |
sum(N_old) | |
# abundance in the "new" area | |
N_new <- predict(dsm_nb_xy_ms, pred[(nrow(predgrid)+1):nrow(pred), ]) | |
sum(N_new) | |
# average abundance | |
mean(c(sum(N_old), sum(N_new))) | |
# average abundance a different way | |
N_avg <- predict(dsm_nb_xy_ms, pred) | |
N_avg <- mean(c(sum(N_avg[1:nrow(predgrid)]), | |
sum(N_avg[(nrow(predgrid)+1):nrow(pred)]))) | |
# same answer! :) | |
## uncertainty | |
# abundance in the "old" area | |
var_old <- dsm.var.gam(dsm_nb_xy_ms, pred[1:nrow(predgrid), ], | |
pred[1:nrow(predgrid), ]$off.set) | |
summary(var_old) | |
# abundance in the "new" area | |
var_new <- dsm.var.gam(dsm_nb_xy_ms, pred[(nrow(predgrid)+1):nrow(pred), ], | |
pred[(nrow(predgrid)+1):nrow(pred), ]$off.set) | |
summary(var_new) | |
### posterior simulation | |
# number of posterior samples to take (larger better but only to a point) | |
Nb <- 200 # e.g. Marra et al (2011) | |
# matrix to map coefficients to linear predictor | |
Xp <- predict(dsm_nb_xy_ms, pred, type="lpmatrix") | |
# generate some posterior coefficients, do this all at once into one big matrix | |
betas <- rmvn(n=Nb, coef(dsm_nb_xy_ms), vcov(dsm_nb_xy_ms)) | |
# storage | |
avg_N_sim <- rep(NA, Nb) | |
# now calculate the averages | |
for(i in 1:Nb){ | |
# get the linear predictor (i.e., prediction on the log-link scale) | |
lp <- Xp %*% betas[i, ] | |
# need to put onto the response scale and take care of the grid areas | |
Ns <- exp(lp)*predgrid$off.set | |
# calculate our summary statistic | |
avg_N_sim[i] <- mean(c(sum(Ns[1:nrow(predgrid)]), | |
sum(Ns[(nrow(predgrid)+1):nrow(pred)]))) | |
} | |
# plot the posterior distribution of the average abundance | |
hist(avg_N_sim) | |
# overplot the average we calculate above | |
abline(v=N_avg, col="red", lwd=2) | |
# if everthing worked the red line should be at the mode of the histogram | |
# calculate the standard error | |
avg_N_sim_se <- sqrt(var(avg_N_sim)) | |
# now add in the detection function CV, assuming independence | |
# picking this values out of the summary() output from df_hn | |
total_CV_sq <- (avg_N_sim_se/N_avg)^2 + 0.06670757^2 | |
total_CV <- sqrt(total_CV_sq) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment