Commit 90d5f0a8 authored by Irvine, Kathryn M's avatar Irvine, Kathryn M

Initial Commit

parents
File added
library(rjags)
sum(380,149,700,595)
191*4
sum(380,149,763.44,595)
46*5
150*5
435+750+230+700
library(ggmap)
exp^.22
exp(.22)
1775.00*6
4035*6
24210 + 10650
.27*400
citation(package = "unmarked")
devtools::document("~/Documents/nabat-postdoc/power-analysis/manuscript /MEE-17-11-927_supporting/DataS1/OCacoustic/.")
? sim_process_gen
# explore single season results #
source("~/Documents/nabat-phase2/east-coast_analysis/code/mapping-funs.R")#
source("~/Documents/nabat-phase2/east-coast_analysis/code/static-datasim.R")#
#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_3pct.Rdata")#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_5pct.Rdata")#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_7pct.Rdata")#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_100.Rdata")#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_60.Rdata")#
load("~/Documents/nabat-phase2/east-coast_analysis/cache/sim_30.Rdata")
sim_list <- list("3%: N = 200" = sim_3pct, #
"5%: N = 340" = sim_5pct, #
"7%: N = 475" = sim_7pct, #
"N = 100" = sim_100, #
"N = 60" = sim_60, #
"N = 30" = sim_30)
names(sim_list)
sc <- sim_compare(sim_list, sim_order= names(sim_list[c(6:4,1:3)]))
sc <- sim_compare(sim_list[4:6], sim_order= names(sim_list[6:4]))
A <- c(4,4,4,5,5,5,5,5,6,6,3,3); A <- data.frame(A)#
p <- ggplot(data = A, aes(x = A)) + geom_dotplot(binaxis = "x")#
p#
p + xlim(0,9)#
#
p + xlim(0,9) + ylim(0,0.25) + theme_bw()
library(ggplo2)
library(ggplot2)
A <- c(4,4,4,5,5,5,5,5,6,6,3,3); A <- data.frame(A)#
p <- ggplot(data = A, aes(x = A)) + geom_dotplot(binaxis = "x")#
p#
p + xlim(0,9)#
#
p + xlim(0,9) + ylim(0,0.25) + theme_bw()
p + xlim(0,9) + ylim(0,0.25) + theme_bw() + ylim("n")
p + xlim(0,9) + ylim(0,0.25) + theme_bw() + ylim(0.25)
p + xlim(0,9) + ylim(0,0.25) + theme_bw() + ylab("")
p + xlim(0,9) + ylim(0,0.25) + theme_minimal() + ylab("")
p + xlim(0,9) + ylim(0,0.25) + theme_classic() + ylab("")
p + xlim(0,9) + ylim(0,0.25) + theme_classic() + ylab("") + geom_text(aes(x = 0.1, y = 0.1), "MAD = 0.819 \n  SD = 0.996 \n Sample Mean = 4.583") +
p + xlim(0,9) + ylim(0,0.25) + theme_classic() + ylab("") + geom_text(aes(x = 0.1, y = 0.1), "MAD = 0.819 \n  SD = 0.996 \n Sample Mean = 4.583")
p + xlim(0,9) + ylim(0,0.25) + theme_classic() + ylab("") + geom_text("MAD = 0.819 \n  SD = 0.996 \n Sample Mean = 4.583", mapping = aes(x = 0.1, y = 0.1))
2^4
setwd("~/Documents/nabat-phase2/paper-dynamic_power/DataS1/NABatpow/")
devtools::load_all(".")
data("ex_datR1")
data("ex_datR1")
data("ex_datR9")
r9_s <- ex_datR9[,c("GRTS_ID", "Percent__2")]
r9_s$perFor <- scale(r9_s$Percent__2)#
r9_s <- r9_s[order(r9_s$GRTS_ID), ]
for_idx <- which(ex_datR9$FORESTNAME != " ")#
for_s <- ex_datR9[for_idx, ]#
head(for_s)
for_dat <- for_s[, c("GRTS_ID", "Percent__2")]#
for_dat$perFor <- scale(for_s$Percent__2)
idx_prov <- which(ex_datR9$MAP_UNIT_2 == "Laurentian Mixed Forest Province")#
prov_s <- ex_datR9[idx_prov, ]#
prov_dat <- prov_s[, c("GRTS_ID", "Percent__2")]#
prov_dat$perFor <- scale(prov_dat$Percent__2)
dim(prov_dat)#
dim(for_dat)#
dim(r9_s)
n <- seq(from = 30, to = 90, by = 10)#
J <- rep(4, 7)#
pp <- rep(5, 7)#
#
# nv = 4, ny = 5, ns = 60,90, 120, 240#
scenarios <- cbind(ns = n, nv = J, ny = pp)
n <- seq(from = 30, to = 90, by = 10)#
J <- rep(4, 7)#
pp <- rep(10, 7)
scenarios <- cbind(ns = n, nv = J, ny = pp)
scenarios
start <- Sys.time()#
LHHten <- dyn_decline_sim(n_sim = 100, decline = 0.9, phi_s = 0.8, #
gam_s = rep(0.01, 9), scenarios = scenarios,#
theta_prc = c(1.4,1.4),#
psi1_form = ~perFor, #
covs_psi1 = r9_s, #
path_model = "../models/sim-model-fixed-detection.txt", #
niter = 10000, nburn = 5000, thin = 3, nchains = 3, #
GRTS = TRUE, #
IDvar = "GRTS_ID")#
save(LHHten, file = "../../simulation-code/cache/LHHten.Rdata")#
end3 <- Sys.time() - start#
end3#
#
########################
## Low initial occu ###
########################
#
start <- Sys.time()#
LHLten <- dyn_decline_sim(n_sim = 100, decline = 0.9, phi_s = 0.8, #
gam_s = rep(0.01, 9), scenarios = scenarios,#
theta_prc = c(0,1.4),#
psi1_form = ~perFor, #
covs_psi1 = r9_s, #
path_model = "../models/sim-model-fixed-detection.txt", #
niter = 10000, nburn = 5000, thin = 3, nchains = 3, #
GRTS = TRUE, #
IDvar = "GRTS_ID")#
save(LHLten, file = "../../simulation-code/cache/LHLten.Rdata")#
end4 <- Sys.time() - start
^.*\.Rproj$
^\.Rproj\.user$
om_prc <- omega_prc[1:ny,]
} else {
om_prc <- omega_prc
}
if(!is.null(alpha_prc)){
al_prc <- alpha_prc[1:ny, ]
} else {
al_prc <- alpha_prc
}
if(!is.null(beta_prc)){
bet_prc <- beta_prc[1:ny, ]
} else {
bet_prc <- beta_prc
}
all_results[[i]] <- simple_covs_sim(n_sim = n_sim, ns = ns, nv = nv,
ny = ny, ab_prc = ab_prc,
path_model = path_model,
pars = pars,
theta_prc = theta_prc,
covs_psi1 = covs_psi1,
psi1_form = psi1_form,
beta_prc = bet_prc,
alpha_prc = al_prc,
alpha_covs = alpha_covs,
beta_covs = beta_covs,
form_alpha = form_alpha,
form_beta = form_beta,
omega_prc = om_prc,
p_covs = p_covs,
form_p = form_p,
conf = conf,
nchains = nchains,
niter = niter,
nburn = nburn,
thin = thin)
sr <- all_results[[i]]$sim_df
sr$ns <- ns
sr$nv <- nv
sr$ny <- ny
cri_results <- rbind(cri_results, sr)
}
i <- 1
ns <- scenarios[i,1]
nv <- scenarios[i,2]
ny <- scenarios[i,3]
ns
nv
ny
ab_prc <- prc_ab(gam = gam_s[1:(ny-1)], phi1 = phi_s,
decline = decline)
if(!is.null(omega_prc)){
om_prc <- omega_prc[1:ny,]
} else {
om_prc <- omega_prc
}
if(!is.null(alpha_prc)){
al_prc <- alpha_prc[1:ny, ]
} else {
al_prc <- alpha_prc
}
if(!is.null(beta_prc)){
bet_prc <- beta_prc[1:ny, ]
} else {
bet_prc <- beta_prc
}
simple_covs_sim(n_sim = n_sim, ns = ns, nv = nv,
ny = ny, ab_prc = ab_prc,
path_model = path_model,
pars = pars,
theta_prc = theta_prc,
covs_psi1 = covs_psi1,
psi1_form = psi1_form,
beta_prc = bet_prc,
alpha_prc = al_prc,
alpha_covs = alpha_covs,
beta_covs = beta_covs,
form_alpha = form_alpha,
form_beta = form_beta,
omega_prc = om_prc,
p_covs = p_covs,
form_p = form_p,
conf = conf,
nchains = nchains,
niter = niter,
nburn = nburn,
thin = thin)
ql <- (1-conf)/2
qu <- conf + ql
df1 <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
true_ab <- c(df1$assumed$pi_prc[,1:2])
names(true_ab) <- c(paste0("a[", 1:(ny-1), "]"),
paste0("b[", 1:(ny-1), "]"))
true_gamma <- expit(df1$assumed$pi_prc[,1])
names(true_gamma) <- paste0("gamma[", 1:(ny-1), "]")
true_phi <- expit(df1$assumed$pi_prc[,1] +
df1$assumed$pi_prc[,2])
names(true_phi) <- paste0("phi[", 1:(ny-1), "]")
true_theta <- df1$assumed[[2]]
names(true_theta) <- paste0("theta[", 1:length(true_theta), "]")
omega <- df1$assumed[[3]]
true_omega <- c(omega)
names(true_omega) <- paste0(paste0("omega[",1:ny, ","),
rep(paste0(1:ncol(omega), "]"), ny))
psi1 <- expit(df1$jags_dat$theta_x %*% df1$assumed[[2]])
true_psi <- matrix(NA, nrow = ns, ncol = ny)
colnames(true_psi) <- paste0("psi[", 1:ny, "]")
true_lambda <- true_psi[, -1]
colnames(true_lambda) <- paste0("lambda[", 1:(ny-1), "]")
true_psi[,1] <- psi1
psi1
true_psi[,1]
df1$jags_dat$theta_x
df1$assumed[[2]]
df1 <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
df1
df1$assumed$theta
df1$jags_dat$theta_x
psi1_form
psi1_form <- ~ perFor
df1 <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
n_sim = 2; decline = 0.9; phi_s = 0.8;
gam_s = rep(0.2, 5); scenarios = scenarios;
theta_prc = c(0,1.4);
psi1_form = ~perFor;
covs_psi1 = r9_s;
path_model = "../models/sim-model-fixed-detection.txt";
pars = c("gamma", "phi", "omega",
"a", "b", "theta",
"lambda_T", "avg_lambda", "avg_psi", "lam_T", "lam_avg");
niter = 10000; nburn = 5000; thin = 3; nchains = 3
n_designs <- nrow(scenarios)
all_results <- as.list(1:n_designs)
cri_results <- data.frame()
cri_summ <- data.frame()
i
ns <- scenarios[i,1]
nv <- scenarios[i,2]
ny <- scenarios[i,3]
ab_prc <- prc_ab(gam = gam_s[1:(ny-1)], phi1 = phi_s,
decline = decline)
if(!is.null(omega_prc)){
om_prc <- omega_prc[1:ny,]
} else {
om_prc <- omega_prc
}
if(!is.null(alpha_prc)){
al_prc <- alpha_prc[1:ny, ]
} else {
al_prc <- alpha_prc
}
if(!is.null(beta_prc)){
bet_prc <- beta_prc[1:ny, ]
} else {
bet_prc <- beta_prc
}
bet_prc
all_results[[i]] <- simple_covs_sim(n_sim = n_sim, ns = ns, nv = nv,
ny = ny, ab_prc = ab_prc,
path_model = path_model,
pars = pars,
theta_prc = theta_prc,
covs_psi1 = covs_psi1,
psi1_form = psi1_form,
beta_prc = bet_prc,
alpha_prc = al_prc,
alpha_covs = alpha_covs,
beta_covs = beta_covs,
form_alpha = form_alpha,
form_beta = form_beta,
omega_prc = om_prc,
p_covs = p_covs,
form_p = form_p,
conf = conf,
nchains = nchains,
niter = niter,
nburn = nburn,
thin = thin)
ql <- (1-conf)/2
qu <- conf + ql
df1 <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
df1
df1$jags_dat$theta_x
true_ab <- c(df1$assumed$pi_prc[,1:2])
names(true_ab) <- c(paste0("a[", 1:(ny-1), "]"),
paste0("b[", 1:(ny-1), "]"))
true_gamma <- expit(df1$assumed$pi_prc[,1])
names(true_gamma) <- paste0("gamma[", 1:(ny-1), "]")
true_phi <- expit(df1$assumed$pi_prc[,1] +
df1$assumed$pi_prc[,2])
names(true_phi) <- paste0("phi[", 1:(ny-1), "]")
true_theta <- df1$assumed[[2]]
names(true_theta) <- paste0("theta[", 1:length(true_theta), "]")
omega <- df1$assumed[[3]]
true_omega <- c(omega)
names(true_omega) <- paste0(paste0("omega[",1:ny, ","),
rep(paste0(1:ncol(omega), "]"), ny))
psi1 <- expit(df1$jags_dat$theta_x %*% df1$assumed[[2]])
true_psi <- matrix(NA, nrow = ns, ncol = ny)
colnames(true_psi) <- paste0("psi[", 1:ny, "]")
true_lambda <- true_psi[, -1]
colnames(true_lambda) <- paste0("lambda[", 1:(ny-1), "]")
true_psi[,1] <- psi1
for(t in 1:(ny-1)){
true_psi[,t+1] <- true_psi[,t]*true_phi[t] +
(1-true_psi[,t])*true_gamma[t]
}
true_lambda <- true_psi[,2:ny]/true_psi[,1:(ny-1)]
t_avg_lambda <- apply(true_lambda, 2, mean)
true_lambda_T <- true_psi[,ny]/true_psi[,1]
t_der <- cbind(true_psi, "lambda_T" = true_lambda_T)
t_avg_psi <- apply(true_psi, 2, mean)
truth_der <- c(mean(true_lambda_T), mean(t_avg_lambda),
apply(true_psi, 2, mean))
a_lam <- mean(t_avg_psi[2:ny]/t_avg_psi[1:(ny-1)])
l_T <- t_avg_psi[ny]/t_avg_psi[1]
overall_lam <- c(l_T, a_lam)
names(overall_lam) <- c("lam_T", "lam_avg")
names(truth_der) <- c("lambda_T", "avg_lambda",
paste0("avg_psi[", 1:ny, "]"))
assumed <- c(true_ab[1:(ny-1)], truth_der[2], truth_der[3:length(truth_der)],
true_ab[ny:(2*ny-2)],
true_gamma, overall_lam, truth_der[1],
true_omega, true_phi, true_theta)
n_gen <- length(assumed)
n_gen
all_draws <- as.list(1:n_sim)
all_draws
par_df <- data.frame()
pb <- txtProgressBar(min = 0, max = n_sim, style = 3)
i
test <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
if(strsplit(path_model, split = "/")[[1]][length(strsplit(path_model, split = "/")[[1]])] ==
"sim-model-fixed-detection.txt") {
dat <- test$jags_dat
} else {
idx_out <- which(names(test$jags_dat) %in% c("omega_x"))
dat <- test$jags_dat[-idx_out]
}
dat
z_init <- apply(test$jags_dat$y, 1, sum)
z_init[z_init > 0] <- 1
z <- matrix(rep(z_init, ny), ncol = ny)
inits <- list("z" = z)
jm <- rjags::jags.model(path_model,
data = dat,
inits = inits,
n.chains = nchains, n.adapt = nburn)
path_model
strsplit(path_model, split = "/")[[1]][length(strsplit(path_model, split = "/")[[1]])
]
if(strsplit(path_model, split = "/")[[1]][length(strsplit(path_model, split = "/")[[1]])] ==
"sim-model_general.txt") {
dat <- test$jags_dat
} else {
idx_out <- which(names(test$jags_dat) %in% c("omega_x"))
dat <- test$jags_dat[-idx_out]
}
z_init <- apply(test$jags_dat$y, 1, sum)
z_init[z_init > 0] <- 1
z <- matrix(rep(z_init, ny), ncol = ny)
inits <- list("z" = z)
jm <- rjags::jags.model(path_model,
data = dat,
inits = inits,
n.chains = nchains, n.adapt = nburn)
samps <- rjags::coda.samples(jm, variable.names = pars,
n.iter = niter, thin = thin)
all_draws[[i]] <- samps
all_samps <- do.call(rbind, samps)
medians <- apply(all_samps, 2, median)
means <- apply(all_samps, 2, mean)
cri <- t(apply(all_samps, 2, function(x){
c(quantile(x, c(ql, qu)), mean(x))
}))
dimnames(cri)[[2]] <- c(paste0("lower", ql),
paste0("upper", qu),
"est")
ses <- apply(all_samps, 2, sd)
est_df <- data.frame(cri)
est_df$par <- rownames(cri)
est_df$est_med <- medians
est_df$post_se <- ses
est_df$capture <- (assumed >= cri[,1] & assumed <= cri[,2])
est_df$truth <- assumed
est_df$iter <- i
rownames(est_df) <- NULL
est_df
est_df <- est_df[,c(4,1:3,5:9)]
est_df
par_df <- rbind(par_df, est_df)
par_df
setTxtProgressBar(pb, i)
print(paste("Completed iteration", i, "of", n_sim))
ql <- (1-conf)/2
qu <- conf + ql
df1 <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
# compute true values for derived parameters from
# assumed parameters
# a, b
true_ab <- c(df1$assumed$pi_prc[,1:2])
names(true_ab) <- c(paste0("a[", 1:(ny-1), "]"),
paste0("b[", 1:(ny-1), "]"))
# gamma/phi
true_gamma <- expit(df1$assumed$pi_prc[,1])
names(true_gamma) <- paste0("gamma[", 1:(ny-1), "]")
true_phi <- expit(df1$assumed$pi_prc[,1] +
df1$assumed$pi_prc[,2])
names(true_phi) <- paste0("phi[", 1:(ny-1), "]")
true_theta <- df1$assumed[[2]]
names(true_theta) <- paste0("theta[", 1:length(true_theta), "]")
omega <- df1$assumed[[3]]
true_omega <- c(omega)
names(true_omega) <- paste0(paste0("omega[",1:ny, ","),
rep(paste0(1:ncol(omega), "]"), ny))
# true values for derived parameters
psi1 <- expit(df1$jags_dat$theta_x %*% df1$assumed[[2]])
true_psi <- matrix(NA, nrow = ns, ncol = ny)
colnames(true_psi) <- paste0("psi[", 1:ny, "]")
true_lambda <- true_psi[, -1]
colnames(true_lambda) <- paste0("lambda[", 1:(ny-1), "]")
true_psi[,1] <- psi1
for(t in 1:(ny-1)){
true_psi[,t+1] <- true_psi[,t]*true_phi[t] +
(1-true_psi[,t])*true_gamma[t]
}
true_lambda <- true_psi[,2:ny]/true_psi[,1:(ny-1)]
t_avg_lambda <- apply(true_lambda, 2, mean)
true_lambda_T <- true_psi[,ny]/true_psi[,1]
t_der <- cbind(true_psi, "lambda_T" = true_lambda_T)
t_avg_psi <- apply(true_psi, 2, mean)
truth_der <- c(mean(true_lambda_T), mean(t_avg_lambda),
apply(true_psi, 2, mean))
a_lam <- mean(t_avg_psi[2:ny]/t_avg_psi[1:(ny-1)])
l_T <- t_avg_psi[ny]/t_avg_psi[1]
overall_lam <- c(l_T, a_lam)
names(overall_lam) <- c("lam_T", "lam_avg")
names(truth_der) <- c("lambda_T", "avg_lambda",
paste0("avg_psi[", 1:ny, "]"))
# not includeing alpha and beta in this simulation, but could...
assumed <- c(true_ab[1:(ny-1)], truth_der[2], truth_der[3:length(truth_der)],
true_ab[ny:(2*ny-2)],
true_gamma, overall_lam, truth_der[1],
true_omega, true_phi, true_theta)
n_gen <- length(assumed)
# initialize storage
all_draws <- as.list(1:n_sim)
par_df <- data.frame()
pb <- txtProgressBar(min = 0, max = n_sim, style = 3)
for(i in 1:n_sim){
# generate data
test <- dynoccu_datsim(ns = ns, nv = nv, ab_prc = ab_prc,
ny = ny, theta_prc = theta_prc,
covs_psi1 = covs_psi1, psi1_form = psi1_form,
beta_prc = beta_prc, alpha_prc = alpha_prc,
alpha_covs = alpha_covs, beta_covs = beta_covs,
form_alpha = form_alpha, form_beta = form_beta,
omega_prc = omega_prc,
form_p = form_p,
p_covs = p_covs)
# set up jags data and pars to track
if(strsplit(path_model, split = "/")[[1]][length(strsplit(path_model, split = "/")[[1]])] ==
"sim-model_general.txt") {
dat <- test$jags_dat
} else {
idx_out <- which(names(test$jags_dat) %in% c("omega_x"))
dat <- test$jags_dat[-idx_out]
}
#inits
z_init <- apply(test$jags_dat$y, 1, sum)
z_init[z_init > 0] <- 1
z <- matrix(rep(z_init, ny), ncol = ny)
inits <- list("z" = z)
# initialize model
jm <- rjags::jags.model(path_model,
data = dat,
inits = inits,
n.chains = nchains, n.adapt = nburn)
# posterior samples
samps <- rjags::coda.samples(jm, variable.names = pars,
n.iter = niter, thin = thin)
all_draws[[i]] <- samps
# summarize all draws
all_samps