Bayesian GLM Part8
1 Preparations
Load the necessary libraries
2 Scenario
To investigate the effects of low light and a herbicide (diuron) on algae, King et al. (2022) set up an experiment in which algae were exposed to different combinations of light levels and diuron concentration. Over a 72 hour period, chlorophyll-a fluorescence and cell density were measured.
The data are in the file king.csv in the data folder.
| Time | Light | Diuron | Block | t0 | celld | lncelld | PI |
|---|---|---|---|---|---|---|---|
| 72 | 20 | 0.11 | 2 | 35.24612 | 125.94344 | 4.835832917 | 1.012390192 |
| 72 | 20 | 0.11 | 3 | 39.6347 | 114.65852 | 4.74195832 | 0.929296694 |
| 72 | 20 | 0.11 | 4 | 34.61918 | 167.32148 | 5.119916992 | 0.977363859 |
| 72 | 20 | 0.33 | 2 | 35.4551 | 126.77936 | 4.842448253 | 1.02771364 |
| 72 | 20 | 0.33 | 3 | 45.69512 | 155.20064 | 5.044718731 | 0.95016989 |
| 72 | 20 | 0.33 | 4 | 46.53104 | 162.72392 | 5.092055022 | 0.972060034 |
| 72 | 20 | 1 | 2 | 48.62084 | 165.85862 | 5.111135739 | 0.914524138 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| Time: | Hour of measurements |
| Light: | Light level treatment (5, 20, 80μmol photons m-2 s-1) - Predictor variable |
| Diuron: | Diuron herbicide concentration (0, 0.11, 0.33, 1, 3μg/L) - Predictor variable |
| Block: | Block - we will ignore this for now |
| t0: | Cell density at time 0 |
| celld | Cell density at treatment level - Response variable |
| PI | Photosynthetic inhibition (%), photosynthetic yield (as a % of controls) - Response variable |
3 Read in the data
Rows: 45 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (8): Time, Light, Diuron, Block, t0, celld, lncelld, PI
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 45
Columns: 8
$ Time <dbl> 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72…
$ Light <dbl> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 5,…
$ Diuron <dbl> 0.11, 0.11, 0.11, 0.33, 0.33, 0.33, 1.00, 1.00, 1.00, 3.00, 3.…
$ Block <dbl> 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4,…
$ t0 <dbl> 35.24612, 39.63470, 34.61918, 35.45510, 45.69512, 46.53104, 48…
$ celld <dbl> 125.94344, 114.65852, 167.32148, 126.77936, 155.20064, 162.723…
$ lncelld <dbl> 4.835833, 4.741958, 5.119917, 4.842448, 5.044719, 5.092055, 5.…
$ PI <dbl> 1.0123902, 0.9292967, 0.9773639, 1.0277136, 0.9501699, 0.97206…# A tibble: 6 × 8
Time Light Diuron Block t0 celld lncelld PI
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 72 20 0.11 2 35.2 126. 4.84 1.01
2 72 20 0.11 3 39.6 115. 4.74 0.929
3 72 20 0.11 4 34.6 167. 5.12 0.977
4 72 20 0.33 2 35.5 127. 4.84 1.03
5 72 20 0.33 3 45.7 155. 5.04 0.950
6 72 20 0.33 4 46.5 163. 5.09 0.972spc_tbl_ [45 × 8] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ Time : num [1:45] 72 72 72 72 72 72 72 72 72 72 ...
$ Light : num [1:45] 20 20 20 20 20 20 20 20 20 20 ...
$ Diuron : num [1:45] 0.11 0.11 0.11 0.33 0.33 0.33 1 1 1 3 ...
$ Block : num [1:45] 2 3 4 2 3 4 2 3 4 2 ...
$ t0 : num [1:45] 35.2 39.6 34.6 35.5 45.7 ...
$ celld : num [1:45] 126 115 167 127 155 ...
$ lncelld: num [1:45] 4.84 4.74 5.12 4.84 5.04 ...
$ PI : num [1:45] 1.012 0.929 0.977 1.028 0.95 ...
- attr(*, "spec")=
.. cols(
.. Time = col_double(),
.. Light = col_double(),
.. Diuron = col_double(),
.. Block = col_double(),
.. t0 = col_double(),
.. celld = col_double(),
.. lncelld = col_double(),
.. PI = col_double()
.. )
- attr(*, "problems")=<externalptr> king (45 rows and 8 variables, 8 shown)
ID | Name | Type | Missings | Values | N
---+---------+---------+----------+-----------------+------------
1 | Time | numeric | 0 (0.0%) | 72 | 45 (100.0%)
---+---------+---------+----------+-----------------+------------
2 | Light | numeric | 0 (0.0%) | 5 | 15 ( 33.3%)
| | | | 20 | 15 ( 33.3%)
| | | | 80 | 15 ( 33.3%)
---+---------+---------+----------+-----------------+------------
3 | Diuron | numeric | 0 (0.0%) | 0 | 9 ( 20.0%)
| | | | 0.11 | 9 ( 20.0%)
| | | | 0.33 | 9 ( 20.0%)
| | | | 1 | 9 ( 20.0%)
| | | | 3 | 9 ( 20.0%)
---+---------+---------+----------+-----------------+------------
4 | Block | numeric | 0 (0.0%) | 2 | 15 ( 33.3%)
| | | | 3 | 15 ( 33.3%)
| | | | 4 | 15 ( 33.3%)
---+---------+---------+----------+-----------------+------------
5 | t0 | numeric | 0 (0.0%) | [29.19, 68.89] | 45
---+---------+---------+----------+-----------------+------------
6 | celld | numeric | 0 (0.0%) | [37.96, 270.56] | 45
---+---------+---------+----------+-----------------+------------
7 | lncelld | numeric | 0 (0.0%) | [3.64, 5.6] | 45
---+---------+---------+----------+-----------------+------------
8 | PI | numeric | 0 (0.0%) | [0.63, 1.53] | 45
-----------------------------------------------------------------4 Process data
5 Research questions
- Q1: estimate the effect of light on photosynthetic inhibition
- Q2: estimate the effect of photosynthetic inhibition on algal cell density
- Q3: estimate the effect of light on algal cell density
- Q4: estimate the effect of diuron (herbicide) on photosynthetic inhibition
- Q5: estimate the effect of diuron (herbicide) on algal cell density
6 Naive approach
g1 <-
king |>
ggplot(aes(y = PI, x = Diuron, colour = factor(Light))) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic()
g2 <-
king |>
ggplot(aes(y = celld, x = PI, colour = factor(Light))) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic()
g3 <-
king |>
ggplot(aes(y = celld, x = Diuron, colour = factor(Light))) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic()
g1 + g2 + g3 + plot_layout(guides = "collect")`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
form <- bf(PI ~ scale(Diuron) * factor(Light))
## summarise data to help inform priors
king |>
group_by(Light) |>
summarise(PI_mu = median(PI), PI_sd = sd(PI), PI_mad = mad(PI))# A tibble: 3 × 4
Light PI_mu PI_sd PI_mad
<dbl> <dbl> <dbl> <dbl>
1 5 1.36 0.241 0.256
2 20 0.950 0.112 0.0739
3 80 0.915 0.114 0.125 ## define priors
priors <- prior(normal(1.4, 0.3), class = "Intercept") +
prior(normal(0, 0.5), class = "b") +
prior(student_t(3, 0, 0.3), class = "sigma")
## fit the mode
king_brms0a <- brm(
form,
prior = priors,
data = king,
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: PI ~ scale(Diuron) * factor(Light)
Data: king (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 1.26 0.03 1.19 1.32 1.00 2946
scaleDiuron -0.11 0.04 -0.18 -0.05 1.00 2117
factorLight20 -0.34 0.05 -0.44 -0.24 1.00 3231
factorLight80 -0.39 0.05 -0.48 -0.29 1.00 3198
scaleDiuron:factorLight20 0.01 0.05 -0.09 0.11 1.00 2211
scaleDiuron:factorLight80 0.01 0.05 -0.09 0.11 1.00 2390
Tail_ESS
Intercept 2599
scaleDiuron 2340
factorLight20 3042
factorLight80 2715
scaleDiuron:factorLight20 2687
scaleDiuron:factorLight80 2367
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.13 0.02 0.11 0.17 1.00 3162 3160
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- increasing diuron concentration associated with decreasing photosynthetic inhibition ()
- increasing light is associated with decreasing photosynthetic inhibition
form <- bf(celld ~ scale(Diuron) + scale(PI) + Light)
## summarise data to help inform priors
king |>
group_by(Light) |>
summarise(celld_mu = median(celld), celld_sd = sd(celld), celld_mad = mad(celld))# A tibble: 3 × 4
Light celld_mu celld_sd celld_mad
<dbl> <dbl> <dbl> <dbl>
1 5 146. 38.0 33.8
2 20 163. 43.4 23.9
3 80 155. 30.1 20.8## define priors
priors <- prior(normal(150, 50), class = "Intercept") +
prior(normal(0, 20), class = "b") +
prior(student_t(3, 0, 50), class = "sigma")
## fit the mode
king_brms0b <- brm(
form,
prior = priors,
data = king |> mutate(Light = factor(Light)),
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: celld ~ scale(Diuron) + scale(PI) + Light
Data: mutate(king, Light = factor(Light)) (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 148.05 9.65 129.02 166.89 1.00 2798 2451
scaleDiuron -1.34 6.57 -14.31 11.55 1.00 3116 2770
scalePI -8.23 7.70 -23.70 6.52 1.00 2791 2748
Light20 12.93 13.26 -13.33 38.84 1.00 3161 3017
Light80 3.76 13.59 -22.81 29.61 1.00 2837 2809
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 38.89 4.42 31.56 48.48 1.00 3922 2861
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- no effect of diuron concentration on algal cell density
- no effect of photosynthetic inhibition on algal cell density
- no effect of light on algal cell density
7 Structural Causal Modelling
Create a conceptual model of the system of interest
Choose a statistical test
Test the consistency of your conceptual model with the data
Identify biases (confounding, overcontrol, collider) that you need to adjust for when testing the cause of interest
Test for causality with the appropriate statistical test
Repeat steps 3-5 for each cause of interest
Repeat steps 1-6 for other conceptual models
8 Step 1 - conceptual model
King et al. (2022) hypothesised that photosynthetic inhibition would be effected by light levels and that the diuron herbicide would effect photosynthesis by blocking certain chemical reactions. In tern, changes in photosynthesis could have an effect on algal cell density. Furthermore, diuron might also have a direct effect on cell density by causing cell damage.
Lets start by creating a DAG that reflects this conceptual model.
9 Step 2 - choose a statistical test
For the main inference tests we will fit Bayesian generalized linear models, however, for the purpose of testing the consistency of our conceptual model (DAG) with data, we will just use more simple frequentist analyses.
10 Step 3 - test the consistency of the DAG with observed data
Lets now determine what all the implied assumptions (conditional independencies) of our DAG are.
Conclusions:
- Diuron is independent of Light
- Light is independent of celld when conditioned on both Diuron and PI
And challenge those assumptions with the data.
Warning in cov2cor(cov(data)): diag(V) had non-positive or NA entries; the
non-finite result may be dubious estimate p.value 2.5% 97.5%
Dirn _||_ Lght 0.00000000 1.000000 -0.2936838 0.2936838
Lght _||_ clld | Dirn, PI -0.04685251 0.766819 -0.3427060 0.2571860Conclusions:
- it seems that our conceptual model is consistent with the data.
11 Step 4 - identify biases that need to be adjusted for
Our goal is to be able to investigate the effects of light and diuron on both algal photosynthesis and cell density. So lets explore each one separately, starting with photosynthesis.
ggdag_paths(king_dag, from = "Light", to = "PI", text_col = "black", shadow = TRUE) +
theme_dag_blank(panel.border = element_rect(fill = NA))
Conclusions:
- there is a single direct path
{}We can also ask whether there are any other optional covariates that we could include without introducing biases.
{ Diuron }- we could also include diuron
ggdag_paths(king_dag, from = "Light", to = "celld", text_col = "black", shadow = TRUE) +
theme_dag_blank(panel.border = element_rect(fill = NA))
Conclusions:
- there is a single path
- a pipe (
Light -> PI -> celld)
- a pipe (
adjustmentSets(king_dag,
exposure = "Light",
outcome = "celld",
type = "minimal",
effect = "total"
) {}- to explore the total effects of Light on celld, we do not condition on anything
We can also ask whether there are any other optional covariates that we could include without introducing biases.
adjustmentSets(king_dag,
exposure = "Light",
outcome = "celld",
type = "canonical",
effect = "total"
){ Diuron }- we could also include diuron
The conceptual model does not include a direct pathway from Light to cell density, but if we wanted to examine this, what would we need to condition on?
adjustmentSets(king_dag,
exposure = "Light",
outcome = "celld",
type = "minimal",
effect = "direct"
){ Diuron, PI }- we would need to include both Diuron and PI
ggdag_paths(king_dag, from = "PI", to = "celld", text_col = "black", shadow = TRUE) +
theme_dag_blank(panel.border = element_rect(fill = NA))
Conclusions:
- there is are two paths
- direct (
PI -> celld) - a forl (
PI <- Diuron -> celld)- this is a backdoor
- we want to block this path
- direct (
{ Diuron }- to explore the total effects of photosynthesis inhibition on cell density, we must condition on diuron
We can also ask whether there are any other optional covariates that we could include without introducing biases.
{ Diuron, Light }- we could also include light, but we already know that we must include diuron
The conceptual model also includes a direct pathway from PI to cell density, so if we wanted to examine this, what would we need to condition on?
{ Diuron }- again, as a minimum, we must condition on diuron
ggdag_paths(king_dag, from = "Diuron", to = "celld", text_col = "black", shadow = TRUE) +
theme_dag_blank(panel.border = element_rect(fill = NA))
Conclusions:
- there are two paths
- one is a pipe (
Diron -> PI -> celld) - one is direct (
Diuron -> celld)
- one is a pipe (
adjustmentSets(king_dag,
exposure = "Diuron",
outcome = "celld",
type = "minimal",
effect = "total"
) {}- to explore the total effects of Diuron on cell density, we do not condition on any other covariate
We can also ask whether there are any other optional covariates that we could include without introducing biases.
adjustmentSets(king_dag,
exposure = "Diuron",
outcome = "celld",
type = "canonical",
effect = "total"
){ Light }- we could also include Light
If we wanted to examine the direct effects of Diuron on cell density, what would we need to condition on?
adjustmentSets(king_dag,
exposure = "Diuron",
outcome = "celld",
type = "minimal",
effect = "direct"
){ PI }- we would need to include PI
ggdag_paths(king_dag, from = "Diuron", to = "PI", text_col = "black", shadow = TRUE) +
theme_dag_blank(panel.border = element_rect(fill = NA))
Conclusions:
- there are two paths
- one is a pipe (
Diron -> PI -> celld) - one is direct (
Diuron -> celld)
- one is a pipe (
{}- to explore the total effects of Diuron on cell density, we do not condition on any other covariate
We can also ask whether there are any other optional covariates that we could include without introducing biases.
adjustmentSets(king_dag,
exposure = "Diuron",
outcome = "PI",
type = "canonical",
effect = "total"
){ Light }- we could also include Light
If we wanted to examine the direct effects of Diuron on cell density, what would we need to condition on?
{}- we would need to include PI
Conclusions:
- Response: PI
- Light
- can include diuron
- Diuron
- can also have light
- Light
- Response: celld
- Light
- total: can include diuron
- direct: must include PI and diuron
- PI
- total: must include diuron
- total: can also include Light
- direct: must include diuron
- Diuron
- total: can include Light
- direct: must include PI
- Light
PI: ~ Light + Diuron Celld (total): ~ Light + Diuron
12 Step 5 - test for causality
12.0.1 Exploratory data analysis
form <- bf(PI ~ fLight * scale(Diuron))
## summarise data to help inform priors
king |>
group_by(fLight) |>
summarise(PI_mu = median(PI), PI_sd = sd(PI), PI_mad = mad(PI))# A tibble: 3 × 4
fLight PI_mu PI_sd PI_mad
<fct> <dbl> <dbl> <dbl>
1 5 1.36 0.241 0.256
2 20 0.950 0.112 0.0739
3 80 0.915 0.114 0.125 ## define priors
priors <- prior(normal(1.4, 0.3), class = "Intercept") +
prior(normal(0, 0.5), class = "b") +
prior(student_t(3, 0, 0.3), class = "sigma")
## fit the mode
king_brms1a <- brm(
form,
prior = priors,
data = king,
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: PI ~ fLight * scale(Diuron)
Data: king (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 1.26 0.03 1.18 1.32 1.00 3291
fLight20 -0.34 0.05 -0.43 -0.24 1.00 3413
fLight80 -0.38 0.05 -0.48 -0.28 1.00 3758
scaleDiuron -0.11 0.04 -0.18 -0.04 1.00 2458
fLight20:scaleDiuron 0.01 0.05 -0.09 0.11 1.00 2846
fLight80:scaleDiuron 0.01 0.05 -0.09 0.11 1.00 2652
Tail_ESS
Intercept 2804
fLight20 2576
fLight80 3132
scaleDiuron 2670
fLight20:scaleDiuron 2987
fLight80:scaleDiuron 3093
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.13 0.02 0.11 0.17 1.00 3574 2791
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- increasing light reduces photosynthetic inhibition
- increasing diuron reduces photosynthetic inhibition
12.0.2 Exploratory data analysis
form <- bf(celld ~ fLight * scale(Diuron))
## summarise data to help inform priors
king |>
group_by(fLight) |>
summarise(celld_mu = median(celld), celld_sd = sd(celld), celld_mad = mad(celld))# A tibble: 3 × 4
fLight celld_mu celld_sd celld_mad
<fct> <dbl> <dbl> <dbl>
1 5 146. 38.0 33.8
2 20 163. 43.4 23.9
3 80 155. 30.1 20.8## define priors
priors <- prior(normal(150, 50), class = "Intercept") +
prior(normal(0, 30), class = "b") +
prior(student_t(3, 0, 50), class = "sigma")
## fit the mode
king_brms2a <- brm(
form,
prior = priors,
data = king,
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: celld ~ fLight * scale(Diuron)
Data: king (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 139.18 7.82 124.41 154.62 1.00 3658
fLight20 26.06 11.04 3.96 47.55 1.00 3892
fLight80 17.08 10.85 -4.63 37.80 1.00 4003
scaleDiuron -19.02 7.80 -33.91 -3.77 1.00 2221
fLight20:scaleDiuron 32.30 11.07 10.05 53.06 1.00 2551
fLight80:scaleDiuron 33.94 10.97 11.58 55.41 1.00 2587
Tail_ESS
Intercept 3047
fLight20 2938
fLight80 3266
scaleDiuron 2654
fLight20:scaleDiuron 2935
fLight80:scaleDiuron 2600
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.39 3.89 26.67 41.73 1.00 3645 2917
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- increasing light increases algal cell density (although only at 20 units)
- increasing diuron concentration reduces algal cell density
- however, there is an interaction between diuron and light levels
12.0.3 Exploratory data analysis
form <- bf(celld ~ scale(PI) * fLight * scale(Diuron))
## summarise data to help inform priors
king |>
group_by(fLight) |>
summarise(celld_mu = median(celld), celld_sd = sd(celld), celld_mad = mad(celld))# A tibble: 3 × 4
fLight celld_mu celld_sd celld_mad
<fct> <dbl> <dbl> <dbl>
1 5 146. 38.0 33.8
2 20 163. 43.4 23.9
3 80 155. 30.1 20.8## define priors
priors <- prior(normal(150, 50), class = "Intercept") +
prior(normal(0, 30), class = "b") +
prior(student_t(3, 0, 50), class = "sigma")
## fit the mode
king_brms3a <- brm(
form,
prior = priors,
data = king,
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: celld ~ scale(PI) * fLight * scale(Diuron)
Data: king (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 143.48 10.75 122.50 165.35 1.00 4029
scalePI -6.71 8.84 -24.40 10.48 1.00 3343
fLight20 16.14 16.74 -17.59 49.26 1.00 3137
fLight80 2.87 18.23 -32.74 39.26 1.00 2977
scaleDiuron -22.76 8.97 -39.76 -4.43 1.00 4077
scalePI:fLight20 -18.44 25.39 -67.52 31.47 1.00 4272
scalePI:fLight80 0.22 22.32 -45.19 44.48 1.00 3956
scalePI:scaleDiuron -1.01 10.97 -23.25 20.66 1.00 3692
fLight20:scaleDiuron 33.06 19.61 -5.98 71.11 1.00 3311
fLight80:scaleDiuron 18.48 20.20 -21.10 57.90 1.00 3269
scalePI:fLight20:scaleDiuron 10.00 19.37 -28.09 47.66 1.00 3600
scalePI:fLight80:scaleDiuron -16.35 18.56 -51.42 20.67 1.00 3535
Tail_ESS
Intercept 3018
scalePI 2946
fLight20 2289
fLight80 2409
scaleDiuron 3029
scalePI:fLight20 3012
scalePI:fLight80 2757
scalePI:scaleDiuron 2723
fLight20:scaleDiuron 2958
fLight80:scaleDiuron 2739
scalePI:fLight20:scaleDiuron 2854
scalePI:fLight80:scaleDiuron 3145
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.76 4.11 26.81 42.84 1.00 3647 3169
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- no evidence that photosynthetic inhibition has an effect on algal cell density
- but ignore diuron and light
12.0.4 Exploratory data analysis
form <- bf(celld ~ fLight * scale(Diuron))
## summarise data to help inform priors
king |>
group_by(fLight) |>
summarise(celld_mu = median(celld), celld_sd = sd(celld), celld_mad = mad(celld))# A tibble: 3 × 4
fLight celld_mu celld_sd celld_mad
<fct> <dbl> <dbl> <dbl>
1 5 146. 38.0 33.8
2 20 163. 43.4 23.9
3 80 155. 30.1 20.8## define priors
priors <- prior(normal(150, 50), class = "Intercept") +
prior(normal(0, 30), class = "b") +
prior(student_t(3, 0, 50), class = "sigma")
## fit the mode
king_brms4a <- brm(
form,
prior = priors,
data = king,
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling Family: gaussian
Links: mu = identity; sigma = identity
Formula: celld ~ fLight * scale(Diuron)
Data: king (Number of observations: 45)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 139.23 8.16 123.67 155.23 1.00 3742
fLight20 26.08 11.35 3.29 48.46 1.00 3797
fLight80 17.05 11.18 -4.85 38.58 1.00 3406
scaleDiuron -18.89 8.12 -34.43 -2.58 1.00 2358
fLight20:scaleDiuron 31.95 11.49 8.91 53.64 1.00 2597
fLight80:scaleDiuron 33.69 11.20 11.07 55.12 1.00 2975
Tail_ESS
Intercept 2949
fLight20 3137
fLight80 3113
scaleDiuron 2315
fLight20:scaleDiuron 2489
fLight80:scaleDiuron 2899
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.44 4.03 26.60 42.62 1.00 3517 2702
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- increasing light levels increase algal cell density
- increasing diuron concentration decreases algal cell density
- there is evidence of an interaction between light and diuron concentration
References
King, Oliver C., Jason P. van de Merwe, Christopher J. Brown, Michael St J. Warne, and A. Smith Rachael. 2022. “Individual and Combined Effects of Diuron and Light Reduction on Marine Microalgae.” Ecotoxicology and Environmental Safety. https://doi.org/https://doi.org/10.1016/j.ecoenv.2022.113729.