Piecewise Linear Fit of FIO Data - pooled parameters
Leighton Pritchard and Nicholas Waters
2018-06-14
Introduction
To analyse the Fecal Indicator Organism (FIO) levels under the effect of various treatments (10L/50mL vessels; pasteurised/unpasteurised; pre-pasteurisation; post-pasteurisation) we use piecewise linear regression to estimate the rate of FIO loss at two stages of fermentation: “early” and “late”. We pool the estimated parameters, assuming that each is drawn from a single distribution that applies to the bacterial family.
Piecewise linear regression
It would be natural to assume a smooth (rather than piecewise) exponential decay rate for FIO as fermentation progresses. Arguably, the choice of piecewise linear regression is mechanistically inappropriate, as it is unlikely that population change occurs in two distinct, linear phases. However, several features of the data lead us to use piecewise linear regression as a reasonable approximation to the measured change in population size for the purpose of identifying treatment effects in this study, despite it being an unlikely representation of the true biological mechanism.
- There are a small number of widely-spaced measured timepoints. These do not describe the early stage reduction of bacteria in sufficient resolution to fit a single detailed mechanism of bacterial population change. This is particularly so in the early stages of population reduction.
- The timepoints from 7d onwards describe a regime with little overall change in bacterial population. Datapoints from 7d onwards are the bulk of the dataset. A linear fit to this regime permits estimates of:
- slope: which we might expect to be zero if the bacterial population is at equilibrium, and non-zero if there continues to be growth or loss
- intercept: which we may interpret as a proxy for the equilibrium population level, if the credibility interval for the slope contains zero
- The timepoints up to 7d cover the majority of bacterial population reduction, but not at a detailed resolution capable of discriminating between alternative models of population reduction. In several cases we have only data at 0d and 7d. Here, a linear fit to this regime provides estimates of:
- slope: which we may interpret as a proxy for average rate of reduction over seven days, regardless of the total or initial rate of loss for that treatment
- intercept: which we may interpret as the initial bacterial load
By using piecewise linear regression with two regimes we approximate an early ‘loss’ regime, and late ‘equilibrium’ regime for the FIO population. While mechanistically far from ideal, we remain able to assess differences between: rate of FIO reduction (taking into account initial bacterial load); final bacterial load; and to what extent the bacterial population is at an equilibrium, for each treatment.
Location of transition point
It would be possible to assume that the transition point between ‘early’ (loss) and ‘late’ (equilibrium) could occur at any timepoint. However, we assert a transition point at seven days, for the following reasons:
- For all treatments, there appears to be no effective change in FIO load after the 7d timepoint. The transition between ‘early’ and ‘late’ stages is therefore between 0d and 7d.
- For some treatments, we do not have timepoints between 0d and 7d. In order to estimate a rate of FIO loss, we require a minimum of two timepoints, and would be unable to infer ‘early’ stage FIO loss if we took the transition to be any earlier than 7d. For these treatments, we must therefore take 0d and 7d to be the start and end of the ‘early’ period.
The model
We fit the following piecewise linear model to the cleaned dataset for each bacterial group:
\[\hat{y} = \bigg \{ \begin{array}{cccl} m1_i x & + & c1_i & , \textrm{if } x <= 7 \\ m2_i x & + & c2_i & , \textrm{if } x > 7\end{array} \bigg \} + \epsilon \]
where:
- \(\hat{y}\) is the FIO population level
- \(x\) is the measurement timepoint
- \(i\) is an index representing the treatment
- \(i = 1\): 10L vessel, no pasteurisation
- \(i = 2\): 50mL vessel, no pasteurisation (P0)
- \(i = 3\): 50mL vessel, pre-pasteurisation (P1)
- \(i = 4\): 50mL vessel, pre-pasteurisation (P2)
- \(i = 5\): 50mL vessel, post-pasteurisation (P3)
- \(i = 6\): 50mL vessel, post-pasteurisation (P4)
- \(m1\) is the early stage slope
- \(c1\) is the early stage intercept
- \(m2\) is the late stage slope
- \(c2\) is the late stage intercept
- \(\epsilon\) is irreducible (measurement) error
and \(m1_i\), \(c1_i\), \(m2_i\), and \(c2_i\) are values drawn from the following distributions:
- \(m1_i ~ \textrm{Cauchy}(\mu_{m1}, \sigma_{m1})\)
- \(c1_i ~ \textrm{Cauchy}(\mu_{c1}, \sigma_{c1})\)
- \(m2_i ~ \textrm{Cauchy}(\mu_{m2}, \sigma_{m2})\)
- \(c2_i ~ \textrm{Cauchy}(\mu_{c2}, \sigma_{c2})\)
Data Processing
We use the cleaned dataset in ../data/processed/aggregated_clean_summarized.csv, as described elsewhere in the Supplementary Information and manuscript, and make it suitable for analysis.
- First the data is loaded, and columns converted to the appropriate datatype.
# Load the experimental data
data = read.table(datafile, header=TRUE, sep=",")
# Replicate labels should be factors, not integers
data$biorep = as.factor(data$biorep)
data$techrep = as.factor(data$techrep)- The feedstock (
FS) data are removed, along with the non-indicator measurements:
# Remove the feedstock data
data = data %>%
filter(treatment != "FS") %>%
filter(measure %in% c("Coliforms", "E. coli", "Enterococci")) %>%
droplevels()- Columns are renamed for ease of interpretation:
day: measurement timepointtreatment: string describing the treatmentreplicate: the replicate identifier for the measurement at that timepoint for that treatmentmeasurement: identifier for the measurement (meaning is dependent on the vessel being measured - see manuscript for details)measure: string describing what was measuredvalue: the measured value
# Rename techrep column
names(data) = c("day", "treatment", "replicate", "measurement", "measure", "value")- Create a new column:
vessel: a unique combination of timepoint, replicate and treatment to clarify the distinction between repeated measures and destructive sampling; this contains a unique value for each physically-measured vessel
# Map the biorep column to a new column called vessel, where each vessel is a
# unique combination of timepoint, biorep, and treatment. This helps reflects which values are repeated measures over time, and which are the result of destructive sampling.
data = data %>%
mutate(vessel = paste(treatment, replicate, day, sep="_"))
# Set the chamber for the 10L treatments to be the biorep value
data[data$treatment == "10L",] = data %>%
filter(treatment == "10L") %>%
mutate(vessel=replicate)
# Make the chamber column a factor
data$vessel = as.factor(data$vessel)- Add columns containing indices for each treatment. These are useful as indices for tracking treatments, and calculating contrasts.
# Add "small" index (50ml vessel, applies to all P0-P4 treatments, not 10L)
# I'm going to rename this "small", just to avoid the confusion I experienced at first
data$small = !data$treatment %in% c("10L")
# Add index for each pasteurisation treatment
data$P1 = data$treatment == "P1"
data$P2 = data$treatment == "P2"
data$P3 = data$treatment == "P3"
data$P4 = data$treatment == "P4"- Change FIO names, for easier manipulation in
R.
# Replace bacterial measurement names
data$measure = data$measure %>%
str_replace("E. coli", "ecoli") %>%
str_replace("Enterococci", "enterococci") %>%
str_replace("Coliforms", "coliforms")- Create a new column (
logvalue) with log-transformed bacterial count data. We expect the limit of detection to be approximately a count of 100, so we take any “undetected” value and raise it to half of the expected detection limit, before taking the log.
# Convert bacterial counts to logs, and set zeroes to counts of '50' (below detection limit)
#data$logvalue = log10(pmax(data$value, 50))
data$logvalue = log10(ifelse(data$value < 100, 50, data$value))A summary of the dataset is included below:
# Columnwise summary of data
summary(data)## day treatment replicate measurement measure
## Min. : 0.00 10L:76 1:112 1:294 Length:329
## 1st Qu.: 0.00 P0 :51 2:108 2: 35 Class :character
## Median : 9.00 P1 :49 3:109 Mode :character
## Mean :11.85 P2 :51
## 3rd Qu.:21.00 P3 :51
## Max. :28.00 P4 :51
##
## value vessel small P1
## Min. : 0 1 : 26 Mode :logical Mode :logical
## 1st Qu.: 200 3 : 26 FALSE:76 FALSE:280
## Median : 700 2 : 24 TRUE :253 TRUE :49
## Mean : 337097 P0_1_0 : 6
## 3rd Qu.: 86500 P0_2_0 : 5
## Max. :3850000 P1_1_0 : 5
## (Other):237
## P2 P3 P4 logvalue
## Mode :logical Mode :logical Mode :logical Min. :1.699
## FALSE:278 FALSE:278 FALSE:278 1st Qu.:2.301
## TRUE :51 TRUE :51 TRUE :51 Median :2.845
## Mean :3.408
## 3rd Qu.:4.937
## Max. :6.585
## FIO datasets
We generate three subsets of the main dataset: one for each of the FIO groups:
coliform_data = data %>% filter(measure == "coliforms")
enterococci_data = data %>% filter(measure == "enterococci")
ecoli_data = data %>% filter(measure == "ecoli")Fitting bacterial datasets
We fit each dataset to the model described above using STAN. We define a single model in STAN code that will be used for each dataset. This describes the following data to be input for each fit:
N: the total number of datapoints (bacterial counts)Ntrt: the total number of treatments we are consideringy: a list of the bacterial counts at each timepointt: the corresponding timepoints for each measurementytrt: the corresponding treatments for each measurementy
and describes the model stated above in STAN code (see ../models/pooled_indicator_LP.stan).
Coliform fit
We first convert the coliform data into the appropriate format for the STAN model as coliform_trt_stan_data:
coliform_data$treatment = droplevels(coliform_data$treatment, "FS")
coliform_trt_stan_data = list(N = nrow(coliform_data),
y = coliform_data$logvalue,
t = coliform_data$day,
trt = as.integer(coliform_data$treatment),
Ntrt = nlevels(coliform_data$treatment))Then call STAN to fit the model, producing output in coliform_fit:
coliform_fit = stan(file = "../models/pooled_indicator_LP.stan",
data = coliform_trt_stan_data,
iter = 4000,
chains = 4)Once the model has been fit, we can inspect the fitted values:
summary_coliform = summary(coliform_fit, pars=c("m1", "c1", "m2", "c2",
"mu_m1", "sigma_m1", "mu_c1", "sigma_c1",
"mu_m2", "sigma_m2", "mu_c2", "sigma_c2"), probs = c(0.25, 0.75))
summary_coliform$summary## mean se_mean sd 25% 75%
## m1[1] -0.538622140 0.0006183026 0.022063101 -0.553269029 -0.523840099
## m1[2] -0.515936250 0.0005853137 0.022844778 -0.531284659 -0.501064332
## m1[3] -0.533492709 0.0007208585 0.023034567 -0.548538951 -0.518652030
## m1[4] -0.490334148 0.0008258817 0.025023223 -0.507571946 -0.472543991
## m1[5] -0.656299042 0.0007241618 0.026872436 -0.674603791 -0.638388272
## m1[6] -0.583525277 0.0017685899 0.030064980 -0.604013846 -0.562320620
## c1[1] 6.223247492 0.0020869455 0.087609826 6.168818751 6.281555700
## c1[2] 6.355757484 0.0055615954 0.090515065 6.292821590 6.416084253
## c1[3] 6.342024061 0.0047784095 0.090902929 6.278646944 6.400121311
## c1[4] 6.336775484 0.0033019160 0.087131132 6.275629549 6.391310401
## c1[5] 6.338616063 0.0047971239 0.091982880 6.274668936 6.398363293
## c1[6] 5.944903751 0.0120587583 0.133946384 5.849581168 6.030296390
## m2[1] -0.014451074 0.0004614339 0.008006746 -0.019415766 -0.009100185
## m2[2] -0.006785323 0.0002423941 0.010001312 -0.013321667 -0.001687031
## m2[3] -0.013044534 0.0006140444 0.009389064 -0.018394323 -0.007208506
## m2[4] -0.011475389 0.0002227714 0.008719108 -0.016664344 -0.006037301
## m2[5] -0.012949579 0.0002980180 0.009619629 -0.018344674 -0.006737304
## m2[6] -0.010972930 0.0002452669 0.009416411 -0.016646351 -0.005288013
## c2[1] 2.821306810 0.0096744969 0.163234623 2.706420122 2.927747847
## c2[2] 2.931418858 0.0052977181 0.222987618 2.810933703 3.077158227
## c2[3] 3.024382750 0.0155277174 0.215839048 2.890937763 3.151608734
## c2[4] 3.175547683 0.0049754659 0.201157882 3.046480710 3.296725431
## c2[5] 2.052220611 0.0055590110 0.221460732 1.910989777 2.183394736
## c2[6] 1.947915905 0.0047035275 0.214763002 1.819533350 2.082668463
## mu_m1 -0.529690941 0.0011684313 0.031319192 -0.549081430 -0.512191457
## sigma_m1 0.052100130 0.0008316895 0.041039395 0.026843605 0.064788384
## mu_c1 6.303017616 0.0034811185 0.102451034 6.244493910 6.360178162
## sigma_c1 0.134830415 0.0071247559 0.152836003 0.051395496 0.169780956
## mu_m2 -0.011746162 0.0002574586 0.007468915 -0.016274278 -0.006967018
## sigma_m2 0.005859133 0.0001945276 0.007074991 0.001544283 0.007522596
## mu_c2 2.783067099 0.0089114979 0.340935206 2.598203344 2.996411169
## sigma_c2 0.540091565 0.0131965860 0.400705374 0.279914498 0.696353001
## n_eff Rhat
## m1[1] 1273.3005 1.006983
## m1[2] 1523.3390 1.000791
## m1[3] 1021.0808 1.002212
## m1[4] 918.0172 1.004187
## m1[5] 1377.0283 1.003325
## m1[6] 288.9799 1.013845
## c1[1] 1762.3146 1.003142
## c1[2] 264.8762 1.013695
## c1[3] 361.9004 1.009782
## c1[4] 696.3293 1.006548
## c1[5] 367.6649 1.011238
## c1[6] 123.3834 1.032306
## m2[1] 301.0877 1.015906
## m2[2] 1702.4330 1.004788
## m2[3] 233.8003 1.019394
## m2[4] 1531.8821 1.003512
## m2[5] 1041.9136 1.003686
## m2[6] 1473.9844 1.005445
## c2[1] 284.6871 1.017336
## c2[2] 1771.6740 1.004331
## c2[3] 193.2168 1.022675
## c2[4] 1634.5816 1.000503
## c2[5] 1587.0784 1.001387
## c2[6] 2084.8344 1.001820
## mu_m1 718.4801 1.004462
## sigma_m1 2434.8907 1.000955
## mu_c1 866.1540 1.005788
## sigma_c1 460.1626 1.001627
## mu_m2 841.5894 1.011672
## sigma_m2 1322.7847 1.003240
## mu_c2 1463.6670 1.003346
## sigma_c2 921.9920 1.003666Plots
This can be visualised graphically, showing 80% (red, filled area) and 95% (distribution ends) CIs:
# Plot the early stage slopes
plot(coliform_fit, pars=c("m1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the early stage intercepts
plot(coliform_fit, pars=c("c1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage slopes
plot(coliform_fit, pars=c("m2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage intercepts
plot(coliform_fit, pars=c("c2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
E. coli fit
We first convert the E. coli data into the appropriate format for the STAN model as ecoli_trt_stan_data:
ecoli_data$treatment = droplevels(ecoli_data$treatment, "FS")
ecoli_trt_stan_data = list(N = nrow(ecoli_data),
y = ecoli_data$logvalue,
t = ecoli_data$day,
trt = as.integer(ecoli_data$treatment),
Ntrt = nlevels(ecoli_data$treatment))We call STAN to fit the model, storing the result in ecoli_fit:
ecoli_fit = stan(file = "../models/pooled_indicator_LP.stan",
data = ecoli_trt_stan_data,
iter = 4000,
chains = 4)We again inspect the fitted values:
summary_ecoli = summary(ecoli_fit,
pars=c("m1", "c1", "m2", "c2",
"mu_m1", "sigma_m1", "mu_c1", "sigma_c1",
"mu_m2", "sigma_m2", "mu_c2", "sigma_c2"),
probs = c(0.25, 0.75))
summary_ecoli$summary## mean se_mean sd 25% 75%
## m1[1] -0.510556575 0.0003913658 0.02532992 -0.527440165 -0.493557566
## m1[2] -0.481374091 0.0004521828 0.02673520 -0.499580423 -0.463577819
## m1[3] -0.445109444 0.0004531668 0.02893628 -0.464981074 -0.425482610
## m1[4] -0.361636735 0.0005170796 0.02884568 -0.381014240 -0.342541736
## m1[5] -0.635653166 0.0006735025 0.02880730 -0.655275031 -0.615654888
## m1[6] -0.566647097 0.0004328783 0.02807280 -0.584919635 -0.548460933
## c1[1] 5.922972071 0.0018573980 0.11012005 5.852246365 5.999022948
## c1[2] 5.947582457 0.0022891363 0.11569067 5.869246084 6.025352230
## c1[3] 4.993686617 0.0022403884 0.13651749 4.902291563 5.085807286
## c1[4] 4.906299543 0.0017780219 0.11339313 4.828473655 4.982087912
## c1[5] 6.187606151 0.0017237364 0.11611116 6.108017470 6.269568730
## c1[6] 5.804690390 0.0021306217 0.11229092 5.727609242 5.879881284
## m2[1] -0.037026169 0.0004748342 0.01340107 -0.046464179 -0.027451279
## m2[2] -0.013194676 0.0002138076 0.01109790 -0.020207550 -0.005868108
## m2[3] -0.014130327 0.0002428642 0.01062554 -0.021104315 -0.007349750
## m2[4] -0.008767309 0.0002802369 0.01154569 -0.016738311 -0.001862780
## m2[5] -0.013481098 0.0002516803 0.01121188 -0.021199589 -0.006149645
## m2[6] -0.013360509 0.0002444195 0.01174745 -0.021085859 -0.005952481
## c2[1] 2.875956877 0.0071838981 0.25265949 2.699834557 3.052010430
## c2[2] 2.701025397 0.0046603255 0.24748509 2.536052698 2.860705457
## c2[3] 2.394293394 0.0051202847 0.23446902 2.239980504 2.545730643
## c2[4] 2.355873510 0.0060901366 0.25295326 2.194696493 2.534237761
## c2[5] 2.002688914 0.0055458224 0.25337415 1.837422386 2.177743915
## c2[6] 1.999529908 0.0054226880 0.26612157 1.827559972 2.174140525
## mu_m1 -0.479412207 0.0007494363 0.04815333 -0.510505920 -0.450926257
## sigma_m1 0.097053819 0.0016145893 0.07284897 0.053114086 0.117666351
## mu_c1 5.777840671 0.0067144857 0.32844875 5.622759530 5.966049784
## sigma_c1 0.531515389 0.0093659891 0.41765738 0.270062969 0.658513431
## mu_m2 -0.014410829 0.0002499299 0.01001374 -0.020861351 -0.008001434
## sigma_m2 0.010206578 0.0002455893 0.01018128 0.003638571 0.013221948
## mu_c2 2.377245669 0.0053283061 0.29005722 2.198674222 2.559641083
## sigma_c2 0.404703275 0.0060810470 0.34503914 0.195805504 0.513957789
## n_eff Rhat
## m1[1] 4188.9174 1.0008363
## m1[2] 3495.7377 1.0007217
## m1[3] 4077.2693 1.0004043
## m1[4] 3112.0510 1.0015605
## m1[5] 1829.4758 1.0017150
## m1[6] 4205.7152 1.0005841
## c1[1] 3514.9798 1.0015875
## c1[2] 2554.1908 1.0021964
## c1[3] 3713.0427 1.0011949
## c1[4] 4067.2354 1.0004656
## c1[5] 4537.3910 1.0007945
## c1[6] 2777.6443 1.0002073
## m2[1] 796.5163 1.0068766
## m2[2] 2694.2342 1.0023375
## m2[3] 1914.1445 1.0014926
## m2[4] 1697.4180 1.0019176
## m2[5] 1984.5350 1.0020436
## m2[6] 2310.0181 1.0010603
## c2[1] 1236.9464 1.0048877
## c2[2] 2820.1069 1.0016184
## c2[3] 2096.9243 1.0016087
## c2[4] 1725.1486 1.0012178
## c2[5] 2087.3375 1.0015346
## c2[6] 2408.4128 1.0010976
## mu_m1 4128.4129 1.0002544
## sigma_m1 2035.7419 1.0017665
## mu_c1 2392.8183 1.0003689
## sigma_c1 1988.5340 1.0009751
## mu_m2 1605.2997 1.0021439
## sigma_m2 1718.6449 1.0040638
## mu_c2 2963.3914 1.0012429
## sigma_c2 3219.4375 0.9999462basic_fit_plots(fit=ecoli_fit, condition_names = ecoli_data$treatment, ignored_params = c("lp__", "y_hat", grep("mu", ecoli_fit@model_pars, value = T), grep("sigma", ecoli_fit@model_pars, value = T)))## [1] "making plot 1 of 4"## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.## [1] "making plot 2 of 4"## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## [1] "making plot 3 of 4"## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## [1] "making plot 4 of 4"## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Plots
# Plot the early stage slopes
plot(ecoli_fit, pars=c("m1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the early stage intercepts
plot(ecoli_fit, pars=c("c1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage slopes
plot(ecoli_fit, pars=c("m2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage intercepts
plot(ecoli_fit, pars=c("c2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
Enterococci fit
We convert the Enterococci data for the STAN model as enterococci_trt_stan_data:
enterococci_data$treatment = droplevels(enterococci_data$treatment, "FS")
enterococci_trt_stan_data = list(N = nrow(enterococci_data),
y = enterococci_data$logvalue,
t = enterococci_data$day,
trt = as.integer(enterococci_data$treatment),
Ntrt = nlevels(enterococci_data$treatment))We call STAN to fit the model, storing the result in ecoli_fit:
enterococci_fit = stan(file = "../models/pooled_indicator_LP.stan",
data = enterococci_trt_stan_data,
iter = 4000,
chains = 4)We again inspect the fitted values:
summary_enterococci = summary(enterococci_fit, pars=c("m1", "c1", "m2", "c2",
"mu_m1", "sigma_m1", "mu_c1", "sigma_c1",
"mu_m2", "sigma_m2", "mu_c2", "sigma_c2"), probs = c(0.25, 0.75))
summary_enterococci$summary## mean se_mean sd 25% 75%
## m1[1] -0.281344262 0.0005728545 0.023538817 -0.297203935 -0.265354282
## m1[2] -0.270136431 0.0007756025 0.028068005 -0.289579435 -0.251512035
## m1[3] -0.287942826 0.0007777311 0.028056960 -0.305986947 -0.270834705
## m1[4] -0.278660120 0.0006257326 0.026985185 -0.297307035 -0.261465893
## m1[5] -0.361698485 0.0007724236 0.033754585 -0.384615845 -0.338840314
## m1[6] -0.321057506 0.0006325509 0.028959412 -0.339979906 -0.301063735
## c1[1] 5.086841137 0.0021152516 0.087022030 5.028607759 5.142883286
## c1[2] 5.118927661 0.0027262618 0.110062692 5.044940413 5.184305885
## c1[3] 4.991043121 0.0036229274 0.136035714 4.917407370 5.083471737
## c1[4] 5.068472139 0.0027486149 0.106672084 5.004770690 5.134062726
## c1[5] 5.069769634 0.0027022986 0.110973615 5.004636373 5.138153293
## c1[6] 5.093735484 0.0027695099 0.110399728 5.021393133 5.157349186
## m2[1] -0.052135015 0.0003124614 0.010355567 -0.058309880 -0.045319309
## m2[2] -0.053904946 0.0003629130 0.011691357 -0.060122459 -0.046208362
## m2[3] -0.052175207 0.0002831769 0.009717227 -0.058384361 -0.045912647
## m2[4] -0.050178851 0.0002941799 0.009764337 -0.056545548 -0.044107964
## m2[5] -0.051774252 0.0002737027 0.010451366 -0.058261597 -0.045621893
## m2[6] -0.048988818 0.0003314735 0.011570556 -0.056555336 -0.043081949
## c2[1] 3.911131575 0.0072128932 0.238877535 3.752137261 4.053524880
## c2[2] 3.989394404 0.0085153116 0.272057513 3.807359357 4.137205762
## c2[3] 3.736027367 0.0063023549 0.217914938 3.596411218 3.876363825
## c2[4] 3.815496989 0.0066956031 0.219532563 3.672494385 3.960022742
## c2[5] 3.508666383 0.0062436442 0.243374586 3.365736302 3.667866112
## c2[6] 3.544508711 0.0079051413 0.270125916 3.400039993 3.722971626
## mu_m1 -0.298833755 0.0008366874 0.027840963 -0.313444839 -0.281313108
## sigma_m1 0.036671185 0.0007722946 0.034717437 0.015022104 0.047287810
## mu_c1 5.075148920 0.0024130327 0.089189334 5.016838567 5.131553590
## sigma_c1 0.072607996 0.0022674952 0.080673112 0.020918117 0.093644191
## mu_m2 -0.051451051 0.0002684489 0.008892604 -0.057297518 -0.045742911
## sigma_m2 0.005288959 0.0001804567 0.006556395 0.001440616 0.006557658
## mu_c2 3.759578339 0.0063630018 0.227262506 3.616276648 3.901997848
## sigma_c2 0.227832877 0.0068178952 0.230712721 0.088555057 0.288552293
## n_eff Rhat
## m1[1] 1688.420 1.004024
## m1[2] 1309.618 1.002537
## m1[3] 1301.434 1.002709
## m1[4] 1859.830 1.001608
## m1[5] 1909.653 1.001452
## m1[6] 2095.987 1.000927
## c1[1] 1692.522 1.002753
## c1[2] 1629.841 1.003384
## c1[3] 1409.895 1.003340
## c1[4] 1506.169 1.003729
## c1[5] 1686.447 1.001258
## c1[6] 1589.023 1.003158
## m2[1] 1098.386 1.001541
## m2[2] 1037.827 1.000816
## m2[3] 1177.522 1.000260
## m2[4] 1101.691 1.000174
## m2[5] 1458.102 1.000115
## m2[6] 1218.459 1.000715
## c2[1] 1096.811 1.001677
## c2[2] 1020.752 1.000997
## c2[3] 1195.552 1.000815
## c2[4] 1075.025 1.000257
## c2[5] 1519.407 1.000819
## c2[6] 1167.651 1.001539
## mu_m1 1107.241 1.002433
## sigma_m1 2020.828 1.001438
## mu_c1 1366.153 1.003081
## sigma_c1 1265.799 1.002471
## mu_m2 1097.324 1.000144
## sigma_m2 1320.031 1.001940
## mu_c2 1275.650 1.001105
## sigma_c2 1145.097 1.001738Plots
# Plot the early stage slopes
plot(enterococci_fit, pars=c("m1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the early stage intercepts
plot(enterococci_fit, pars=c("c1"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage slopes
plot(enterococci_fit, pars=c("m2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
# Plot the late stage intercepts
plot(enterococci_fit, pars=c("c2"), show_density=TRUE)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
Interpretation
As a reminder, we have two main research questions:
- Are die-off rates similar between the 10L and 50ml reactors?
- Do any of the 4 pasteurization schemes have an effect on die-off?
Days 0 – 7
Lets start by looking at the combined results for the first phase (day 0 to day 7):
for (param in c("m1", "c1")){
plot_all_fits(fits = c("E. coli" = ecoli_fit,
"Enterococci" = enterococci_fit,
"Coliforms" = coliform_fit),
param=param,
condition_names = ecoli_data$treatment)
}
From these figures, we can make a few assertions about this initial phase.
- There is no difference in the rate or the intercept between the 10L condition and the 50ml (unpasteurized) from days 0-7.
- There is no difference between the rate of die-off in either pre-pasteurization conditions
- Both pre-pasteurization conditions show about 1 log lower E. coli counts than the other conditions, but no change in the other two FIO.
- There appears to be a slight decrease in starting coliforms in the P4 post-pasteurized condition, but as the 95% CIs overlap with the 10L reactors, we cannot confirm that this difference is relevant.
- The P3 post-pasteurization treatment appears to have a decreased rate of die-off of E. coli and coliforms.
- There is no rate difference for P4.
From this, we can confirm that there is no difference in die-off between the 50ml batches and the 10L lab-scale bioreactors. Pre and Post-pasteurization does not effect all FIO equally; pre-pasteurization appears to effect initial E. coli counts negatively, though the rate of decline is ok; Enterococci and Coliforms are unaffected by pre-pasteurization. Post-pasteurization condition P3 shows statistically decreased rate of E. coli and Coliform die-off.
Days 7 – 28
Now lets plot the results for the late stage, from day 7 - day 28
for (param in c("m2", "c2")){
plot_all_fits(fits = c("E. coli" = ecoli_fit,
"Enterococci" = enterococci_fit,
"Coliforms" = coliform_fit),
param=param,
condition_names = ecoli_data$treatment)
}
We can conclude the following:
- There is no difference in the rate or the intercept between the 10L condition and the 50ml (unpasteurized) from days 7-28.
- For E. coli and the coliforms, the slope’s CIs include 0. Thus, we can conclude there is essentially no change after day 7. The 10L E. coli rate does not include 0, indicating that there is a slight reduction of E. coli over time.
- Enterococci decreases slightly over time, equally across all the conditions.
- Post pasteurization shows a significantly lower intercept comparing coliform unpasteurized conditions to the post-pasteurized treatments. This trend was similar for the other FIO, but not significant.
- Pre-pasteurization appears to have no effect on the intercept.