Howdy! I’d be grateful if someone would take a look at the below serial dilution model. I keep getting initialization errors…I guess I have looked at it for too long and can’t find anything… or maybe the model is just not very good.
This is a small part of a much larger and complicated model of a laboratory experiment involving several variations and different serial dilutions. I’m trying to get the simplest part to work first.
In this part, the concept is simple: you have a stock of spores (unknown volume and number of spores), where each 10ul sample from stock contains a latent number of spores lambda_mu. We can think of lambda_mu as a normal distribution with some standard deviation, that represent a population of 10ul samples from stock. Some samples from stock are taken and the goal is to find the number of spores in each sample, lambda_con. Each sample is serially diluted and then samples are put on growth plates and the number of spores are counted. The seral dilution has to be done, because the number of spores in the original sample is quite large with maybe mean 5e9. Some of the early serial dilutions have censored values because the spores on the plate are too numerous to count.
At every stage of dilution you get pipetting error. I am treating the error as coming from a single normal distribution with some mean and sd, and I scale the error appropriately at each step of the dilution.
Non-centered parameterization is used.
The data and model are below, which can be run and fully reproduce the problem (apologies for using an old rstan version with different array syntax, but that is all I have on this computer, and can’t currently update):
library(rstan)
#define data used in Stan model in a list
stan_data <- list(
N1_sp = 48,
s = c(330,330,49,48,23,23,0,0,330,330,51,74,11,10,0,8,330,330,57,49,8,2,1,0,330,330,34,39,6,2,0,0,330,330,49,49,14,2,1,2,330,330,47,52,2,5,0,1),
Ex = 6,
Ex_1 = rep(1:6, each=8),
censored_sp = rep(c(1,1,0,0,0,0,0,0), times=6),
U = 330,
dilution_plate_sp = rep(c(6,6,7,7,8,8,9,9), times=6)
)
#define Stan model
stan_model <- "
data {
int N1_sp; //number of observations (including censored values)
int s[N1_sp]; //observed colony counts on plates
int Ex; //number of straight plate experiments
int Ex_1[N1_sp]; //index for the straight plate experiments
int censored_sp[N1_sp]; //indicator for censored observations
int U; //upper bound for censoring
int dilution_plate_sp[N1_sp]; //the dilution plate indicator for the if else statements for straight plate serial dilutions
}
parameters {
vector<lower=0>[Ex] lambda_con_std;
real<lower=0> lambda_mu;
real<lower=0> tau;
vector[Ex] error_6_std;
vector[Ex] error_7_std;
vector[Ex] error_8_std;
real<lower=0> tau_er;
}
transformed parameters {
vector<lower=0>[Ex] lambda_con = lambda_mu + tau * lambda_con_std; //the mean number of spores in the 10ul sample from stock (for straight plate data; same distribution parameters lambda_mu and tau as for the coupons)
vector[Ex] error_6 = tau_er * error_6_std;
vector[Ex] error_7 = tau_er * error_7_std;
vector[Ex] error_8 = tau_er * error_8_std;
print(lambda_con);
}
model {
// likelihood for straight plating data, including censoring
for (n in 1:N1_sp) {
if (dilution_plate_sp[n]==6) {
if (censored_sp[n]==0) {
target += poisson_lpmf( s[n] | 0.000001 * lambda_con[Ex_1[n]] );
}
else if (censored_sp[n]==1) {
target += poisson_lccdf(U | 0.000001 * lambda_con[Ex_1[n]] );
}
}
if (dilution_plate_sp[n]==7) {
if (censored_sp[n]==0) {
target += poisson_lpmf( s[n] | 0.0000001 * lambda_con[Ex_1[n]] + 0.000001 * error_6[Ex_1[n]] );
}
else if (censored_sp[n]==1) {
target += poisson_lccdf(U | 0.0000001 * lambda_con[Ex_1[n]] + 0.000001 * error_6[Ex_1[n]] );
}
}
if (dilution_plate_sp[n]==8) {
if (censored_sp[n]==0) {
target += poisson_lpmf( s[n] | 0.00000001 * lambda_con[Ex_1[n]] + 0.0000001 * error_6[Ex_1[n]] + 0.0000001 * error_7[Ex_1[n]] );
}
else if (censored_sp[n]==1) {
target += poisson_lccdf(U | 0.00000001 * lambda_con[Ex_1[n]] + 0.0000001 * error_6[Ex_1[n]] + 0.0000001 * error_7[Ex_1[n]] );
}
}
if (dilution_plate_sp[n]==9) {
if (censored_sp[n]==0) {
target += poisson_lpmf( s[n] | 0.000000001 * lambda_con[Ex_1[n]] + 0.00000001 * error_6[Ex_1[n]] + 0.00000001 * error_7[Ex_1[n]] + 0.00000001 * error_8[Ex_1[n]] );
}
else if (censored_sp[n]==1) {
target += poisson_lccdf(U | 0.000000001 * lambda_con[Ex_1[n]] + 0.00000001 * error_6[Ex_1[n]] + 0.00000001 * error_7[Ex_1[n]] + 0.00000001 * error_8[Ex_1[n]] );
}
}
}
//priors
target += normal_lpdf(lambda_mu | 5e9, 1e9);
target += normal_lpdf(tau | 1.5e9, 5e8) - 1 * normal_lccdf(0 | 1.5e9, 5e8);
target += normal_lpdf(tau_er | 0, 7.5e7) - 1 * normal_lccdf(0 | 0, 7.5e7);
target += std_normal_lpdf(lambda_con_std);
target += std_normal_lpdf(error_6_std);
target += std_normal_lpdf(error_7_std);
target += std_normal_lpdf(error_8_std);
}
generated quantities {
vector[N1_sp] s_preds;
// predictions for straight plating
for (n in 1:N1_sp) {
if (dilution_plate_sp[n]==6) {
s_preds[n] = poisson_rng( 0.000001 * lambda_con[Ex_1[n]] );
}
if (dilution_plate_sp[n]==7) {
s_preds[n] = poisson_rng( 0.0000001 * lambda_con[Ex_1[n]] + 0.000001 * error_6[Ex_1[n]] );
}
if (dilution_plate_sp[n]==8) {
s_preds[n] = poisson_rng( 0.00000001 * lambda_con[Ex_1[n]] + 0.0000001 * error_6[Ex_1[n]] + 0.0000001 * error_7[Ex_1[n]] );
}
if (dilution_plate_sp[n]==9) {
s_preds[n] = poisson_rng( 0.000000001 * lambda_con[Ex_1[n]] + 0.00000001 * error_6[Ex_1[n]] + 0.00000001 * error_7[Ex_1[n]] + 0.00000001 * error_8[Ex_1[n]] );
}
}
}
"
#fit Stan model
un20_sp <- stan(model_code = stan_model, data = stan_data,
chains = 1, iter = 10, init="0",
thin = 1, cores = 4)
Perhaps it’s too flexible with all the error terms? Priors reflect the manufacturer specs for accuracy error from the pipette.
Thanks!
