Skip to content

Instantly share code, notes, and snippets.

@dill
Created September 8, 2018 07:40
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save dill/1f601a4ad43a62cb327456f6bbbb49a4 to your computer and use it in GitHub Desktop.
Save dill/1f601a4ad43a62cb327456f6bbbb49a4 to your computer and use it in GitHub Desktop.
Calculate the average abundance estimates and uncertainty over multiple time periods
# 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