Skip to content

Commit

Permalink
lots of testing
Browse files Browse the repository at this point in the history
  • Loading branch information
@doserjef
doserjef committed Oct 12, 2023
1 parent fcaf9dc commit ac45869
Show file tree
Hide file tree
Showing 52 changed files with 1,777 additions and 2,109 deletions.
35 changes: 22 additions & 13 deletions R/NMix.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,8 @@ NMix <- function(abund.formula, det.formula, data, inits, priors, tuning,
# Make both covariates a data frame. Unlist is necessary for when factors
# are supplied.
data$det.covs <- data.frame(lapply(data$det.covs, function(a) unlist(c(a))))
# Indicator of whether all det.covs are site level or not
site.level.ind <- ifelse(nrow(data$det.covs) == nrow(y), TRUE, FALSE)
data$abund.covs <- as.data.frame(data$abund.covs)

# Check whether random effects are sent in as numeric, and
Expand Down Expand Up @@ -159,21 +161,28 @@ NMix <- function(abund.formula, det.formula, data, inits, priors, tuning,
stop("error: missing values in abund.covs. Please remove these sites from all objects in data or somehow replace the NA values with non-missing values (e.g., mean imputation).")
}
# det.covs ------------------------
for (i in 1:ncol(data$det.covs)) {
if (sum(is.na(data$det.covs[, i])) > sum(is.na(y))) {
stop("error: some elements in det.covs have missing values where there is an observed data value in y. Please either replace the NA values in det.covs with non-missing values (e.g., mean imputation) or set the corresponding values in y to NA where the covariate is missing.")
if (!site.level.ind) {
for (i in 1:ncol(data$det.covs)) {
if (sum(is.na(data$det.covs[, i])) > sum(is.na(y))) {
stop("error: some elements in det.covs have missing values where there is an observed data value in y. Please either replace the NA values in det.covs with non-missing values (e.g., mean imputation) or set the corresponding values in y to NA where the covariate is missing.")
}
}
# Misalignment between y and det.covs
y.missing <- which(is.na(data$y))
det.covs.missing <- lapply(data$det.covs, function(a) which(is.na(a)))
for (i in 1:length(det.covs.missing)) {
tmp.indx <- !(y.missing %in% det.covs.missing[[i]])
if (sum(tmp.indx) > 0) {
if (i == 1 & verbose) {
message("There are missing values in data$y with corresponding non-missing values in data$det.covs.\nRemoving these site/replicate combinations for fitting the model.\n")
}
data$det.covs[y.missing, i] <- NA
}
}
}
# Misalignment between y and det.covs
y.missing <- which(is.na(data$y))
det.covs.missing <- lapply(data$det.covs, function(a) which(is.na(a)))
for (i in 1:length(det.covs.missing)) {
tmp.indx <- !(y.missing %in% det.covs.missing[[i]])
if (sum(tmp.indx) > 0) {
if (i == 1 & verbose) {
message("There are missing values in data$y with corresponding non-missing values in data$det.covs.\nRemoving these site/replicate combinations for fitting the model.\n")
}
data$det.covs[y.missing, i] <- NA
if (site.level.ind) {
if (sum(is.na(data$det.covs)) != 0) {
stop("missing values in site-level det.covs. Please remove these sites from all objects in data or somehow replace the NA values with non-missing values (e.g., mean imputation).")
}
}

Expand Down
4 changes: 2 additions & 2 deletions R/abund.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ abund <- function(formula, data, inits, priors, tuning,
if (family %in% c('Gaussian', 'zi-Gaussian')) {
abundGaussian(formula, data, inits, priors, tuning,
n.batch, batch.length, accept.rate, family,
n.omp.threads, verbose, n.report, n.burn, n.thin, n.chains)
n.omp.threads, verbose, n.report, n.burn, n.thin, n.chains, save.fitted)
} else {

# Make it look nice
Expand Down Expand Up @@ -114,7 +114,7 @@ abund <- function(formula, data, inits, priors, tuning,
data$covs <- data.frame(lapply(data$covs, function(a) unlist(c(a))))
# Check if only site-level covariates are included
if (nrow(data$covs) == dim(y)[1]) {
data$covs <- as.data.frame(mapply(rep, data$covs, dim(y)[2]))
data$covs <- as.data.frame(lapply(data$covs, rep, dim(y)[2]))
}

# Check whether random effects are sent in as numeric, and
Expand Down
37 changes: 22 additions & 15 deletions R/abundGaussian.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ abundGaussian <- function(formula, data, inits, priors, tuning, n.batch,
batch.length, accept.rate = 0.43, family = 'Gaussian',
n.omp.threads = 1, verbose = TRUE, n.report = 100,
n.burn = round(.10 * n.batch * batch.length),
n.thin = 1, n.chains = 1, ...){
n.thin = 1, n.chains = 1, save.fitted = TRUE, ...){

ptm <- proc.time()

Expand Down Expand Up @@ -136,6 +136,11 @@ abundGaussian <- function(formula, data, inits, priors, tuning, n.batch,
}
}

# Check save.fitted ---------------------------------------------------
if (!(save.fitted %in% c(TRUE, FALSE))) {
stop("save.fitted must be either TRUE or FALSE")
}

# Formula -------------------------------------------------------------
# Occupancy -----------------------
if (is(formula, 'formula')) {
Expand Down Expand Up @@ -393,7 +398,7 @@ abundGaussian <- function(formula, data, inits, priors, tuning, n.batch,
# Set storage for all variables ---------------------------------------
storage.mode(y) <- "double"
storage.mode(X) <- "double"
consts <- c(J.est, p, p.re, n.re, J.zero)
consts <- c(J.est, p, p.re, n.re, J.zero, save.fitted)
storage.mode(consts) <- "integer"
storage.mode(beta.inits) <- "double"
storage.mode(tau.sq.inits) <- "double"
Expand Down Expand Up @@ -485,22 +490,24 @@ abundGaussian <- function(formula, data, inits, priors, tuning, n.batch,
out$tau.sq.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$tau.sq.samples))))
colnames(out$tau.sq.samples) <- 'tau.sq'
# Get everything back in the original order
if (!two.stage) {
out$y.rep.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.samples))))
out$like.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$like.samples))))
} else {
y.rep.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.samples))))
y.rep.zero.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.zero.samples))))
out$y.rep.samples <- matrix(NA, n.post.samples * n.chains, J.est + J.zero)
out$y.rep.samples[, z.indx] <- y.rep.samples
out$y.rep.samples[, -z.indx] <- y.rep.zero.samples
out$y.rep.samples <- mcmc(out$y.rep.samples)
out$like.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$like.samples))))
if (save.fitted) {
if (!two.stage) {
out$y.rep.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.samples))))
out$like.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$like.samples))))
} else {
y.rep.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.samples))))
y.rep.zero.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$y.rep.zero.samples))))
out$y.rep.samples <- matrix(NA, n.post.samples * n.chains, J.est + J.zero)
out$y.rep.samples[, z.indx] <- y.rep.samples
out$y.rep.samples[, -z.indx] <- y.rep.zero.samples
out$y.rep.samples <- mcmc(out$y.rep.samples)
out$like.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$like.samples))))
}
out$mu.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$mu.samples))))
out$mu.samples <- mcmc(out$mu.samples)
}
out$X <- X
out$X.re <- X.re
out$mu.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$mu.samples))))
out$mu.samples <- mcmc(out$mu.samples)
out$y <- y.orig
if (p.re > 0) {
out$sigma.sq.mu.samples <- mcmc(
Expand Down
Loading

0 comments on commit ac45869

Please sign in to comment.