diff --git a/NEWS.md b/NEWS.md
index 6d503e0..0c1fc03 100644
--- a/NEWS.md
+++ b/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.
diff --git a/R/msNMix.R b/R/msNMix.R
index 611c7a4..15d8ffd 100644
--- a/R/msNMix.R
+++ b/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 ---------------
diff --git a/R/simNMix.R b/R/simNMix.R
index 0718ecb..cb04f12 100644
--- a/R/simNMix.R
+++ b/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))
diff --git a/R/simTNMix.R b/R/simTNMix.R
new file mode 100644
index 0000000..3d138d7
--- /dev/null
+++ b/R/simTNMix.R
@@ -0,0 +1,378 @@
+simTNMix <- function(J.x, J.y, n.time, n.rep, n.rep.max, beta, alpha, sp.only = 0,
+                     trend = TRUE, kappa, mu.RE = list(), p.RE = list(),
+                     offset = 1, sp = FALSE, cov.model, sigma.sq, phi,
+                     nu, family = 'Poisson', ar1 = FALSE, rho, sigma.sq.t, ...) {
+
+  # Check for unused arguments ------------------------------------------
+  formal.args <- names(formals(sys.function(sys.parent())))
+  elip.args <- names(list(...))
+  for(i in elip.args){
+      if(! i %in% formal.args)
+          warning("'",i, "' is not an argument")
+  }
+
+  # Check function inputs -------------------------------------------------
+  # J.x -------------------------------
+  if (missing(J.x)) {
+    stop("error: J.x must be specified")
+  }
+  if (length(J.x) != 1) {
+    stop("error: J.x must be a single numeric value.")
+  }
+  # J.y -------------------------------
+  if (missing(J.y)) {
+    stop("error: J.y must be specified")
+  }
+  if (length(J.y) != 1) {
+    stop("error: J.y must be a single numeric value.")
+  }
+  J <- J.x * J.y
+  # n.time ---------------------------
+  if (missing(n.time)) {
+    stop("error: n.time must be specified.")
+  }
+  # n.rep -----------------------------
+  if (missing(n.rep)) {
+    stop("error: n.rep must be specified.")
+  }
+  if (!is.matrix(n.rep)) {
+    stop(paste("error: n.rep must be a matrix with ", J, " rows and ", max(n.time), " columns", sep = ''))
+  }
+  if (missing(n.rep.max)) {
+    n.rep.max <- max(n.rep, na.rm = TRUE)
+  }
+  if (nrow(n.rep) != J | ncol(n.rep) != max(n.time)) {
+    stop(paste("error: n.rep must be a matrix with ", J, " rows and ", max(n.time), " columns", sep = ''))
+  }
+  # beta ------------------------------
+  if (missing(beta)) {
+    stop("error: beta must be specified.")
+    if (length(beta) <= 1) {
+      stop("error: beta must have at least two elements (intercept and trend)")
+    }
+  }
+  # alpha -----------------------------
+  if (missing(alpha)) {
+    stop("error: alpha must be specified.")
+  }
+  # family ------------------------------
+  if (! (family %in% c('NB', 'Poisson'))) {
+    stop("error: family must be either NB (negative binomial) or Poisson")
+  }
+  # kappa -----------------------------
+  if (family == 'NB') {
+    if (missing(kappa)) {
+      stop("error: kappa (overdispersion parameter) must be specified when family = 'NB'.")
+    }
+  }
+  if (family == 'Poisson' & !missing(kappa)) {
+    message("overdispersion parameter (kappa) is ignored when family == 'Poisson'")
+  }
+  # mu.RE ----------------------------
+  names(mu.RE) <- tolower(names(mu.RE))
+  if (!is.list(mu.RE)) {
+    stop("error: if specified, mu.RE must be a list with tags 'levels' and 'sigma.sq.mu'")
+  }
+  if (length(names(mu.RE)) > 0) {
+    if (!'sigma.sq.mu' %in% names(mu.RE)) {
+      stop("error: sigma.sq.mu must be a tag in mu.RE with values for the abundance random effect variances")
+    }
+    if (!'levels' %in% names(mu.RE)) {
+      stop("error: levels must be a tag in mu.RE with the number of random effect levels for each abundance random intercept.")
+    }
+    if (!'beta.indx' %in% names(mu.RE)) {
+      mu.RE$beta.indx <- list()
+      for (i in 1:length(mu.RE$sigma.sq.mu)) {
+        mu.RE$beta.indx[[i]] <- 1
+      }
+    }
+  }
+  # p.RE ----------------------------
+  names(p.RE) <- tolower(names(p.RE))
+  if (!is.list(p.RE)) {
+    stop("error: if specified, p.RE must be a list with tags 'levels' and 'sigma.sq.p'")
+  }
+  if (length(names(p.RE)) > 0) {
+    if (!'sigma.sq.p' %in% names(p.RE)) {
+      stop("error: sigma.sq.p must be a tag in p.RE with values for the detection random effect variances")
+    }
+    if (!'levels' %in% names(p.RE)) {
+      stop("error: levels must be a tag in p.RE with the number of random effect levels for each detection random intercept.")
+    }
+    if (!'alpha.indx' %in% names(p.RE)) {
+      p.RE$alpha.indx <- list()
+      for (i in 1:length(p.RE$sigma.sq.p)) {
+        p.RE$alpha.indx[[i]] <- 1
+      }
+    }
+  }
+  # Spatial parameters ----------------
+  if (sp) {
+    if(missing(sigma.sq)) {
+      stop("error: sigma.sq must be specified when sp = TRUE")
+    }
+    if(missing(phi)) {
+      stop("error: phi must be specified when sp = TRUE")
+    }
+    if(missing(cov.model)) {
+      stop("error: cov.model must be specified when sp = TRUE")
+    }
+    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 ",
+           paste(cov.model.names, collapse=", ", sep="") ,".")
+    }
+    if (cov.model == 'matern' & missing(nu)) {
+      stop("error: nu must be specified when cov.model = 'matern'")
+    }
+  }
+  # AR1 -------------------------------
+  if (ar1) {
+    if (missing(rho)) {
+      stop("error: rho must be specified when ar1 = TRUE")
+    }
+    if (missing(sigma.sq.t)) {
+      stop("error: sigma.sq.t must be specified when ar1 = TRUE")
+    }
+  }
+
+  # Subroutines -----------------------------------------------------------
+  # MVN
+  rmvn <- function(n, mu=0, V = matrix(1)) {
+    p <- length(mu)
+    if(any(is.na(match(dim(V),p))))
+      stop("Dimension problem!")
+    D <- chol(V)
+    t(matrix(rnorm(n*p), ncol=p)%*%D + rep(mu,rep(n,p)))
+  }
+  logit <- function(theta, a = 0, b = 1){log((theta-a)/(b-theta))}
+  logit.inv <- function(z, a = 0, b = 1){b-(b-a)/(1+exp(z))}
+
+  # Form abundance covariates if any --------------------------------------
+  n.beta <- length(beta)
+  n.time.max <- max(n.time, na.rm = TRUE)
+  time.indx <- list()
+  for (j in 1:J) {
+    time.indx[[j]] <- sample(which(!is.na(n.rep[j, ])), n.time[j], replace = FALSE)
+  }
+  X <- array(NA, dim = c(J, n.time.max, n.beta))
+  X[, , 1] <- 1
+  if (n.beta > 1) {
+    if (trend) { # If simulating data with a trend
+      # By default the second simulated covariate is a standardized trend
+      X[, , 2] <- scale(c(matrix(rep(1:n.time.max, each = J), nrow = J, ncol = n.time.max)))
+      if (n.beta > 2) {
+        for (i in 3:n.beta) {
+          if (i %in% sp.only) {
+            X[, , i] <- rep(rnorm(J), n.time.max)
+          } else {
+            X[, , i] <- rnorm(J * n.time.max)
+          }
+        }
+      }
+    } else { # If not simulating data with a trend
+      if (n.beta > 1) {
+        for (i in 2:n.beta) {
+          if (i %in% sp.only) {
+            X[, , i] <- rep(rnorm(J), n.time.max)
+          } else {
+            X[, , i] <- rnorm(J * n.time.max)
+          }
+        }
+      }
+    }
+  }
+  # Form detection covariates (if any) ------------------------------------
+  # Time dependent --------------------
+  rep.indx <- list()
+  for (j in 1:J) {
+    rep.indx[[j]] <- list()
+    for (t in time.indx[[j]]) {
+      rep.indx[[j]][[t]] <- sample(1:n.rep.max, n.rep[j, t], replace = FALSE)
+    }
+  }
+  n.alpha <- length(alpha)
+  X.p <- array(NA, dim = c(J, n.time.max, n.rep.max, n.alpha))
+  X.p[, , , 1] <- 1
+  if (n.alpha > 1) {
+    for (j in 1:J) {
+      for (t in time.indx[[j]]) {
+        for (k in rep.indx[[j]][[t]]) {
+          X.p[j, t, k, 2:n.alpha] <- rnorm(n.alpha - 1)
+        } # k
+      } # t
+    } # j
+  }
+
+  # Simulate spatial random effect ----------------------------------------
+  # Matrix of spatial locations
+  s.x <- seq(0, 1, length.out = J.x)
+  s.y <- seq(0, 1, length.out = J.y)
+  coords <- as.matrix(expand.grid(s.x, s.y))
+  if (sp) {
+    if (cov.model == 'matern') {
+      theta <- c(phi, nu)
+    } else {
+      theta <- phi
+    }
+    Sigma <- mkSpCov(coords, as.matrix(sigma.sq), as.matrix(0), theta, cov.model)
+    # Random spatial process
+    w <- rmvn(1, rep(0, J), Sigma)
+  } else {
+    w <- NA
+  }
+
+  # Random effects --------------------------------------------------------
+  # Abundance -------------------------
+  if (length(mu.RE) > 0) {
+    p.nmix.re <- length(unlist(mu.RE$beta.indx))
+    tmp <- sapply(mu.RE$beta.indx, length)
+    re.col.indx <- unlist(lapply(1:length(mu.RE$beta.indx), function(a) rep(a, tmp[a])))
+    sigma.sq.mu <- mu.RE$sigma.sq.mu[re.col.indx]
+    n.nmix.re.long <- mu.RE$levels[re.col.indx]
+    n.nmix.re <- sum(n.nmix.re.long)
+    beta.star.indx <- rep(1:p.nmix.re, n.nmix.re.long)
+    beta.star <- rep(0, n.nmix.re)
+    X.random <- X[, , unlist(mu.RE$beta.indx), drop = FALSE]
+    n.random <- dim(X.random)[3]
+    X.re <- array(NA, dim = c(J, n.time.max, length(mu.RE$levels)))
+    for (i in 1:length(mu.RE$levels)) {
+      X.re[, , i] <- sample(1:mu.RE$levels[i], J * n.time.max, 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,
+                                                     sqrt(sigma.sq.mu[i]))
+    }
+    if (length(mu.RE$levels) > 1) {
+      for (j in 2:length(mu.RE$levels)) {
+        X.re[, , j] <- X.re[, , j] + max(X.re[, , j - 1], na.rm = TRUE)
+      }
+    }
+    if (p.nmix.re > 1) {
+      for (j in 2:p.nmix.re) {
+        indx.mat[, , j] <- indx.mat[, , j] + max(indx.mat[, , j - 1], na.rm = TRUE)
+      }
+    }
+
+    beta.star.sites <- matrix(NA, J, n.time.max)
+    for (j in 1:J) {
+      for (t in 1:n.time.max) {
+        beta.star.sites[j, t] <- beta.star[indx.mat[j, t, ]] %*% as.matrix(X.random[j, t, ])
+      } # k
+    } # j
+  } else {
+    X.re <- NA
+    beta.star <- NA
+  }
+  # Detection -------------------------
+  if (length(p.RE) > 0) {
+    p.det.re <- length(unlist(p.RE$alpha.indx))
+    tmp <- sapply(p.RE$alpha.indx, length)
+    p.re.col.indx <- unlist(lapply(1:length(p.RE$alpha.indx), function(a) rep(a, tmp[a])))
+    sigma.sq.p <- p.RE$sigma.sq.p[p.re.col.indx]
+    n.det.re.long <- p.RE$levels[p.re.col.indx]
+    n.det.re <- sum(n.det.re.long)
+    alpha.star.indx <- rep(1:p.det.re, n.det.re.long)
+    alpha.star <- rep(0, n.det.re)
+    X.p.random <- X.p[, , , unlist(p.RE$alpha.indx), drop = FALSE]
+    X.p.re <- array(NA, dim = c(J, n.time.max, n.rep.max, length(p.RE$levels)))
+    for (i in 1:length(p.RE$levels)) {
+      X.p.re[, , , i] <- array(sample(1:p.RE$levels[i], J * n.rep.max * n.time.max, replace = TRUE),
+                               dim = c(J, n.time.max, 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]))
+    }
+    for (j in 1:J) {
+      for (t in time.indx[[j]]) {
+        X.p.re[j, t, -rep.indx[[j]][[t]], ] <- NA
+      }
+    }
+    indx.mat <- X.p.re[, , , p.re.col.indx, drop = FALSE]
+    if (length(p.RE$levels) > 1) {
+      for (j in 2:length(p.RE$levels)) {
+        X.p.re[, , , j] <- X.p.re[, , , j] + max(X.p.re[, , , j - 1], na.rm = TRUE)
+      }
+    }
+    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)
+      }
+    }
+    alpha.star.sites <- array(NA, dim = c(J, n.time.max, n.rep.max))
+    for (j in 1:J) {
+      for (t in time.indx[[j]]) {
+        for (k in rep.indx[[j]][[t]]) {
+          alpha.star.sites[j, t, k] <- alpha.star[indx.mat[j, t, k, ]] %*% X.p.random[j, t, k, ]
+        }
+      }
+    }
+  } else {
+    X.p.re <- NA
+    alpha.star <- NA
+  }
+
+  # Simulate temporal (AR1) random effect ---------------------------------
+  if (ar1) {
+    exponent <- abs(matrix(1:n.time.max - 1, nrow = n.time.max,
+                           ncol = n.time.max, byrow = TRUE) - (1:n.time.max - 1))
+    Sigma.eta <- sigma.sq.t * rho^exponent
+    eta <- rmvn(1, rep(0, n.time.max), Sigma.eta)
+  } else {
+    eta <- matrix(rep(0, n.time.max))
+  }
+
+  # Latent abundance process ----------------------------------------------
+  mu <- matrix(NA, J, n.time.max)
+  N <- matrix(NA, J, n.time.max)
+  for (j in 1:J) {
+    for (t in 1:n.time.max) {
+      if (sp) {
+        if (length(mu.RE) > 0) {
+          mu[j, t] <- exp(X[j, t, ] %*% as.matrix(beta) + w[j] +
+                          beta.star.sites[j, t] + eta[t])
+        } else {
+          mu[j, t] <- exp(X[j, t, ] %*% as.matrix(beta) + w[j] + eta[t])
+        }
+      } else {
+        if (length(mu.RE) > 0) {
+          mu[j, t] <- exp(X[j, t, ] %*% as.matrix(beta) +
+                          beta.star.sites[j, t] + eta[t])
+        } else {
+          mu[j, t] <- exp(X[j, t, ] %*% as.matrix(beta) + eta[t])
+        }
+      }
+      if (family == 'NB') {
+        # Get mean and overdispersion parameter
+        N[j, t] <- rnbinom(1, size = kappa, mu = mu[j, t] * offset)
+      } else if (family == 'Poisson') {
+        N[j, t] <- rpois(1, lambda = mu[j, t] * offset)
+      }
+    } # t
+  } # j
+
+  # Data Formation --------------------------------------------------------
+  p <- array(NA, dim = c(J, n.time.max, n.rep.max))
+  y <- array(NA, dim = c(J, n.time.max, n.rep.max))
+  for (j in 1:J) {
+    for (t in time.indx[[j]]) {
+      if (length(p.RE) > 0) {
+        p[j, t, rep.indx[[j]][[t]]] <- logit.inv(X.p[j, t, rep.indx[[j]][[t]], ]
+                                                 %*% as.matrix(alpha) +
+                                                 alpha.star.sites[j, t, rep.indx[[j]][[t]]])
+      } else {
+        p[j, t, rep.indx[[j]][[t]]] <- logit.inv(X.p[j, t, rep.indx[[j]][[t]], ]
+                                                 %*% as.matrix(alpha))
+      }
+      y[j, t, rep.indx[[j]][[t]]] <- rbinom(n.rep[j, t], N[j, t], p[j, t, rep.indx[[j]][[t]]])
+    } # t
+  } # j
+
+  # Return list -----------------------------------------------------------
+  return(
+    list(X = X, X.p = X.p, coords = coords, w = w, mu = mu, N = N,
+         y = y, X.re = X.re, X.p.re = X.p.re, beta.star = beta.star,  p = p,
+         alpha.star = alpha.star, eta = eta)
+  )
+}
diff --git a/src/abundGaussian.cpp b/src/abundGaussian.cpp
index 35c7521..2e0dcd6 100644
--- a/src/abundGaussian.cpp
+++ b/src/abundGaussian.cpp
@@ -298,9 +298,9 @@ extern "C" {
             for (j = 0; j < J; j++) {
               if (XRE[betaStarIndx[l] * J + j] == betaLevelIndx[l]) {
                 tmp_02 = 0.0;
-		for (ll = 0; ll < pRE; ll++) {
+                for (ll = 0; ll < pRE; ll++) {
                   tmp_02 += betaStar[betaStarLongIndx[ll * J + j]] * XRandom[ll * J + j];
-		}
+                }
                 tmp_one[0] += XRandom[betaStarIndx[l] * J + j] * (y[j] - F77_NAME(ddot)(&p, &X[j], &J, beta, &inc) -
           		    tmp_02 + (betaStar[l] * XRandom[betaStarIndx[l] * J + j])) / tauSq;
                 tmp_0 += XRandom[betaStarIndx[l] * J + j] * XRandom[betaStarIndx[l] * J + j] / tauSq;
diff --git a/src/msNMix.cpp b/src/msNMix.cpp
index 6bf8875..41f85b7 100644
--- a/src/msNMix.cpp
+++ b/src/msNMix.cpp
@@ -71,6 +71,8 @@ extern "C" {
     int pDet = INTEGER(consts_r)[6];
     int pDetRE = INTEGER(consts_r)[7];
     int nDetRE = INTEGER(consts_r)[8];
+    int indBetas = INTEGER(consts_r)[9];
+    int indAlphas = INTEGER(consts_r)[10];
     int ppDet = pDet * pDet;
     int ppAbund = pAbund * pAbund;
     double *muBetaComm = REAL(muBetaComm_r);
@@ -451,95 +453,103 @@ extern "C" {
         /********************************************************************
          Update Community level Abundance Coefficients
          *******************************************************************/
-        /********************************
-         Compute b.beta.comm
-         *******************************/
-        zeros(tmp_pAbund, pAbund);
-        for (i = 0; i < nSp; i++) {
-          F77_NAME(dgemv)(ytran, &pAbund, &pAbund, &one, TauBetaInv, &pAbund, &beta[i], &nSp, &one, tmp_pAbund, &inc FCONE);
-        } // i
-        for (q = 0; q < pAbund; q++) {
-          tmp_pAbund[q] += SigmaBetaCommInvMuBeta[q];
-        } // j
-
-        /********************************
-         Compute A.beta.comm
-         *******************************/
-        for (q = 0; q < ppAbund; q++) {
-          tmp_ppAbund[q] = SigmaBetaCommInv[q] + nSp * TauBetaInv[q];
+        if (indBetas == 0) {
+          /********************************
+           Compute b.beta.comm
+           *******************************/
+          zeros(tmp_pAbund, pAbund);
+          for (i = 0; i < nSp; i++) {
+            F77_NAME(dgemv)(ytran, &pAbund, &pAbund, &one, TauBetaInv, &pAbund, &beta[i], &nSp, &one, tmp_pAbund, &inc FCONE);
+          } // i
+          for (q = 0; q < pAbund; q++) {
+            tmp_pAbund[q] += SigmaBetaCommInvMuBeta[q];
+          } // j
+
+          /********************************
+           Compute A.beta.comm
+           *******************************/
+          for (q = 0; q < ppAbund; q++) {
+            tmp_ppAbund[q] = SigmaBetaCommInv[q] + nSp * TauBetaInv[q];
+          }
+          F77_NAME(dpotrf)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf ABetaComm failed\n");}
+          F77_NAME(dpotri)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotri ABetaComm failed\n");}
+          F77_NAME(dsymv)(lower, &pAbund, &one, tmp_ppAbund, &pAbund, tmp_pAbund, &inc, &zero, tmp_pAbund2, &inc FCONE);
+          F77_NAME(dpotrf)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf ABetaComm failed\n");}
+          mvrnorm(betaComm, tmp_pAbund2, tmp_ppAbund, pAbund);
         }
-        F77_NAME(dpotrf)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf ABetaComm failed\n");}
-        F77_NAME(dpotri)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotri ABetaComm failed\n");}
-        F77_NAME(dsymv)(lower, &pAbund, &one, tmp_ppAbund, &pAbund, tmp_pAbund, &inc, &zero, tmp_pAbund2, &inc FCONE);
-        F77_NAME(dpotrf)(lower, &pAbund, tmp_ppAbund, &pAbund, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf ABetaComm failed\n");}
-        mvrnorm(betaComm, tmp_pAbund2, tmp_ppAbund, pAbund);
         /********************************************************************
          Update Community level Detection Coefficients
          *******************************************************************/
-        /********************************
-         * Compute b.alpha.comm
-         *******************************/
-         zeros(tmp_pDet, pDet);
-         for (i = 0; i < nSp; i++) {
-           F77_NAME(dgemv)(ytran, &pDet, &pDet, &one, TauAlphaInv, &pDet, &alpha[i], &nSp, &one, tmp_pDet, &inc FCONE);
-         } // i
-         for (q = 0; q < pDet; q++) {
-           tmp_pDet[q] += SigmaAlphaCommInvMuAlpha[q];
-         } // j
-        /********************************
-         * Compute A.alpha.comm
-         *******************************/
-        for (q = 0; q < ppDet; q++) {
-          tmp_ppDet[q] = SigmaAlphaCommInv[q] + nSp * TauAlphaInv[q];
+        if (indAlphas == 0) {
+          /********************************
+           * Compute b.alpha.comm
+           *******************************/
+           zeros(tmp_pDet, pDet);
+           for (i = 0; i < nSp; i++) {
+             F77_NAME(dgemv)(ytran, &pDet, &pDet, &one, TauAlphaInv, &pDet, &alpha[i], &nSp, &one, tmp_pDet, &inc FCONE);
+           } // i
+           for (q = 0; q < pDet; q++) {
+             tmp_pDet[q] += SigmaAlphaCommInvMuAlpha[q];
+           } // j
+          /********************************
+           * Compute A.alpha.comm
+           *******************************/
+          for (q = 0; q < ppDet; q++) {
+            tmp_ppDet[q] = SigmaAlphaCommInv[q] + nSp * TauAlphaInv[q];
+          }
+          F77_NAME(dpotrf)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf AAlphaComm failed\n");}
+          F77_NAME(dpotri)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotri AAlphaComm failed\n");}
+          F77_NAME(dsymv)(lower, &pDet, &one, tmp_ppDet, &pDet, tmp_pDet, &inc, &zero, tmp_pDet2, &inc FCONE);
+          F77_NAME(dpotrf)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf AAlphaComm failed\n");}
+          mvrnorm(alphaComm, tmp_pDet2, tmp_ppDet, pDet);
         }
-        F77_NAME(dpotrf)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf AAlphaComm failed\n");}
-        F77_NAME(dpotri)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotri AAlphaComm failed\n");}
-        F77_NAME(dsymv)(lower, &pDet, &one, tmp_ppDet, &pDet, tmp_pDet, &inc, &zero, tmp_pDet2, &inc FCONE);
-        F77_NAME(dpotrf)(lower, &pDet, tmp_ppDet, &pDet, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf AAlphaComm failed\n");}
-        mvrnorm(alphaComm, tmp_pDet2, tmp_ppDet, pDet);
 
         /********************************************************************
          Update Community Abundance Variance Parameter
         ********************************************************************/
-        for (q = 0; q < pAbund; q++) {
-          tmp_0 = 0.0;
-          for (i = 0; i < nSp; i++) {
-            tmp_0 += (beta[q * nSp + i] - betaComm[q]) * (beta[q * nSp + i] - betaComm[q]);
-          } // i
-          tmp_0 *= 0.5;
-          tauSqBeta[q] = rigamma(tauSqBetaA[q] + nSp / 2.0, tauSqBetaB[q] + tmp_0);
-        } // q
-        for (q = 0; q < pAbund; q++) {
-          TauBetaInv[q * pAbund + q] = tauSqBeta[q];
-        } // q
-        F77_NAME(dpotrf)(lower, &pAbund, TauBetaInv, &pAbund, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf TauBetaInv failed\n");}
-        F77_NAME(dpotri)(lower, &pAbund, TauBetaInv, &pAbund, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotri TauBetaInv failed\n");}
+        if (indBetas == 0) {
+          for (q = 0; q < pAbund; q++) {
+            tmp_0 = 0.0;
+            for (i = 0; i < nSp; i++) {
+              tmp_0 += (beta[q * nSp + i] - betaComm[q]) * (beta[q * nSp + i] - betaComm[q]);
+            } // i
+            tmp_0 *= 0.5;
+            tauSqBeta[q] = rigamma(tauSqBetaA[q] + nSp / 2.0, tauSqBetaB[q] + tmp_0);
+          } // q
+          for (q = 0; q < pAbund; q++) {
+            TauBetaInv[q * pAbund + q] = tauSqBeta[q];
+          } // q
+          F77_NAME(dpotrf)(lower, &pAbund, TauBetaInv, &pAbund, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf TauBetaInv failed\n");}
+          F77_NAME(dpotri)(lower, &pAbund, TauBetaInv, &pAbund, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotri TauBetaInv failed\n");}
+        }
         /********************************************************************
          Update Community Detection Variance Parameter
         ********************************************************************/
-        for (q = 0; q < pDet; q++) {
-          tmp_0 = 0.0;
-          for (i = 0; i < nSp; i++) {
-            tmp_0 += (alpha[q * nSp + i] - alphaComm[q]) * (alpha[q * nSp + i] - alphaComm[q]);
-          } // i
-          tmp_0 *= 0.5;
-          tauSqAlpha[q] = rigamma(tauSqAlphaA[q] + nSp / 2.0, tauSqAlphaB[q] + tmp_0);
-        } // q
-        for (q = 0; q < pDet; q++) {
-          TauAlphaInv[q * pDet + q] = tauSqAlpha[q];
-        } // q
-        F77_NAME(dpotrf)(lower, &pDet, TauAlphaInv, &pDet, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotrf TauAlphaInv failed\n");}
-        F77_NAME(dpotri)(lower, &pDet, TauAlphaInv, &pDet, &info FCONE);
-        if(info != 0){Rf_error("c++ error: dpotri TauAlphaInv failed\n");}
+        if (indAlphas == 0) {
+          for (q = 0; q < pDet; q++) {
+            tmp_0 = 0.0;
+            for (i = 0; i < nSp; i++) {
+              tmp_0 += (alpha[q * nSp + i] - alphaComm[q]) * (alpha[q * nSp + i] - alphaComm[q]);
+            } // i
+            tmp_0 *= 0.5;
+            tauSqAlpha[q] = rigamma(tauSqAlphaA[q] + nSp / 2.0, tauSqAlphaB[q] + tmp_0);
+          } // q
+          for (q = 0; q < pDet; q++) {
+            TauAlphaInv[q * pDet + q] = tauSqAlpha[q];
+          } // q
+          F77_NAME(dpotrf)(lower, &pDet, TauAlphaInv, &pDet, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotrf TauAlphaInv failed\n");}
+          F77_NAME(dpotri)(lower, &pDet, TauAlphaInv, &pDet, &info FCONE);
+          if(info != 0){Rf_error("c++ error: dpotri TauAlphaInv failed\n");}
+        }
 
         /********************************************************************
          *Update Abundance random effects variance