Dichotomous endpoint case study
Step 1: Planning phase
This binary case study follows the EFSA stepwise framework for evaluating and
using historical control data (HCD). The endpoint is dichotomous: the presence
or absence of a pathological finding.
Step 2: Data acceptability
The HCD come from the same laboratory and animal strain as the index study.
Following the case-study protocol, HCD are first restricted to the 5-year period
centred around the index study date.
The data meet the statistical requirements for a dichotomous endpoint because,
for each control group, both the number of animals examined and the number with
the finding are available.
## Load packages used in this vignette
library(efsaHCD)
library(openxlsx)
## Locate the binary case-study files shipped with the package.
## These files are stored in inst/extdata/ in the source package.
hcd_file <- system.file(
"extdata",
"appendix_binary_hcd.xlsx",
package = "efsaHCD"
)
index_file <- system.file(
"extdata",
"appendix_binary_index-study.xlsx",
package = "efsaHCD"
)
## Read the historical control data and the index-study data
library(openxlsx)
hcd_raw <- as.data.frame(openxlsx::read.xlsx(hcd_file, check.names = FALSE))
index_raw <- as.data.frame(openxlsx::read.xlsx(index_file, check.names = FALSE))
## Inspect the imported data
str(hcd_raw)
## 'data.frame': 23 obs. of 9 variables:
## $ Date : num 32842 33604 33117 33116 33360 ...
## $ Route : chr "dt" "sc" "dt" "dt" ...
## $ Type : chr "car" "car" "car" "car" ...
## $ Diet : chr "lad2" "lad2" "lad2" "lad2" ...
## $ No./cage : num 5 5 5 5 5 5 5 5 5 5 ...
## $ Duration : num 104 104 104 104 104 104 92 104 104 104 ...
## $ No..of.animals : num 50 50 50 50 55 50 55 60 55 50 ...
## $ No..examined : num 50 49 50 49 54 49 55 56 55 50 ...
## $ No..with.finding: num 2 2 1 5 3 5 2 2 2 4 ...
str(index_raw)
## 'data.frame': 5 obs. of 10 variables:
## $ Date : num 33603 33603 33603 33603 33603
## $ Route : chr "dt" "dt" "dt" "dt" ...
## $ Type : chr "ct" "ct" "ct" "ct" ...
## $ Diet : chr "lad2" "lad2" "lad2" "lad2" ...
## $ No./cage : num 5 5 5 5 5
## $ Duration : num 104 104 104 104 104
## $ No..of.animals : num 50 50 50 50 50
## $ No..examined : num 47 40 44 49 49
## $ No..with.finding: num 3 4 5 7 8
## $ Dose : num 0 10 50 1500 10000
head(hcd_raw)
## Date Route Type Diet No./cage Duration No..of.animals No..examined
## 1 32842 dt car lad2 5 104 50 50
## 2 33604 sc car lad2 5 104 50 49
## 3 33117 dt car lad2 5 104 50 50
## 4 33116 dt car lad2 5 104 50 49
## 5 33360 dt ct lad1 5 104 55 54
## 6 33908 dt ct lad2 5 104 50 49
## No..with.finding
## 1 2
## 2 2
## 3 1
## 4 5
## 5 3
## 6 5
head(index_raw)
## Date Route Type Diet No./cage Duration No..of.animals No..examined
## 1 33603 dt ct lad2 5 104 50 47
## 2 33603 dt ct lad2 5 104 50 40
## 3 33603 dt ct lad2 5 104 50 44
## 4 33603 dt ct lad2 5 104 50 49
## 5 33603 dt ct lad2 5 104 50 49
## No..with.finding Dose
## 1 3 0
## 2 4 10
## 3 5 50
## 4 7 1500
## 5 8 10000
dim(hcd_raw)
## [1] 23 9
dim(index_raw)
## [1] 5 10
## Create short analysis variable names used throughout the vignette.
## k = number of animals with the finding
## n = number of animals examined
hcd <- hcd_raw
index_study <- index_raw
hcd$k <- hcd[["No..with.finding"]]
hcd$n <- hcd[["No..examined"]]
index_study$k <- index_study[["No..with.finding"]]
index_study$n <- index_study[["No..examined"]]
## Convert the study date to Date format if needed.
## Excel dates may already be imported as Date objects, but this keeps the code robust.
hcd$Date <- as.Date(hcd$Date) - lubridate::years(70)
index_study$Date <- as.Date(index_study$Date) - lubridate::years(70)
## Step 2a: restrict HCD to the 5-year period centred around the index-study date.
## The index-study date is 31/12/1991, so the 5-year window is approximately
## 30/06/1989 to 30/06/1994.
hcd_5y <- hcd[
hcd$Date >= as.Date("1989-06-30") &
hcd$Date <= as.Date("1994-06-30"),
]
## Check the number of HCD groups before and after the 5-year restriction
nrow(hcd)
## [1] 23
nrow(hcd_5y)
## [1] 17
## Check the index-study groups
index_study
## Date Route Type Diet No./cage Duration No..of.animals No..examined
## 1 1992-01-01 dt ct lad2 5 104 50 47
## 2 1992-01-01 dt ct lad2 5 104 50 40
## 3 1992-01-01 dt ct lad2 5 104 50 44
## 4 1992-01-01 dt ct lad2 5 104 50 49
## 5 1992-01-01 dt ct lad2 5 104 50 49
## No..with.finding Dose k n
## 1 3 0 3 47
## 2 4 10 4 40
## 3 5 50 5 44
## 4 7 1500 7 49
## 5 8 10000 8 49
Step 3: Preliminary investigation of HCD
In Step 3, the selected historical control groups are investigated statistically
before modelling HCD variability. The objectives are:
- to explore the distribution of sample incidences,
- to identify potential outlying studies,
- and to assess whether the endpoint is associated with study-level covariates
listed in the protocol.
For dichotomous endpoints, these investigations are performed using exploratory
plots and generalized linear models (GLMs) based on binomial or quasi-binomial
distributions, following the EFSA case-study approach.
The outcome of Step 3 is the final selection of HCD groups that will be used for
modelling between-study variability in Step 4.
The following boxplot summarises the distribution of sample incidences across
the 17 historical control groups included in the 5-year subset. Most studies
show incidences between approximately 2% and 7%, although some lower and higher
incidences are also observed. No historical control group appears to be an
extreme outlier based on the conventional 1.5 × IQR criterion. The symbols
indicate the number of animals with the finding in each study.
## Step 3a: exploratory investigation of HCD
##
## We first inspect the distribution of sample incidence k / n.
## The boxplot includes the conventional 1.5 * IQR fences, and the points
## show individual HCD groups. Point shape indicates the number of animals
## with the finding.
binary_step3a_plots <- make_binary_hcd_plots(
dat = hcd_5y,
k_col = "k",
n_col = "n",
covariate_cols = c(
"Route",
"Type",
"Diet",
"No./cage",
"Duration",
"Date"
),
dataset_label = "5-year HCD subset"
)
binary_step3a_plots$boxplot
The exploratory covariate plots do not suggest a clear association between
sample incidence and route of exposure, study type, diet, number of animals per
cage, or study start date. In contrast, study duration appears to be associated
with the observed incidence, with shorter studies tending to show lower
incidences than the 104-week studies.
## Combine the covariate plots in one figure, similar to the EFSA report
library(patchwork)
p_covariates <-
(
binary_step3a_plots$scatterplots$Route +
binary_step3a_plots$scatterplots$Type
) /
(
binary_step3a_plots$scatterplots$Diet +
binary_step3a_plots$scatterplots$`No./cage`
) /
(
binary_step3a_plots$scatterplots$Duration +
binary_step3a_plots$scatterplots$Date
)
p_covariates
The formal screening uses one univariate GLM per covariate. For each model, the
binary outcome is represented by the number of animals with the finding (k)
and the number without the finding (n - k).
First, a binomial GLM is fitted. The dispersion statistic is then used to assess
whether the variability in the data is larger than expected under a binomial
model. When there is evidence of overdispersion, a quasi-binomial GLM is also
useful because it allows the variance to be inflated relative to the binomial
variance.
The screening table reports the association between sample incidence and each
study-level covariate. These results are used to decide whether any covariate
should influence the selection of HCD before Step 4. The following output shows
a formatted summary table. The underlying raw results can be accessed using
binary_step3a_screen$summary_table.
binary_step3a_screen <- screen_binary_hcd_covariates(
dat = hcd_5y,
k_col = "k",
n_col = "n",
covariate_cols = c(
"Route",
"Type",
"Diet",
"No./cage",
"Duration",
"Date"
),
grouping_var = NULL
)
binary_step3a_screen$formatted_summary_table
## Covariate Preferred p-value
## 1 Route 0.839
## 2 Type 0.885
## 3 Diet 0.763
## 4 No./cage 0.902
## 5 Duration 0.008
## 6 Date 0.450
The coefficient table can be accessed using
binary_step3a_screen$coefficient_table. The results show no evidence
of association for route of exposure, study type, diet, number of animals per
cage, or study start date. In contrast, study duration shows evidence of
association with incidence. The estimated coefficient for Duration is
positive, indicating that longer studies tend to have higher observed incidence.
This supports the case-study decision to exclude shorter-duration studies before
fitting the HCD variability model. The following shows a formatted version of
this table.
binary_step3a_screen$formatted_coefficient_table
## covariate model term beta sd p_value
## 1 Route binomial (Intercept) -3.004444e+00 0.1619823348 8.468349e-77
## 2 Route binomial Routesc -1.525562e-01 0.7399424750 8.366557e-01
## 3 Route quasibinomial (Intercept) -3.004444e+00 0.1647882796 1.200757e-11
## 4 Route quasibinomial Routesc -1.525562e-01 0.7527601547 8.421222e-01
## 5 Type binomial (Intercept) -3.032546e+00 0.2089823909 1.031850e-47
## 6 Type binomial Typect 4.798727e-02 0.3194310076 8.805852e-01
## 7 Type quasibinomial (Intercept) -3.032546e+00 0.2133202675 4.132147e-10
## 8 Type quasibinomial Typect 4.798727e-02 0.3260614815 8.849558e-01
## 9 Diet binomial (Intercept) -3.178054e+00 0.2830692544 3.000853e-29
## 10 Diet binomial Dietlad2 2.797769e-01 0.3610705226 4.384260e-01
## 11 Diet binomial Dietrm1 1.760910e-01 0.4597166333 7.016883e-01
## 12 Diet quasibinomial (Intercept) -3.178054e+00 0.2978668571 4.162587e-08
## 13 Diet quasibinomial Dietlad2 2.797769e-01 0.3799456850 4.736704e-01
## 14 Diet quasibinomial Dietrm1 1.760910e-01 0.4837485760 7.212867e-01
## 15 No./cage binomial (Intercept) -2.947182e+00 0.5250782715 1.990286e-08
## 16 No./cage binomial `No./cage` -1.464854e-02 0.1130317481 8.968855e-01
## 17 No./cage quasibinomial (Intercept) -2.947182e+00 0.5382237104 6.389405e-05
## 18 No./cage quasibinomial `No./cage` -1.464854e-02 0.1158615203 9.010696e-01
## 19 Duration binomial (Intercept) -1.119484e+01 3.5331305839 1.532100e-03
## 20 Duration binomial Duration 8.091505e-02 0.0345598695 1.921668e-02
## 21 Duration quasibinomial (Intercept) -1.119484e+01 2.7531684946 1.013480e-03
## 22 Duration quasibinomial Duration 8.091505e-02 0.0269305483 8.889402e-03
## 23 Date binomial (Intercept) -5.192254e+00 2.7509766705 5.910360e-02
## 24 Date binomial Date 2.765304e-04 0.0003472537 4.258371e-01
## 25 Date quasibinomial (Intercept) -5.192254e+00 2.8096055137 8.441018e-02
## 26 Date quasibinomial Date 2.765304e-04 0.0003546544 4.476857e-01
Step 4: Modelling HCD variability
In Step 4, the selected HCD are used to quantify the variability in control
incidence across studies. For a dichotomous endpoint, the EFSA case study uses a
beta-binomial model. This model separates the observed variability into:
- within-study binomial variability, due to the finite number of animals
examined in each control group;
- between-study variability, reflecting differences in the underlying control
incidence between historical studies.
Let \(k_i\) be the number of animals with the finding in historical control
group \(i\), and \(n_i\) the number of animals examined. Conditional on the true
control incidence \(p_i\), the outcome is modelled as:
\[
k_i \mid p_i \sim \text{Binomial}(n_i, p_i)
\]
The beta-binomial model assumes that the true study-specific incidences vary
between studies according to a beta distribution:
\[
p_i \sim \text{Beta}(\alpha, \beta)
\]
Together, these assumptions imply a beta-binomial distribution for the observed
counts \(k_i\). The model therefore allows for extra variability beyond a simple
binomial model.
The package function model_binary_hcd_variability() fits this beta-binomial
model by maximum likelihood. It estimates the average historical control
incidence and the between-study variability parameter. The model is
parameterised using:
\[
\pi = \frac{\alpha}{\alpha + \beta}
\]
where \(\pi\) is the mean population incidence, and
\[
\rho = \alpha + \beta
\]
where \(\rho\) is interpreted as an effective sample size. Larger values of
\(\rho\) imply less between-study variability, because the beta distribution is
more concentrated around its mean. Smaller values imply more between-study
variability.
The package also reports the variability parameter
\[
\phi = \frac{1}{\rho + 1}
\]
so that the effective sample size can be written as:
\[
\rho = \frac{1}{\phi} - 1
\]
Following the EFSA case-study approach, uncertainty in the estimated
between-study variability is accounted for using a profile likelihood confidence
limit.
The beta-binomial model estimates a very small amount of between-study
variability in the selected HCD set. The maximum likelihood estimate (MLE) of
the variability parameter \(\phi\) is close to zero, corresponding to a large
effective sample size. However, following the EFSA case-study approach,
uncertainty in the variability estimate is incorporated using a profile
likelihood confidence interval for the between-study variability parameter.
The upper confidence limit for \(\phi\) corresponds to a lower confidence limit
for the effective sample size of 64.78 animals. In the EFSA Annex A report, this
value is reported as 64 animals. In the present analysis, we keep the exact
value returned by model_binary_hcd_variability() and pass the Step 4 object
directly to the subsequent steps.
The formatted summary table is shown below. The underlying raw results can be
accessed using binary_step4$summary_table.
## Step 3b selection: retain 104-week studies for Step 4
hcd_step4 <- hcd_5y[hcd_5y$Duration == 104, ]
nrow(hcd_step4)
## [1] 12
binary_step4 <- model_binary_hcd_variability(
dat = hcd_step4,
k_col = "k",
n_col = "n"
)
binary_step4$formatted_summary_table
## Number of rows used Alpha Phi MLE Phi lower limit Phi upper limit
## 12.000 0.050 0.000 0.000 0.015
## Effective sample size MLE Effective sample size lower limit
## 15126.015 64.782
## Effective sample size upper limit
## Inf
The profile likelihood plot shows the support in the HCD for different values of
the between-study variability parameter \(\phi\). The x-axis is \(\phi\), and
the y-axis is the relative profile log-likelihood, calculated relative to the
maximum profile log-likelihood.
The dashed vertical line marks the maximum likelihood estimate of \(\phi\). In
this example, the MLE is very close to zero, indicating little estimated
between-study variability.
The red horizontal line is the likelihood-ratio cutoff used to obtain the
profile likelihood confidence interval. The upper confidence limit for \(\phi\)
is obtained from the point where the profile log-likelihood curve crosses this
red line on the right-hand side. This upper confidence limit is then converted
to the conservative lower confidence limit for the effective sample size.
binary_step4$profile_plot
The distribution plot shows the fitted beta distribution for the underlying
control incidence \(p\) across historical studies. This distribution represents
the between-study variability estimated from the selected HCD.
binary_step4$distribution_plot
The fitted distribution indicates the estimated between-study variability in the
underlying control incidence across the selected HCD. The grey curve represents
the HCD-based distribution obtained using a very large effective sample size,
corresponding to negligible between-study variability. The black curve
represents the more conservative HCD-based distribution using the lower
confidence limit for the effective sample size, and this black curve is carried
forward to the subsequent steps.
Step 5: Comparison of the concurrent control with HCD
In Step 5, the concurrent control (CC) is compared with the distribution
expected from the historical control data. This step checks whether the CC is
compatible with the HCD before using HCD information in the analysis of the
index study.
The EFSA framework describes two possible approaches for this comparison. The
chosen approach should be specified in the analysis protocol before carrying out
the analysis.
Approach 1: Comparison on the incidence scale
In Approach 1, the comparison is made on the scale of the underlying control
incidence. The HCD model from Step 4 is used to derive a prediction interval for
the underlying control incidence expected in a new comparable study. This
interval is then compared with the confidence interval for the CC incidence.
For the CC group, the observed incidence is:
\[
\hat{p}_{CC} = \frac{k_{CC}}{n_{CC}},
\]
where \(k_{CC}\) is the number of animals with the finding and \(n_{CC}\) is the
number examined.
This approach asks whether the CC incidence, allowing for its sampling
uncertainty, is compatible with the HCD-based prediction interval for the
underlying control incidence.
Approach 2: Comparison on the observed count scale
In Approach 2, the comparison is made on the scale of the observed number of
animals with the finding in the CC group. The Step 4 beta-binomial model is used
to derive the predictive distribution for a new control group with the same
number of examined animals as the CC.
If the CC has \(n_{CC}\) examined animals, the predictive distribution is:
\[
K_{new} \sim \text{Beta-binomial}(n_{CC}, \alpha, \beta),
\]
where \(\alpha\) and \(\beta\) are derived from the fitted HCD model. A
prediction interval is then obtained from this predictive distribution and
compared with the observed number of cases in the concurrent control group,
\(k_{CC}\).
This approach directly asks whether the observed number of cases in the
concurrent control group is plausible under the HCD model.
For the binary compatibility decision, the package also allows a small
rounding tolerance through comparison_digits (3 decimal places by default).
The observed CC incidence is rounded towards the HCD predictive interval,
while the HCD predictive interval expressed on the incidence scale is rounded
outwards: the
lower limit is rounded down and the upper limit is rounded up. This can only
make the predictive interval wider for the decision. The original count and
incidence values are retained unchanged in the reported output. The comparison
uses inclusive >= and <= boundaries.
Approach used in this case study
In the EFSA binary case study, Approach 2 is used. This choice is appropriate
for reproducing the case study because the endpoint is a rare dichotomous
finding and the report compares the observed number of animals with the finding
in the CC group with a beta-binomial prediction interval.
Therefore, in this vignette, Step 5 is implemented using the case-based
comparison: the observed number of cases in the concurrent control group is
compared with the HCD-based predictive distribution for a new control group of
the same size.
The concurrent control group had 3 animals with the finding among 47 examined
animals, corresponding to a sample incidence of approximately 6.4%. Using the
Step 4 distribution based on the upper confidence limit for
\(\phi\), the beta-binomial predictive interval for the number of affected
animals in a new control group of size 47 ranged from 0 to 7 animals, with
coverage probability 96.3%.
The observed concurrent control outcome lies well within this prediction
interval and is therefore considered fully consistent with the HCD. This
conclusion agrees with the EFSA case study, and the analysis proceeds to Step 6.
## Step 5: compare concurrent control with HCD
##
## Approach 2 compares the observed number of CC cases with the beta-binomial
## predictive distribution for a new control group with the same number
## of examined animals as the CC.
cc_binary <- index_study[index_study$Dose == 0, ]
binary_step5 <- compare_binary_cc_to_hcd(
cc_dat = cc_binary,
step4_result = binary_step4,
k_col = "k",
n_col = "n",
approach = "new_group_outcome",
prediction_basis = "upper_limit",
pi_conf_level = 0.95,
comparison_digits = 3
)
binary_step5$formatted_summary_table
## Group Approach CC number with finding CC number examined CC incidence
## <NA> new_group_outcome 3.000 47.000 0.064
## Prediction basis Effective sample size used HCD incidence used for prediction
## upper_limit 64.782 0.063
## HCD PI lower number of cases HCD PI upper number of cases
## 0.000 7.000
## Left-tail probability Right-tail probability Comparison rounding digits
## 0.000 0.037 3.000
## CC within HCD prediction interval before rounding
## TRUE
## CC within HCD prediction interval Compatibility changed after rounding
## TRUE FALSE
## Compatibility assessment
## within before rounding
The observed number of CC cases lies within the HCD-based prediction interval.
This indicates that the concurrent control outcome is compatible with the
selected HCD under the Step 5 comparison. Therefore, in this case study, the
workflow proceeds to Step 6.
The Step 5 plot shows the beta-binomial predictive distribution for the number
of animals with the finding in a new control group of 47 examined animals. Each
point gives the predicted probability of observing a given count. The horizontal
black bar shows the HCD prediction interval, from 0 to 7 animals with the
finding. The solid vertical line shows the observed concurrent control count,
which is 3 animals. The dashed vertical line marks the first count above the
prediction interval, here 8 animals. Since the observed concurrent control count
lies inside the prediction interval, no inconsistency between the concurrent
control and the HCD is identified.
binary_step5$plot$new_group_outcome
Illustration of Approach 1
For illustration purposes only, the comparison can also be performed using
Approach 1, where the confidence interval for the concurrent control incidence
is compared with the HCD prediction interval for the underlying population
incidence.
According to the EFSA framework, the analysis protocol should prespecify which
Step 5 approach will be used. The two approaches are alternative strategies and
should not both be used to choose the most favourable conclusion after
inspecting the results.
In this vignette, Approach 2 is the primary analysis because it is the approach
used in the EFSA binary case study. The following code simply illustrates how
the package can also implement Approach 1.
In the main analysis, we use approach = "new_group_outcome" because this is
the approach prespecified in the analysis of this binary case study. Here, we
additionally illustrate the alternative option
approach = "population_incidence".
This command compares the 95% confidence interval for the CC incidence with the
HCD prediction interval for the underlying population incidence.
The CC incidence is 3/47, or approximately 6.4%. Its 95% confidence interval is
2.2% to 17.2%. The HCD prediction interval is 1.8% to 13.3%.
For the compatibility decision, the package uses directional rounding towards
compatibility. With comparison_digits = 3 by default, the HCD prediction
interval is widened by rounding its lower limit down and its upper limit up,
while the CC confidence interval is narrowed by rounding its lower limit up and
its upper limit down. The inclusive checks >= and <= are then applied. The
raw interval values are retained unchanged in the returned results; the
directional rounding is used only for the compatibility check.
Because the upper part of the CC confidence interval remains outside the HCD
prediction interval even after this three-decimal directional rounding, the
result remains FALSE in this example.
binary_step5_a1 <- compare_binary_cc_to_hcd(
cc_dat = cc_binary,
step4_result = binary_step4,
k_col = "k",
n_col = "n",
approach = "population_incidence",
prediction_basis = "upper_limit",
pi_conf_level = 0.95,
cc_conf_level = 0.95,
comparison_digits = 3
)
binary_step5_a1$formatted_summary_table
## Group Approach CC number with finding CC number examined
## <NA> population_incidence 3.000 47.000
## CC incidence Prediction basis Effective sample size used
## 0.064 upper_limit 64.782
## HCD incidence used for prediction CC CI lower limit CC CI upper limit
## 0.063 0.022 0.172
## HCD PI lower limit HCD PI upper limit Comparison rounding digits
## 0.018 0.133 3.000
## CC within HCD prediction interval before rounding
## FALSE
## CC within HCD prediction interval Compatibility changed after rounding
## FALSE FALSE
## Compatibility assessment
## outside after directional rounding to 3 decimal places
The Approach 1 plot shows the fitted HCD distribution for the underlying
population incidence. The black horizontal bar is the HCD prediction interval.
The grey horizontal bar is the 95% confidence interval for the concurrent
control incidence. Because the grey interval is wider and is not fully contained
within the HCD prediction interval, Approach 1 would flag the CC as not fully
compatible with the HCD.
binary_step5_a1$plots$population_incidence
Step 5b: Investigation of possible inconsistency
Step 5b is relevant when the Step 5a comparison suggests that the concurrent
control is not fully consistent with the HCD. In the main analysis of this case
study, Approach 2 was used and no inconsistency was found. However, Approach 1
flags a possible mismatch because the CC confidence interval is not fully
contained within the HCD prediction interval.
For illustration, we use investigate_binary_cc_hcd_mismatch() to summarise
this possible mismatch and to explore simple sensitivity checks. For Approach 1,
the function compares the gap between the CC confidence interval and the HCD
prediction interval, and can repeat the comparison under alternative confidence
levels, interval methods, or Step 4 prediction bases.
The object binary_step5b_a1$formatted_main_diagnostics gives the main Step 5b
diagnostic summary for the selected settings. It compares the CC confidence
interval with the HCD prediction interval and reports the interval overlap and
gap using the original unrounded values. It also reports whether the interval
was within the HCD prediction interval before rounding, the final compatibility
result after applying the directional rounding controlled by
comparison_digits, and whether rounding changed the decision. The descriptive status is reported as within or outside; it is
not intended to replace the scientific judgement required in Step 5b.
In this case, the CC confidence interval is not fully contained within the HCD
prediction interval. The upper part of the CC interval extends beyond the HCD
interval by about 3.8 percentage points. The descriptive status is therefore
outside. This indicates that further investigation and scientific judgement
are needed before deciding how to proceed in the EFSA workflow.
The formatted diagnostic table is shown below. The raw results can be accessed
using binary_step5b_a1$main_diagnostics.
The object binary_step5b_a1$formatted_sensitivity_table repeats the same
comparison under alternative settings, including different CC confidence
interval methods, confidence levels, HCD prediction interval levels, and Step 4
prediction bases.
In this case, all sensitivity settings still show that the CC interval is not
within the HCD prediction interval, and the descriptive status remains
outside. This suggests that, under Approach 1, the apparent mismatch is not
driven by a single interval method or confidence-level choice. The sensitivity
results are diagnostic and do not themselves determine the scientific action to
be taken in Step 5b.
The formatted sensitivity table is shown below. The raw results can be accessed
using binary_step5b_a1$sensitivity_table.
binary_step5b_a1 <- investigate_binary_cc_hcd_mismatch(
cc_dat = cc_binary,
step4_result = binary_step4,
k_col = "k",
n_col = "n",
approach = "population_incidence",
prediction_basis_main = "upper_limit",
cc_conf_level_main = 0.95,
pi_conf_level_main = 0.95,
comparison_digits = 3
)
binary_step5b_a1$formatted_main_diagnostics
## Group Approach Prediction basis CC confidence level
## <NA> population_incidence upper_limit 0.950
## HCD PI confidence level CC CI method Comparison rounding digits
## 0.950 wilson 3.000
## CC CI lower limit CC CI upper limit HCD PI lower limit HCD PI upper limit
## 0.022 0.172 0.018 0.133
## Lower gap Upper gap Maximum gap CC interval width HCD interval width
## 0.000 0.038 0.038 0.150 0.115
## Overlap width Overlap proportion of CC interval Within HCD PI before rounding
## 0.111 0.743 FALSE
## Within HCD PI Compatibility changed after rounding
## FALSE FALSE
## Compatibility assessment Status
## outside after directional rounding to 3 decimal places outside
binary_step5b_a1$formatted_sensitivity_table
## Group Approach Prediction basis CC CI method
## <NA> population_incidence lower_limit wilson
## <NA> population_incidence lower_limit wilson
## <NA> population_incidence lower_limit wilson
## <NA> population_incidence lower_limit wilson
## <NA> population_incidence lower_limit clopper-pearson
## <NA> population_incidence lower_limit clopper-pearson
## <NA> population_incidence lower_limit clopper-pearson
## <NA> population_incidence lower_limit clopper-pearson
## <NA> population_incidence lower_limit agresti-coull
## <NA> population_incidence lower_limit agresti-coull
## <NA> population_incidence lower_limit agresti-coull
## <NA> population_incidence lower_limit agresti-coull
## <NA> population_incidence point_estimate wilson
## <NA> population_incidence point_estimate wilson
## <NA> population_incidence point_estimate wilson
## <NA> population_incidence point_estimate wilson
## <NA> population_incidence point_estimate clopper-pearson
## <NA> population_incidence point_estimate clopper-pearson
## <NA> population_incidence point_estimate clopper-pearson
## <NA> population_incidence point_estimate clopper-pearson
## <NA> population_incidence point_estimate agresti-coull
## <NA> population_incidence point_estimate agresti-coull
## <NA> population_incidence point_estimate agresti-coull
## <NA> population_incidence point_estimate agresti-coull
## <NA> population_incidence upper_limit wilson
## <NA> population_incidence upper_limit wilson
## <NA> population_incidence upper_limit wilson
## <NA> population_incidence upper_limit wilson
## <NA> population_incidence upper_limit clopper-pearson
## <NA> population_incidence upper_limit clopper-pearson
## <NA> population_incidence upper_limit clopper-pearson
## <NA> population_incidence upper_limit clopper-pearson
## <NA> population_incidence upper_limit agresti-coull
## <NA> population_incidence upper_limit agresti-coull
## <NA> population_incidence upper_limit agresti-coull
## <NA> population_incidence upper_limit agresti-coull
## CC confidence level HCD PI confidence level Comparison rounding digits
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## 0.900 0.900 3.000
## 0.900 0.950 3.000
## 0.950 0.900 3.000
## 0.950 0.950 3.000
## CC CI lower limit CC CI upper limit HCD PI lower limit HCD PI upper limit
## 0.026 0.149 0.058 0.059
## 0.026 0.149 0.058 0.059
## 0.022 0.172 0.058 0.059
## 0.022 0.172 0.058 0.059
## 0.018 0.157 0.058 0.059
## 0.018 0.157 0.058 0.059
## 0.013 0.175 0.058 0.059
## 0.013 0.175 0.058 0.059
## 0.022 0.154 0.058 0.059
## 0.022 0.154 0.058 0.059
## 0.016 0.178 0.058 0.059
## 0.016 0.178 0.058 0.059
## 0.026 0.149 0.055 0.061
## 0.026 0.149 0.055 0.062
## 0.022 0.172 0.055 0.061
## 0.022 0.172 0.055 0.062
## 0.018 0.157 0.055 0.061
## 0.018 0.157 0.055 0.062
## 0.013 0.175 0.055 0.061
## 0.013 0.175 0.055 0.062
## 0.022 0.154 0.055 0.061
## 0.022 0.154 0.055 0.062
## 0.016 0.178 0.055 0.061
## 0.016 0.178 0.055 0.062
## 0.026 0.149 0.022 0.119
## 0.026 0.149 0.018 0.133
## 0.022 0.172 0.022 0.119
## 0.022 0.172 0.018 0.133
## 0.018 0.157 0.022 0.119
## 0.018 0.157 0.018 0.133
## 0.013 0.175 0.022 0.119
## 0.013 0.175 0.018 0.133
## 0.022 0.154 0.022 0.119
## 0.022 0.154 0.018 0.133
## 0.016 0.178 0.022 0.119
## 0.016 0.178 0.018 0.133
## Within HCD PI before rounding Within HCD PI
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## FALSE FALSE
## Compatibility changed after rounding
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## FALSE
## Compatibility assessment Status
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
## outside after directional rounding to 3 decimal places outside
Step 6: Analysis of the index study with and without HCD
In Step 6, the index study is analysed both without HCD and with HCD
incorporated into the control group. The goal is to evaluate whether inclusion
of HCD changes the statistical conclusion regarding a possible dose-related
increase in the finding.
The package function analyze_binary_index_with_hcd() implements this step for
grouped binary data.
Incorporating HCD using pseudo-data
The EFSA framework incorporates HCD by converting the fitted Step 4 beta
distribution into pseudo-data for the control response. Let
\[
p \sim \text{Beta}(\alpha,\beta)
\]
be the fitted historical control distribution from Step 4. The effective sample
size is:
\[
n_{\mathrm{eff}} = \alpha + \beta
\]
and the mean historical incidence is:
\[
\pi = \frac{\alpha}{\alpha+\beta}.
\]
Pseudo-data are then constructed as:
\[
k^{*} \approx \alpha,
\qquad
n^{*} \approx \alpha + \beta,
\]
where \(k^{*}\) is the pseudo-number of animals with the finding and \(n^{*}\)
is the pseudo-number examined.
In this case study, the Step 4 distribution based on the upper
confidence limit for \(\phi\) is used. This gives an effective sample size of
approximately 64 animals and an HCD incidence of approximately 6.3%,
corresponding to pseudo-data of:
\[
4/64.
\]
These pseudo-data are then pooled with the concurrent control group:
\[
k_{\mathrm{pooled}}
=
k_{\mathrm{CC}} + k^{*},
\]
\[
n_{\mathrm{pooled}}
=
n_{\mathrm{CC}} + n^{*}.
\]
Statistical tests
This section shows the statistical tests used to analyse the index study. The
tests are presented both without and with HCD to illustrate how the
corresponding package functions are applied and how the outputs are structured.
Important: This comparison is shown only for illustrative purposes. In an
actual application of the EFSA workflow, the decision to analyse the index study
with or without HCD should be driven by the conclusions from the previous
workflow steps, especially the assessment of compatibility between the
concurrent control and the HCD. Therefore, only the analysis justified by the
preceding steps would normally be carried out.
Statistical tests without HCD
The function first performs analyses without HCD using only the concurrent
control and dose groups from the index study.
The following tests are performed:
Cochran-Armitage (CA) trend test: evaluates dose-related trend across dose
groups;
Exact CA trend test: exact version of the CA trend test;
Fisher’s exact test: compares the highest dose group with the concurrent
control.
Statistical tests with HCD
The function then repeats the analysis after incorporating the HCD pseudo-data
into the control group.
The following tests are performed:
Tarone asymptotic trend test: an extension of the CA trend test incorporating
HCD pseudo-data;
HYT (Hoel-Yanagawa exact test): exact trend test incorporating HCD
pseudo-data;
Exact highest-dose versus control comparison with HCD.
Main outputs of analyze_binary_index_with_hcd()
The function returns several outputs, including:
formatted_summary_table_without_hcd: display-ready summary of the analyses
based only on the index study;
formatted_summary_table_with_hcd: display-ready summary of the analyses
incorporating HCD;
summary_table_without_hcd and summary_table_with_hcd: raw numeric versions
of the two summary tables;
pseudo_data: details of the constructed HCD pseudo-data and pooled control
data;
tests_without_hcd: results of analyses using only the index study;
tests_with_hcd: results of analyses incorporating HCD pseudo-data;
plot: graphical summary of the concurrent control, HCD pseudo-data, pooled
control, and dose groups.
The Step 6 results are returned in two separate summary tables: one for the
analysis without HCD and one for the analysis with HCD. The with-HCD table also
reports the HCD basis used for the analysis. Using the Step 4 distribution
(upper_limit basis), the estimated HCD incidence is approximately 6.3%, with
an effective sample size of 64.8 animals. For the grouped binary analysis, this
corresponds to rounded pseudo-data of 4 affected animals out of 64
pseudo-animals.
The concurrent control group contains 3 affected animals among 47 examined
animals. After incorporation of the HCD pseudo-data, the pooled control consists
of 7 affected animals among 111 animals.
Without HCD, the exact Cochran-Armitage trend test gives a p-value of
approximately 0.104, and Fisher’s exact test comparing the highest dose group
with the concurrent control gives a p-value of approximately 0.113. Therefore,
without HCD, the analysis is not statistically significant at the 5% level.
After incorporation of HCD, the HYT exact trend test gives a p-value of
approximately 0.036, and the exact highest-dose versus control comparison gives
a p-value of approximately 0.031. Therefore, when HCD are incorporated according
to the EFSA framework, the analysis becomes statistically significant at the 5%
level.
binary_step6 <- analyze_binary_index_with_hcd(
index_dat = index_study,
step4_result = binary_step4,
dose_col = "Dose",
k_col = "k",
n_col = "n",
control_dose = 0,
hcd_basis = "upper_limit",
ci_method = "wilson",
ci_conf_level = 0.95,
make_plot = TRUE,
add_pooled_to_plot = TRUE
)
binary_step6$formatted_summary_table_without_hcd
## Group Control dose Highest dose Number of CC cases CC number examined
## <NA> 0.000 10000.000 3.000 47.000
## CA asymptotic p-value CA exact p-value Fisher control vs max p-value
## 0.1073 0.1042 0.1128
## CA exact control vs max p-value
## 0.1128
binary_step6$formatted_summary_table_with_hcd
## Group HCD basis Control dose Highest dose Tarone asymptotic p-value
## <NA> upper_limit 0.000 10000.000 0.0543
## HYT exact p-value HYT exact control vs max p-value
## 0.0358 0.0306
These results show how the index-study analysis is carried out after the
concurrent control has been judged compatible with the selected HCD in Step 5.
In this case study, the HCD information contributes additional control
information to the trend analysis, while the evidence for a treatment-related
increase is still driven by the dose groups in the index study.
Step 6 returns two source-comparison versions of the descriptive plot. Both
views retain all active-dose groups from the index study. The first shows the
index study alone, including its concurrent control. The second shows the same
index-study dose groups together with the HCD pseudo-control and the pooled
index/HCD control. The all-sources plot is the same view that was returned as
binary_step6$plot in earlier versions; that object is retained as a
backward-compatible alias.
binary_step6$plots$index_study_only
binary_step6$plots$all_sources
Step 7: Sensitivity analyses and alternative scenarios
The EFSA framework treats Step 7 as a case-specific sensitivity-analysis stage
rather than as a single fixed procedure. The goal is to evaluate whether the
main conclusions remain stable under alternative scientifically reasonable
assumptions.
The efsaHCD package supports many such sensitivity analyses directly. Examples
include:
- changing the Step 4 prediction basis (lower_limit, point_estimate,
upper_limit);
- varying confidence interval methods and confidence levels;
- exploring Step 5 approaches in different settings;
- comparing analyses with and without HCD.
However, the precise sensitivity analyses to perform depend strongly on the
endpoint, study design, available HCD, and scientific context. Therefore, Step 7
is intentionally flexible and should be tailored case by case rather than
treated as a fully standardised workflow step.
Important: When implementing the EFSA HCD workflow, the key decisions,
assumptions, and methods to be explored in the sensitivity analysis should be
prespecified in the analysis protocol as far as possible.
Continuous endpoint case study
This section illustrates the use of the efsaHCD package for the EFSA HCD
workflow using the continuous endpoint case study from Annex B. The endpoint is
a clinical chemistry measurement, expressed as a continuous outcome. In a real
regulatory assessment, a similar endpoint could be a liver enzyme level, organ
weight, or another quantitative toxicological measurement.
The example uses aggregated historical control data and aggregated index study
data. The historical control database contains 63 historical control studies,
with study sample sizes ranging from 5 to 12 animals. The index study contains
four groups: one concurrent control group and three treated groups receiving
600, 1400, and 5000 ppm of the test item. The data are artificial, but they are
designed to mimic the type of information that may be available in a real
assessment.
The aim of the case study is to assess whether the historical control data are
suitable for use, to model the variability among historical control studies, to
compare the concurrent control with the historical control distribution, and, if
appropriate, to include the historical control information in the dose-response
analysis of the index study.
Data preparation for reproducing the Annex example
The data preparation in this vignette follows the logic of the EFSA Annex B R
script. This step is included to reproduce the published case study as closely
as possible. It should not be interpreted as a general package function or as a
recommended preprocessing step for all analyses.
In the Annex example, the historical control data are provided as aggregated
study-level means and sample sizes. The EFSA script does not read standard
errors directly from the Excel file. Instead, it reconstructs artificial
variability values for the example. More specifically, standard deviations are
generated from a truncated normal distribution, with mean equal to 20% of the
study mean, standard deviation equal to 0.15, and truncation at zero. The
standard error is then computed from the generated standard deviation and the
study sample size.
This reconstruction is only needed because the Annex example is based on
artificial data. In the efsaHCD package workflow, the analysis functions
assume that the required summary statistics, such as the sample mean, standard
error, and sample size, are already available in the analysis dataset.
We generated the prepared dataset using the original EFSA Annex B code provided
for the continuous endpoint case study and included the resulting analysis-ready
dataset in the package for illustration and reproducibility purposes.
cont_file <- system.file(
"extdata",
"appendix_cont_Dataset_prepared.csv",
package = "efsaHCD"
)
cont <- read.csv(cont_file)
hcd_cont <- cont[cont$HCD == 1, ]
hcd_cont_legal <- cont[cont$HCD == 1 & cont$year_2 > 2012, ]
hcd_cc_cont <- cont[cont$HCD < 3, ]
index_cont <- cont[cont$HCD > 1, ]
cc_cont <- cont[cont$HCD == 2, ]
Steps 1 and 2: Planning and data acceptability
The aim of this case study is to assess whether the historical control data are
suitable for use, to model the variability among historical control studies, to
compare the concurrent control with the historical control distribution, and, if
appropriate, to include the historical control information in the dose-response
analysis of the index study.
In Step 1, the intended use of HCD is planned according to the EFSA decision
scheme. This includes defining the endpoint of interest, identifying the
relevant HCD studies, specifying comparability criteria, identifying covariates
to be investigated, planning the statistical analyses for the subsequent steps,
and prespecifying important sensitivity analyses and decision criteria in the
protocol as far as possible.
In Step 2, the available HCD are evaluated with respect to biological,
experimental, legal, and statistical acceptability before proceeding to the
formal statistical workflow. The HCD are initially considered as provided.
However, criteria for inclusion or exclusion of studies, as well as other key
decisions such as which aspects to investigate in sensitivity analyses, should
be decided upfront and stated in the protocol as far as possible.
Step 3: Preliminary investigation of the historical control data
The aim of Step 3 is to explore the historical control data before fitting the
final meta-analytical model. This includes assessing the plausibility of the
normality assumption, identifying potentially influential or outlying studies,
and investigating possible associations between the endpoint and available
study-level covariates.
Following the EFSA Annex example, these exploratory assessments are performed on
both the original scale and the log-transformed scale. The purpose is not to
apply automatic exclusion rules, but rather to identify observations or study
characteristics that may require further investigation before proceeding to the
modelling stage.
Step 3a: Normality diagnostics
Before modelling the historical control data, we first inspect whether the data
are reasonably compatible with the normality assumption used in the
meta-analytical model.
The function assess_continuous_summary_hcd_normality() provides two options
through the diagnostic_method argument. With
diagnostic_method = "raw_values", normality is assessed directly on the
observed study-level means. This is a simple descriptive check. With
diagnostic_method = "model_residuals", a random-effects meta-analysis model is
first fitted, and normality is assessed using internally studentized marginal
residuals. This second option follows more closely the approach used in the EFSA
Annex example.
The assessment is performed on both the original scale and the log-transformed
scale. The log scale is included as a sensitivity check, because transformation
may sometimes improve the agreement with the normality assumption.
The diagnostic plots and summary table are used as screening tools only. They
are not automatic exclusion rules. Any decision to exclude studies should also
consider biological and study-design information.
normality_step3 <- assess_continuous_summary_hcd_normality(
data = hcd_cont,
response_col = "bioch_means",
log_response_col = "log_bioch_means",
se_col = "se",
log_se_col = "log_se",
study_id_col = "StudyID",
endpoint_label = "Clinical chemistry endpoint",
use_log = TRUE,
diagnostic_method = "model_residuals",
ks_method = "standardized"
)
normality_step3$formatted_summary_table
## Scale Diagnostic method Diagnostic quantity N Mean
## raw model_residuals Internally studentized marginal residual 63.000 0.067
## log model_residuals Internally studentized marginal residual 63.000 -0.024
## SD Median Minimum Maximum Skewness Shapiro-Wilk W Shapiro-Wilk p-value
## 1.033 0.117 -2.040 2.522 0.252 0.980 0.4015
## 1.006 0.002 -2.159 2.305 0.255 0.981 0.4376
## Kolmogorov-Smirnov D Kolmogorov-Smirnov p-value No. possible outliers
## 0.072 0.8987 4.000
## 0.065 0.9557 4.000
## Possible outlier IDs
## 6, 19, 23, 26
## 19, 23, 26, 30
## Notes
## Ties detected; Kolmogorov-Smirnov p-value should be interpreted cautiously. KS warning: ties should not be present for the one-sample Kolmogorov-Smirnov test
## Ties detected; Kolmogorov-Smirnov p-value should be interpreted cautiously. KS warning: ties should not be present for the one-sample Kolmogorov-Smirnov test
resid_box_step3 <- plot_continuous_hcd_residual_boxplot(
data = hcd_cont,
response_col = "bioch_means",
log_response_col = "log_bioch_means",
se_col = "se",
log_se_col = "log_se",
study_id_col = "StudyID",
endpoint_label = "Clinical chemistry endpoint",
use_log = TRUE,
residual_type = "rstudent",
model_method = "REML",
model_test = "adhoc"
)
resid_box_step3$plot
resid_box_step3$formatted_outlier_table
## Scale Scale Study ID Residual IQR outlier Lower fence
## raw original scale 23.000 2.543 TRUE -2.347
## raw original scale 26.000 2.543 TRUE -2.347
## log log-transformed scale 23.000 2.401 TRUE -2.544
## log log-transformed scale 26.000 2.332 TRUE -2.544
## Upper fence
## 2.339
## 2.339
## 2.285
## 2.285
Step 3b: Influence diagnostics
After assessing normality, we investigate whether some historical control
studies are unusually influential in the meta-analytical model. Influence
diagnostics help identify studies that strongly affect the estimated mean
response, heterogeneity, or model fit.
The function assess_continuous_hcd_influence() fits a random-effects
meta-analysis model and computes several standard influence diagnostics,
including externally studentized residuals, DFFITS values, Cook’s distances,
covariance ratios, hat values, and leave-one-out heterogeneity estimates.
These diagnostics are intended as exploratory tools and should not be
interpreted as automatic exclusion criteria.
infl_step3 <- assess_continuous_hcd_influence(
dat = hcd_cont,
mean_col = "bioch_means",
se_col = "se",
study_id_col = "StudyID",
scale = "orig",
method = "REML",
test = "adhoc",
make_plot = TRUE
)
infl_step3$formatted_diagnostics
## study_id order rstudent dffits cook.d cov.r tau2.del QE.del hat weight
## 1 1 1 -1.086 -0.148 0.022 1.016 0.018 105.673 0.018 1.820
## 2 2 2 -1.726 -0.208 0.042 0.986 0.016 102.230 0.019 1.947
## 3 3 3 0.732 0.085 0.007 1.023 0.018 106.508 0.014 1.447
## 4 4 4 -0.926 -0.131 0.017 1.021 0.018 106.259 0.018 1.793
## 5 5 5 0.202 0.010 0.000 1.033 0.018 107.475 0.016 1.557
## 6 6 6 2.003 0.266 0.068 0.962 0.016 101.375 0.012 1.187
## 7 7 7 -0.665 -0.101 0.010 1.028 0.018 106.987 0.017 1.740
## 8 8 8 0.689 0.079 0.006 1.024 0.018 106.615 0.015 1.458
## 9 9 9 0.247 0.017 0.000 1.032 0.018 107.427 0.015 1.546
## 10 10 10 1.193 0.152 0.023 1.006 0.017 105.038 0.014 1.356
## 11 11 11 1.390 0.180 0.032 0.997 0.017 104.267 0.013 1.309
## 12 12 12 0.064 -0.008 0.000 1.034 0.018 107.581 0.016 1.592
## 13 13 13 0.689 0.079 0.006 1.024 0.018 106.615 0.015 1.458
## 14 14 14 0.816 0.097 0.009 1.020 0.018 106.282 0.014 1.426
## 15 15 15 -1.256 -0.166 0.027 1.009 0.017 104.923 0.019 1.861
## 16 16 16 -0.267 -0.052 0.003 1.034 0.018 107.563 0.017 1.652
## 17 17 17 0.517 0.055 0.003 1.028 0.018 107.000 0.015 1.501
## 18 18 18 0.427 0.042 0.002 1.030 0.018 107.168 0.015 1.512
## 19 19 19 -2.048 -0.231 0.050 0.966 0.015 99.868 0.020 2.007
## 20 20 20 0.471 0.048 0.002 1.029 0.018 107.090 0.015 1.501
## 21 21 21 1.894 0.251 0.061 0.969 0.016 101.929 0.012 1.212
## 22 22 22 -0.464 -0.077 0.006 1.032 0.018 107.356 0.017 1.702
## 23 23 23 2.543 0.337 0.106 0.924 0.015 98.333 0.011 1.086
## 24 24 24 -0.219 -0.046 0.002 1.034 0.019 107.590 0.017 1.652
## 25 25 25 1.193 0.152 0.023 1.006 0.017 105.038 0.014 1.356
## 26 26 26 2.543 0.337 0.106 0.924 0.015 98.333 0.011 1.086
## 27 27 27 0.247 0.017 0.000 1.032 0.018 107.427 0.015 1.546
## 28 28 28 -0.219 -0.046 0.002 1.034 0.019 107.590 0.017 1.652
## 29 29 29 -1.142 -0.154 0.024 1.014 0.018 105.440 0.018 1.834
## 30 30 30 -1.921 -0.223 0.047 0.975 0.015 100.819 0.020 1.992
## 31 31 31 1.820 0.241 0.057 0.973 0.016 102.293 0.012 1.229
## 32 32 32 -0.365 -0.065 0.004 1.033 0.018 107.479 0.017 1.676
## 33 33 33 -0.365 -0.065 0.004 1.033 0.018 107.479 0.017 1.676
## 34 34 34 -0.819 -0.119 0.014 1.024 0.018 106.594 0.018 1.766
## 35 35 35 0.293 0.023 0.001 1.032 0.018 107.369 0.015 1.546
## 36 36 36 0.202 0.010 0.000 1.033 0.018 107.475 0.016 1.557
## 37 37 37 0.517 0.055 0.003 1.028 0.018 107.000 0.015 1.501
## 38 38 38 0.111 -0.002 0.000 1.034 0.018 107.552 0.016 1.580
## 39 39 39 -1.142 -0.154 0.024 1.014 0.018 105.440 0.018 1.834
## 40 40 40 0.064 -0.008 0.000 1.034 0.018 107.581 0.016 1.592
## 41 41 41 0.111 -0.002 0.000 1.034 0.018 107.552 0.016 1.580
## 42 42 42 -0.980 -0.137 0.019 1.019 0.018 106.070 0.018 1.806
## 43 43 43 -0.819 -0.119 0.014 1.024 0.018 106.594 0.018 1.766
## 44 44 44 0.858 0.103 0.011 1.019 0.018 106.163 0.014 1.416
## 45 45 45 1.857 0.246 0.059 0.971 0.016 102.112 0.012 1.220
## 46 46 46 0.858 0.103 0.011 1.019 0.018 106.163 0.014 1.416
## 47 47 47 -0.514 -0.083 0.007 1.031 0.018 107.279 0.017 1.714
## 48 48 48 0.064 -0.008 0.000 1.034 0.018 107.581 0.016 1.592
## 49 49 49 0.202 0.010 0.000 1.033 0.018 107.475 0.016 1.557
## 50 50 50 0.202 0.010 0.000 1.033 0.018 107.475 0.016 1.557
## 51 51 51 -0.076 -0.027 0.001 1.034 0.018 107.622 0.016 1.615
## 52 52 52 0.560 0.061 0.004 1.027 0.018 106.910 0.015 1.490
## 53 53 53 -0.980 -0.137 0.019 1.019 0.018 106.070 0.018 1.806
## 54 54 54 -0.414 -0.071 0.005 1.032 0.018 107.423 0.017 1.689
## 55 55 55 -0.171 -0.040 0.002 1.034 0.019 107.609 0.016 1.640
## 56 56 56 0.427 0.042 0.002 1.030 0.018 107.168 0.015 1.512
## 57 57 57 0.111 -0.002 0.000 1.034 0.018 107.552 0.016 1.580
## 58 58 58 -0.926 -0.131 0.017 1.021 0.018 106.259 0.018 1.793
## 59 59 59 -0.219 -0.046 0.002 1.034 0.019 107.590 0.017 1.652
## 60 60 60 -1.426 -0.182 0.032 1.001 0.017 104.061 0.019 1.890
## 61 61 61 -1.485 -0.188 0.034 0.999 0.017 103.722 0.019 1.904
## 62 62 62 0.111 -0.002 0.000 1.034 0.018 107.552 0.016 1.580
## 63 63 63 0.604 0.067 0.005 1.026 0.018 106.817 0.015 1.479
infl_step3$plot
Step 3b: Outlier diagnostics
After the normality assessment, we screen the historical control studies for
potential outliers. The function assess_continuous_summary_hcd_outliers()
provides two options through the outlier_method argument.
With outlier_method = "raw_values", the outlier checks are applied directly to
the observed study-level means. This gives a descriptive assessment of unusual
endpoint values and includes several screening approaches, such as empirical
quantiles, IQR-based boxplots, z-scores, and other outlier tests.
With outlier_method = "model_residuals", a random-effects meta-analysis model
is first fitted, and the outlier checks are applied to the model residuals. This
option is closer to the approach used in the published EFSA Annex case study,
where externally studentized residuals were used to identify potentially unusual
studies.
In this vignette, we use outlier_method = "model_residuals" with externally
studentized residuals to reproduce the boxplot-style diagnostics from the
published continuous endpoint case study. These results are interpreted as
screening diagnostics only, not as automatic exclusion rules.
outlier_resid_orig <- assess_continuous_summary_hcd_outliers(
data = hcd_cont,
response_col = "bioch_means",
se_col = "se",
study_id_col = "StudyID",
endpoint_label = "Clinical chemistry endpoint",
scale_label = "original scale",
outlier_method = "model_residuals",
residual_type = "rstudent",
model_method = "REML",
model_test = "adhoc",
histogram_binwidth = 0.25
)
outlier_resid_log <- assess_continuous_summary_hcd_outliers(
data = hcd_cont,
response_col = "log_bioch_means",
se_col = "log_se",
study_id_col = "StudyID",
endpoint_label = "Clinical chemistry endpoint",
scale_label = "log-transformed scale",
outlier_method = "model_residuals",
residual_type = "rstudent",
model_method = "REML",
model_test = "adhoc",
histogram_binwidth = 0.25
)
patchwork::wrap_plots(
outlier_resid_orig$plots$iqr_boxplot,
outlier_resid_log$plots$iqr_boxplot,
ncol = 2
)
Step 3d: Investigation of study-level covariates
After evaluating normality and influence diagnostics, we investigate whether the
endpoint is associated with available study-level covariates. Such associations
may indicate systematic differences between studies and can help determine
whether restrictions of the historical control dataset are needed before the
final variability modelling step.
Following the EFSA Annex example, we explore potential associations between the
endpoint and selected study characteristics, including route of exposure and
study year. These analyses are exploratory and are intended to identify
important sources of heterogeneity rather than to perform formal confirmatory
inference.
The assessment includes subgroup summaries, meta-regression models, and
graphical exploration of study-level effects. The results are then used to
decide whether subsets of the historical control data should be excluded or
analysed separately in later steps of the workflow.
The first covariate considered is route of exposure. The function
analyze_continuous_hcd_meta_covariate() fits subgroup meta-analyses for each
route category and provides the pooled mean, 95% confidence interval,
between-study variance, within-subgroup heterogeneity, and the test for subgroup
differences.
route_forest_step3 <- analyze_continuous_hcd_meta_covariate(
data = hcd_cont,
covariate_col = "route",
covariate_type = "categorical",
mean_col = "bioch_means",
sd_col = "sds",
n_col = "n_anim",
se_col = "se",
study_id_col = "StudyID",
scale = "orig",
covariate_label = "Route of exposure",
endpoint_label = "Clinical chemistry endpoint",
method_tau = "REML",
method_random_ci = "HK",
tau_common = TRUE,
prediction = TRUE,
show_plot = TRUE,
digits = 3
)
route_forest_step3$formatted_table
## Covariate Scale Subgroup No. studies Pooled mean 95% CI
## 1 Route of exposure orig dermal 5 1.436 [0.909; 1.962]
## 2 Route of exposure orig diet 26 1.309 [1.242; 1.376]
## 3 Route of exposure orig drinking 1 1.370 [0.948; 1.792]
## 4 Route of exposure orig gavage 23 1.307 [1.228; 1.385]
## 5 Route of exposure orig i.v. 1 1.780 [1.303; 2.257]
## 6 Route of exposure orig inhalation 6 1.326 [0.979; 1.673]
## 7 Route of exposure orig water 1 1.380 [0.956; 1.804]
## tau^2 Q within I^2 (%) Q subgroup df subgroup p subgroup
## 1 0.02 24.346 83.6 4.297 6 0.637
## 2 0.02 26.994 7.4 4.297 6 0.637
## 3 0.02 0.000 NA 4.297 6 0.637
## 4 0.02 28.562 23.0 4.297 6 0.637
## 5 0.02 0.000 NA 4.297 6 0.637
## 6 0.02 20.407 75.5 4.297 6 0.637
## 7 0.02 0.000 NA 4.297 6 0.637
The second covariate considered is study start year. This analysis checks
whether the historical control means differ systematically across calendar
years. The output reproduces the subgroup meta-analysis summary shown in Table 3
of the EFSA Annex example, including the pooled mean by year, 95% confidence
interval, heterogeneity measures, and the test for subgroup differences.
year_forest_step3 <- analyze_continuous_hcd_meta_covariate(
data = hcd_cont,
covariate_col = "year_2",
covariate_type = "categorical",
mean_col = "bioch_means",
sd_col = "sds",
n_col = "n_anim",
se_col = "se",
study_id_col = "StudyID",
scale = "orig",
covariate_label = "Study start year",
endpoint_label = "Clinical chemistry endpoint",
method_tau = "REML",
method_random_ci = "HK",
tau_common = TRUE,
prediction = TRUE,
show_plot = TRUE,
digits = 3
)
year_forest_step3$formatted_table
## Covariate Scale Subgroup No. studies Pooled mean 95% CI tau^2
## 1 Study start year orig 2012 1 1.110 [0.718; 1.502] 0.02
## 2 Study start year orig 2013 17 1.367 [1.259; 1.475] 0.02
## 3 Study start year orig 2014 13 1.373 [1.168; 1.577] 0.02
## 4 Study start year orig 2015 18 1.321 [1.242; 1.399] 0.02
## 5 Study start year orig 2016 7 1.310 [1.207; 1.413] 0.02
## 6 Study start year orig 2017 7 1.226 [1.074; 1.377] 0.02
## Q within I^2 (%) Q subgroup df subgroup p subgroup
## 1 0.000 NA 4.61 5 0.465
## 2 26.738 40.2 4.61 5 0.465
## 3 49.037 75.5 4.61 5 0.465
## 4 15.967 0.0 4.61 5 0.465
## 5 3.020 0.0 4.61 5 0.465
## 6 6.948 13.6 4.61 5 0.465
As an additional package-based screening step, we also apply
assess_continuous_summary_hcd_covariates(). This function performs simple
exploratory tests, such as ANOVA and Kruskal-Wallis tests, directly on the
study-level means.
covariate_screen_step3 <- assess_continuous_summary_hcd_covariates(
data = hcd_cont,
response_col = "bioch_means",
categorical_covariates = c("route", "year_2"),
covariate_labels = c(
route = "Route of exposure",
year_2 = "Study start year"
)
)
covariate_screen_step3$formatted_table
## covariate covariate_label method n n_groups
## 1 route Route of exposure Kruskal-Wallis 63 7
## 2 route Route of exposure ANOVA 63 7
## 3 year_2 Study start year Kruskal-Wallis 63 6
## 4 year_2 Study start year ANOVA 63 6
## statistic_name statistic df1 df2 p_value note
## 1 Kruskal-Wallis chi-squared 3.667 6.000 <NA> 0.7217
## 2 F 1.247 6.000 56.000 0.2969
## 3 Kruskal-Wallis chi-squared 5.238 5.000 <NA> 0.3875
## 4 F 1.128 5.000 57.000 0.3560
Step 3: Selection of the final HCD subset
The decision to restrict the historical control dataset is not based only on
statistical significance. Although the subgroup analyses did not show strong
statistical evidence of an association between the endpoint and route of
exposure or study year, the final selection of HCD also requires scientific and
regulatory judgement.
Following the EFSA Annex example, the statistician and domain experts considered
several exposure routes insufficiently comparable with the index study route.
Therefore, studies with inhalation, water, drinking, dermal, and intravenous
exposure were excluded. Studies with gavage exposure were retained, because
their pooled estimate was very similar to that for diet and gavage represented a
substantial part of the HCD.
The single study conducted in 2012 was also excluded, because it was the oldest
study and furthest away from the index study period. After these restrictions,
the final HCD subset contains studies with diet or gavage exposure only,
excluding studies from 2012. This leads to a final set of 48 historical control
studies, which is used in Step 4.
After applying these selection criteria, it is good practice to repeat the main
exploratory checks, including normality diagnostics, outlier assessment, and
covariate investigation, on the reduced HCD dataset. In the EFSA Annex example,
these repeated checks are mentioned but not shown in detail, and no further
issues were reported.
hcd_cont_step4 <- hcd_cont[
hcd_cont$route %in% c("diet", "gavage") &
hcd_cont$year_2 > 2012,
]
## Quick checks
nrow(hcd_cont_step4)
## [1] 48
table(hcd_cont_step4$route)
##
## diet gavage
## 25 23
table(hcd_cont_step4$year_2)
##
## 2013 2014 2015 2016 2017
## 14 6 14 7 7
## Studies excluded by this restriction
excluded_step3 <- hcd_cont[
!(hcd_cont$StudyID %in% hcd_cont_step4$StudyID),
]
excluded_step3[, c("StudyID", "route", "year_2", "bioch_means", "se")]
## StudyID route year_2 bioch_means se
## 1 1 diet 2012 1.11 0.143
## 4 4 inhalation 2013 1.14 0.145
## 6 6 inhalation 2013 1.81 0.202
## 9 9 water 2013 1.38 0.165
## 19 19 inhalation 2014 0.95 0.130
## 21 21 i.v. 2014 1.78 0.199
## 23 23 dermal 2014 1.96 0.215
## 24 24 inhalation 2014 1.28 0.156
## 26 26 dermal 2014 1.96 0.215
## 30 30 dermal 2014 0.97 0.131
## 31 31 inhalation 2014 1.76 0.197
## 36 36 dermal 2015 1.37 0.164
## 38 38 inhalation 2015 1.35 0.162
## 41 41 dermal 2015 1.35 0.162
## 49 49 drinking 2015 1.37 0.164
Overall, the Step 3 investigation suggests that some study-level differences may
need to be considered before modelling HCD variability. The diagnostic and
covariate checks are used here to support the selection of the HCD set carried
forward to Step 4. Therefore, the subsequent modelling is based on the studies
retained after these exploratory and statistical checks.
Step 4: Modelling HCD variability
After selecting the final HCD subset, the next step is to characterise the
between-study variability of the historical control data using a random-effects
meta-analysis model.
The objective of this step is not only to estimate the average historical
control response, but also to quantify the variability that can reasonably be
expected between future studies. This variability estimate is later used in Step
5 to evaluate whether the concurrent control of the index study is compatible
with the historical control distribution.
In the EFSA Annex example, the analysis is performed using random-effects
meta-analysis models fitted on the study-level summary statistics. The
function model_continuous_hcd_variability() reproduces this workflow and
provides:
- pooled HCD estimates,
- confidence intervals,
- prediction intervals,
- heterogeneity statistics,
- between-study variance estimates,
- optional subgroup-adjusted models,
- and forest plots.
The function can analyse both the original scale and the log-transformed scale.
The function combines implementations based on both the meta and metafor
packages, allowing users to compare results across modelling frameworks while
keeping a unified output structure within the package.
step4_cont_orig <- model_continuous_hcd_variability(
data = hcd_cont_step4,
mean_col = "bioch_means",
sd_col = "sds",
n_col = "n_anim",
se_col = "se",
study_id_col = "StudyID",
hcd_indicator_col = NULL,
endpoint_label = "Clinical chemistry endpoint",
scale = "orig",
method_tau = "REML",
method_random_ci = "HK",
prediction = TRUE,
make_forest = TRUE,
show_forest = TRUE
)
step4_cont_orig$formatted_summary_table
## model fitting_package adjusted covariates prediction_at
## 1 meta_unadjusted meta FALSE none overall HCD mean
## 2 metafor_unadjusted metafor FALSE none overall HCD mean
## estimate ci_lower ci_upper pi_lower pi_upper tau2 tau Q Q_df
## 1 1.3037 1.2543 1.3531 1.1691 1.4383 0.0039 0.0622 53.8425 47
## 2 1.3037 1.2541 1.3533 1.1686 1.4388 0.0039 0.0625 53.8896 47
## Q_pvalue
## 1 0.2290
## 2 0.2276
The package gives a HCD 95% prediction interval of [1.169, 1.438], while the
EFSA Annex B report gives [1.162, 1.446]. The pooled HCD mean, tau² estimate,
and Q statistic are reproduced very closely: pooled mean = 1.304, tau² = 0.004,
and Q = 53.84 with 47 degrees of freedom.
The small difference between the package prediction interval and the published
Annex value is due to the Annex code rounding intermediate quantities before
computing the PI. In particular, the Annex calculation uses the rounded pooled
mean, rounded squared standard error, and rounded tau² estimate. When the same
intermediate rounding is applied, the published PI [1.162, 1.446] is reproduced
exactly.
## -------------------------------------------------------------------------
## Ad hoc reproduction of the Annex B Table 4 PI calculation
## using the same intermediate rounding logic as the EFSA code
## -------------------------------------------------------------------------
m_step4 <- step4_cont_orig$models$meta_unadjusted
k_hcd <- nrow(hcd_cont_step4)
annex_step4_pi <- data.frame(
tau2_estimate_label = "Central 95%CI",
pooled_hcd = round(m_step4$TE.random, 3),
pooled_hcd_se2 = round(m_step4$seTE.random^2, 3),
tau2 = round(m_step4$tau2, 3),
df = k_hcd - 2
)
annex_step4_pi$pi_variance <- round(
annex_step4_pi$pooled_hcd_se2 + annex_step4_pi$tau2,
3
)
annex_step4_pi$pi_lower <- round(
annex_step4_pi$pooled_hcd -
qt(0.975, annex_step4_pi$df) * sqrt(annex_step4_pi$pi_variance),
3
)
annex_step4_pi$pi_upper <- round(
annex_step4_pi$pooled_hcd +
qt(0.975, annex_step4_pi$df) * sqrt(annex_step4_pi$pi_variance),
3
)
annex_step4_pi
## tau2_estimate_label pooled_hcd pooled_hcd_se2 tau2 df pi_variance pi_lower
## 1 Central 95%CI 1.304 0.001 0.004 46 0.005 1.162
## pi_upper
## 1 1.446
Step 5: Comparing the concurrent control with HCD
In Step 5, the concurrent control (CC) from the index study is compared with the
prediction interval obtained from the selected HCD model in Step 4. The aim is
to assess whether the CC appears compatible with the historical control
distribution before deciding whether HCD can be integrated into the index study
analysis.
For continuous endpoints, the package function compare_continuous_cc_to_hcd()
implements this comparison. Two approaches are available through the approach
argument:
approach = "population_mean" compares the CC confidence interval with the
HCD prediction interval for a future control population mean.approach = "new_group_mean" compares the observed CC sample mean with the
HCD prediction interval for a future control group sample mean.
In the EFSA Annex example, Approach 1 is used. In package terms, this
corresponds to approach = "population_mean". The CC is summarised using a 90%
confidence interval, computed using a Student t distribution, matching the Annex
code.
The function returns the CC mean and confidence interval, the HCD estimate and
prediction interval, a compatibility indicator, and an optional plot showing the
CC estimate relative to the HCD prediction interval.
step5a_cont <- compare_continuous_cc_to_hcd(
step4_result = step4_cont_orig,
approach = "population_mean",
cc_data = cc_cont,
mean_col = "bioch_means",
se_col = "se",
n_col = "n_anim",
alpha = 0.05,
use_manual_pi = TRUE,
make_plots = TRUE,
digits = 3
)
step5a_cont$formatted_summary_table
## Model Adjusted Covariates Covariate values used
## 1 meta_unadjusted_manual_HK_style FALSE none None
## CC mean CC 90% CI lower CC 90% CI upper HCD 95% PI lower HCD 95% PI upper
## 1 1.05 0.756 1.344 1.169 1.438
## Compatible
## 1 FALSE
step5a_cont$plots$population_mean_meta_unadjusted_manual_HK_style
The concurrent control mean was 1.05 with a 90% confidence interval of
[0.756, 1.344], while the pooled HCD estimate from Step 4 was 1.304 with a 95%
prediction interval of [1.169, 1.438]. Although the intervals partially overlap,
the lower part of the CC confidence interval extends below the lower bound of
the HCD prediction interval. Accordingly, the compatibility assessment concluded
that similarity between the CC and HCD could not be confirmed
(compatible = FALSE). This reproduces the main conclusion of the EFSA Annex
example and motivates the additional investigations considered in Step 5b.
According to the EFSA workflow, an apparent inconsistency between the concurrent
control and the HCD should not automatically be addressed only through
sensitivity analyses. Step 5b may also involve reconsidering the comparability
of the HCD, reviewing possible biological or experimental differences between
studies, re-evaluating covariates included in the analysis, investigating
potential outlying studies, or assessing whether alternative HCD subsets are
more appropriate. In this case study, the additional investigations mainly focus
on sensitivity analyses similar to those discussed in the EFSA Annex example.
For the original selected HCD dataset, compatibility between the CC and HCD
could not be confirmed. To illustrate how the workflow proceeds when the
apparent mismatch is investigated and a compatible HCD set is obtained, the EFSA
Annex uses an additional artificial dataset. This dataset is based on the second
sheet of the provided Excel file and is used only for illustration in Step 5b
and the following steps.
To keep the vignette reproducible, we stored the generated version of this
artificial dataset in the package. In the next chunk, we read the reduced Step
5b HCD dataset directly from inst/extdata/ and then rerun Steps 4 and 5 using
this dataset.
## Reduced HCD dataset used in Annex Step 5b
## This dataset was generated from DATA_forSO_v2 using the Annex code.
hcd_cont_step5b_file <- system.file(
"extdata",
"appendix_cont_step5b_HCD_final_v2.csv",
package = "efsaHCD"
)
hcd_cont_step5b <- read.csv(hcd_cont_step5b_file)
## Quick checks
nrow(hcd_cont_step5b)
## [1] 25
table(hcd_cont_step5b$route)
##
## diet
## 25
table(hcd_cont_step5b$year_2)
##
## 2013 2014 2015 2016 2017
## 8 2 6 4 5
hcd_cont_step5b[, c("StudyID", "route", "year_2", "bioch_means", "se")]
## StudyID route year_2 bioch_means se
## 1 2 diet 2013 1.00 0.134
## 2 5 diet 2013 1.77 0.198
## 3 7 diet 2013 1.19 0.149
## 4 12 diet 2013 1.20 0.150
## 5 13 diet 2013 1.33 0.161
## 6 14 diet 2013 1.02 0.136
## 7 15 diet 2013 1.08 0.140
## 8 17 diet 2013 0.88 0.125
## 9 22 diet 2014 1.23 0.152
## 10 27 diet 2014 0.89 0.126
## 11 37 diet 2015 1.77 0.198
## 12 39 diet 2015 1.10 0.142
## 13 45 diet 2015 1.45 0.170
## 14 46 diet 2015 0.98 0.132
## 15 47 diet 2015 1.21 0.151
## 16 48 diet 2015 1.16 0.147
## 17 50 diet 2016 1.33 0.161
## 18 51 diet 2016 0.78 0.118
## 19 54 diet 2016 1.05 0.138
## 20 56 diet 2016 0.87 0.124
## 21 57 diet 2017 1.01 0.135
## 22 58 diet 2017 1.16 0.147
## 23 59 diet 2017 0.95 0.130
## 24 62 diet 2017 1.77 0.198
## 25 63 diet 2017 0.77 0.117
Because the output of Step 4 is used directly as the input for Step 5, the
variability model must first be refitted using the reduced Step 5b HCD dataset.
We therefore rerun Step 4 on the reduced dataset and then repeat the Step 5
comparison.
The function investigate_continuous_cc_hcd_mismatch() provides additional
sensitivity analyses for cases where the initial CC-HCD compatibility assessment
is not confirmed. In addition to varying the prediction interval source and
confidence level, the function can also investigate the impact of different
assumptions for the between-study variability (tau²). Here, we evaluate the
prediction interval using the central estimate of tau², as well as the lower
and upper bounds of its 95% confidence interval, reproducing the sensitivity
investigation reported in the EFSA Annex example.
step4_cont_step5b <- model_continuous_hcd_variability(
data = hcd_cont_step5b,
mean_col = "bioch_means",
sd_col = "sds",
n_col = "n_anim",
se_col = "se",
study_id_col = "StudyID",
hcd_indicator_col = NULL,
endpoint_label = "Clinical chemistry endpoint",
scale = "orig",
method_tau = "REML",
method_random_ci = "HK",
prediction = TRUE,
make_forest = TRUE,
show_forest = TRUE
)
step4_cont_step5b$formatted_summary_table
## model fitting_package adjusted covariates prediction_at
## 1 meta_unadjusted meta FALSE none overall HCD mean
## 2 metafor_unadjusted metafor FALSE none overall HCD mean
## estimate ci_lower ci_upper pi_lower pi_upper tau2 tau Q Q_df
## 1 1.1311 1.019 1.2432 0.6747 1.5875 0.0459 0.2143 73.0055 24
## 2 1.1312 1.024 1.2384 0.6757 1.5867 0.0460 0.2145 73.0956 24
## Q_pvalue
## 1 0
## 2 0
step5b_cont <- investigate_continuous_cc_hcd_mismatch(
step4_result = step4_cont_step5b,
approach = "population_mean",
alpha_values = 0.05,
pi_source = "manual",
tau2_scenarios = c("central", "lower", "upper"),
cc_ci_method = "t",
cc_data = cc_cont,
mean_col = "bioch_means",
se_col = "se",
n_col = "n_anim",
make_plot = TRUE,
digits = 3
)
step5b_cont$formatted_sensitivity_table
## Approach CC interval level HCD PI level Tau2 scenario Adjusted
## 1 population_mean 0.9 0.95 central FALSE
## 2 population_mean 0.9 0.95 lower FALSE
## 3 population_mean 0.9 0.95 upper FALSE
## Covariates CC mean CC CI lower CC CI upper HCD PI lower HCD PI upper
## 1 none 1.05 0.756 1.344 0.674 1.589
## 2 none 1.05 0.756 1.344 0.807 1.455
## 3 none 1.05 0.756 1.344 0.379 1.883
## Compatible
## 1 TRUE
## 2 FALSE
## 3 TRUE
Using the reduced HCD dataset, the pooled HCD estimate decreased to 1.131 and
the prediction intervals became wider. Under the central and upper tau²
scenarios, the CC 90% confidence interval was fully contained within the HCD
prediction interval, and compatibility between the CC and HCD was therefore
supported (compatible = TRUE). Under the lower-bound tau² scenario, the
prediction interval became narrower and the lower part of the CC confidence
interval extended slightly below the HCD interval, leading to
compatible = FALSE. However, the discrepancy was small
(gap_outside = 0.051), and the overall conclusions reproduce the
interpretation reported in the EFSA Annex example: after reconsidering the HCD
selection and exploring plausible variability assumptions, the CC and HCD were
considered sufficiently similar for subsequent analyses.
Step 6: Analysis of the index study with and without HCD
After establishing sufficient compatibility between the concurrent control and
HCD in Step 5, the selected HCD information can be incorporated into the
analysis of the index study. For continuous endpoints, the EFSA workflow
combines the selected HCD information with the concurrent control by
constructing a pseudo-control row derived from the pooled HCD estimate obtained
in Step 4.
In the package, the function create_continuous_hcd_pseudo_control() generates
this pseudo-control dataset. The pseudo-control output contains a raw
pseudo-control dataset and a formatted summary table. The formatted table is
intentionally simplified and displays only the original-scale mean and SD, the
log-scale mean and SD, and the effective sample size used for the HCD
pseudo-control. In the EFSA Annex example, the pseudo-control is constructed
from the Step 4 prediction interval used in the compatibility assessment.
The resulting pseudo-control is then combined with the original index study data
and analyzed using the function analyze_continuous_index_with_hcd(). This
function fits meta-regression models both:
- without HCD (concurrent control only), and
- with the HCD pseudo-control included.
Important: The analyses shown with and without HCD are included only for
illustrative purposes to demonstrate the package workflow and outputs. In a real
application, the decision to include HCD should be based on the conclusions from
the previous workflow steps.
Two dose-response parameterizations are evaluated:
- a continuous dose model (
dose_model = "continuous"), corresponding to a
linear dose-response relationship, and
- a categorical dose model (
dose_model = "categorical"), where each dose level
is compared separately with the control group.
The models are fitted using restricted maximum likelihood (REML), reproducing
the approach used in the EFSA Annex example. For the HCD-integrated analysis, a
random-effect structure is included for study-level variability, while the
concurrent-control-only analysis reduces to a standard fixed-effect regression
because only a single study is analyzed.
The function returns descriptive summaries, coefficient tables, overall
dose-effect comparison results, fitted model objects, analysis datasets, and
graphical displays of the dose-response relationship.
## -------------------------------------------------------------------------
## Step 6: index study analysis with and without HCD pseudo-control
## -------------------------------------------------------------------------
## Prepare the index study rows from the generated DATA_forSO_v2 dataset
cont_step5b_file <- system.file(
"extdata",
"appendix_cont_generated_DATA_forSO_v2.csv",
package = "efsaHCD"
)
cont_step6 <- read.csv(cont_step5b_file)
index_cont_step6_raw <- cont_step6[cont_step6$HCD > 1, ]
index_cont_step6 <- prepare_continuous_index_study(
data = index_cont_step6_raw,
endpoint_mean_col = "bioch_means",
endpoint_se_col = "se",
n_col = "n_anim",
dose_col = "dose",
study_id_col = "StudyID",
control_dose = 0,
endpoint_name = "Clinical chemistry endpoint",
print_summary = FALSE
)
## Create the HCD pseudo-control row from the Step 4 result used in Step 5b.
## Use the manual EFSA-style Step 4 estimate, matching the Annex workflow.
hcd_pseudo_step6 <- create_continuous_hcd_pseudo_control(
step4_result = step4_cont_step5b,
index_data = index_cont_step6,
pi_source = "manual",
models = "meta_unadjusted_manual_HK_style",
control_dose = 0,
hcd_study_id = "HCD",
hcd_group_type = "historical_control_pooled",
hcd_index_code = 1,
n_value_method = "mean_hcd",
digits = 4
)
hcd_pseudo_step6$formatted_pseudo_control
## Original mean Original SD Log mean Log SD Effective n
## 1 1.1311 0.5416 0.1232 0.4789 6
## Run Step 6 models:
## - continuous dose model: reproduces Table 6
## - categorical dose model: reproduces Table 7
step6_cont <- analyze_continuous_index_with_hcd(
index_data = index_cont_step6,
hcd_pseudo_control = hcd_pseudo_step6,
dose_model = "both",
covariate_strategy = "none",
include_hcd = TRUE,
include_cc_only = TRUE,
control_dose = 0,
method = "REML",
cc_only_random = FALSE,
hcd_random = TRUE,
make_plots = TRUE,
digits = 4
)
## Descriptive table / Figure 9 data
step6_cont$formatted_descriptive_table
## analysis_type hcd_model covariate_strategy study_ID
## 1 cc_only none none 64
## 2 cc_only none none 65
## 3 cc_only none none 66
## 4 cc_only none none 67
## 9 with_hcd meta_unadjusted_manual_HK_style none HCD_1
## 10 with_hcd meta_unadjusted_manual_HK_style none 64
## 11 with_hcd meta_unadjusted_manual_HK_style none 65
## 12 with_hcd meta_unadjusted_manual_HK_style none 66
## 13 with_hcd meta_unadjusted_manual_HK_style none 67
## dose mean_value se_value ci_lower ci_upper n_value group_type
## 1 0 1.0500 0.1380 0.7795 1.3205 5 concurrent_control
## 2 600 1.0700 0.1400 0.7956 1.3444 5 treated
## 3 1400 1.3700 0.1640 1.0486 1.6914 5 treated
## 4 5000 1.8400 0.2040 1.4402 2.2398 5 treated
## 9 0 1.1311 0.2211 0.6977 1.5645 6 historical_control_pooled
## 10 0 1.0500 0.1380 0.7795 1.3205 5 concurrent_control
## 11 600 1.0700 0.1400 0.7956 1.3444 5 treated
## 12 1400 1.3700 0.1640 1.0486 1.6914 5 treated
## 13 5000 1.8400 0.2040 1.4402 2.2398 5 treated
## hcd_adjusted hcd_covariates
## 1 NA <NA>
## 2 NA <NA>
## 3 NA <NA>
## 4 NA <NA>
## 9 FALSE none
## 10 NA <NA>
## 11 NA <NA>
## 12 NA <NA>
## 13 NA <NA>
## Table 6: continuous dose meta-regression
step6_cont$formatted_coefficient_table[
step6_cont$formatted_coefficient_table$`Dose model` == "continuous",
]
## Analysis HCD model Dose model Model term Estimate
## 1 cc_only none continuous intrcpt 1.0407
## 2 cc_only none continuous .dose_cont 0.0002
## 7 with_hcd meta_unadjusted_manual_HK_style continuous intrcpt 1.0554
## 8 with_hcd meta_unadjusted_manual_HK_style continuous .dose_cont 0.0002
## SE 95% CI lower 95% CI upper Term p-value Overall dose-effect p-value
## 1 0.0975 0.8497 1.2317 < 0.001 < 0.001
## 2 0.0000 0.0001 0.0003 < 0.001 < 0.001
## 7 0.0892 0.8806 1.2302 < 0.001 < 0.001
## 8 0.0000 0.0001 0.0003 < 0.001 < 0.001
## Table 7: dose as categorical factor
step6_cont$formatted_coefficient_table[
step6_cont$formatted_coefficient_table$`Dose model` == "categorical",
]
## Analysis HCD model Dose model Model term Estimate
## 3 cc_only none categorical intrcpt 1.0500
## 4 cc_only none categorical .dose_fac600 0.0200
## 5 cc_only none categorical .dose_fac1400 0.3200
## 6 cc_only none categorical .dose_fac5000 0.7900
## 9 with_hcd meta_unadjusted_manual_HK_style categorical intrcpt 1.0727
## 10 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac600 -0.0027
## 11 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac1400 0.2973
## 12 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac5000 0.7673
## SE 95% CI lower 95% CI upper Term p-value Overall dose-effect p-value
## 3 0.1380 0.7795 1.3205 < 0.001 0.00547
## 4 0.1966 -0.3653 0.4053 0.91896 0.00547
## 5 0.2143 -0.1001 0.7401 0.13544 0.00547
## 6 0.2463 0.3073 1.2727 0.00134 0.00547
## 9 0.1171 0.8433 1.3022 < 0.001 0.00514
## 10 0.1825 -0.3604 0.3550 0.98805 0.00514
## 11 0.2015 -0.0977 0.6922 0.14014 0.00514
## 12 0.2352 0.3063 1.2283 0.00111 0.00514
## Model-level tests: residual heterogeneity and moderator tests
step6_cont$formatted_comparison_table
## Analysis HCD model Dose model Model term Estimate
## 2 cc_only none continuous .dose_cont 0.0002
## 4 cc_only none categorical .dose_fac600 0.0200
## 5 cc_only none categorical .dose_fac1400 0.3200
## 6 cc_only none categorical .dose_fac5000 0.7900
## 8 with_hcd meta_unadjusted_manual_HK_style continuous .dose_cont 0.0002
## 10 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac600 -0.0027
## 11 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac1400 0.2973
## 12 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac5000 0.7673
## SE 95% CI lower 95% CI upper Term p-value Overall dose-effect p-value
## 2 0.0000 0.0001 0.0003 < 0.001 < 0.001
## 4 0.1966 -0.3653 0.4053 0.91896 0.00547
## 5 0.2143 -0.1001 0.7401 0.13544 0.00547
## 6 0.2463 0.3073 1.2727 0.00134 0.00547
## 8 0.0000 0.0001 0.0003 < 0.001 < 0.001
## 10 0.1825 -0.3604 0.3550 0.98805 0.00514
## 11 0.2015 -0.0977 0.6922 0.14014 0.00514
## 12 0.2352 0.3063 1.2283 0.00111 0.00514
## Plots
## Figure 9-style forest plots
step6_cont$plots$forest__with_hcd__none__meta_unadjusted_manual_HK_style__continuous
## Figure 10 left: linear dose-response without HCD
step6_cont$plots$dose_response__cc_only__continuous
## Figure 10 right: linear dose-response with HCD
step6_cont$plots$dose_response__with_hcd__none__meta_unadjusted_manual_HK_style__continuous
The continuous Step 6 output retains the existing model-specific forest and
dose-response plots. No additional source-comparison plotting layer is introduced
for continuous endpoints.
## -------------------------------------------------------------------------
## Ad hoc report-style Step 6 tables
## -------------------------------------------------------------------------
coef_tab <- step6_cont$raw_coefficient_table
make_step6_panel_text <- function(coef_tab, analysis_type, dose_model) {
z <- coef_tab[
coef_tab$analysis_type == analysis_type &
coef_tab$dose_model == dose_model,
]
tau2 <- unique(z$tau2_or_sigma2)[1]
qm <- unique(z$QM)[1]
qm_df <- unique(z$QM_df)[1]
qm_p <- unique(z$QM_p)[1]
out <- data.frame(
section = c(
"Variance components",
"Test for moderators",
"Model results"
),
text = c(
paste0("Tau2: ", round(tau2, 2)),
paste0(
"QM(df = ", qm_df, ") = ",
round(qm, 2),
", p = ",
signif(qm_p, 3)
),
""
)
)
out
}
make_step6_model_results <- function(coef_tab, analysis_type, dose_model) {
z <- coef_tab[
coef_tab$analysis_type == analysis_type &
coef_tab$dose_model == dose_model,
]
z$term_clean <- z$term
z$term_clean[z$term == "intrcpt"] <- "Intercept"
z$term_clean[z$term == ".dose_cont"] <- "Dose"
z$term_clean[z$term == ".dose_fac600"] <- "Dose600"
z$term_clean[z$term == ".dose_fac1400"] <- "Dose1400"
z$term_clean[z$term == ".dose_fac5000"] <- "Dose5000"
data.frame(
term = z$term_clean,
Est = round(z$estimate, 4),
SE = round(z$se, 3),
p = ifelse(
z$p_value < 0.001,
"<0.001",
as.character(round(z$p_value, 4))
),
CILB = round(z$ci_lower, 4),
CIUB = round(z$ci_upper, 4),
row.names = NULL
)
}
## -------------------------------------------------------------------------
## Table 6: continuous dose-response model
## -------------------------------------------------------------------------
table6_with_hcd_header <- make_step6_panel_text(
coef_tab,
analysis_type = "with_hcd",
dose_model = "continuous"
)
table6_cc_only_header <- make_step6_panel_text(
coef_tab,
analysis_type = "cc_only",
dose_model = "continuous"
)
table6_with_hcd_results <- make_step6_model_results(
coef_tab,
analysis_type = "with_hcd",
dose_model = "continuous"
)
table6_cc_only_results <- make_step6_model_results(
coef_tab,
analysis_type = "cc_only",
dose_model = "continuous"
)
table6_with_hcd_header
## section text
## 1 Variance components Tau2: 0
## 2 Test for moderators QM(df = 1) = 12.11, p = 0.000502
## 3 Model results
table6_with_hcd_results
## term Est SE p CILB CIUB
## 1 Intercept 1.0554 0.089 <0.001 0.8806 1.2302
## 2 Dose 0.0002 0.000 <0.001 0.0001 0.0003
table6_cc_only_header
## section text
## 1 Variance components Tau2: 0
## 2 Test for moderators QM(df = 1) = 12.02, p = 0.000527
## 3 Model results
table6_cc_only_results
## term Est SE p CILB CIUB
## 1 Intercept 1.0407 0.097 <0.001 0.8497 1.2317
## 2 Dose 0.0002 0.000 <0.001 0.0001 0.0003
## -------------------------------------------------------------------------
## Table 7: categorical dose-response model
## -------------------------------------------------------------------------
table7_with_hcd_header <- make_step6_panel_text(
coef_tab,
analysis_type = "with_hcd",
dose_model = "categorical"
)
table7_cc_only_header <- make_step6_panel_text(
coef_tab,
analysis_type = "cc_only",
dose_model = "categorical"
)
table7_with_hcd_results <- make_step6_model_results(
coef_tab,
analysis_type = "with_hcd",
dose_model = "categorical"
)
table7_cc_only_results <- make_step6_model_results(
coef_tab,
analysis_type = "cc_only",
dose_model = "categorical"
)
table7_with_hcd_header
## section text
## 1 Variance components Tau2: 0
## 2 Test for moderators QM(df = 3) = 12.78, p = 0.00514
## 3 Model results
table7_with_hcd_results
## term Est SE p CILB CIUB
## 1 Intercept 1.0727 0.117 <0.001 0.8433 1.3022
## 2 Dose600 -0.0027 0.182 0.988 -0.3604 0.3550
## 3 Dose1400 0.2973 0.201 0.1401 -0.0977 0.6922
## 4 Dose5000 0.7673 0.235 0.0011 0.3063 1.2283
table7_cc_only_header
## section text
## 1 Variance components Tau2: 0
## 2 Test for moderators QM(df = 3) = 12.65, p = 0.00547
## 3 Model results
table7_cc_only_results
## term Est SE p CILB CIUB
## 1 Intercept 1.05 0.138 <0.001 0.7795 1.3205
## 2 Dose600 0.02 0.197 0.919 -0.3653 0.4053
## 3 Dose1400 0.32 0.214 0.1354 -0.1001 0.7401
## 4 Dose5000 0.79 0.246 0.0013 0.3073 1.2727
Step 7: Sensitivity analysis
The final step of the EFSA workflow evaluates the robustness of the conclusions
to modelling assumptions and HCD selection choices. In the continuous Annex
example, the main sensitivity investigations were already performed in Step 5b
by varying the between-study variability (tau²) and reconsidering the selected
HCD subset. These analyses showed that the compatibility conclusions were
generally stable across plausible assumptions. Therefore, no additional formal
analyses were reported in the published example.
Important: When implementing the EFSA HCD workflow, the key decisions,
assumptions, and methods to be explored in the sensitivity analysis should be
prespecified in the analysis protocol as far as possible.
Count data case study
This case study is based on the count-data example presented in Annex C of the
EFSA Opinion on the use and reporting of Historical Control Data (HCD). The
endpoint of interest is the number of eggs laid by a mallard duck hen during a
fixed 10-week period.
Before applying the workflow, it is important to distinguish this count endpoint
from the dichotomous endpoint used earlier. In the dichotomous case study, the
outcome is the number of animals with a finding out of a fixed number of animals
examined, and the statistical model is based on binomial-type data. In the count
endpoint case study, the outcome is a count measured for each experimental unit,
here the number of eggs laid by each bird during a fixed period. This is not a
number of animals with a finding, and therefore requires a different modelling
approach.
The individual-level HCD used in this case study were artificially generated
from summary data reported by Valverde-Garcia P, Springer T, Kramer V,
Foudoulakis M and Wheeler JR (2018), An avian reproduction study historical
control database: a tool for data interpretation, Regulatory Toxicology and
Pharmacology, 92, 295–302. The generation used the published summary information
for mallard duck reproduction studies and assumed 12 hens per study and a
10-week egg-laying period for each study.
An artificial index study was also generated for illustration purposes. The
index study includes one concurrent control group and three treatment groups
with increasing dose levels. The assumed study date for the index study is
February 2014.
The package includes the following datasets:
appendix_count_hcd.csv: historical control dataappendix_count_index-study.csv: index study data
Steps 1 and 2: Planning and data acceptability
In Step 1, the report specifies the intended workflow for the count endpoint.
The endpoint is a count variable, namely the number of eggs laid by each hen
during a fixed 10-week period. The report also states that individual-level data
are preferred for this endpoint, because they allow the subsequent analyses to
separate within-study and between-study variability.
In Step 2, the report discusses whether the HCD meet the legal, biological, and
statistical requirements for use. The HCD are based on studies from a single
laboratory and are treated as suitable for illustrating the workflow, although
the report notes that in a real regulatory assessment more detailed study
information would usually be requested. From a statistical perspective, the HCD
meet the minimum requirements because individual egg counts are available for
each hen in each historical control study.
Because these steps are mainly about planning, documentation, and data
acceptability, they are not implemented as separate analysis functions in the
package. The package workflow starts from the point where the data have already
been accepted for statistical evaluation.
Reading and preparing the count data
This section describes a data-management step used to put the HCD and
index-study data into a common analysis format. This is good practice for
reproducible analysis, but it is not a separate step of the EFSA decision
scheme.
The package contains two CSV files for the count-data case study: one for the
historical control data and one for the index study. We first read both files
and harmonise the column names before combining them into one dataset.
library(efsaHCD)
hcd_file <- system.file(
"extdata",
"appendix_count_hcd.csv",
package = "efsaHCD"
)
index_file <- system.file(
"extdata",
"appendix_count_index-study.csv",
package = "efsaHCD"
)
hcd_raw <- read.csv(hcd_file)
index_raw <- read.csv(index_file)
head(hcd_raw)
## study date mdelta bird neggs
## 1 1 2016-08-23 78.9 1 66
## 2 1 2016-08-23 78.9 2 31
## 3 1 2016-08-23 78.9 3 38
## 4 1 2016-08-23 78.9 4 32
## 5 1 2016-08-23 78.9 5 42
## 6 1 2016-08-23 78.9 6 65
head(index_raw)
## dose bird neggsb
## 1 0 1 47
## 2 0 2 45
## 3 0 3 52
## 4 0 4 24
## 5 0 5 59
## 6 0 6 44
The historical control data contain one row per bird, with the study identifier,
study date, bird identifier, and number of eggs. The index study also contains
one row per bird, with the dose group and number of eggs.
For the package workflow, we create a common structure across the two datasets.
The HCD are assigned index_study = "No" and dose 0, while the index-study
records are assigned index_study = "Yes".
hcd_dat <- hcd_raw
hcd_dat$dose <- 0
hcd_dat$index_study <- "No"
index_dat <- index_raw
index_dat$study <- "Index"
index_dat$date <- "2014-02-01"
index_dat$mdelta <- NA
index_dat$index_study <- "Yes"
index_dat$neggs <- index_dat$neggsb
count_dat <- rbind(
hcd_dat[, c("study", "date", "mdelta", "bird", "dose", "index_study", "neggs")],
index_dat[, c("study", "date", "mdelta", "bird", "dose", "index_study", "neggs")]
)
count_dat$study <- as.character(count_dat$study)
count_dat$index_study <- factor(count_dat$index_study, levels = c("No", "Yes"))
count_dat$dose <- as.numeric(count_dat$dose)
count_dat$date <- as.Date(count_dat$date)
count_dat$n_date <- as.numeric(count_dat$date)
count_dat$year <- as.numeric(format(count_dat$date, "%Y"))
head(count_dat)
## study date mdelta bird dose index_study neggs n_date year
## 1 1 2016-08-23 78.9 1 0 No 66 17036 2016
## 2 1 2016-08-23 78.9 2 0 No 31 17036 2016
## 3 1 2016-08-23 78.9 3 0 No 38 17036 2016
## 4 1 2016-08-23 78.9 4 0 No 32 17036 2016
## 5 1 2016-08-23 78.9 5 0 No 42 17036 2016
## 6 1 2016-08-23 78.9 6 0 No 65 17036 2016
We then prepare the data for the EFSA-style HCD workflow using
prepare_count_hcd_data(). This creates the main analysis dataset, the HCD
subset, the concurrent control subset, the full index-study subset, and the
combined control data.
count_prep <- prepare_count_hcd_data(
data = count_dat,
response_col = "neggs",
study_id_col = "study",
dose_col = "dose",
control_dose = 0,
index_study_col = "index_study",
index_study_value = "Yes",
subject_id_cols = c("study", "dose", "bird"),
subject_id_new = "bird_id",
date_col = "date",
date_format = "%Y-%m-%d",
factor_cols = c("index_study"),
covariate_cols = c("n_date", "year")
)
count_hcd <- count_prep$hcd
count_cc <- count_prep$cc
count_hcd_cc <- count_prep$hcd_cc
count_index <- count_prep$index_study
str(count_prep$settings)
## List of 23
## $ response_col : chr "neggs"
## $ study_id_col : chr "study"
## $ dose_col : chr "dose"
## $ control_dose : num 0
## $ index_study_col : chr "index_study"
## $ index_study_value : chr "Yes"
## $ subject_id_cols : chr [1:3] "study" "dose" "bird"
## $ subject_id_new : chr "bird_id"
## $ deduplicate : logi TRUE
## $ date_col : chr "date"
## $ date_format : chr "%Y-%m-%d"
## $ factor_cols : chr "index_study"
## $ relevel_factor_cols : NULL
## $ covariate_cols : chr [1:2] "n_date" "year"
## $ exclude_study_ids : NULL
## $ restrict_expr : NULL
## $ keep_missing_response: logi FALSE
## $ n_rows : int 1320
## $ n_studies : int 107
## $ n_hcd_rows : int 1272
## $ n_cc_rows : int 12
## $ n_index_rows : int 48
## $ n_hcd_cc_rows : int 1284
Step 3: Preliminary investigation of the HCD
Step 3 investigates the historical control data before they are used in later
parts of the workflow. For the count-data case study, this step has two main
aims.
First, we check whether any historical control studies have unusually high or
low mean values for the endpoint. The endpoint is the number of eggs laid per
hen during the 10-week observation period, so the study-level mean is used as a
simple summary for identifying potentially unusual studies.
Second, we investigate whether the endpoint is associated with available
covariates. In this case study, the relevant available covariate is the study
date. If there is evidence that egg counts change over time, this may affect
which historical studies should be carried forward or whether the covariate
should be included in later modelling.
The prepared HCD object count_hcd, created in the previous section, is used as
the input for this step.
Investigation of potential outlying studies
We first investigate whether any historical control studies show unusually high
or low mean egg counts compared with the remaining HCD studies. The function
explore_count_hcd_outliers() calculates study-level summary statistics,
evaluates the approximate normality of the study means using a Shapiro-Wilk
test, and identifies studies outside the boxplot fences based on the
interquartile range (IQR) rule.
The function also produces a graphical summary of the study means. In the plot,
the dashed vertical lines represent the lower and upper boxplot fences, while
studies outside these limits are highlighted.
count_outliers <- explore_count_hcd_outliers(
data = count_hcd,
response_col = "neggs",
study_id_col = "study",
aggregate_by = "study",
make_plots = TRUE
)
count_outliers$study_mean_shapiro_test
##
## Shapiro-Wilk normality test
##
## data: response_mean
## W = 0.99294, p-value = 0.8629
count_outliers$study_fence_table
## lower_fence upper_fence box_coef
## 1 24.875 55.875 1.5
count_outliers$study_outlier_table
## study_id response_mean outside_fence
## 47 45 24.66667 TRUE
## 77 72 23.91667 TRUE
count_outliers$plots$study_mean_outlier
Investigation of association with study date
We next investigate whether the historical control endpoint is associated with
study date. In this case study, the date covariate is represented by mdelta, a
centred version of the study date used to improve numerical stability of the
model fitting.
The function screen_count_hcd_covariate() fits a count-data mixed model with
and without the date covariate and compares the two models. The function
supports several model families and flexible random-effect structures, allowing
the modelling strategy to be adapted to the characteristics of the endpoint and
study design. The choice of model and random effects should therefore be guided
by the structure of the data, the amount of dispersion, and the scientific
context.
Here, we use a Poisson GLMM with identity link, a study-level random effect, and
an additional bird-level random effect to account for within-study variability.
This follows the model structure used in the original count-data case study.
The function returns the model comparison, a descriptive summary plot of study
means over time, a model-based association plot, and a simple scatter plot of
the observed individual counts.
count_date_screen <- screen_count_hcd_covariate(
data = count_hcd,
response_col = "neggs",
covariate_col = "mdelta",
covariate_type = "continuous",
study_id_col = "study",
random_effects = c("study", "bird_id"),
model_type = "glmer_poisson",
link = "identity",
standardize = FALSE,
make_plot = TRUE,
summary_x_col = "date",
summary_scale = 70,
summary_y_label = "No. eggs laid / hen / day"
)
count_date_screen$model_comparison
## Data: dat
## Models:
## reduced_model: neggs ~ 1 + (1 | study) + (1 | bird_id)
## model: neggs ~ 1 + mdelta + (1 | study) + (1 | bird_id)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## reduced_model 3 9888.2 9903.7 -4941.1 9882.2
## model 4 9890.2 9910.8 -4941.1 9882.2 0.0079 1 0.9291
count_date_screen$plots$summary
count_date_screen$plots$association
count_date_screen$plots$scatter
The graphical exploration does not suggest any clear association between the
endpoint and study date. The study-level summaries remain relatively stable over
time, with no obvious trend or substantial change in variability. The
scatterplot of individual observations similarly does not indicate any
systematic pattern with study date. Formal testing using the GLMM also found no
statistically significant association between the endpoint and study date (p =
0.93). Therefore, there was no evidence that study date needed to be accounted
for in the subsequent steps of the workflow.
Step 4: Modelling HCD variability
The goal of Step 4 is to estimate the variability represented by the historical
control data. For this count endpoint, the variability has two main components:
variability between historical control studies and variability between
individual hens within the same study.
The function model_count_hcd_variability() fits a mixed model to the
individual HCD observations and extracts the estimated between-study
distribution. This distribution summarises the expected control population mean
for a new comparable study and is used in the later comparison between the HCD
and the concurrent control.
Here, we use a Poisson GLMM with identity link, a study-level random effect, and
an additional bird-level random effect. The argument
between_study_re = "study" specifies which random-effect component should be
interpreted as the between-study variability.
count_hcd_model <- model_count_hcd_variability(
data = count_hcd,
response_col = "neggs",
study_id_col = "study",
random_effects = c("study", "bird_id"),
between_study_re = "study",
model_type = "glmer_poisson",
link = "identity",
make_plots = TRUE
)
count_hcd_model$hcd_distribution
## component mean sd lower upper link model_type
## 1 between_study 39.82849 5.455057 29.13678 50.52021 identity glmer_poisson
## between_study_re
## 1 study
count_hcd_model$random_effects
## group sd_intercept variance_intercept
## 1 bird_id 9.230139 85.19547
## 2 study 5.455057 29.75765
count_hcd_model$plots$between_study_distribution
The fitted GLMM estimates the average number of eggs per hen across the selected
HCD together with the between-study variability. The estimated between-study
distribution is approximately centred around 40 eggs per hen, with moderate
variability across studies. This distribution is subsequently used in Step 5 to
derive prediction intervals for future comparable control groups and to assess
compatibility between the HCD and the concurrent control.
Step 5 - Comparison of Concurrent Control to HCD
In Step 5, the concurrent control group from the index study is compared with
the prediction distribution derived from the HCD model fitted in Step 4.
For this case study, we use Approach 2, implemented as
"new_control_group_total". This approach compares the observed total number of
eggs in the concurrent-control group with the predictive distribution for the
total number of eggs in a new comparable control group of the same size.
The Step 4 model contains two random-effect components: a study-level random
effect and a bird-level random effect. Both sources of variability are included
in the predictive distribution. The study-level effect is sampled once for each
simulated control group and is shared by all birds in that group, whereas the
bird-level effect is sampled separately for each bird.
Because this model contains two random-effect levels, these settings are
specified explicitly. This overrides the simpler default used for models that
contain only a between-study random effect.
count_step5 <- compare_count_cc_to_hcd(
cc_data = count_cc,
step4_result = count_hcd_model,
response_col = "neggs",
approach = "new_control_group_total",
prediction_random_effects = c("study", "bird_id"),
group_level_random_effects = "study",
N_MC = 100000,
seed = 123,
make_plot = TRUE
)
count_step5$result$prediction_interval
## lower upper
## 330 628
count_step5$result$prediction_interval_raw
## [1] 329.3231 627.9802
count_step5$result$observed_total
## [1] 546
count_step5$result$inside_interval
## [1] TRUE
count_step5$result$mass_inside_interval
## [1] 0.9502684
count_step5$plots$predictive_distribution
The observed concurrent-control total is 546 eggs, which lies within the
HCD-based prediction interval of 330 to 628 eggs. This indicates that the
concurrent control is compatible with the selected HCD under the pre-specified
Step 5 comparison. Therefore, the workflow can proceed to Step 6.
Optional illustration: Approach 1
The EFSA workflow requires the Step 5 approach to be selected before the
analysis is carried out. The count-data case study uses Approach 2, implemented
above as "new_control_group_total". The code below is included only to
illustrate how the package can also run Approach 1.
Approach 1 compares the concurrent-control mean and its confidence interval with
the HCD prediction interval for the expected control mean in a future comparable
study.
count_step5_approach1 <- compare_count_cc_to_hcd(
cc_data = count_cc,
step4_result = count_hcd_model,
response_col = "neggs",
approach = "population_mean",
sigma_source = "estimate",
make_plot = TRUE
)
count_step5_approach1$result$hcd_prediction_interval
## [1] 29.13678 50.52021
count_step5_approach1$result$cc_mean
## [1] 45.5
count_step5_approach1$result$cc_confidence_interval
## [1] 39.37871 51.62129
count_step5_approach1$plots$comparison
Optional Step 5 sensitivity check
The function investigate_count_cc_hcd_mismatch() provides a descriptive
sensitivity check within one selected Step 5 approach. It does not compare
different Step 5 approaches and does not make regulatory or biological
conclusions.
For approach = "population_mean", the function varies the HCD prediction
interval level and the concurrent-control confidence interval level. It reports
the two intervals, their overlap, any gap between them, and their factual
relationship.
Here, the function is illustrated using Approach 1 only. This is not part of the
main EFSA count-data case-study decision path, but demonstrates how sensitivity
settings can be examined for a selected approach.
count_step5b_approach1 <- investigate_count_cc_hcd_mismatch(
cc_data = count_cc,
step4_result = count_hcd_model,
response_col = "neggs",
approach = "population_mean",
pi_levels = c(0.90, 0.95, 0.99),
cc_ci_levels = c(0.90, 0.95),
sigma_sources = "estimate",
make_plot = TRUE
)
count_step5b_approach1$sensitivity_table
## pi_level cc_ci_level sigma_source sigma_study observed_mean hcd_lower
## 1 0.90 0.90 estimate 5.455057 45.5 30.85572
## 2 0.90 0.95 estimate 5.455057 45.5 30.85572
## 3 0.95 0.90 estimate 5.455057 45.5 29.13678
## 4 0.95 0.95 estimate 5.455057 45.5 29.13678
## 5 0.99 0.90 estimate 5.455057 45.5 25.77719
## 6 0.99 0.95 estimate 5.455057 45.5 25.77719
## hcd_upper cc_lower cc_upper interval_overlap gap_outside interval_relation
## 1 48.80126 40.50536 50.49464 8.295904 0 intervals_overlap
## 2 48.80126 39.37871 51.62129 9.422553 0 intervals_overlap
## 3 50.52021 40.50536 50.49464 9.989284 0 cc_ci_inside_hcd_pi
## 4 50.52021 39.37871 51.62129 11.141498 0 intervals_overlap
## 5 53.87979 40.50536 50.49464 9.989284 0 cc_ci_inside_hcd_pi
## 6 53.87979 39.37871 51.62129 12.242582 0 cc_ci_inside_hcd_pi
count_step5b_approach1$plot
Step 6: Analysis of the current study including HCD
Step 6 is carried out because the concurrent control was found to be compatible
with the HCD in Step 5. The aim is to analyse the index study while
incorporating information from the HCD about the expected control response.
For the count-data case study, this is done with
analyze_count_index_with_hcd(). The function fits one or more user-specified
models to the index study alone and/or to the combined index-study and HCD data.
The model specification is supplied through model_specs, where the user can
choose the data source, fixed-effect predictors, random-effect structure, model
type, link function, and model-control options.
In this case study, the main model combines HCD with the index study and fits a
Poisson GLMM with identity link. The HCD records are treated as historical
control observations with dose zero. The HCD contribute information about the
expected control response and its variability, while the treatment effect is
driven by the contrast between the control and treated groups in the index
study.
The dose is scaled as dose100 = dose / 100 for the combined model, following
the original case-study code, to improve numerical stability. The coefficient
for dose100 can be interpreted as the change in the expected number of eggs
per hen for each 100 mg/day increase in dose.
For illustrative purposes, we also fit the index-study model without HCD. This
is included to reproduce and document the published case-study results and to
show the corresponding package output. In an actual application of the EFSA
workflow, the decision to analyse the index study with or without HCD should be
based on the conclusions from the previous steps, rather than on a comparison of
both analyses.
count_index_all <- count_prep$data[
count_prep$data$index_study == "Yes",
]
count_index_all$dose100 <- count_index_all$dose / 100
count_hcd_step6 <- count_hcd
count_hcd_step6$dose100 <- 0
count_step6 <- analyze_count_index_with_hcd(
index_data = count_index_all,
hcd_data = count_hcd_step6,
step4_result = count_hcd_model,
response_col = "neggs",
dose_col = "dose",
study_id_col = "study",
control_dose = 0,
model_specs = list(
with_hcd = list(
name = "with_hcd",
data_source = "combined",
predictor_cols = "dose100",
random_effect_cols = c("study", "study:dose:bird"),
model_type = "glmer_poisson",
link = "identity",
glmer_control = lme4::glmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)
)
),
without_hcd = list(
name = "without_hcd",
data_source = "index",
predictor_cols = "dose",
random_effect_cols = "dose:bird",
model_type = "glmer_poisson",
link = "identity"
)
),
make_plots = TRUE
)
count_step6$formatted_coefficient_tables$with_hcd
## Term Estimate Standard error p-value
## 1 (Intercept) 39.8796 0.6086 <0.001
## 2 dose100 -1.8486 1.0191 0.070
count_step6$formatted_coefficient_tables$without_hcd
## Term Estimate Standard error p-value
## 1 (Intercept) 41.8248 2.3301 <0.001
## 2 dose -0.0197 0.0100 0.049
## Source-comparison views
count_step6$plots$comparison_index_study_only
count_step6$plots$comparison_all_sources
Both count source-comparison views retain all active-dose groups from the index
study. The index-study-only view shows the concurrent control without HCD. The
all-sources view adds the HCD and pooled-control information. The previous
count_step6$plots$comparison_dose_response object is retained as an alias for
comparison_all_sources, so existing code continues to work.
The fitted model including HCD gives a negative dose coefficient and a trend
test p-value of 0.070, indicating a decreasing trend in the number of eggs with
increasing dose. This reproduces the conclusion of the EFSA Annex example: after
the concurrent control was found to be consistent with the HCD, the index study
was analysed including HCD, and the evidence for a dose-related trend was not
statistically significant at the conventional 5% level.
