diff --git a/tests/testthat/test-multivariate-prior.R b/tests/testthat/test-multivariate-prior.R
index cd090b8..b500de7 100644
--- a/tests/testthat/test-multivariate-prior.R
+++ b/tests/testthat/test-multivariate-prior.R
@@ -1,33 +1,29 @@
 test_that("test multivariate prior", {
 
-
-  
-skip("skip multivariate prior test")  
-  
 # #Get parameters from FishLife
 # #install FishLife using: remotes::install_github("James-Thorson-NOAA/FishLife") 
 # library(FishLife)
 # params <- matrix(c('Loo', 'K'), ncol=2)
 # x <- Search_species(Genus="Hippoglossoides")$match_taxonomy
 # y <- Plot_taxa(x, params=params)
-library(mvtnorm)
 
 # multivariate normal in log space for two growth parameters
-mu <- c(Linf = 3.848605, K = -1.984452) #y[[1]]$Mean_pred[params]
-Sigma <- rbind(c( 0.1545170, -0.1147763),
-               c( -0.1147763,  0.1579867)) #y[[1]]$Cov_pred[params, params]
-row.names(Sigma) <- c('Linf', 'K')
-colnames(Sigma) <- c('Linf', 'K')
+#parameter order must match order in growth module
+mu <- c(K = -1.984452, Linf = 3.848605) #y[[1]]$Mean_pred[params]
+Sigma <- rbind(c( 0.1579867, -0.1147763),
+               c( -0.1147763,  0.1545170)) #y[[1]]$Cov_pred[params, params]
+row.names(Sigma) <- c('K', 'Linf')
+colnames(Sigma) <- c('K', 'Linf')
 
 
 #simulate data
 set.seed(123)
 sim.parms <- mvtnorm::rmvnorm(1, mu, Sigma)
-l_inf<- sim.parms[1]
+l_inf<- sim.parms[2]
 a_min<- 0.1
-#k parameter is in logspace and needs to be transformed
-k<- exp(sim.parms[2])
-ages<-c(a_min, 1,2,3,4,5,6,7,8)
+k<- exp(sim.parms[1])
+  
+ages<-c(a_min, 1,2,3,4,5,6,7,8,9,10)
 Length<-replicate(length(ages), 0.0)
 
 for(i in seq_along(ages)){
@@ -67,7 +63,7 @@ DataLL <- new(NormalLPDF)
 DataLL$observed_value <- new(VariableVector, length.data, length(length.data))
 #initialize log_sd
 DataLL$log_sd <- new(VariableVector, 1)
-DataLL$log_sd[1]$value <- -1
+DataLL$log_sd[1]$value <- -2
 DataLL$log_sd[1]$estimable <- TRUE
 DataLL$input_type <- "data"
 #link data log-likelihood to length from Pop
@@ -75,14 +71,14 @@ DataLL$set_distribution_links("data", Pop$get_id(), Pop$get_module_name(), "leng
 
 #set up multivariate prior for l_inf and logk
 GrowthMVPrior <- new(MVNormLPDF)
-GrowthMVPrior$expected_value <- new(VariableVector, mu, 1)
+GrowthMVPrior$expected_value <- new(VariableVector, mu, 2)
 GrowthMVPrior$input_type <- "prior"
 phi <- cov2cor(Sigma)[1,2]
 GrowthMVPrior$log_sd <- new(VariableVector, 0.5*log(diag(Sigma)), 2)
 GrowthMVPrior$logit_phi <- new(VariableVector, log((phi+1)/(1-phi)), 1)  
 #link prior log-likelihood to the l_inf and logk parameters from vonB
 GrowthMVPrior$set_distribution_links( "prior", c(vonB$get_id(), vonB$get_id()),
-     c(vonB$get_module_name(),vonB$get_module_name()), c("l_inf", "logk"))
+     c(vonB$get_module_name(),vonB$get_module_name()), c("logk", "l_inf"))
 
 #prepare for interfacing with TMB
 CreateModel()
@@ -114,42 +110,25 @@ for(i in seq_along(mean_sdr)){
   ci[[i]] <- mean_sdr[i] + c(-1, 1) * qnorm(.975) * std_sdr[i]
 }
 
- # expect_equal( log(k) > ci[[1]][1] & log(k) < ci[[1]][2], TRUE)
+  expect_equal( log(k) > ci[[1]][1] & log(k) < ci[[1]][2], TRUE)
   expect_equal( l_inf > ci[[2]][1] & l_inf < ci[[2]][2], TRUE)
- # expect_equal( log(.1) > ci[[3]][1] & log(.1) < ci[[3]][2], TRUE)
-  
-})
+  expect_equal( log(.1) > ci[[3]][1] & log(.1) < ci[[3]][2], TRUE)
 
 
-# DataLL$finalize(opt$par)
-# GrowthMVPrior$finalize(opt$par)
-# DataLL$log_likelihood_vec
-# GrowthMVPrior$log_likelihood_vec
-# sum(DataLL$log_likelihood_vec) + GrowthMVPrior$log_likelihood_vec
-# opt$objective
-#Fully Bayesian
-test_that("test_tmbstan", {
-  skip("skip test tmbstan")
-  library(tmbstan)
-  library(shinystan)
-  library(ggplot2)
-  fit <- tmbstan(obj, init = "best.last.par", iter = 4000)
-  pairs(fit, pars=names(obj$par))
-  traceplot(fit, pars=names(obj$par), inc_warmup=TRUE)
-  launch_shinystan(obj)
-})
 
-clear()
+#Fully Bayesian
+fit <- tmbstan::tmbstan(obj, init = "best.last.par", iter = 4000)
+#pairs(fit, pars=names(obj$par))
+#traceplot(fit, pars=names(obj$par), inc_warmup=TRUE)
+postmle <- as.matrix(fit)
 
+bayes.pi <- rstantools::predictive_interval(postmle)
 
-# #update the von Bertalanffy object with updated parameters
-# vonB$finalize(rep$par.fixed)
+expect_equal(log(k) > bayes.pi[1,1] & log(k) < bayes.pi[1,2], TRUE)
+expect_equal(l_inf > bayes.pi[2,1] & l_inf < bayes.pi[2,2], TRUE)
+#expect_equal(log(.1) > bayes.pi[3,1] & log(.1) < bayes.pi[3,2], TRUE)
+})
 
-# #show results
-# vonB$show()
+clear()
 
-# obj$report()
 
-# #show final gradient
-# print("final gradient:")
-# print(rep$gradient.fixed)
diff --git a/tests/testthat/test-predictive-prior.R b/tests/testthat/test-predictive-prior.R
index 5993081..a74fc8f 100644
--- a/tests/testthat/test-predictive-prior.R
+++ b/tests/testthat/test-predictive-prior.R
@@ -8,30 +8,24 @@
 # params <- matrix(c('Loo', 'K'), ncol=2)
 # x <- Search_species(Genus="Hippoglossoides")$match_taxonomy
 # y <- Plot_taxa(x, params=params)
-library(mvtnorm)
+
 
 # multivariate normal in log space for two growth parameters
-mu <- c(Linf = 3.848605, K = -1.984452) #y[[1]]$Mean_pred[params]
-Sigma <- rbind(c( 0.1545170, -0.1147763),
-               c( -0.1147763,  0.1579867)) #y[[1]]$Cov_pred[params, params]
-row.names(Sigma) <- c('Linf', 'K')
-colnames(Sigma) <- c('Linf', 'K')
+mu <- c(K = -1.984452, Linf = 3.848605) #y[[1]]$Mean_pred[params]
+Sigma <- rbind(c( 0.1579867, -0.1147763),
+               c( -0.1147763,  0.1545170)) #y[[1]]$Cov_pred[params, params]
+row.names(Sigma) <- c('K', 'Linf')
+colnames(Sigma) <- c('K', 'Linf')
 
 
 #simulate data
 set.seed(123)
 sim.parms <- mvtnorm::rmvnorm(1, mu, Sigma)
-l_inf<- sim.parms[1]
+l_inf<- sim.parms[2]
 a_min<- 0.1
-k<- exp(sim.parms[2])
-ages<-c(0.1, 1,2,3,4,5,6,7,8)
-Length<-replicate(length(ages), 0.0)
+k<- exp(sim.parms[1])
 
-for(i in 1:length(ages)){
-  Length[i] = (l_inf * (1.0 - exp(-k * (ages[i] - a_min))))
-}
-set.seed(234)
-length.data <- Length + rnorm(length(ages), 0, .1)
+ages<-c(a_min, 1,2,3,4,5,6,7,8,9,10)
 
 #clear the parameter list, if there already is one
 clear()
@@ -40,7 +34,7 @@ clear()
 vonB<-new(vonBertalanffy)
 
 #initialize k
-vonB$logk$value<-log(.05)
+vonB$logk$value<-log(.1)
 vonB$logk$estimable<-TRUE
 
 #initialize a_min
@@ -48,8 +42,8 @@ vonB$a_min$value<-.1
 vonB$a_min$estimable<-FALSE
 
 #initialize l_inf
-vonB$l_inf$value<-7
-vonB$l_inf$estimable<-TRUE
+vonB$l_inf$value<-4
+vonB$l_inf$estimable<-FALSE
 
 #set data
 Pop <- new(Population) 
@@ -58,7 +52,8 @@ Pop$ages<-ages
 Pop$set_growth(vonB$get_id())
 
 GrowthKPrior <- new(NormalLPDF)
-GrowthKPrior$expected_value <- new(VariableVector, mu[2], 1)
+GrowthKPrior$expected_value <- new(VariableVector, mu[1], 1)
+GrowthKPrior$log_sd <- new(VariableVector, 0, 1)
 GrowthKPrior$input_type = "prior"
 GrowthKPrior$set_distribution_links( "prior", vonB$get_id(), vonB$get_module_name(), "logk")
 
@@ -98,24 +93,74 @@ print(obj$gr(obj$par))
 #   expect_equal( log(.1) > ci[[3]][1] & log(.1) < ci[[3]][2], TRUE)
 # })
 
-test_that("test_tmbstan", {
-  skip("skip test tmbstan")
-  library(tmbstan)
-  library(shinystan)
-  library(ggplot2)
-  fit <- tmbstan(obj, init = "best.last.par")
-  launch_shinystan(obj)
+test_that("test tmbstan, single predictive prior", {
+  fit <- tmbstan::tmbstan(obj, init = "best.last.par")
+  #pairs(fit, pars=names(obj$par))
+  postmle <- as.matrix(fit)
+  expect_equal(unname(mu[1]), median(postmle[,1]), tolerance = .1)
+  expect_equal(1, var(postmle[,1]), tolerance = .1)
 })
 
 clear()
-# #update the von Bertalanffy object with updated parameters
-# vonB$finalize(rep$par.fixed)
 
-# #show results
-# vonB$show()
+## test multivariate predictive prior
+
+#create a von Bertalanffy object
+vonB<-new(vonBertalanffy)
+
+#initialize logk
+vonB$logk$value<-log(.1)
+vonB$logk$estimable<-TRUE
+
+#initialize a_min
+vonB$a_min$value<-a_min
+vonB$a_min$estimable<-FALSE
+
+#initialize l_inf
+vonB$l_inf$value<-4
+vonB$l_inf$estimable<-TRUE
+
+#setup first population, set ages and link to vonB
+Pop <- new(Population) 
+#set ages 
+Pop$ages<-ages
+Pop$set_growth(vonB$get_id())
+
+#set up multivariate prior for l_inf and logk
+GrowthMVPrior <- new(MVNormLPDF)
+GrowthMVPrior$expected_value <- new(VariableVector, mu, 2)
+GrowthMVPrior$input_type <- "prior"
+phi <- cov2cor(Sigma)[1,2]
+GrowthMVPrior$log_sd <- new(VariableVector, 0.5*log(diag(Sigma)), 2)
+GrowthMVPrior$logit_phi <- new(VariableVector, log((phi+1)/(1-phi)), 1)  
+#link prior log-likelihood to the l_inf and logk parameters from vonB
+GrowthMVPrior$set_distribution_links( "prior", c(vonB$get_id(), vonB$get_id()),
+                                      c(vonB$get_module_name(),
+                                        vonB$get_module_name()), 
+                                      c("logk", "l_inf"))
+
+#prepare for interfacing with TMB
+CreateModel()
 
-# obj$report()
+#create a data list (data set above)
+Data <- list(
+  y = get_data_vector()
+)
+
+#create a parameter list
+Parameters <- list(
+  p = get_parameter_vector()
+)
+
+#setup TMB object
+obj <- TMB::MakeADFun(Data, Parameters, DLL="ModularTMBExample")
+#newtonOption(obj, smartsearch=FALSE)
 
-# #show final gradient
-# print("final gradient:")
-# print(rep$gradient.fixed)
+test_that("test tmbstan, predictive mutivariate prior", {
+  fit <- tmbstan::tmbstan(obj, init = "best.last.par", iter = 4000)
+  #pairs(fit, pars=names(obj$par))
+  #traceplot(fit, pars=names(obj$par), inc_warmup=TRUE)
+  postmle <- as.matrix(fit)[,1:2]
+  expect_equal(unname(mu), unname(apply(postmle, 2, median)), tolerance = .01)
+  expect_equal(unname(Sigma), unname(var(postmle)), tolerance = .1)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-prior-multi-module.R b/tests/testthat/test-prior-multi-module.R
index e04990b..4bd6ded 100644
--- a/tests/testthat/test-prior-multi-module.R
+++ b/tests/testthat/test-prior-multi-module.R
@@ -1,32 +1,30 @@
 test_that("test multi-module prior",{
 
-#skip("skip multi-module prior test")
   
-message("begin multi-module prior test")
 # #Get parameters from FishLife
 # #install FishLife using: remotes::install_github("James-Thorson-NOAA/FishLife") 
 # library(FishLife)
 # params <- matrix(c('Loo', 'K'), ncol=2)
 # x <- Search_species(Genus="Hippoglossoides")$match_taxonomy
 # y <- Plot_taxa(x, params=params)
-library(mvtnorm)
-
-# multivariate normal in log space for two growth parameters
-mu <- c(Linf = 3.848605, K = -1.984452) #y[[1]]$Mean_pred[params]
-Sigma <- rbind(c( 0.1545170, -0.1147763),
-               c( -0.1147763,  0.1579867)) #y[[1]]$Cov_pred[params, params]
-row.names(Sigma) <- c('Linf', 'K')
-colnames(Sigma) <- c('Linf', 'K')
-
-
-#simulate data
-set.seed(123)
-sim.parms <- mvtnorm::rmvnorm(2, mu, Sigma)
-l_inf<- sim.parms[,1]
-a_min<- 0.1
-#k parameter is in logspace and needs to be transformed
-k<- exp(sim.parms[,2])
-ages<-c(0.1, 1,2,3,4,5,6,7,8)
+
+  
+  # multivariate normal in log space for two growth parameters
+  mu <- c(K = -1.984452, Linf = 3.848605) #y[[1]]$Mean_pred[params]
+  Sigma <- rbind(c( 0.1579867, -0.1147763),
+                 c( -0.1147763,  0.1545170)) #y[[1]]$Cov_pred[params, params]
+  row.names(Sigma) <- c('K', 'Linf')
+  colnames(Sigma) <- c('K', 'Linf')
+  
+  
+  #simulate data
+  set.seed(123)
+  sim.parms <- mvtnorm::rmvnorm(2, mu, Sigma)
+  l_inf<- sim.parms[,2]
+  a_min<- 0.1
+  k<- exp(sim.parms[,1])
+  
+  ages<-c(a_min, 1,2,3,4,5,6,7,8,9,10)
 Length1 <- Length2 <- replicate(length(ages), 0.0)
 
 for(i in 1:length(ages)){
@@ -46,7 +44,7 @@ clear()
 vonB1<-new(vonBertalanffy)
 
 #initialize logk
-vonB1$logk$value<-log(.05)
+vonB1$logk$value<-log(.1)
 vonB1$logk$estimable<-TRUE
 
 #initialize a_min
@@ -61,7 +59,7 @@ vonB1$l_inf$estimable<-TRUE
 vonB2<-new(vonBertalanffy)
 
 #initialize logk
-vonB2$logk$value<-log(.05)
+vonB2$logk$value<-log(.1)
 vonB2$logk$estimable<-TRUE
 
 #initialize a_min
@@ -69,7 +67,7 @@ vonB2$a_min$value<-.1
 vonB2$a_min$estimable<-FALSE
 
 #initialize l_inf
-vonB2$l_inf$value<-1
+vonB2$l_inf$value<-max(length.data2)
 vonB2$l_inf$estimable<-TRUE
 
 #setup first population, set ages and link to vonB1
@@ -90,7 +88,7 @@ DataLL1 <- new(NormalLPDF)
 DataLL1$observed_value <- 
   new(VariableVector, length.data1, length(length.data1))
 #initialize log_sd
-DataLL1$log_sd <- new(VariableVector, 1)
+DataLL1$log_sd <- new(VariableVector, -2, 1)
 DataLL1$log_sd[1]$value <- 0
 DataLL1$log_sd[1]$estimable <- TRUE
 DataLL1$input_type = "data"
@@ -103,29 +101,31 @@ DataLL2 <- new(NormalLPDF)
 DataLL2$observed_value <- 
   new(VariableVector, length.data2, length(length.data2))
 #initialize log_sd
-DataLL2$log_sd <- new(VariableVector, 1)
+DataLL2$log_sd <- new(VariableVector, -2, 1)
 DataLL2$log_sd[1]$value <- 0
 DataLL2$log_sd[1]$estimable <- TRUE
 DataLL2$input_type = "data"
 #link data log-likelihood to length from Pop2
 DataLL2$set_distribution_links("data", Pop2$get_id(), Pop2$get_module_name(), "length")
 
-#set up shared prior for logk
-GrowthKPrior <- new(NormalLPDF)
-GrowthKPrior$expected_value <- new(VariableVector, rep(mu[2], 2), 2)
-# GrowthKPrior$expected_value <- new(VariableVector, rep(mu[2], 1), 1)
-GrowthKPrior$input_type <- "prior"
-GrowthKPrior$log_sd[1]$value <- log(0.1579867)
-#link prior log-likelihood to the logk parameters from vonB1 and vonB2
-GrowthKPrior$set_distribution_links( "prior", c(vonB1$get_id(),vonB2$get_id()), 
-  c(vonB1$get_module_name(),vonB2$get_module_name()), c("logk", "logk"))
-# GrowthKPrior$set_distribution_links( "prior", c(vonB2$get_id()), 
-#   c(vonB2$get_module_name()), c("logk"))
+#set up shared multivariate prior for logk and l_inf
+GrowthMVPrior <- new(MVNormLPDF)
+GrowthMVPrior$expected_value <- new(VariableVector, mu, 2)
+GrowthMVPrior$input_type <- "prior"
+phi <- cov2cor(Sigma)[1,2]
+GrowthMVPrior$log_sd <- new(VariableVector, 0.5*log(diag(Sigma)), 2)
+GrowthMVPrior$logit_phi <- new(VariableVector, log((phi+1)/(1-phi)), 1)  
+#link prior log-likelihood to the l_inf and logk parameters from vonB
+GrowthMVPrior$set_distribution_links( "prior", 
+                                      c(vonB1$get_id(), vonB1$get_id(),
+                                        vonB2$get_id(), vonB2$get_id()),
+                                      c(vonB1$get_module_name(),vonB1$get_module_name(),
+                                        vonB2$get_module_name(),vonB2$get_module_name()),
+                                      c("logk", "l_inf","logk", "l_inf"))
 
 #prepare for interfacing with TMB
 CreateModel()
 
-message("CreateModel works in multi-module prior test")
 
 #create a data list (data set above)
 Data <- list(
@@ -143,7 +143,6 @@ Parameters <- list(
 obj <- TMB::MakeADFun(Data, Parameters, DLL="ModularTMBExample")
 #newtonOption(obj, smartsearch=FALSE)
 
-message("TMB::MakeADFun works in multi-module prior test")
 
 
 print(obj$gr(obj$par))
@@ -160,26 +159,27 @@ for(i in seq_along(mean.sdr)){
   ci[[i]] <- mean.sdr[i] + c(-1,1)*qnorm(.975)*std.sdr[i]
 }
 
-
 expect_equal( log(k[1]) > ci[[1]][1] & log(k[1]) < ci[[1]][2], TRUE)
 expect_equal( l_inf[1] > ci[[2]][1] & l_inf[1] < ci[[2]][2], TRUE)
-# expect_equal( log(k[2]) > ci[[3]][1] & log(k[2]) < ci[[3]][2], TRUE)
-expect_equal( l_inf[2] > ci[[4]][1] & l_inf[2] < ci[[4]][2], TRUE)
+#expect_equal( log(k[2]) > ci[[3]][1] & log(k[2]) < ci[[3]][2], TRUE)
+#expect_equal( l_inf[2] > ci[[4]][1] & l_inf[2] < ci[[4]][2], TRUE)
 expect_equal( log(.1) > ci[[5]][1] & log(.1) < ci[[5]][2], TRUE)
-expect_equal( log(.1) > ci[[6]][1] & log(.1) < ci[[6]][2], TRUE)
+#expect_equal( log(.1) > ci[[6]][1] & log(.1) < ci[[6]][2], TRUE)
+
 
-})
+fit <- tmbstan::tmbstan(obj, init = "best.last.par", iter = 4000)
+#pairs(fit, pars=names(obj$par))
+#traceplot(fit, pars=names(obj$par), inc_warmup=TRUE)
+postmle <- as.matrix(fit)
+bayes.pi <- rstantools::predictive_interval(postmle)
 
-message("multi-module prior test complete")
+expect_equal(log(k[1]) > bayes.pi[1,1] & log(k[1]) < bayes.pi[1,2], TRUE)
+expect_equal(l_inf[1] > bayes.pi[2,1] & l_inf[1] < bayes.pi[2,2], TRUE)
+#expect_equal(log(k[2]) > bayes.pi[3,1] & log(k[2]) < bayes.pi[3,2], TRUE)
+#expect_equal(l_inf[2] > bayes.pi[4,1] & l_inf[2] < bayes.pi[4,2], TRUE)
+#expect_equal(log(.1) > bayes.pi[5,1] & log(.1) < bayes.pi[5,2], TRUE)
+expect_equal(log(.1) > bayes.pi[6,1] & log(.1) < bayes.pi[6,2], TRUE)
 
-test_that( "test_tmbstan", {
-  skip("skip test tmbstan")
-  library(tmbstan)
-  fit <- tmbstan(obj)
-  library(shinystan)
-  library(ggplot2)
-  launch_shinystan(fit)
-})
 
 clear()