Skip to content

Commit

Permalink
Added functionality to use Negative Binomial distribution for observa…
Browse files Browse the repository at this point in the history 
…tion likelihood. Flag for whether to use the NB (0= Poisson, 1=NB) is set in the data. A prior needs to be added for the phi (precision) term in the NB.
  • Loading branch information
@marksorel8
marksorel8 committed Sep 9, 2019
1 parent c8d1a84 commit 6d9a311
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 72 deletions.
19 changes: 15 additions & 4 deletions src/Stan_demo/juv_trap.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ data {
matrix[max(trap_day),NX_M] X_M; # design matrix of outmigrant covariates (first column is 1)
int<lower=0> C[N_trap]; # daily trap catch observations
int<lower=1> elapsed_time[N_trap]; # time (days) of sampling for each trap catch obs
int<lower=0, upper=1> Use_NB; //Flag for whether to use the the Negative Binomial distribution in the observation model
}

transformed data {
Expand All @@ -30,6 +31,7 @@ parameters {
real<lower=-1,upper=1> phi_M; # AR(1) coefficient for log-mean daily outmigrants
real<lower=0> sigma_M; # AR(1) process error SD for log-mean daily outmigrants
vector[max_day] log_M_hat_z; # log-means of daily outmigrants (z-scores)
real<lower=0> phi_obs[Use_NB]; // dispersion paramater for the Negative Binomial observation model for catch.
}

transformed parameters {
Expand Down Expand Up @@ -84,16 +86,25 @@ model {
recap ~ binomial(mark, p[MR_week]);

# Trap catch observations
if(Use_NB){
C ~ neg_binomial_2(C_hat,phi_obs[1]);
}
else{
# Note that a Poisson RV thinned by binomial sampling is Poisson
C ~ poisson(C_hat);
C ~ poisson(C_hat);}
}

generated quantities {
int M[max_day]; # daily outmigrants
int M_tot; # total outmigrants

for(t in 1:max_day)
M[t] = poisson_log_rng(log_M_hat[t]);

for(t in 1:max_day){

if(Use_NB){
M[t] = neg_binomial_2_rng(M_hat[t],phi_obs[1]);}
else{
M[t] = poisson_log_rng(log_M_hat[t]);}
}

M_tot = sum(M);
}
17 changes: 12 additions & 5 deletions src/Stan_demo/juv_trap_Stan_script.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@ stan_dat <- list(N_MR = nrow(MR),
NX_M = ncol(X_M),
X_M = X_M,
C = trap_catch$catch,
elapsed_time = round(trap_catch$elapsed_time))
elapsed_time = round(trap_catch$elapsed_time),
Use_NB=1) #Flag to use Negative Binomial (1) or Poisson (0)

# Stan data for multi-year model
trap_catch_my <- na.omit(trap_catch_all[is.element(trap_catch_all$brood_year, MR_all$brood_year),])
Expand Down Expand Up @@ -120,22 +121,26 @@ stan_init <- function(data, chains)
sigma_p = runif(1,0.1,2),
beta_M = array(rnorm(NX_M, c(log(20), rep(0, NX_M - 1)), 0.5), dim = NX_M),
phi_M = runif(1,0,0.9),
sigma_M = runif(1,0.5,2))))
sigma_M = runif(1,0.5,2),
phi_obs=array(runif(stan_dat$Use_NB,.2,2000)))))
})
}


# Call Stan to fit model
juv_trap_fit <- stan(file = here("src","Stan_demo","juv_trap.stan"),
data = stan_dat,
init = stan_init(stan_dat,3),
pars = c("beta_M","phi_M","sigma_M",
"beta_p","sigma_p","p",
"M_hat","M","M_tot","C_hat"),
"M_hat","M","M_tot","C_hat","phi_obs"),
chains = 3, iter = 1500, warmup = 500, thin = 1, cores = 3,
control = list(adapt_delta = 0.99, max_treedepth = 13))

# Print fitted model
print(juv_trap_fit, pars = c("M_hat","M","p","C_hat"), include = F, probs = c(0.05,0.5,0.95))
print(juv_trap_fit, pars = c("phi_M","sigma_M",
"beta_p","sigma_p"
,"M_tot","phi_obs"), include = T, probs = c(0.05,0.5,0.95))

# Check it out in Shinystan
launch_shinystan(juv_trap_fit)
Expand Down Expand Up @@ -168,7 +173,9 @@ stan_init_my <- function(data, chains)
# beta_M = array(rnorm(NX_M, c(log(20), rep(0, NX_M - 1)), 0.5), dim = NX_M),
mu_M = array(rnorm(N_year, log(20), 0.5), dim = N_year),
phi_M = runif(N_year,0,0.9),
sigma_M = runif(N_year,0.5,2))))
sigma_M = runif(N_year,0.5,2),
phi_obs = runif(1,0.1,2)
)))
})
}

Expand Down
111 changes: 48 additions & 63 deletions src/Stan_demo/juv_trap_Stan_script_Wenatchee.R
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,5 @@
if(.Platform$OS.type == "windows") options(device=windows)



install.packages("dataRetrieval")
install.packages("matrixStats")
install.packages("lubridate")
install.packages("rstan")
install.packages("loo")
install.packages("shinystan")
install.packages("yarrr")
install.packages("corrplot")
install.packages("here")


library(dataRetrieval)
library(matrixStats)
library(lubridate)
Expand All @@ -23,8 +10,6 @@ library(yarrr)
library(corrplot)
library(here)



if(file.exists(here("src","Stan_demo","juv_trap_fit.RData")))
load(here("src","Stan_demo","juv_trap_fit.RData"))
if(file.exists(here("src","Stan_demo","juv_trap_my_fit.RData")))
Expand Down Expand Up @@ -52,22 +37,30 @@ source(here("src","Load Screw Trap Data.R"))
screw_trap_dat<-load_dat()
screw_trap_dat$chiw$chiw_catch<-droplevels(screw_trap_dat$chiw$chiw_catch[!is.na(screw_trap_dat$chiw$chiw_catch$allSubs),])

# Read in Tucannon trap catch and mark-recapture datasets
trap_catch_all <- read.table(here("src","Stan_demo","trap_catch_all.txt"), sep = "\t", header = T)
screw_trap_dat$chiw$chiw_catch$week<-ceiling( (as.numeric(screw_trap_dat$chiw$chiw_catch$DOY)-50)/7)
#as.numeric(format(screw_trap_dat$chiw$chiw_catch$EndDat,"%W"))#add week


min_DOY<-min(as.numeric(screw_trap_dat$chiw$chiw_catch$DOY))
screw_trap_dat$chiw$chiw_catch$DOY<-as.numeric(screw_trap_dat$chiw$chiw_catch$DOY)-min_DOY+1

# add week
screw_trap_dat$chiw$chiw_catch$week<-ceiling(as.numeric(screw_trap_dat$chiw$chiw_catch$DOY)/7)
#min_week<-min(screw_trap_dat$chiw$chiw_catch$week)
#screw_trap_dat$chiw$chiw_catch$week<-screw_trap_dat$chiw$chiw_catch$week
#add elapsed time
screw_trap_dat$chiw$chiw_catch$elapsed_time<-1

trap_catch_all <-screw_trap_dat$chiw$chiw_catch[,c("year","DOY","week","allSubs","elapsed_time")]
trap_catch_all$year<-as.integer(trap_catch_all$year)
trap_catch_all$DOY<-as.integer(trap_catch_all$DOY)-50
trap_catch_all$DOY<-as.integer(trap_catch_all$DOY)

MR_all <- read.table(here("src","Stan_demo","MR_all.txt"), sep = "\t", header = T)
head(MR_all)

#Mark Recapture data
screw_trap_dat$chiw$chiw_effic$Date2<-as.Date(c(as.Date(screw_trap_dat$chiw$chiw_effic$Date[1:88],format="%d-%b-%y"),as.Date(screw_trap_dat$chiw$chiw_effic$Date[89:143],format="%m/%d/%y"),rep(NA, times=length(144:329))))
screw_trap_dat$chiw$chiw_effic$year<-as.integer(format(screw_trap_dat$chiw$chiw_effic$Date2,"%Y"))
screw_trap_dat$chiw$chiw_effic$week<-ceiling( (as.integer(format(screw_trap_dat$chiw$chiw_effic$Date2,"%j"))-50)/7)
screw_trap_dat$chiw$chiw_effic$year<-as.integer(format(screw_trap_dat$chiw$chiw_effic$Date2,"%Y"))#add year
screw_trap_dat$chiw$chiw_effic$week<-ceiling((as.numeric(format(screw_trap_dat$chiw$chiw_effic$Date2,"%j"))-min_DOY+1)/7)#add week




rel<-reshape::melt( tapply(screw_trap_dat$chiw$chiw_effic$rel, screw_trap_dat$chiw$chiw_effic[,c("year","week")], sum))
rec<-reshape::melt( tapply(screw_trap_dat$chiw$chiw_effic$rec, screw_trap_dat$chiw$chiw_effic[,c("year","week")], sum))
Expand All @@ -76,41 +69,30 @@ MR_all <-cbind(rel[,1:3],rec[,3])[!is.na(rec[,3]),]
colnames(MR_all)[3:4]<-c("mark","recap")

# Download USGS daily mean discharge data

Chiwawa<-"12456500"

flow_all <- renameNWISColumns(readNWISdv(siteNumbers = Chiwawa, parameterCd = "00060", statCd = "00003",
startDate = paste(min(trap_catch_all$year) , "01", "01", sep = "-"),
endDate = paste(max(trap_catch_all$year), "12", "31", sep = "-")))
names(flow_all)[3:4] <- c("date","flow")
write.csv(flow_all, here("src","Stan_demo","flow_all.csv"), row.names = FALSE)




# Align flow data to trapping dates and merge with trap and MR data
dates <- data.frame(year = rep(unique(trap_catch_all$year),
tapply(trap_catch_all$DOY, trap_catch_all$year, max)),
day = unlist(sapply(tapply(trap_catch_all$DOY, trap_catch_all$year, max),
function(x) 1:x)))
row.names(dates) <- NULL

<<<<<<< HEAD
dates$date <- as.Date(paste0(dates$year , "-01-01"), format = "%Y-%m-%d") + dates$day-1
dates$week <- as.integer(format(dates$date,"%W"))
dates$date <- as.Date(paste0(dates$year , "-02-22"), format = "%Y-%m-%d") + dates$day-1
dates$week <- ceiling( dates$day/7)

flow_daily <- merge(dates, flow_all[,c("date","flow")], all.x = TRUE)

=======
dates$date <- as.Date(paste0(dates$year , "-02-19"), format = "%Y-%m-%d") + dates$day-1
dates$week <- ceiling(dates$day/7)
#dates$week <- as.integer(format(dates$date,"%W"))
flow_daily <- merge(dates, flow_all[,c("date","flow")], all.x = TRUE)
>>>>>>> 6d93b8738a73bf75edd34acf1e2462a2b5beda23
flow_weekly <- aggregate(flow ~ week + year, data = flow_daily, mean)

flow_weekly<-merge(trap_catch_all,flow_weekly,by.x = c(1,3),by.y = c(2,1),all.x = T,all.y = F)
flow_weekly<-flow_weekly[!duplicated(flow_weekly[,c(1,2)]),]


# Stan data for single-year model
year <- 2011
trap_catch <- na.omit(trap_catch_all[trap_catch_all$year==year,])
Expand All @@ -131,10 +113,11 @@ stan_dat <- list(N_MR = nrow(MR),
NX_M = ncol(X_M),
X_M = X_M,
C = trap_catch$allSubs,
elapsed_time = round(trap_catch$elapsed_time))
elapsed_time = round(trap_catch$elapsed_time),
Use_NB=)

# Stan data for multi-year model
trap_catch_my <-trap_catch_all #na.omit(trap_catch_all[is.element(trap_catch_all$year, MR_all$year),])
trap_catch_my <- na.omit(trap_catch_all[is.element(trap_catch_all$year, MR_all$year),])
trap_catch_my <- trap_catch_my[trap_catch_my$elapsed_time != 0,]
trap_catch_my$year_f <- as.numeric(factor(trap_catch_my$year))
MR_my <- na.omit(MR_all[is.element(MR_all$year, trap_catch_all$year),])
Expand All @@ -150,12 +133,12 @@ stan_dat_my <- list(N_MR = nrow(MR_my),
NX_p = ncol(X_p_my),
X_p = X_p_my,
N_trap = nrow(trap_catch_my),
trap_year = trap_catch_my$year_f,
trap_year = trap_catch_my$brood_year_f,
trap_week = trap_catch_my$week,
trap_day = trap_catch_my$DOY,
trap_day = trap_catch_my$day,
NX_M = ncol(X_M_my),
X_M = X_M_my,
C = trap_catch_my$allSubs,
C = trap_catch_my$catch,
elapsed_time = round(trap_catch_my$elapsed_time))


Expand All @@ -177,7 +160,8 @@ stan_init <- function(data, chains)
sigma_p = runif(1,0.1,2),
beta_M = array(rnorm(NX_M, c(log(20), rep(0, NX_M - 1)), 0.5), dim = NX_M),
phi_M = runif(1,0,0.9),
sigma_M = runif(1,0.5,2))))
sigma_M = runif(1,0.5,2),
phi_obs=array(runif(stan_dat$Use_NB,.2,2000)))))
})
}

Expand All @@ -187,12 +171,16 @@ juv_trap_fit <- stan(file = here::here("src","Stan_demo","juv_trap.stan"),
init = stan_init(stan_dat,3),
pars = c("beta_M","phi_M","sigma_M",
"beta_p","sigma_p","p",
"M_hat","M","M_tot","C_hat"),
"M_hat","M","M_tot","C_hat","phi_obs"),
chains = 3, iter = 1500, warmup = 500, thin = 1, cores = 3,
control = list(adapt_delta = 0.99, max_treedepth = 13))



# Print fitted model
print(juv_trap_fit, pars = c("M_hat","M","p","C_hat"), include = F, probs = c(0.05,0.5,0.95))
print(juv_trap_fit, pars = c("phi_M","sigma_M",
"beta_p","sigma_p"
,"M_tot","phi_obs"), include =T, probs = c(0.05,0.5,0.95))

# Check it out in Shinystan
launch_shinystan(juv_trap_fit)
Expand Down Expand Up @@ -230,7 +218,7 @@ stan_init_my <- function(data, chains)
}

# Call Stan to fit model
juv_trap_my_fit_2 <- stan(file = here("src","Stan_demo","juv_trap_multiyear.stan"),
juv_trap_my_fit <- stan(file = here("src","Stan_demo","juv_trap_multiyear.stan"),
data = stan_dat_my,
init = stan_init_my(stan_dat_my,3),
pars = c("mu_M","phi_M","sigma_M","Q_M",
Expand Down Expand Up @@ -296,10 +284,11 @@ with(stan_dat, {
legend("topleft",
paste("total = ", round(mean(M_tot), 1),
" (", round(quantile(M_tot, 0.025), 1), ", ",
round(quantile(M_tot, 0.975), 1), ") \n",
"proportion fall = ", round(quantile(p_fall, 0.975), 2),
" (", round(quantile(p_fall, 0.025), 2), ", ",
round(quantile(p_fall, 0.975), 2), ")", sep = ""),
round(quantile(M_tot, 0.975), 1), ")"#,
#"proportion fall = ", round(quantile(p_fall, 0.975), 2),
# " (", round(quantile(p_fall, 0.025), 2), ", ",
#round(quantile(p_fall, 0.975), 2), ")", sep = ""
),
bty = "n")
})

Expand All @@ -312,12 +301,7 @@ with(stan_dat, {
C_hat <- extract1(juv_trap_my_fit_2,"C_hat")

dev.new(width = 10, height = 14)
<<<<<<< HEAD
par(mfcol = c(5,2), mar = c(3,4.5,1,0.5), oma = c(1,0,1,0))
=======
pdf("catch.pdf")
par(mfcol = c(7,2), mar = c(3,4.5,1,0.5), oma = c(1,0,1,0))
>>>>>>> 6d93b8738a73bf75edd34acf1e2462a2b5beda23

with(stan_dat_my, {
c1 <- transparent("blue", 0.3)
Expand All @@ -335,10 +319,10 @@ with(stan_dat_my, {
col = c1, border = NA)
if(par("mfg")[1]==par("mfg")[3]) mtext(side = 1, line = 2.5, "Date")
if(par("mfg")[2]==1) mtext(side = 2, line = 3, "Catch")
mtext(side = 3, line = 0.1, sort(unique(trap_catch_my$brood_year))[i])
mtext(side = 3, line = 0.1, sort(unique(trap_catch_my$year))[i])
}
})
dev.off()

# Annual time series of predicted true outmigrants
M <- extract1(juv_trap_my_fit_2,"M")
M_tot <- extract1(juv_trap_my_fit_2,"M_tot")
Expand All @@ -360,7 +344,7 @@ with(stan_dat_my, {
cex.axis = 1)
if(par("mfg")[1]==par("mfg")[3]) mtext(side = 1, line = 2.5, "Date")
if(par("mfg")[2]==1) mtext(side = 2, line = 3.5, "Outmigrants")
mtext(side = 3, line = 0.1, sort(unique(trap_catch_my$brood_year))[i])
mtext(side = 3, line = 0.1, sort(unique(trap_catch_my$year))[i])
p_fall <- rowSums(M[,i,1:125])/rowSums(M[,i,])
legend("topleft",
paste("total = ", round(mean(M_tot[,i]), 1),
Expand All @@ -380,7 +364,8 @@ log_M_hat_z <- matrix(stan_mean(juv_trap_my_fit_2,"log_M_hat_z"),
byrow = TRUE)

dev.new(width = 10, height = 14)
par(mfcol = c(3,2), mar = c(3,2,1,0.5), oma = c(1,1.5,1,0))
pdf("Stan_chiwawa_subyearling_process_error.pdf" )
par(mfcol = c(7,2), mar = c(3,2,1,0.5), oma = c(1,1.5,1,0))

with(stan_dat_my, {
c1 <- transparent("blue", 0.3)
Expand All @@ -397,12 +382,12 @@ with(stan_dat_my, {
mtext("Process error", side = 2, outer = TRUE)
})


dev.off()
# Time series of total predicted true outmigrants and proportion fall by brood year
M <- extract1(juv_trap_my_fit_2,"M")/1000
M_tot <- extract1(juv_trap_my_fit_2,"M_tot")/1000
p_fall <- apply(M[,,1:125], 1:2, sum)/M_tot
y <- sort(unique(trap_catch_my$year))
y <- sort(unique(trap_catch_my$brood_year))
c1 <- transparent("blue", 0.7)

dev.new()
Expand Down

0 comments on commit 6d9a311

Please sign in to comment.