#' GetInits
#'
#' Estimate model coefficients in Presence; save output for mu/se values & inits
#' @param alpha four letter alpha code for species of interest; if NULL, runs Presence model for all species in parallel
#' @export
GetInits <- function(alpha = NULL){
if(is.null(alpha)){
### Read species list
spp_list <- read.csv("inst/spp_list.csv")
### Register cores
cores <- parallel::detectCores()
if(length(spp_list$spp) < cores) cores <- length(spp_list$spp)
doParallel::registerDoParallel(cores = cores)
### Run models in parallel
inits <- foreach::foreach(i = 1:length(spp_list$spp), .combine = c,
.packages = c("dplyr", "BayesCorrOcc")) %dopar%{
## Read count data
dat <- readRDS(paste0("inst/output/", spp_list$spp[i], "/bbs_data.rds"))
## Write pao file
BayesCorrOcc::write_pao(alpha = spp_list$spp[i])
## Read pao file
pao <- RPresence::read.pao(paste0("inst/output/", spp_list$spp[i], "/pres_in.pao"))
## Create design matrices
start_yr <- dat$start_year
n_surv <- pao$nsurveys # total num of surveys
n_seas <- pao$nseasons # num seasons
n_count <- pao$nsurveyseason[1]
years <- seq(from = start_yr, to = start_yr + n_seas - 1)
num.betas <- c(11, # psi
1, #th0
1, #th1
11, #gamma
11, #epsilon
5) #p
# 1st design matrix, 1 row for psi, K rows for theta0, K rows for theta1
dm1 <- matrix('0', n_surv * 2 + 1, sum(num.betas[1:3]))
# psi
dm1[1, 1:num.betas[1]] <- c(1, colnames(pao$unitcov)[grepl(start_yr, colnames(pao$unitcov))])
# th0
dm1[2:(n_surv + 1),(1 + num.betas[1]):sum(num.betas[1:2])] <- rep(1, n_surv)
# th1
dm1[(n_surv + 2):(2 * n_surv + 1), (1 + sum(num.betas[1:2])):sum(num.betas[1:3])] <- rep(1, n_surv)
rownames(dm1) <- c('psi', paste0('th0(', 1:n_surv,')'), paste0('th1(', 1:n_surv,')'))
colnames(dm1) <- paste0('a',1:sum(num.betas[1:3]))
# gam
gam.dm <- c(1, colnames(pao$unitcov)[grepl((start_yr + 1), colnames(pao$unitcov))])
for (ii in 2:(length(years) - 1)) {
gam.dm <- c(gam.dm, c(1, colnames(pao$unitcov)[grepl((start_yr + ii), colnames(pao$unitcov))]))
}
dm2 <- matrix(gam.dm, n_seas - 1, num.betas[4], byrow = T,
dimnames = list(paste0('gam', 1:(n_seas-1)), paste0('b', 1:num.betas[4])))
# eps
dm3 <- matrix(gam.dm, n_seas - 1, num.betas[5], byrow = T,
dimnames = list(paste0('eps', 1:(n_seas-1)), paste0('c', 1:num.betas[4])))
# p
scale.stop <- c("Stop1", "sq_Stop1")
for (ii in 2:n_count) {
scale.stop <- c(scale.stop, paste0("Stop", ii), paste0("sq_Stop", ii))
}
scale.stop <- matrix(scale.stop, nrow = n_surv, ncol = 2, byrow = T)
coord.mat <- matrix(rep(c("Lat", "Lon"), n_surv), nrow = n_surv, byrow = TRUE)
p1.intercept <- rep(1, n_surv)
dm4 <- cbind(p1.intercept, scale.stop, coord.mat)
rownames(dm4) <- c(paste0('p1(',1:n_surv,')')) # could get rid of parentheses
colnames(dm4) <- paste0('d', 1:dim(dm4)[2])
# theta.pi
dm5 <- matrix(0, n_seas, 1, dimnames = list(paste0('thta0pi', 1:n_seas), NULL)) # NOTE THE ZERO when fixing, note no colname
dm_list <- list(dm1 = dm1, dm2 = dm2, dm3 = dm3, dm4 = dm4, dm5 = dm5, dm6 = NULL)
fixedpars <- matrix(rep("eq", pao$nseasons), pao$nseasons, 1)
r1 <- dim(dm_list$dm1)[1] + dim(dm_list$dm2)[1] + dim(dm_list$dm3)[1] + dim(dm_list$dm4)[1]
rownames(fixedpars) <- (r1 + 1):(r1 + pao$nseasons)
## Run model
RPresence::write_dm_and_run(paoname = pao$paoname, fixed = fixedpars,
dms = dm_list, model = "Inits",
noderived = TRUE, limit.real = TRUE,
modname = paste0("inits_", spp_list$spp[i]))
suppressMessages(file.rename(from = paste0("pres_inits_", tolower(spp_list$spp[i]), ".out"),
to = paste0("inst/output/", toupper(spp_list$spp[i]), "/inits.out")))
inits <- suppressWarnings(GetBetas(alpha = spp_list$spp[i]))
saveRDS(inits, paste0("inst/output/", spp_list$spp[i], "/inits.rds"))
return(as.character(spp_list$spp[i]))
}
return(inits)
}else{
## Read count data
dat <- readRDS(paste0("inst/output/", alpha, "/bbs_data.rds"))
## Write pao file
BayesCorrOcc::write_pao(alpha = alpha)
## Read pao file
pao <- RPresence::read.pao(paste0("inst/output/", alpha, "/pres_in.pao"))
## Create design matrices
start_yr <- dat$start_year
n_surv <- pao$nsurveys # total num of surveys
n_seas <- pao$nseasons # num seasons
n_count <- pao$nsurveyseason[1]
years <- seq(from = start_yr, to = start_yr + n_seas - 1)
num.betas <- c(11, # psi
1, #th0
1, #th1
11, #gamma
11, #epsilon
5) #p
# 1st design matrix, 1 row for psi, K rows for theta0, K rows for theta1
dm1 <- matrix('0', n_surv * 2 + 1, sum(num.betas[1:3]))
# psi
dm1[1, 1:num.betas[1]] <- c(1, colnames(pao$unitcov)[grepl(start_yr, colnames(pao$unitcov))])
# th0
dm1[2:(n_surv + 1),(1 + num.betas[1]):sum(num.betas[1:2])] <- rep(1, n_surv)
# th1
dm1[(n_surv + 2):(2 * n_surv + 1), (1 + sum(num.betas[1:2])):sum(num.betas[1:3])] <- rep(1, n_surv)
rownames(dm1) <- c('psi', paste0('th0(', 1:n_surv,')'), paste0('th1(', 1:n_surv,')'))
colnames(dm1) <- paste0('a',1:sum(num.betas[1:3]))
# gam
gam.dm <- c(1, colnames(pao$unitcov)[grepl((start_yr + 1), colnames(pao$unitcov))])
for (ii in 2:(length(years) - 1)) {
gam.dm <- c(gam.dm, c(1, colnames(pao$unitcov)[grepl((start_yr + ii), colnames(pao$unitcov))]))
}
dm2 <- matrix(gam.dm, n_seas - 1, num.betas[4], byrow = T,
dimnames = list(paste0('gam', 1:(n_seas-1)), paste0('b', 1:num.betas[4])))
# eps
dm3 <- matrix(gam.dm, n_seas - 1, num.betas[5], byrow = T,
dimnames = list(paste0('eps', 1:(n_seas-1)), paste0('c', 1:num.betas[4])))
# p
scale.stop <- c("Stop1", "sq_Stop1")
for (ii in 2:n_count) {
scale.stop <- c(scale.stop, paste0("Stop", ii), paste0("sq_Stop", ii))
}
scale.stop <- matrix(scale.stop, nrow = n_surv, ncol = 2, byrow = T)
coord.mat <- matrix(rep(c("Lat", "Lon"), n_surv), nrow = n_surv, byrow = TRUE)
p1.intercept <- rep(1, n_surv)
dm4 <- cbind(p1.intercept, scale.stop, coord.mat)
rownames(dm4) <- c(paste0('p1(',1:n_surv,')')) # could get rid of parentheses
colnames(dm4) <- paste0('d', 1:dim(dm4)[2])
# theta.pi
dm5 <- matrix(0, n_seas, 1, dimnames = list(paste0('thta0pi', 1:n_seas), NULL)) # NOTE THE ZERO when fixing, note no colname
dm_list <- list(dm1 = dm1, dm2 = dm2, dm3 = dm3, dm4 = dm4, dm5 = dm5, dm6 = NULL)
fixedpars <- matrix(rep("eq", pao$nseasons), pao$nseasons, 1)
r1 <- dim(dm_list$dm1)[1] + dim(dm_list$dm2)[1] + dim(dm_list$dm3)[1] + dim(dm_list$dm4)[1]
rownames(fixedpars) <- (r1 + 1):(r1 + pao$nseasons)
## Run model
RPresence::write_dm_and_run(paoname = pao$paoname, fixed = fixedpars,
dms = dm_list, model = "Inits",
noderived = TRUE, limit.real = TRUE,
modname = "inits")
suppressMessages(file.rename(from = "pres_inits.out",
to = paste0("inst/output/", alpha, "/inits.out")))
inits <- suppressWarnings(GetBetas(alpha = alpha))
saveRDS(inits, paste0("inst/output/", alpha, "/inits.rds"))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.