library(jagsUI) library(tidyverse) setwd("c:\\users\\s429217\\onedrive\\data\\koala") ################################################################################ ################################################################################ # JAGS Exponential model mod1 <- "model { for(i in 1:N) { # Occupancy as Bernoulli variable z[i] ~ dbern(psi) # Detection time follows exponential distribution. ttd[i] ~ dexp(lambda) # Model for censoring isCensored[i] ~ dbern(theta[i]) theta[i] <- z[i] * step(ttd[i] - Tmax) + (1 - z[i]) # Log-likelihood ll[i] <- ifelse(isCensored[i], 0, log(lambda * exp(-lambda * ttd[i]))) } # Uniformative priors for probability of occupancy and detection rate psi ~ dunif(0, 1) lambda ~ dgamma(0.01, 0.01) # Estimated detection probability within search period D.hat <- 1 - exp(-lambda * Tmax) # Total estimated occupancies N_occ <- sum(z[]) # log likelihood log.like <- sum(ll[]) }" # write model write(mod1, "exp.model.txt") ################################################################################ # JAGS Weibull model mod2 <- "model { for(i in 1:N) { # Occupancy as Bernoulli variable z[i] ~ dbern(psi) # Detection times follow Weibull distribution ttd[i] ~ dweib(shape, lambda) # Model for censoring isCensored[i] ~ dbern(theta[i]) theta[i] <- z[i] * step(ttd[i] - Tmax) + (1 - z[i]) # Log-likelihood ll[i] <- ifelse(isCensored[i], 0, log(shape * lambda * ttd[i] ^ (shape - 1) * exp(-lambda * ttd[i] ^ shape))) } # Uninformative priors psi ~ dunif(0, 1) lambda ~ dgamma(0.01, 0.01) shape ~ dgamma(0.01, 0.01) # Estimated detection probability within search period D.hat <- 1 - exp(-lambda * Tmax ^ shape) # Estimated number of occupancies N_occ <- sum(z[]) # log likelihood log.like <- sum(ll[]) }" # write model write(mod2, "weib.model.txt") ################################################################################ # JAGS Exponential mixture model mod3 <- "model { for(i in 1:N) { # Occupancy as Bernoulli variable z[i] ~ dbern(psi) # Detection time follows exponential distribution. ttd[i] ~ dexp(lambda[i]) lambda[i] ~ dgamma(alpha, beta) # Model for censoring isCensored[i] ~ dbern(theta[i]) theta[i] <- z[i] * step(ttd[i] - Tmax) + (1 - z[i]) # Log-likelihood ll[i] <- ifelse(isCensored[i], 0, log(lambda[i] * exp(-lambda[i] * ttd[i]))) } # Uniformative priors for probability of occupancy and detection rate psi ~ dunif(0, 1) alpha ~ dunif(0, 100) beta ~ dunif(0, 100) mu <- alpha / beta # Estimated detection probability within search period D.hat <- 1 - exp(-mu * Tmax) # Total estimated occupancies N_occ <- sum(z[]) # log likelihood log.like <- sum(ll[]) }" # write model write(mod3, "exp.mixture.model.txt") ################################################################################ ################################################################################ # function to fit exponential and weibull models to the data # takes as imput a list generated by sim.nohet or sim.het fit <- function(pars) { y <- pars$y yt <- pars$yt n.tree <- pars$n.tree n.occ <- pars$n.occ sd.lambda <- pars$sd.lambda # run the exponential model mod1.jags <- jags(model = "exp.model.txt", data = list(ttd = pars$ttd, isCensored = pars$isCensored, Tmax = pars$Tmax, N = pars$N, z = pars$z), inits = function() list(ttd = pars$ttdst), param = c("psi", "lambda", "N_occ", "D.hat", "log.like"), n.chains = 3, n.iter = 10000, n.burnin = 2000, parallel = TRUE) mod1.sum <- mod1.jags$summary exp.lambda <- mod1.sum[2, 1] # run the Weibull model mod2.jags <- jags(model = "weib.model.txt", data = list(ttd = pars$ttd, isCensored = pars$isCensored, Tmax = pars$Tmax, N = pars$N, z = pars$z), inits = function() list(ttd = pars$ttdst), param = c("psi", "lambda", "shape", "N_occ", "D.hat", "log.like"), n.chains = 3, n.iter = 10000, n.burnin = 2000, parallel = TRUE) mod2.sum <- mod2.jags$summary weib.lambda <- mod2.sum[2, 1] shape <- mod2.sum[3, 1] scale <- 1 / (weib.lambda ^ (1/shape)) # run the Exponential mixture model mod3.jags <- jags(model = "exp.mixture.model.txt", data = list(ttd = pars$ttd, isCensored = pars$isCensored, Tmax = pars$Tmax, N = pars$N, z = pars$z), inits = function() list(ttd = pars$ttdst), param = c("psi", "mu", "alpha", "beta", "N_occ", "D.hat", "log.like"), n.chains = 3, n.iter = 10000, n.burnin = 2000, parallel = TRUE) mod3.sum <- mod3.jags$summary # plot full distribution of ttd and ttd censored at Tmax par(mfrow = c(2, 3)) hist(yt, breaks = 50, freq = F) curve(dexp(x, exp.lambda), 0, max(yt, na.rm = T), add = T, col = "red", lwd = 2) curve(dweibull(x, shape, scale), 0.01, max(yt, na.rm = T), add = T, col = "blue", lwd = 2) hist(y, breaks = 50, freq = F) curve(dexp(x, exp.lambda), 0, pars$Tmax, add = T, col = "red", lwd = 2) curve(dweibull(x, shape, scale), 0.01, pars$Tmax, add = T, col = "blue", lwd = 2) # plot posteriors hist(mod1.jags$sims.list$N, xlim = c(0, n.tree)) abline(v = n.occ, col = "red") hist(mod2.jags$sims.list$N, xlim = c(0, n.tree)) abline(v = n.occ, col = "red") hist(mod3.jags$sims.list$N, xlim = c(0, n.tree)) abline(v = n.occ, col = "red") return(list(sd.lambda = sd.lambda, exp = mod1.jags$sims.list$N, weib = mod2.jags$sims.list$N, exp.mix = mod3.jags$sims.list$N)) } ################################################################################ # function to simulate a population with heterogeneity among trees in hazard sim.het <- function(sd.lambda, n.tree = 319, n.occ = 258, Tmax = 10, mean.lambda = 0.35) { if(sd.lambda >0) { # Convert to shape and scale parameters alpha = (mean.lambda^2) / (sd.lambda^2) beta = (sd.lambda^2) / mean.lambda # generate rate values from a gamma distribution r <- rgamma(n.tree, shape = alpha, rate = 1 / beta) } else r = mean.lambda # hist(r) # Full ttd yt <- rexp(n = n.tree, rate = r) # Set unoccupied trees to NA nt <- 1:length(yt) unocc <- sample(nt, (n.tree - n.occ)) yt[unocc] <- NA # ttd with censoring at Tmax y <- ifelse(yt > Tmax, NA, yt) ttd <- y isCensored <- ifelse(is.na(ttd) == T, 1, 0) z <- ifelse(is.na(ttd) == F, 1, NA) N <- length(ttd) ttdst <- rep(Tmax + 1, N) ttdst[isCensored == 0] <- NA return(list(sd.lambda = sd.lambda, yt = yt, y = y, ttd = ttd, isCensored = isCensored, z = z, N = N, ttdst = ttdst, Tmax = Tmax, n.tree = n.tree, n.occ = n.occ)) } # a <- fit(sim.het(0.4)) ################################################################################ # levels of standard deviation to examine sd_level <- c(0, 0.1, 0.2, 0.3, 0.4) n_level <- length(sd_level) # number of simulations per level of standard deviation n_sim <- 100 out <- list() pb <- winProgressBar(title = "progress bar", min = 0, max = (n_level*n_sim), width = 300) count <- 0 for(i in 1:length(sd_level)) { for(j in 1:n_sim) { count <- count + 1 b <- sim.het(sd_level[i]) a <- fit(b) sm <- length(a$exp) mat <- matrix(nrow = sm, ncol = 5) mat[, 1] <- rep(a$sd.lambda, sm) mat[, 2] <- rep(count, sm) mat[, 3] <- a$exp mat[, 4] <- a$exp.mix mat[, 5] <- a$weib out[[count]] <- mat Sys.sleep(0.1) setWinProgressBar(pb, i, title=paste( round(count/(n_level*n_sim)*100, 0), "% done")) } } close(pb) save.image("c:\\users\\s429217\\onedrive\\data\\koala\\simulation output gamma.RData")