some stuff

NEWS.md

@@ -1,12 +1,13 @@
 # spAbundance 0.1.4
 
++ Added functionality for independent priors on the species-specific effects to allow species-effects to be treated as fixed effects as opposed to random effects from a common community-level distribution for the following model types: `msNMix()`.  
 + Added in a check at the top of all model fitting functions to return an error when the number of posterior samples saved based on the MCMC criteria (`n.batch`, `batch.length`, `n.samples`, `n.burn`, `n.thin`, `n.chains`) are specified in a way that leads to a non-integer value. In such situations, models would previously run and return without an error, but sometimes the last posterior sample in any given chain could have widely inaccurate values, or values that prevented subsequent functions from working. Thanks to Wendy Leuenberger for bringing this to my attention. 
 + Updated C++ code to adhere to the new lack of re-mapping of functions in Rinternals.h and R_ext/Error.h when building packages on CRAN. 
 
 
 # spAbundance 0.1.3
 
-+ Added in the `independent.betas` tag into the `priors` list for certain multi-species model types to allow for specifying an independent prior on the species-specific random effects as opposed to treating species-specific effects as random effects. This can be useful under certain circumstances when the distribution of effects across species may not be adequately represented by a Gaussian distribution. This tag is available for the following functions: `lfMsAbund` (Gaussian only), `sfMsAbund` (Gaussian only), and `svcMsAbund`.
++ Added in the `independent.betas` tag into the `priors` list for certain multi-species model types to allow for specifying an independent prior on the species-specific effects as opposed to treating species-specific effects as random effects. This can be useful under certain circumstances when the distribution of effects across species may not be adequately represented by a Gaussian distribution. This tag is available for the following functions: `lfMsAbund` (Gaussian only), `sfMsAbund` (Gaussian only), and `svcMsAbund`.
 + Fixed a bug in zero-inflated Gaussian latent factor abundance models (`lfMsAbund`). 
 + Fixed a bug in `waicAbund()` that prevented it from working with `svcMsAbund` models.
 + Fixed minor issues in C++ to pass CRAN additional checks.

R/msNMix.R

@@ -299,6 +299,29 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
   }
   names(priors) <- tolower(names(priors))
 
+  # Independent beta parameters -----
+  if ('independent.betas' %in% names(priors)) {
+    if (priors$independent.betas == TRUE) {
+      message("Beta parameters will be estimated independently\n")
+      ind.betas <- TRUE
+    } else if (priors$independent.betas == FALSE) {
+      ind.betas <- FALSE
+    }
+  } else {
+    ind.betas <- FALSE
+  }
+  # Independent alpha parameters -----
+  if ('independent.alphas' %in% names(priors)) {
+    if (priors$independent.alphas == TRUE) {
+      message("Alpha parameters will be estimated independently\n")
+      ind.alphas <- TRUE
+    } else if (priors$independent.alphas == FALSE) {
+      ind.alphas <- FALSE
+    }
+  } else {
+    ind.alphas <- FALSE
+  }
+
   # beta.comm -----------------------
   if ("beta.comm.normal" %in% names(priors)) {
     if (!is.list(priors$beta.comm.normal) | length(priors$beta.comm.normal) != 2) {
@@ -332,7 +355,7 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
     }
     Sigma.beta.comm <- sigma.beta.comm * diag(p.abund)
   } else {
-    if (verbose) {
+    if (verbose & !ind.betas) {
       message("No prior specified for beta.comm.normal.\nSetting prior mean to 0 and prior variance to 100\n")
     }
     mu.beta.comm <- rep(0, p.abund)
@@ -373,7 +396,7 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
     }
     Sigma.alpha.comm <- sigma.alpha.comm * diag(p.det)
   } else {
-    if (verbose) {
+    if (verbose & !ind.alphas) {
       message("No prior specified for alpha.comm.normal.\nSetting prior mean to 0 and prior variance to 2.72\n")
     }
     mu.alpha.comm <- rep(0, p.det)
@@ -413,7 +436,7 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
       tau.sq.beta.b <- rep(tau.sq.beta.b, p.abund)
     }
   } else {
-    if (verbose) {
+    if (verbose & !ind.betas) {
       message("No prior specified for tau.sq.beta.ig.\nSetting prior shape to 0.1 and prior scale to 0.1\n")
     }
     tau.sq.beta.a <- rep(0.1, p.abund)
@@ -452,7 +475,7 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
       tau.sq.alpha.b <- rep(tau.sq.alpha.b, p.det)
     }
   } else {
-    if (verbose) {
+    if (verbose & !ind.alphas) {
       message("No prior specified for tau.sq.alpha.ig.\nSetting prior shape to 0.1 and prior scale to 0.1\n")
     }
     tau.sq.alpha.a <- rep(0.1, p.det)
@@ -939,7 +962,8 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
   storage.mode(K) <- "double"
   storage.mode(offset) <- "double"
   storage.mode(y.max) <- "double"
-  consts <- c(n.sp, J, n.obs, p.abund, p.abund.re, n.abund.re, p.det, p.det.re, n.det.re)
+  consts <- c(n.sp, J, n.obs, p.abund, p.abund.re, n.abund.re, p.det, p.det.re, n.det.re,
+              ind.betas, ind.alphas)
   storage.mode(consts) <- "integer"
   storage.mode(beta.inits) <- "double"
   storage.mode(alpha.inits) <- "double"
@@ -1005,10 +1029,14 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
   for (i in 1:n.chains) {
     # Change initial values if i > 1
     if ((i > 1) & (!fix.inits)) {
-      beta.comm.inits <- rnorm(p.abund, 0, 1)
-      alpha.comm.inits <- rnorm(p.det, 0, 1)
-      tau.sq.beta.inits <- runif(p.abund, 0.05, 1)
-      tau.sq.alpha.inits <- runif(p.det, 0.05, 1)
+      if (!ind.betas) {
+        beta.comm.inits <- rnorm(p.abund, 0, 1)
+        tau.sq.beta.inits <- runif(p.abund, 0.05, 1)
+      }
+      if (!ind.alphas) {
+        alpha.comm.inits <- rnorm(p.det, 0, 1)
+        tau.sq.alpha.inits <- runif(p.det, 0.05, 1)
+      }
       beta.inits <- matrix(rnorm(n.sp * p.abund, beta.comm.inits,
             		     sqrt(tau.sq.beta.inits)), n.sp, p.abund)
       alpha.inits <- matrix(rnorm(n.sp * p.det, alpha.comm.inits,
@@ -1030,20 +1058,20 @@ msNMix <- function(abund.formula, det.formula, data, inits, priors,
 
     storage.mode(chain.info) <- "integer"
     out.tmp[[i]] <- .Call("msNMix", y, X, X.p, X.re, X.p.re,
-        		  X.random, X.p.random, y.max, consts,
-    	                  n.abund.re.long, n.det.re.long, beta.inits,
-        		  alpha.inits, kappa.inits, N.inits, beta.comm.inits,
+                          X.random, X.p.random, y.max, consts,
+                          n.abund.re.long, n.det.re.long, beta.inits,
+                          alpha.inits, kappa.inits, N.inits, beta.comm.inits,
       	                  alpha.comm.inits, tau.sq.beta.inits, tau.sq.alpha.inits,
-      		          sigma.sq.mu.inits, sigma.sq.p.inits,
-    	                  beta.star.inits, alpha.star.inits, N.long.indx,
-      		          beta.star.indx, beta.level.indx, alpha.star.indx,
-      		          alpha.level.indx, mu.beta.comm, mu.alpha.comm,
-      		          Sigma.beta.comm, Sigma.alpha.comm, kappa.a,
-        		  kappa.b, tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a,
+                          sigma.sq.mu.inits, sigma.sq.p.inits,
+    	                    beta.star.inits, alpha.star.inits, N.long.indx,
+                          beta.star.indx, beta.level.indx, alpha.star.indx,
+                          alpha.level.indx, mu.beta.comm, mu.alpha.comm,
+                          Sigma.beta.comm, Sigma.alpha.comm, kappa.a,
+                          kappa.b, tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a,
       	                  tau.sq.alpha.b, sigma.sq.mu.a, sigma.sq.mu.b,
-      		          sigma.sq.p.a, sigma.sq.p.b,
-    	                  tuning.c, n.batch, batch.length, accept.rate, n.omp.threads,
-    	                  verbose, n.report, samples.info, chain.info, family.c, offset)
+                          sigma.sq.p.a, sigma.sq.p.b,
+                          tuning.c, n.batch, batch.length, accept.rate, n.omp.threads,
+                          verbose, n.report, samples.info, chain.info, family.c, offset)
     chain.info[1] <- chain.info[1] + 1
   }
   # Calculate R-Hat ---------------

R/simNMix.R

@@ -1,6 +1,6 @@
-simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha, 
-		    kappa, mu.RE = list(), p.RE = list(), 
-		    offset = 1, sp = FALSE, cov.model, sigma.sq, 
+simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
+		    kappa, mu.RE = list(), p.RE = list(),
+		    offset = 1, sp = FALSE, cov.model, sigma.sq,
 		    phi, nu, family = 'Poisson', ...) {
 
   # Check for unused arguments ------------------------------------------
@@ -95,7 +95,7 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
     }
     cov.model.names <- c("exponential", "spherical", "matern", "gaussian")
     if(! cov.model %in% cov.model.names){
-      stop("error: specified cov.model '",cov.model,"' is not a valid option; choose from ", 
+      stop("error: specified cov.model '",cov.model,"' is not a valid option; choose from ",
            paste(cov.model.names, collapse=", ", sep="") ,".")
     }
     if (cov.model == 'matern' & missing(nu)) {
@@ -128,7 +128,7 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
   n.alpha <- length(alpha)
   X.p <- array(NA, dim = c(J, n.rep.max, n.alpha))
   X.p[, , 1] <- 1
-  # Get index of surveyed replicates for each site. 
+  # Get index of surveyed replicates for each site.
   rep.indx <- list()
   for (j in 1:J) {
     rep.indx[[j]] <- sample(1:n.rep.max, n.rep[j], replace = FALSE)
@@ -174,11 +174,11 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
     n.random <- ncol(X.random)
     X.re <- matrix(NA, J, length(mu.RE$levels))
     for (i in 1:length(mu.RE$levels)) {
-      X.re[, i] <- sample(1:mu.RE$levels[i], J, replace = TRUE)  
+      X.re[, i] <- sample(1:mu.RE$levels[i], J, replace = TRUE)
     }
     indx.mat <- X.re[, re.col.indx, drop = FALSE]
     for (i in 1:p.nmix.re) {
-      beta.star[which(beta.star.indx == i)] <- rnorm(n.nmix.re.long[i], 0, 
+      beta.star[which(beta.star.indx == i)] <- rnorm(n.nmix.re.long[i], 0,
 						     sqrt(sigma.sq.mu[i]))
     }
     if (length(mu.RE$levels) > 1) {
@@ -212,8 +212,8 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
     X.p.random <- X.p[, , unlist(p.RE$alpha.indx), drop = FALSE]
     X.p.re <- array(NA, dim = c(J, n.rep.max, length(p.RE$levels)))
     for (i in 1:length(p.RE$levels)) {
-      X.p.re[, , i] <- matrix(sample(1:p.RE$levels[i], J * n.rep.max, replace = TRUE), 
-		              J, n.rep.max)	      
+      X.p.re[, , i] <- matrix(sample(1:p.RE$levels[i], J * n.rep.max, replace = TRUE),
+		              J, n.rep.max)
     }
     for (i in 1:p.det.re) {
       alpha.star[which(alpha.star.indx == i)] <- rnorm(n.det.re.long[i], 0, sqrt(sigma.sq.p[i]))
@@ -229,7 +229,7 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
     }
     if (p.det.re > 1) {
       for (j in 2:p.det.re) {
-        indx.mat[, , j] <- indx.mat[, , j] + max(indx.mat[, , j - 1], na.rm = TRUE) 
+        indx.mat[, , j] <- indx.mat[, , j] + max(indx.mat[, , j - 1], na.rm = TRUE)
       }
     }
     alpha.star.sites <- matrix(NA, J, n.rep.max)
@@ -258,7 +258,7 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
         mu <- exp(X %*% as.matrix(beta))
       }
     }
-    # Get mean and overdispersion parameter 
+    # Get mean and overdispersion parameter
     N <- rnbinom(J, size = kappa, mu = mu * offset)
   } else if (family == 'Poisson') {
     if (sp) {
@@ -282,7 +282,7 @@ simNMix <- function(J.x, J.y, n.rep, n.rep.max, beta, alpha,
   y <- matrix(NA, nrow = J, ncol = n.rep.max)
   for (j in 1:J) {
     if (length(p.RE) > 0) {
-      p[j, rep.indx[[j]]] <- logit.inv(X.p[j, rep.indx[[j]], ] %*% as.matrix(alpha) + 
+      p[j, rep.indx[[j]]] <- logit.inv(X.p[j, rep.indx[[j]], ] %*% as.matrix(alpha) +
 				    alpha.star.sites[j, rep.indx[[j]]])
     } else {
       p[j, rep.indx[[j]]] <- logit.inv(X.p[j, rep.indx[[j]], ] %*% as.matrix(alpha))