WT Control Series Equivalence

Two wild-type/WT variant datasets were collected, one as part of (i) the solvent-exchange experiment series, and (ii) for the mutant experiment series.

Our question can be phrased as:

Are there systematic differences between these two WT control series that would prevent us pooling the data from the two series in any analysis?

Essentially, the question asks “can we ascribe a plausible non-zero difference in WT maximum amplitude or decay rate to the effect of a Series/batch parameter for the experiment, conditioned on LPR (which we anticipate affects both measurements)?”

Neither answer, yes or no, would prevent downstream statistical analysis, but it would inform how we use the control data.

The two WT datasets will always be reasonable controls for their corresponding experimental series, but if there is no effect ascribable to Series, we would have some justification for pooling the data. If there is a discernible effect due to Series we can incorporate Series as an explanatory variable and estimate its effect separately from that of the treatments we are interested in, should we require to pool the data.

The Dataset

We process the NeRh50_Table_S1_raw_data_used_in_figures.csv dataset, to generate two dataframes - one describing maximum amplitude values, the other describing rate data, with a column indicating experimental series (Table 1).

# Use helper script to load dataframes
source("assets/r/load_dataframes.R")
dfms <- load_dataframes()

# Show interactive tables of data
DT::datatable(dfms$dfm_wt_amp)
DT::datatable(dfms$dfm_wt_k)
# Write datasets to file (commented out for website)
#readr::write_csv(dfms$dfm_wt_amp, "assets/data/00_WT_controls_amplitude.csv")
#readr::write_csv(dfms$dfm_wt_k, "assets/data/00_WT_controls_rate.csv")

Table 1: Measurements for WT across two series. Values are measured at two values of LPR: LPR5 and LPR10

(a) Maximum amplitude

(b) Decay rate

We visualise the distributions of datapoints between series, for each LPR value (Figure 1).

p1_amp <- ggplot2::ggplot(dfms$dfm_wt_amp, ggplot2::aes(x=Series, y=Value, color=Series)) +
  ggplot2::geom_boxplot(outlier.shape = NA) +
  ggplot2::geom_jitter() +
  ggplot2::facet_wrap(~LPR) +
  ggplot2::labs(y="Maximum amplitude (nA)", x="") +
  ggplot2::theme_minimal() + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 0.5, hjust=1))

# Visualise datasets as box/whisker with superimposed datapoints
p1_k <- ggplot2::ggplot(dfms$dfm_wt_k, ggplot2::aes(x=Series, y=Value, color=Series)) +
  ggplot2::geom_boxplot(outlier.shape = NA) +
  ggplot2::geom_jitter() +
  ggplot2::facet_wrap(~LPR) +
  ggplot2::labs(y="Decay rate constant k (/s)", x="") +
  ggplot2::theme_minimal() + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 0.5, hjust=1))

ggpubr::ggarrange(p1_amp, p1_k)

Figure 1: WT measurements across LPR5 and LPR10 for two batches/series

Visually we might expect from Figure 1 that the means of maximum amplitude between Series do not differ. We might have similar expectations for decay rate at LPR10 (though there is substantial variance), but certainly not at LPR5.

An apparent difference in decay rate at LPR5 may induce an interaction effect in any linear model fit. There also appears, visually, to be a difference in variance (heteroscedasticity) between the two series: Series 2 looks to have greater variance than Series 1 in general.

Downloads

• Download 00_WT_controls_amplitude.csv
• Download 00_WT_controls_rate.csv

The Question

Our central question is whether there is a nonzero effect on maximum amplitude or decay rate constant in WT, conditioned on LPR value, due to the choice of Series. If the inferred coefficient corresponding to Series has a confidence interval plausibly including zero, we will conclude that there is no effect on the measured WT properties due to the experimental series.

The Model

We address the question using a linear model, treating LPR and Series as main effects. As we suspect from Figure 1 that the effects of Series may not be independent of LPR, we use a linear model with interactions.

model1_amp_lm <- lm(Value ~ LPR * Series, data=dfms$dfm_wt_amp)
model1_k_lm <- lm(Value ~ LPR * Series, data=dfms$dfm_wt_k)

gtsummary::tbl_regression(model1_amp_lm)
gtsummary::tbl_regression(model1_k_lm)

Table 2: Linear model fits of Value ~ LPR * Series

(a) Maximum amplitude

Characteristic	Beta	95% CI	p-value
LPR
LPR5	—	—
LPR10	-0.91	-1.1, -0.71	<0.001
Series
Series 1	—	—
Series 2	-0.20	-0.43, 0.04	0.094
LPR * Series
LPR10 * Series 2	0.15	-0.19, 0.49	0.4
Abbreviation: CI = Confidence Interval

(b) Decay rate

Characteristic	Beta	95% CI	p-value
LPR
LPR5	—	—
LPR10	-11	-20, -2.2	0.016
Series
Series 1	—	—
Series 2	15	5.4, 25	0.004
LPR * Series
LPR10 * Series 2	-12	-27, 1.9	0.086
Abbreviation: CI = Confidence Interval

We interpret these model results (Table 2) to indicate that there is a reduction of maximum amplitude (-0.91nA, 95%CI[-1.1, -0.71]) and decay rate constant (-11/s, 95%CI[-20, -2.2]) on increasing LPR from LPR5 to LPR10.

The model fit indicates no effect on maximum amplitude due to the Series (-0.20nA [-0.43, -0.04]). However, Series does appear to be associated with a difference in decay rate constant (15/s, 95%CI[5.4, 25]), at a similar magnitude to the effect of increasing LPR.

Note

Modifying LPR appears to affect amplitude and rate constant.

Modifying Series appears to affect decay rate constant, but not maximum amplitude.

There is, for both models, no strong evidence of an effect due to the way Series interacts with LPR.

Checking Assumptions

We have so far only assumed that our data are suitable for linear modelling. We check whether the data meet requirements in Figure 2.

performance::check_model(model1_amp_lm, size_title=6)

performance::check_model(model1_k_lm, size_title=6)

(a) Maximum amplitude

(b) Decay rate

Figure 2: Diagnostic plots of Value ~ LPR + Dataset_use fit

While both models are acceptably linear with a normal distribution of residuals (Figure 2), there is a suggestion of heteroscedasticity, particularly for decay rate data (though the curve fit is suspect). To be cautious we therefore calculate heteroscedasticity-consistent standard errors, and recalculate effect sizes for the model parameters in Table 3.

parameters::model_parameters(model1_amp_lm, vcov="HC3") |>
  kableExtra::kbl(digits=3) |>
  kableExtra::kable_styling(
    bootstrap_options = c("hover", "responsive")
    )
parameters::model_parameters(model1_k_lm, vcov="HC3") |>
  kableExtra::kbl(digits=3) |>
  kableExtra::kable_styling(
    bootstrap_options = c("hover", "responsive")
    )

Table 3: Revised parameters for each model, taking into account heteroscedasticity

(a) Maximum amplitude

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	2.027	0.053	0.95	1.918	2.136	38.381	25	0.000
LPRLPR10	-0.914	0.083	0.95	-1.085	-0.742	-10.993	25	0.000
SeriesSeries 2	-0.196	0.162	0.95	-0.530	0.138	-1.211	25	0.237
LPRLPR10:SeriesSeries 2	0.151	0.194	0.95	-0.248	0.551	0.779	25	0.443

(b) Decay rate

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	44.332	1.024	0.95	42.223	46.442	43.285	25	0.000
LPRLPR10	-11.023	1.176	0.95	-13.444	-8.601	-9.375	25	0.000
SeriesSeries 2	15.269	6.086	0.95	2.735	27.804	2.509	25	0.019
LPRLPR10:SeriesSeries 2	-12.459	10.224	0.95	-33.516	8.598	-1.219	25	0.234

Taking into account heteroscedasticity, the model output (Table 3) now indicates that, for both maximum amplitude (-0.914nA, 95%CI[-1.08,-0.742]) and decay rate constant (-11.0/s, 95%CI[-13.4, -8.60]) there is a decrease with increasing LPR.

For maximum amplitude the model indicates no difference associated with Series, but there is a difference between Series for decay rate constant (15.3/s [2.74, 27.8]).

Note

Taking heteroscedasticity into account there is no strong evidence of interaction between terms for maximum amplitude or decay rate.

The contrasts of interest at each LPR are presented in Table 4.

EMM <- emmeans::emmeans(model1_amp_lm, ~ Series | LPR, vcov=sandwich::vcovHC(model1_amp_lm))
confint(emmeans::contrast(EMM, "pairwise", adjust="sidak")) |>
  kableExtra::kbl() |>
  kableExtra::kable_styling(
    bootstrap_options = c("hover", "responsive")
  )
EMM <- emmeans::emmeans(model1_k_lm, ~ Series | LPR, vcov=sandwich::vcovHC(model1_k_lm))
confint(emmeans::contrast(EMM, "pairwise", adjust="sidak")) |>
  kableExtra::kbl() |>
  kableExtra::kable_styling(
    bootstrap_options = c("hover", "responsive")
  )

Table 4: Contrasts for Series, accounting for heteroscedasticity, with Sidak adjustments.

(a) Maximum amplitude

contrast	LPR	estimate	SE	df	lower.CL	upper.CL
Series 1 - Series 2	LPR5	0.1963889	0.1622254	25	-0.1377207	0.5304984
Series 1 - Series 2	LPR10	0.0452222	0.1063565	25	-0.1738232	0.2642676

(b) Decay rate

contrast	LPR	estimate	SE	df	lower.CL	upper.CL
Series 1 - Series 2	LPR5	-15.269444	6.086282	25	-27.80438	-2.734511
Series 1 - Series 2	LPR10	-2.810444	8.215022	25	-19.72960	14.108711

The contrasts in Table 4 indicate that there is no difference in maximum amplitude associated with Series at either LPR. There is also no difference in decay rate between the two series at LPR10, but at LPR5 an difference is seen (-15.3/s [-27.8,-2.73]).

Figure

We present these relationships as in the manuscript, in Figure 3.

# Get p-values, accounting for heteroscedasticity and multiple testing
# Also generate group1 and group2 columns
EMM_amp <- as.data.frame(pairs(emmeans::emmeans(model1_amp_lm, ~ Series | LPR, 
                                                vcov=sandwich::vcovHC(model1_amp_lm)),
                               adjust="sidak")) |>
  tidyr::separate(contrast, into=c("group1", "group2"), sep=" - ") |>
  dplyr::mutate(y.position=3) |>
  rstatix::add_significance()
EMM_k <- as.data.frame(pairs(emmeans::emmeans(model1_k_lm, ~ Series | LPR, 
                                              vcov=sandwich::vcovHC(model1_k_lm)),
                               adjust="sidak")) |>
  tidyr::separate(contrast, into=c("group1", "group2"), sep=" - ") |>
  dplyr::mutate(y.position=90) |>
  rstatix::add_significance()

p1_amp_annot <- p1_amp +
  ggpubr::stat_pvalue_manual(EMM_amp, label="p = {scales::pvalue(p.value)}", size=3)

p1_k_annot <- p1_k +
  ggpubr::stat_pvalue_manual(EMM_k, label="p = {scales::pvalue(p.value)}", size=3)

p_comp <- ggpubr::ggarrange(p1_amp_annot, p1_k_annot)
p_comp

# Write figures to file (commented out for webpage)
#ggplot2::ggsave("assets/images/fig00_WT_measurements_for_two_series_with_stats.pdf", p_comp)
#ggplot2::ggsave("assets/images/fig00_WT_measurements_for_two_series_with_stats.png", p_comp)

Figure 3: WT measurements across LPR5 and LPR10 for two batches/series, Sidak-adjusted p-values indicated

Downloads

• Download as .pdf
• Download as .png