efsaHCD: Implementing the EFSA Workflow for Historical Control Data
Extra Case Studies – efsaHCD version 0.0.7
Open Analytics
15 July 2026
1 Introduction
This vignette presents three additional case studies used to illustrate the
application of the efsaHCD package for the analysis of historical control data
(HCD).
The aim is to show how the package functions can be used across the main EFSA workflow steps, from the preliminary investigation of HCD to the comparison with the concurrent control and, where appropriate, the analysis of the index study including HCD.
The case studies are based on the BfR report on testing the EFSA approach for
the use and reporting of HCD in regulatory studies. In this vignette, the
analyses are revisited using the current version of the efsaHCD package.
2 Dichotomous endpoint case study: C-cell Carcinomas in F344 Rats
This case study investigates the occurrence of thyroid C-cell carcinomas in male and female F344 rats from a combined chronic toxicity and carcinogenicity study. The endpoint represents dichotomous incidence data and was selected as an example of a relatively common tumour finding in long-term rodent studies.
The objective is to illustrate how the efsaHCD package can be used to support
the EFSA workflow for historical control data (HCD), including exploratory
investigation of HCD, modelling of between-study variability, comparison of the
concurrent control with the HCD, and statistical analysis of the index study
incorporating HCD.
2.1 Step 1: Planning Phase
Step 1 defines the overall strategy for the HCD analysis. This includes specifying the endpoint of interest, identifying potentially relevant covariates, defining comparability criteria for selecting HCD studies, and planning the statistical analyses to be performed in later workflow steps.
For this case study, the endpoint of primary interest was thyroid C-cell carcinoma in F344 rats from 2-year carcinogenicity studies. Considerations included rat strain, study type, study period, laboratory, breeding facility, and possible covariates such as study year, diet, and study facility.
The efsaHCD package is mainly intended for the statistical workflow beginning
at Step 3, after suitable HCD have been collected and prepared.
2.2 Step 2: Data Acceptability
Step 2 evaluates whether the available HCD are sufficiently comparable to the index study and whether they are suitable for statistical analysis.
In the original BfR case study, the HCD initially submitted with the study were considered insufficient because important information was missing or not comparable to the index study, including differences in laboratory, study period, and animal source. Consequently, additional HCD studies were identified and collected before proceeding to the statistical workflow.
At this stage, the focus is primarily on scientific and regulatory assessment of HCD relevance and quality rather than formal statistical modelling.
2.3 Step 3 - Preliminary Investigation of HCD
Step 3 is the first step where the efsaHCD package is used directly. The aim
is to inspect the historical control data before modelling, identify possible
outliers, and assess whether the endpoint is associated with relevant
covariates.
For this binary case study, the data are grouped incidence data. Each row represents one study group and contains the number of animals examined and the number of animals with the finding. The index study and historical control data are stored in separate Excel files included with the package.
library(readxl)
file_hcd <- system.file(
"extdata",
"input_hcd_case study 1.xlsx",
package = "efsaHCD"
)
file_index <- system.file(
"extdata",
"input_index study_case study 1.xlsx",
package = "efsaHCD"
)
full_hcd <- readxl::read_xlsx(file_hcd)
index_study <- readxl::read_xlsx(file_index)
full_hcd$n <- full_hcd$`No. examined`
full_hcd$k <- full_hcd$`No. with finding`
index_study$n <- index_study$`No. examined`
index_study$k <- index_study$`No. with finding_carcinoma`
full_hcd$year <- format(
as.Date(full_hcd$study_period_end_date),
"%Y"
)
full_hcd$date <- as.Date(full_hcd$date, format = "%d.%m.%Y")
binary_hcd <- full_hcd[
full_hcd$index_study == "no" &
full_hcd$parameter != "C-cell adenoma",
]
# making sure variables are numeric
binary_hcd$k <- as.numeric(binary_hcd$k)
binary_hcd$n <- as.numeric(binary_hcd$n)
binary_hcd$incidence <- binary_hcd$k / binary_hcd$n2.3.1 Step 3a - Exploratory Visual Inspection
We first inspect the historical control incidences graphically. The function
make_binary_hcd_plots() creates a boxplot of the observed incidences and
scatterplots against the selected covariates.
Because grouping_var = "sex" is specified, the function returns the results
separately for females and males. This is consistent with the case-study
analysis, where both sexes are analysed separately throughout the workflow.
plots_all <- make_binary_hcd_plots(
dat = binary_hcd,
k_col = "k",
n_col = "n",
covariate_cols = c("diet", "facility", "year"),
study_id_col = NULL,
grouping_var = "sex",
dataset_label = "HCD",
point_shape_by_k = TRUE,
save_plots = FALSE
)
plots_all$female$boxplot
plots_all$male$boxplot
plots_all$female$scatterplots$year
plots_all$male$scatterplots$year
The boxplots summarise the distribution of observed historical control incidences. The plots against year are shown here because study year was an important covariate in the original case study and was later used to restrict the HCD dataset for the main analysis.
2.3.2 Step 3a - Model-based Covariate Screening
The visual inspection is followed by a model-based screening of the main
covariates. The function screen_binary_hcd_covariates() fits separate
univariate models for each covariate. Since grouping_var = "sex" is used, the
screening is performed separately for females and males.
For each covariate, both binomial and quasi-binomial models are fitted. The quasi-binomial model is useful when there is evidence of overdispersion.
step3a_res <- screen_binary_hcd_covariates(
dat = binary_hcd,
k_col = "k",
n_col = "n",
covariate_cols = c("diet", "facility", "year"),
grouping_var = "sex",
dispersion_threshold = 1
)
step3a_res$female$formatted_summary_table## Covariate Preferred p-value
## 1 diet 0.201
## 2 facility 0.125
## 3 year 0.001step3a_res$male$formatted_summary_table## Covariate Preferred p-value
## 1 diet 0.240
## 2 facility 0.243
## 3 year 0.0002.3.3 Step 3b - Selection of HCD for Further Analysis
Based on the exploratory plots and the model-based covariate screening, study year appears to be an important source of heterogeneity in the full HCD dataset. The plots show that more recent studies tend to have higher observed incidences of C-cell carcinomas than earlier studies.
This pattern is also supported by the covariate screening results. For females, the association between incidence and year was statistically significant (binomial GLM p-value = 0.0005). By contrast, diet and facility were not considered statistically significant after accounting for overdispersion using the quasi-binomial models.
Therefore, before proceeding to Step 4, the HCD dataset is restricted to studies closer to the index-study period. This defines the selected HCD dataset used for the main modelling analysis.
2.4 Step 4 - Modelling HCD Variability
In Step 4, the selected HCD dataset from Step 3b is used to estimate the between-study variability in the control incidence. For binary incidence data, this is done using a beta-binomial model.
Based on the Step 3 results, the main analysis is restricted to historical control studies with study end year up to 2001. This keeps studies closer to the index-study period and avoids combining earlier and later studies that appear to come from different background-incidence distributions.
hcd_step4 <- binary_hcd[binary_hcd$year <= 2001, ]
step4_res <- model_binary_hcd_variability(
dat = hcd_step4,
k_col = "k",
n_col = "n",
grouping_var = "sex",
covariate_cols = NULL,
reference_phi = 0,
alpha = 0.05,
root_interval = NULL,
plot_profile = TRUE,
plot_distributions = TRUE,
save_plots = FALSE
)
step4_res$female$summary_table## n_rows_used alpha phi_mle phi_lcl phi_ucl neff_mle neff_lcl neff_ucl
## 1 8 0.05 6.61069e-05 0 0.04534557 15126.02 21.05287 Infstep4_res$male$summary_table## n_rows_used alpha phi_mle phi_lcl phi_ucl neff_mle neff_lcl neff_ucl
## 1 8 0.05 6.61069e-05 0 0.02054816 15126.02 47.66616 InfThe model estimates the average background incidence and the between-study variability parameter. The upper confidence limit for the between-study variability is then used to obtain a conservative estimate of the effective sample size and to define the beta distribution used in later steps.
The grey curve represents the HCD-based distribution obtained using a very large effective sample size, corresponding to negligible between-study variability, whereas the black curve represents the estimated between-study variability distribution from the selected HCD.
step4_res$female$distribution_plot
step4_res$male$distribution_plot
The profile likelihood plots show how the likelihood changes across possible values of the between-study variability parameter. These plots are useful for checking whether the upper confidence limit was obtained in a stable way.
step4_res$female$profile_plot
step4_res$male$profile_plot
2.5 Step 5 - Comparison of Concurrent Control to HCD
In Step 5, the concurrent control group from the index study is compared with the selected HCD model from Step 4. The aim is to assess whether the concurrent control is compatible with the historical control distribution before using HCD in the analysis of the index study.
In practice, the Step 5 approach should be selected in the analysis protocol before seeing the results. In this vignette, we use the new-group-outcome approach (approach 2 in the opinion paper). This approach compares the observed number of animals with the finding in the concurrent control group with the predictive distribution for a new control group of the same size.
if (!("n" %in% names(index_study))) {
index_study$n <- index_study[["No. examined"]]
}
if (!("k" %in% names(index_study))) {
index_study$k <- index_study[["No. with finding_carcinoma"]]
}
cc_dat <- index_study[index_study$Dose == 0, ]
cc_dat$sex <- ifelse(
toupper(cc_dat$Sex) == "F", "female",
ifelse(toupper(cc_dat$Sex) == "M", "male", NA)
)
cc_dat <- cc_dat[, c("sex", "k", "n", "Dose")]The group labels in the concurrent control data must match the group labels in the Step 4 result. Here, the sex labels are therefore harmonised to “female” and “male” before running the comparison.
step5_res <- compare_binary_cc_to_hcd(
cc_dat = cc_dat,
step4_result = step4_res,
k_col = "k",
n_col = "n",
grouping_var = "sex",
approach = "new_group_outcome",
prediction_basis = "upper_limit",
cc_conf_level = 0.95,
pi_conf_level = 0.95,
ci_method = "wilson",
plot_population_incidence = FALSE,
plot_new_group_outcome = TRUE,
save_plots = FALSE
)
step5_res$female$formatted_summary_table## Group Approach CC number with finding CC number examined
## female new_group_outcome 1.000 50.000
## CC incidence Prediction basis Effective sample size used
## 0.020 upper_limit 21.053
## HCD incidence used for prediction HCD PI lower number of cases
## 0.013 0.000
## HCD PI upper number of cases Left-tail probability Right-tail probability
## 4.000 0.000 0.031
## CC within HCD prediction interval
## TRUEstep5_res$male$formatted_summary_table## Group Approach CC number with finding CC number examined CC incidence
## male new_group_outcome 0.000 50.000 0.000
## Prediction basis Effective sample size used HCD incidence used for prediction
## upper_limit 47.666 0.027
## HCD PI lower number of cases HCD PI upper number of cases
## 0.000 5.000
## Left-tail probability Right-tail probability CC within HCD prediction interval
## 0.000 0.027 TRUEFor females, the observed concurrent-control count was 1 out of 50 animals. This falls within the HCD predictive interval of 0 to 4 animals. For males, the observed concurrent-control count was 0 out of 50 animals. This falls within the HCD predictive interval of 0 to 5 animals. Based on the selected new-group-outcome approach, the concurrent control groups are therefore compatible with the selected HCD dataset.
The same comparison can also be inspected visually.
In this plot, the dots show the predictive probabilities for different possible numbers of animals with the finding in a new control group of the same size as the concurrent control group.
The horizontal black bar shows the HCD prediction interval. The solid vertical line shows the observed number of animals with the finding in the concurrent control group.
The dashed vertical line marks the boundary beyond the upper end of the prediction interval. In this example, counts from 0 to 4 animals are considered within the HCD prediction interval, whereas counts of 5 or more animals would fall outside the interval.
Since the observed concurrent-control count lies within the prediction interval, the concurrent control is visually compatible with the selected HCD dataset.
2.5.1 Step 5: Illustration Using the Population-Incidence Approach
The Step 5 approach should normally be selected in the analysis protocol before the results are inspected. In this vignette, the new-group-outcome approach (Approach 2 in the opinion paper) was selected for the main analysis. However, for illustration, we also apply the population-incidence approach (Approach 1 in the opinion paper).
This approach compares the confidence interval for the concurrent-control incidence with the prediction interval for the HCD population incidence. We use a 95% confidence interval for the concurrent control and a 95% prediction interval for the HCD.
step5_res_pop <- compare_binary_cc_to_hcd(
cc_dat = cc_dat,
step4_result = step4_res,
k_col = "k",
n_col = "n",
grouping_var = "sex",
approach = "population_incidence",
prediction_basis = "upper_limit",
cc_conf_level = 0.95,
pi_conf_level = 0.95,
ci_method = "wilson",
plot_population_incidence = TRUE,
plot_new_group_outcome = FALSE,
save_plots = FALSE
)
step5_res_pop$female$formatted_summary_table## Group Approach CC number with finding CC number examined
## female population_incidence 1.000 50.000
## CC incidence Prediction basis Effective sample size used
## 0.020 upper_limit 21.053
## HCD incidence used for prediction CC CI lower limit CC CI upper limit
## 0.013 0.004 0.105
## HCD PI lower limit HCD PI upper limit CC within HCD prediction interval
## 0.000 0.085 FALSEstep5_res_pop$male$formatted_summary_table## Group Approach CC number with finding CC number examined
## male population_incidence 0.000 50.000
## CC incidence Prediction basis Effective sample size used
## 0.000 upper_limit 47.666
## HCD incidence used for prediction CC CI lower limit CC CI upper limit
## 0.027 0.000 0.071
## HCD PI lower limit HCD PI upper limit CC within HCD prediction interval
## 0.001 0.087 FALSEThe corresponding plots show the HCD prediction interval for the population incidence together with the confidence interval for the concurrent-control incidence.
For females, the 95% confidence interval for the concurrent-control incidence was 0.004 to 0.105, while the 95% prediction interval from the HCD model ranged from approximately 0.000 to 0.085. Thus, the upper part of the concurrent-control confidence interval extended slightly beyond the HCD prediction interval.
For males, the 95% confidence interval for the concurrent-control incidence was 0.000 to 0.071, whereas the HCD prediction interval ranged from approximately 0.001 to 0.087. Here, the intervals showed only minimal non-overlap near zero.
Overall, compared with the new-group-outcome approach, the population-incidence approach leads to a stricter assessment of compatibility between the concurrent control and the HCD distribution.
2.5.2 Step 5b: Investigation of Apparent Inconsistencies Between HCD and Concurrent Control
Under the population-incidence approach, the confidence intervals for the concurrent controls and the prediction intervals from the HCD model did not fully overlap. According to the EFSA workflow, this motivates a more detailed investigation of the apparent inconsistency between the concurrent control and the HCD before deciding how to proceed further.
Step 5b is intended as an investigative step. Possible outcomes include: revisiting the HCD selection, analysing the index study without HCD, excluding the index study from further consideration, or accepting the mismatch and proceeding to Step 6 with appropriate justification.
Here, we use the exploratory Step 5b investigation function to examine the source and magnitude of the apparent mismatch between the concurrent control and the HCD prediction intervals.
step5b_res <- investigate_binary_cc_hcd_mismatch(
cc_dat = cc_dat,
step4_result = step4_res,
k_col = "k",
n_col = "n",
grouping_var = "sex",
approach = "population_incidence",
prediction_basis_main = "upper_limit",
cc_conf_level_main = 0.95,
pi_conf_level_main = 0.95,
ci_method_main = "wilson",
make_plot = TRUE,
save_plots = FALSE
)
step5b_res$female$formatted_main_diagnostics## Group Approach Prediction basis CC confidence level
## female population_incidence upper_limit 0.950
## HCD PI confidence level CC CI method CC CI lower limit CC CI upper limit
## 0.950 wilson 0.004 0.105
## HCD PI lower limit HCD PI upper limit Lower gap Upper gap Maximum gap
## 0.000 0.085 0.000 0.020 0.020
## CC interval width HCD interval width Overlap width
## 0.101 0.085 0.082
## Overlap proportion of CC interval Strictly within Borderline flag
## 0.807 FALSE FALSE
## Tolerance used Status
## 0.004 investigate_furtherstep5b_res$male$formatted_main_diagnostics## Group Approach Prediction basis CC confidence level
## male population_incidence upper_limit 0.950
## HCD PI confidence level CC CI method CC CI lower limit CC CI upper limit
## 0.950 wilson 0.000 0.071
## HCD PI lower limit HCD PI upper limit Lower gap Upper gap Maximum gap
## 0.001 0.087 0.001 0.000 0.001
## CC interval width HCD interval width Overlap width
## 0.071 0.086 0.070
## Overlap proportion of CC interval Strictly within Borderline flag
## 0.981 FALSE TRUE
## Tolerance used Status
## 0.004 borderline_reviewThe Step 5b function re-runs the selected Step 5 comparison and adds diagnostic quantities that help interpret the apparent mismatch. The plot has the same structure as the Step 5 population-incidence plot; the additional information is provided mainly by the diagnostic table.
For females, the upper bound of the concurrent-control confidence interval extends beyond the upper bound of the HCD prediction interval by 0.020. The overlap still covers about 81% of the concurrent-control interval, but the result is classified as requiring further investigation.
For males, the mismatch occurs only at the lower boundary: the HCD prediction interval starts at 0.001, while the concurrent-control confidence interval starts at 0.000. The gap is very small and the overlap covers about 98% of the concurrent-control interval. This is therefore classified as a borderline case rather than a clear incompatibility.
The Step 5b function also allows a small sensitivity investigation by repeating the comparison under different methodological choices, including:
- different prediction bases from Step 4,
- different confidence levels,
- and different confidence-interval methods for the concurrent control.
step5b_res$female$formatted_sensitivity_table## Group Approach Prediction basis CC CI method
## female population_incidence lower_limit wilson
## female population_incidence lower_limit wilson
## female population_incidence lower_limit wilson
## female population_incidence lower_limit wilson
## female population_incidence lower_limit clopper-pearson
## female population_incidence lower_limit clopper-pearson
## female population_incidence lower_limit clopper-pearson
## female population_incidence lower_limit clopper-pearson
## female population_incidence lower_limit agresti-coull
## female population_incidence lower_limit agresti-coull
## female population_incidence lower_limit agresti-coull
## female population_incidence lower_limit agresti-coull
## female population_incidence point_estimate wilson
## female population_incidence point_estimate wilson
## female population_incidence point_estimate wilson
## female population_incidence point_estimate wilson
## female population_incidence point_estimate clopper-pearson
## female population_incidence point_estimate clopper-pearson
## female population_incidence point_estimate clopper-pearson
## female population_incidence point_estimate clopper-pearson
## female population_incidence point_estimate agresti-coull
## female population_incidence point_estimate agresti-coull
## female population_incidence point_estimate agresti-coull
## female population_incidence point_estimate agresti-coull
## female population_incidence upper_limit wilson
## female population_incidence upper_limit wilson
## female population_incidence upper_limit wilson
## female population_incidence upper_limit wilson
## female population_incidence upper_limit clopper-pearson
## female population_incidence upper_limit clopper-pearson
## female population_incidence upper_limit clopper-pearson
## female population_incidence upper_limit clopper-pearson
## female population_incidence upper_limit agresti-coull
## female population_incidence upper_limit agresti-coull
## female population_incidence upper_limit agresti-coull
## female population_incidence upper_limit agresti-coull
## CC confidence level HCD PI confidence level CC CI lower limit
## 0.900 0.900 0.004
## 0.900 0.950 0.004
## 0.950 0.900 0.004
## 0.950 0.950 0.004
## 0.900 0.900 0.001
## 0.900 0.950 0.001
## 0.950 0.900 0.001
## 0.950 0.950 0.001
## 0.900 0.900 0.000
## 0.900 0.950 0.000
## 0.950 0.900 0.000
## 0.950 0.950 0.000
## 0.900 0.900 0.004
## 0.900 0.950 0.004
## 0.950 0.900 0.004
## 0.950 0.950 0.004
## 0.900 0.900 0.001
## 0.900 0.950 0.001
## 0.950 0.900 0.001
## 0.950 0.950 0.001
## 0.900 0.900 0.000
## 0.900 0.950 0.000
## 0.950 0.900 0.000
## 0.950 0.950 0.000
## 0.900 0.900 0.004
## 0.900 0.950 0.004
## 0.950 0.900 0.004
## 0.950 0.950 0.004
## 0.900 0.900 0.001
## 0.900 0.950 0.001
## 0.950 0.900 0.001
## 0.950 0.950 0.001
## 0.900 0.900 0.000
## 0.900 0.950 0.000
## 0.950 0.900 0.000
## 0.950 0.950 0.000
## CC CI upper limit HCD PI lower limit HCD PI upper limit Strictly within
## 0.085 0.006 0.007 FALSE
## 0.085 0.006 0.007 FALSE
## 0.105 0.006 0.007 FALSE
## 0.105 0.006 0.007 FALSE
## 0.091 0.006 0.007 FALSE
## 0.091 0.006 0.007 FALSE
## 0.106 0.006 0.007 FALSE
## 0.106 0.006 0.007 FALSE
## 0.091 0.006 0.007 FALSE
## 0.091 0.006 0.007 FALSE
## 0.115 0.006 0.007 FALSE
## 0.115 0.006 0.007 FALSE
## 0.085 0.005 0.008 FALSE
## 0.085 0.005 0.008 FALSE
## 0.105 0.005 0.008 FALSE
## 0.105 0.005 0.008 FALSE
## 0.091 0.005 0.008 FALSE
## 0.091 0.005 0.008 FALSE
## 0.106 0.005 0.008 FALSE
## 0.106 0.005 0.008 FALSE
## 0.091 0.005 0.008 FALSE
## 0.091 0.005 0.008 FALSE
## 0.115 0.005 0.008 FALSE
## 0.115 0.005 0.008 FALSE
## 0.085 0.000 0.062 FALSE
## 0.085 0.000 0.085 TRUE
## 0.105 0.000 0.062 FALSE
## 0.105 0.000 0.085 FALSE
## 0.091 0.000 0.062 FALSE
## 0.091 0.000 0.085 FALSE
## 0.106 0.000 0.062 FALSE
## 0.106 0.000 0.085 FALSE
## 0.091 0.000 0.062 FALSE
## 0.091 0.000 0.085 FALSE
## 0.115 0.000 0.062 FALSE
## 0.115 0.000 0.085 FALSE
## Borderline flag Status
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE within
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_further
## FALSE investigate_furtherThis illustrates that the conclusion from the population-incidence approach can be sensitive to methodological choices and interval definitions, supporting the need for careful interpretation and further investigation rather than a purely automatic decision rule.
2.6 Step 6: Analysis of the Index Study Including HCD
Since the concurrent-control groups were considered sufficiently compatible with the selected HCD dataset under the new-group-outcome approach, the HCD can now be incorporated into the analysis of the index study.
In Step 6, the HCD contribute additional information about the expected control response and its variability. The first table shows the index-study analysis without HCD. The second table shows the corresponding analysis after incorporating the HCD pseudo-control.
Before running the analysis, the grouping variable is harmonised to match the labels used in the Step 4 HCD model.
index_study$sex <- ifelse(
toupper(index_study$Sex) == "F", "female",
ifelse(
toupper(index_study$Sex) == "M", "male",
NA
)
)The function analyze_binary_index_with_hcd() performs the Step 6 analysis
separately for each level of the grouping variable. Here, hcd_basis = "upper_limit" is used, so the HCD contribution is based on the upper confidence
limit for the between-study variability estimated in Step 4.
step6_res <- analyze_binary_index_with_hcd(
index_dat = index_study,
step4_result = step4_res,
dose_col = "Dose",
k_col = "k",
n_col = "n",
grouping_var = "sex",
control_dose = 0,
hcd_basis = "upper_limit",
ci_method = "wilson",
ci_conf_level = 0.95,
make_plot = TRUE,
add_pooled_to_plot = TRUE,
save_plots = FALSE
)
step6_res$female$formatted_summary_table_without_hcd## Group Control dose Highest dose Number of CC cases CC number examined
## female 0.000 1600.000 1.000 50.000
## CA asymptotic p-value CA exact p-value Fisher control vs max p-value
## 0.4555 0.4229 0.7525
## CA exact control vs max p-value
## 0.7525step6_res$female$formatted_summary_table_with_hcd## Group HCD basis Control dose Highest dose Tarone asymptotic p-value
## female upper_limit 0.000 1600.000 0.4359
## HYT exact p-value HYT exact control vs max p-value
## 0.3375 0.4982step6_res$male$formatted_summary_table_without_hcd## Group Control dose Highest dose Number of CC cases CC number examined
## male 0.000 1600.000 0.000 50.000
## CA asymptotic p-value CA exact p-value Fisher control vs max p-value
## 0.1645 0.1779 0.2475
## CA exact control vs max p-value
## 0.2475step6_res$male$formatted_summary_table_with_hcd## Group HCD basis Control dose Highest dose Tarone asymptotic p-value
## male upper_limit 0.000 1600.000 0.2396
## HYT exact p-value HYT exact control vs max p-value
## 0.1538 0.1555The Step 6 plot shows the observed incidence in the index study at each dose, together with the HCD-derived pseudo-control and the pooled control estimate.
Important: The analyses shown with and without HCD are included only for illustrative purposes to demonstrate the corresponding Step 6 workflow and package outputs. In a real application of the EFSA workflow, the decision to include HCD in the analysis should be driven by the conclusions from the previous workflow steps.
2.7 Step 7 - Influence and Sensitivity Analysis
In this case study, Step 7 was used to assess the impact of the main HCD selection decision made in Step 3b. The main analysis used the restricted HCD dataset, keeping studies with end year up to 2001. As a sensitivity analysis, the BfR report repeated the relevant analyses using the full HCD dataset, without this year restriction.
In this vignette, we do not repeat the full Step 7 analysis because it would reuse the same functions already illustrated in Steps 4 to 6, but with the full HCD dataset instead of the restricted one.
More generally, Step 7 should be defined and interpreted case by case. Depending on the dataset and scientific question, sensitivity analyses may consider alternative HCD selections, alternative modelling assumptions, different interval definitions, or other statistical methods.
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.
3 Case Study 2: Continuous Summary Data
This case study illustrates the HCD workflow for a continuous endpoint using summary-level data. The endpoint is a developmental timing outcome. Unlike individual animal-level data, each row in this dataset represents a study or study-group summary.
The main response variables are the group mean, standard deviation, and sample size. Additional summary-level covariates are also available, including offspring body weight and parental body weight. Study-level covariates include sex, supplier, diet, facility, and study date.
The historical control data and the index study are stored in separate files.
The variable index_study identifies whether a record belongs to the HCD or the
index study, with "no" indicating historical control data and "yes"
indicating index-study records.
file_hcd_cont <- system.file(
"extdata",
"input_hcd_case study 2.csv",
package = "efsaHCD"
)
file_index_cont <- system.file(
"extdata",
"input_index study_case study 2.csv",
package = "efsaHCD"
)
orig_hcd <- read.csv(
file = file_hcd_cont,
header = TRUE,
stringsAsFactors = FALSE
)
orig_index <- read.csv(
file = file_index_cont,
header = TRUE,
stringsAsFactors = FALSE
)
str(orig_hcd)## 'data.frame': 42 obs. of 25 variables:
## $ study_ID : chr "TOX2002-248" "ASB2013-9188" "ASB2017-6563" "ASB2019-9537" ...
## $ endpoint : chr "vaginal opening" "vaginal opening" "vaginal opening" "vaginal opening" ...
## $ sex : chr "female" "female" "female" "female" ...
## $ No_pups : int 24 27 20 26 23 20 19 24 30 30 ...
## $ age_offspring_endpoint : num 32.4 32 34 31.7 33.8 34 35 35 31.2 31.6 ...
## $ sd_age_offspring_endpoint: num 1.52 0.9 2 1.2 2.01 2.3 2.2 3 1.9 1.7 ...
## $ bw_offspring_endpoint : num 105 109 118 112 130 ...
## $ sd_bw_offspring_endpoint : num 13 6.3 12 9.8 11.9 ...
## $ bw_offspring_3wPP : num 48.2 50.6 54.4 52.3 54.4 46 44.1 50.1 52.1 54.8 ...
## $ sd_bw_offspring_3wPP : num 4.8 4.7 7.4 5.3 4.97 4.7 5.6 5.3 6.6 6.8 ...
## $ No_parents : int 24 27 24 27 25 25 28 27 30 30 ...
## $ bw_parent_3wPP : num 338 335 299 338 308 ...
## $ sd_bw_parent_3wPP : num 20.6 19 20 19.1 20.9 16.5 27 30 23.5 15.1 ...
## $ index_study : chr "yes" "no" "no" "no" ...
## $ study_type : chr "OECD Guideline 416 (Two-Generation Reproduction Toxicity Study)" "EPA OPPTS 870.3800 (Reproduction and Fertility Effects)" "OECD Guideline 443 (Extended One-Generation Reproductive Toxicity Study)" "EPA OPPTS 870.3800 (Reproduction and Fertility Effects)" ...
## $ study : chr "TG416" "EPA870" "TG443" "EPA870" ...
## $ study_period_start_date : chr "04-Jun-98" "03-Jul-08" "31-Aug-10" "14-Jun-11" ...
## $ date : chr "04.06.1998" "03.07.2008" "31.08.2010" "14.06.2011" ...
## $ testing_facility_location: chr "Isengard Research Division" "Rohan Environmental Sciences" "Bree Laboratory Services" "Rohan Environmental Sciences" ...
## $ facility : chr "ISG" "RHN" "BRE" "RHN" ...
## $ type_of_diet_brand : chr "Purina Mills Rodent Lab Chow 5001-4" "LabDiet Certified Rodent Diet #5002" "Harlan Teklad Certified Rodent Diet 8728C" "LabDiet Certified Rodent Diet #5002" ...
## $ diet : chr "PMRLC" "CRD" "HTCRD" "LCRD" ...
## $ supplier : chr "Charles River Laboratories US" "Charles River Laboratories US" "Taconic Farms (Germantown. NY)" "Charles River Laboratories US" ...
## $ supplier_abb : chr "Crl_US" "Crl_US" "TF_US" "Crl_US" ...
## $ supplier_binned : chr "Crl" "Crl" "Other" "Crl" ...str(orig_index)## 'data.frame': 8 obs. of 22 variables:
## $ study_ID : chr "TOX2002-248" "TOX2002-248" "TOX2002-248" "TOX2002-248" ...
## $ endpoint : chr "preputial separation" "vaginal opening" "preputial separation" "vaginal opening" ...
## $ sex : chr "male" "female" "male" "female" ...
## $ No..of.animals_pups : int 23 24 29 29 27 28 27 27
## $ age_offspring_endpoint : num 41.2 32.4 41.9 32.2 42.5 32.1 47.9 34.7
## $ sd_age_offspring_endpoint: num 1.29 1.52 1.35 1.72 1.56 2.43 2.75 2.34
## $ bw_offspring_endpoint : num 190 105 199 104 192 ...
## $ sd_bw_offspring_endpoint : num 17.6 13 12.8 10.6 19.4 10.6 15.4 14.3
## $ bw_offspring_3wPP : num 49.8 48.2 49.6 46.8 47.2 45.5 37 35.8
## $ sd_bw_offspring_3wPP : num 4.46 4.8 3.55 3.39 5.77 4.39 3.74 4.16
## $ No_parents : int 30 30 30 30 30 30 30 30
## $ bw_parent_3wPP : num 338 338 333 333 330 ...
## $ sd_bw_parent_3wPP : num 20.6 20.6 22.7 22.7 18.2 ...
## $ index_study : chr "yes" "yes" "yes" "yes" ...
## $ study : chr "TG416" "TG416" "TG416" "TG416" ...
## $ study_period_start_date : chr "04-Jun-98" "04-Jun-98" "04-Jun-98" "04-Jun-98" ...
## $ date : chr "04.06.1998" "04.06.1998" "04.06.1998" "04.06.1998" ...
## $ facility : chr "BCAD" "BCAD" "BCAD" "BCAD" ...
## $ diet : chr "PMRLC" "PMRLC" "PMRLC" "PMRLC" ...
## $ supplier_abb : chr "Crl_US" "Crl_US" "Crl_US" "Crl_US" ...
## $ supplier_binned : chr "Crl" "Crl" "Crl" "Crl" ...
## $ dose : int 0 0 150 150 500 500 2500 25003.1 Steps 1 and 2: Planning and Data Acceptability
This case study concerns age at balano-preputial separation in Sprague Dawley rats, with body weight considered as a possible confounder. In the original report, the endpoint was selected as an example of continuous developmental data where the interpretation of HCD may depend on covariates such as body weight, supplier, diet, facility and study year.
In Step 1, the endpoint, target study type, strain, relevant covariates and planned statistical approach are defined. In Step 2, the available HCD are checked for biological and statistical acceptability before formal analysis. For this case study, the originally available HCD were limited, so new HCD were collected using criteria related to strain, study type, laboratory, breeding facility and study period.
The code below prepares the summary-level HCD data in the standard format used
by the continuous-data functions in efsaHCD.
hcd_prep <- prepare_continuous_summary_hcd(
data = orig_hcd,
endpoint_mean_col = "age_offspring_endpoint",
endpoint_sd_col = "sd_age_offspring_endpoint",
n_col = "No_pups",
date_col = "date",
date_format = "%d.%m.%Y",
endpoint_name = "age at preputial separation",
grouping_vars = c("endpoint", "sex"),
categorical_covariates = c(
"study_type",
"testing_facility_location",
"facility",
"type_of_diet_brand",
"diet",
"supplier",
"supplier_abb",
"supplier_binned"
),
continuous_covariates = c(
"bw_offspring_endpoint",
"bw_offspring_3wPP",
"bw_parent_3wPP"
),
hcd_indicator_col = "index_study",
hcd_indicator_value = "no",
create_log_variables = TRUE,
split_by_groups = TRUE
)##
## Prepared continuous summary-level dataset
##
## Rows: 42
## Endpoint mean column: age_offspring_endpoint
## Sample size column: No_pups
## SD column: sd_age_offspring_endpoint
## Grouping variables: endpoint, sex
## HCD rows: 40names(hcd_prep$grouped_data)## [1] "endpoint_vaginal_opening__sex_female"
## [2] "endpoint_preputial_separation__sex_male"head(hcd_prep$grouped_data$endpoint_preputial_separation__sex_male)## study_ID endpoint sex No_pups age_offspring_endpoint
## 22 TOX2002-248 preputial separation male 23 41.2
## 23 ASB2013-9188 preputial separation male 27 43.0
## 24 ASB2017-6563 preputial separation male 20 37.0
## 25 ASB2019-9537 preputial separation male 26 43.5
## 26 ASB2022-21072 preputial separation male 23 45.1
## 27 ASB2023-19175 preputial separation male 20 48.0
## sd_age_offspring_endpoint bw_offspring_endpoint sd_bw_offspring_endpoint
## 22 1.29 190.1 17.60
## 23 1.30 223.2 13.70
## 24 2.00 161.0 15.00
## 25 1.70 245.4 19.40
## 26 2.91 250.5 42.46
## 27 3.30 240.0 24.40
## bw_offspring_3wPP sd_bw_offspring_3wPP No_parents bw_parent_3wPP
## 22 49.8 4.46 24 338.4
## 23 52.5 4.50 27 335.2
## 24 56.2 6.90 24 299.0
## 25 54.1 5.40 27 338.4
## 26 55.5 4.97 25 308.0
## 27 47.7 5.30 25 355.0
## sd_bw_parent_3wPP index_study
## 22 20.6 yes
## 23 19.0 no
## 24 20.0 no
## 25 19.1 no
## 26 20.9 no
## 27 16.5 no
## study_type
## 22 OECD Guideline 416 (Two-Generation Reproduction Toxicity Study)
## 23 EPA OPPTS 870.3800 (Reproduction and Fertility Effects)
## 24 OECD Guideline 443 (Extended One-Generation Reproductive Toxicity Study)
## 25 EPA OPPTS 870.3800 (Reproduction and Fertility Effects)
## 26 OECD Guideline 443 (Extended One-Generation Reproductive Toxicity Study)
## 27 OECD Guideline 443 (Extended One-Generation Reproductive Toxicity Study)
## study study_period_start_date date testing_facility_location
## 22 TG416 04-Jun-98 04.06.1998 Isengard Research Division
## 23 EPA870 03-Jul-08 03.07.2008 Rohan Environmental Sciences
## 24 TG443 31-Aug-10 31.08.2010 Bree Laboratory Services
## 25 EPA870 14-Jun-11 14.06.2011 Rohan Environmental Sciences
## 26 TG443 24-Jun-21 24.06.2021 Shire Preclinical Services
## 27 TG443 05-Jul-16 05.07.2016 Lothlórien Contract Research
## facility type_of_diet_brand
## 22 ISG Purina Mills Rodent Lab Chow 5001-4
## 23 RHN LabDiet Certified Rodent Diet #5002
## 24 BRE Harlan Teklad Certified Rodent Diet 8728C
## 25 RHN LabDiet Certified Rodent Diet #5002
## 26 SHR PMI Nutrition International. LLC Certified Rodent LabDiet 5002
## 27 LTH SDS VRF1 Certified powdered diet
## diet supplier supplier_abb supplier_binned
## 22 PMRLC Charles River Laboratories US Crl_US Crl
## 23 CRD Charles River Laboratories US Crl_US Crl
## 24 HTCRD Taconic Farms (Germantown. NY) TF_US Other
## 25 LCRD Charles River Laboratories US Crl_US Crl
## 26 PMICRD Charles River Laboratories US Crl_US Crl
## 27 SVCPD Charles River Laboratories UK Crl_UK Crl
## mean_value n_value sd_value se_value variance_value study_date study_year
## 22 41.2 23 1.29 0.2689836 1.6641 1998-06-04 1998
## 23 43.0 27 1.30 0.2501851 1.6900 2008-07-03 2008
## 24 37.0 20 2.00 0.4472136 4.0000 2010-08-31 2010
## 25 43.5 26 1.70 0.3333974 2.8900 2011-06-14 2011
## 26 45.1 23 2.91 0.6067769 8.4681 2021-06-24 2021
## 27 48.0 20 3.30 0.7379024 10.8900 2016-07-05 2016
## v_value log_mean_value log_sd_value log_variance_value log_se_value
## 22 1.6641 3.717948 0.03130301 0.0009798784 0.006527129
## 23 1.6900 3.760743 0.03022565 0.0009135901 0.005816930
## 24 4.0000 3.609459 0.05401463 0.0029175805 0.012078039
## 25 2.8900 3.771998 0.03906555 0.0015261172 0.007661385
## 26 8.4681 3.806805 0.06445628 0.0041546115 0.013440063
## 27 10.8900 3.868843 0.06866897 0.0047154274 0.0153548483.2 Step 3: Descriptive Investigation of HCD
In Step 3, the HCD are inspected before formal modelling. For the continuous
summary-data workflow, the function plot_continuous_summary_hcd_descriptives()
is used to explore the distribution of the endpoint across relevant categorical
and continuous covariates.
Here, we focus on males for the endpoint age at preputial separation. The function creates descriptive plots for the response variable against study-level covariates such as diet, supplier, facility and study year, as well as continuous covariates such as offspring and parental body weight.
The first plot shows the relationship between mean age at preputial separation and mean offspring body weight at the endpoint. The second plot provides the same information in grouped form, making it easier to identify whether the response appears to vary systematically across body-weight ranges.
These descriptive plots are used to identify potential covariates that may need to be considered in the later modelling steps.
3.2.1 Step 3: Distributional Check of HCD Summary Means
As part of the descriptive investigation, the distribution of the HCD study-level means is inspected on both the original and log scale. This is a direct distributional check of the summary means, rather than a residual-based diagnostic from a fitted model.
The function assess_continuous_summary_hcd_normality() provides summary tests
and diagnostic plots, including histograms and Q-Q plots. Here, use_log = TRUE
is used so that both the original and log-transformed response are inspected.
## Scale Diagnostic method Diagnostic quantity N Mean SD Median
## raw raw_values Observed study-level mean 20.000 40.675 6.675 43.050
## log raw_values Observed study-level mean 20.000 3.688 0.197 3.762
## Minimum Maximum Skewness Shapiro-Wilk W Shapiro-Wilk p-value
## 23.050 48.000 -1.559 0.754 0.0002
## 3.137 3.869 -1.813 0.688 < 0.0001
## Kolmogorov-Smirnov D Kolmogorov-Smirnov p-value No. possible outliers
## 1.000 < 0.0001 0.000
## 0.999 < 0.0001 0.000
## Possible outlier IDs
##
##
## 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
## 
These plots provide a visual assessment of whether the HCD study-level means are approximately compatible with a normal distribution on the original scale or on the log scale. The log-scale diagnostics are included because later modelling steps may be performed on a transformed scale when this is more appropriate for the endpoint distribution.
3.2.2 Step 3: Outlier Detection
The next part of Step 3 is to check whether any historical control study appears unusual compared with the rest of the HCD distribution. For this summary-level continuous dataset, the outlier assessment is applied to the study-level means.
The function assess_continuous_summary_hcd_outliers() combines several
descriptive outlier checks, including quantile-based inspection, IQR-based
boxplots, z-scores and a gamma-distribution-based view. These outputs are used
as screening tools and should be interpreted together with biological and study
context.
## Method Statistic Value
## Empirical quantile 2.5% quantile 23.691
## Empirical quantile 95% quantile 45.720
## IQR boxplot Outside Q1 - 1.5*IQR or Q3 + 1.5*IQR 23.050
## IQR boxplot Outside Q1 - 1.5*IQR or Q3 + 1.5*IQR 24.400
## Z-score |z| > 2 -2.640
## Z-score |z| > 2 -2.438
## Hampel filter Median +/- 3*MAD 37.000
## Hampel filter Median +/- 3*MAD 23.050
## Hampel filter Median +/- 3*MAD 24.400
## Hampel filter Median +/- 3*MAD 34.040
## Hampel filter Median +/- 3*MAD 36.400
## Grubbs test 1.09744785546779; 0.933274692771698 <NA>
## Grubbs test 2.64036171301509; 0.613766760357009 <NA>
## Dixon test 0.496162528216704 <NA>
## Dixon test 0.200573065902579 <NA>
## Rosner test R.1 23.050
## Medcouple-adjusted boxplot Adjusted boxplot fences 48.000
## Gamma distribution 99% central interval 23.050
## Gamma distribution 99% central interval 24.400
## Gamma distribution 95% central interval 23.050
## Gamma distribution 95% central interval 24.400
## Lower limit Upper limit Observation Study ID p-value
## <NA> <NA> <NA> <NA> <NA>
## <NA> <NA> <NA> <NA> <NA>
## 33.850 50.850 11 ASB2008-1290 <NA>
## 33.850 50.850 12 ASB2009-5540 <NA>
## -2.000 2.000 11 ASB2008-1290 <NA>
## -2.000 2.000 12 ASB2009-5540 <NA>
## 37.950 48.150 2 ASB2017-6563 <NA>
## 37.950 48.150 11 ASB2008-1290 <NA>
## 37.950 48.150 12 ASB2009-5540 <NA>
## 37.950 48.150 14 TOX2001-343 <NA>
## 37.950 48.150 16 ASB2021-25643 <NA>
## <NA> <NA> <NA> <NA> 1.0000
## <NA> <NA> <NA> <NA> 0.0344
## <NA> <NA> <NA> <NA> 0.0451
## <NA> <NA> <NA> <NA> 0.8461
## <NA> <NA> 11 ASB2008-1290 <NA>
## 4.475 45.715 5 ASB2023-19175 <NA>
## 28.306 62.686 11 ASB2008-1290 <NA>
## 28.306 62.686 12 ASB2009-5540 <NA>
## 31.491 57.603 11 ASB2008-1290 <NA>
## 31.491 57.603 12 ASB2009-5540 <NA>
## Notes
##
##
##
##
## Original value = 23.05
## Original value = 24.4
##
##
##
##
##
## highest value 48 is an outlier
## lowest value 23.05 is an outlier
## lowest value 23.05 is an outlier
## highest value 48 is an outlier
## Outlier = FALSE
## Medcouple = -0.494
## shape = 42.626; rate = 0.978
## shape = 42.626; rate = 0.978
## shape = 42.626; rate = 0.978
## shape = 42.626; rate = 0.978
These plots help identify whether the extreme study means are isolated from the rest of the HCD distribution. Any potential outlier identified here should not be excluded automatically; rather, it should trigger further investigation of the study conditions and biological plausibility before deciding whether the study should be retained for Step 4.
3.2.3 Step 3 - Testing Associations with Categorical Covariates
After the descriptive plots, formal exploratory screening is performed for the
categorical covariates. The function
assess_continuous_summary_hcd_covariates() applies the tests separately for
each selected covariate, using only the HCD rows.
For each categorical covariate, the function reports both a Kruskal-Wallis test and a one-way ANOVA. The Kruskal-Wallis test provides a rank-based comparison of the response across covariate groups, while the ANOVA provides the corresponding parametric comparison of group means.
cov_res_male <- assess_continuous_summary_hcd_covariates(
data = hcd_prep$grouped_data$endpoint_preputial_separation__sex_male,
response_col = "age_offspring_endpoint",
categorical_covariates = c(
"diet",
"study",
"facility",
"supplier_abb",
"study_year"
),
hcd_indicator_col = "index_study",
hcd_indicator_value = "no",
covariate_labels = c(
diet = "Diet",
study = "Guideline",
facility = "Facility",
supplier_abb = "Supplier",
study_year = "Start date"
)
)
cov_res_male$formatted_table## covariate covariate_label method n n_groups
## 1 diet Diet Kruskal-Wallis 20 9
## 2 diet Diet ANOVA 20 9
## 3 study Guideline Kruskal-Wallis 20 3
## 4 study Guideline ANOVA 20 3
## 5 facility Facility Kruskal-Wallis 20 12
## 6 facility Facility ANOVA 20 12
## 7 supplier_abb Supplier Kruskal-Wallis 20 6
## 8 supplier_abb Supplier ANOVA 20 6
## 9 study_year Start date Kruskal-Wallis 20 12
## 10 study_year Start date ANOVA 20 12
## statistic_name statistic df1 df2 p_value note
## 1 Kruskal-Wallis chi-squared 13.290 8.000 <NA> 0.1023
## 2 F 7.530 8.000 11.000 0.0016
## 3 Kruskal-Wallis chi-squared 1.908 2.000 <NA> 0.3851
## 4 F 0.919 2.000 17.000 0.4180
## 5 Kruskal-Wallis chi-squared 14.482 11.000 <NA> 0.2074
## 6 F 12.384 11.000 8.000 0.0007
## 7 Kruskal-Wallis chi-squared 11.191 5.000 <NA> 0.0477
## 8 F 26.897 5.000 14.000 < 0.0001
## 9 Kruskal-Wallis chi-squared 10.828 11.000 <NA> 0.4578
## 10 F 1.011 11.000 8.000 0.5073Because this case study is based on aggregated study-level summaries, the association between covariates and the endpoint is evaluated at the study level. The subgroup forest plots should therefore be interpreted cautiously: apparent study-level associations do not necessarily reflect individual-level biological associations.
3.2.4 Step 3: Forest Plot by Subgroup
As an additional descriptive check, a subgroup forest plot is used to inspect the HCD study means within levels of a selected categorical covariate. Here, the studies are grouped by guideline. This helps assess whether the endpoint values appear broadly comparable across subgroups before proceeding to modelling.
3.2.5 Step 3 - Formal Covariate Assessment Using Meta-analysis
The descriptive plots are followed by a model-based assessment of selected
covariates using analyze_continuous_hcd_meta_covariate(). For continuous
summary-level data, this function uses the study-level mean, standard deviation
and sample size to assess whether the endpoint differs across levels of a
categorical covariate or changes with a continuous covariate.
For categorical covariates, such as guideline or animal supplier, the function
performs a subgroup meta-analysis. For continuous covariates, such as study year
or offspring body weight, it performs a meta-regression. In each case, the user
specifies the covariate name, whether it is "categorical" or "continuous",
and a readable label for the output.
Here, the analyses are first performed on the full male HCD dataset. In Step 3b, we will then define the restricted dataset to be carried forward to Step 4.
## Covariate Scale Subgroup No. studies Pooled mean 95% CI tau^2
## 1 Guideline orig EPA870 4 43.500 [42.485; 44.514] 44.907
## 2 Guideline orig TG416 13 39.177 [34.566; 43.787] 44.907
## 3 Guideline orig TG443 3 43.353 [29.183; 57.523] 44.907
## Q within I^2 (%) Q subgroup df subgroup p subgroup
## 1 11.83 74.6 4.081 2 0.13
## 2 14767.64 99.9 4.081 2 0.13
## 3 213.98 99.1 4.081 2 0.13
## Covariate Scale Subgroup No. studies Pooled mean 95% CI tau^2
## 1 Animal supplier orig Crl_G 2 23.724 [15.147; 32.3] 5.508
## 2 Animal supplier orig Crl_J 2 41.601 [37.789; 45.412] 5.508
## 3 Animal supplier orig Crl_UK 4 44.571 [40.908; 48.233] 5.508
## 4 Animal supplier orig Crl_US 10 43.312 [41.424; 45.2] 5.508
## 5 Animal supplier orig HW_G 1 34.040 [29.391; 38.689] 5.508
## 6 Animal supplier orig TF_US 1 37.000 [32.317; 41.683] 5.508
## Q within I^2 (%) Q subgroup df subgroup p subgroup
## 1 29.747 96.6 641.237 5 2.48e-136
## 2 2.153 53.5 641.237 5 2.48e-136
## 3 42.879 93.0 641.237 5 2.48e-136
## 4 611.924 98.5 641.237 5 2.48e-136
## 5 0.000 NA 641.237 5 2.48e-136
## 6 0.000 NA 641.237 5 2.48e-136
## Covariate Scale Term No. studies Estimate SE t value df
## 1 Study year orig intrcpt 20 -482.124 497.842 -0.968 18
## 2 Study year orig study_year 20 0.261 0.248 1.050 18
## 95% CI p value QM QM df QM p value tau^2
## 1 [-1528.05; 563.803] 0.346 1.103 1 0.308 44.286
## 2 [-0.261; 0.783] 0.308 1.103 1 0.308 44.286
## Covariate Scale Term No. studies
## 1 Offspring body weight at endpoint (g) orig intrcpt 20
## 2 Offspring body weight at endpoint (g) orig bw_offspring_endpoint 20
## Estimate SE t value df 95% CI p value QM QM df QM p value
## 1 18.280 1.369 13.357 18 [15.404; 21.155] 8.83e-11 286.286 1 1.69e-12
## 2 0.111 0.007 16.920 18 [0.098; 0.125] 1.69e-12 286.286 1 1.69e-12
## tau^2
## 1 2.69
## 2 2.69
## Covariate Scale Term
## 1 Offspring body weight 21 days post-partum (g) orig intrcpt
## 2 Offspring body weight 21 days post-partum (g) orig bw_offspring_3wPP
## No. studies Estimate SE t value df 95% CI p value QM QM df
## 1 20 13.414 14.853 0.903 18 [-17.791; 44.62] 0.3780 3.396 1
## 2 20 0.518 0.281 1.843 18 [-0.073; 1.109] 0.0819 3.396 1
## QM p value tau^2
## 1 0.0819 39.486
## 2 0.0819 39.486
## Covariate Scale Term No. studies
## 1 Dam body weight 21 days post-partum (g) orig intrcpt 20
## 2 Dam body weight 21 days post-partum (g) orig bw_parent_3wPP 20
## Estimate SE t value df 95% CI p value QM QM df QM p value
## 1 8.384 26.865 0.312 18 [-48.057; 64.826] 0.759 1.448 1 0.244
## 2 0.098 0.081 1.203 18 [-0.073; 0.269] 0.244 1.448 1 0.244
## tau^2
## 1 43.497
## 2 43.4973.3 Step 3b - Selection of HCD for Step 4
Based on the Step 3 investigations, five HCD studies were excluded before proceeding to Step 4.
The first two excluded studies, ASB2008-1290 and ASB2009-5540, had clearly
low study means and were flagged as potential outliers by the gamma-distribution
and boxplot-based outlier checks. These studies were also linked to the animal
supplier group identified in the covariate screening.
After excluding these two studies, the outlier assessment was repeated. Three
additional studies were then excluded: TOX2001-343, ASB2017-6563, and
ASB2021-25643. The report notes that these studies also had contextual
features supporting further investigation, including early study years, unique
animal suppliers, or unique study facilities.
Thus, the Step 3b decision was not based on a single automatic statistical rule. It combined outlier detection, covariate/context investigation, and expert assessment of the study characteristics.
The final excluded studies were: TOX2001-343, ASB2008-1290, ASB2009-5540,
ASB2017-6563, and ASB2021-25643. The Step 4 modelling of HCD variability was
then performed on the restricted HCD dataset with 15 HCD studies.
excluded_studies_male <- c(
"TOX2001-343",
"ASB2008-1290",
"ASB2009-5540",
"ASB2017-6563",
"ASB2021-25643"
)
hcd_male_step4 <- hcd_prep$grouped_data$endpoint_preputial_separation__sex_male[
!(hcd_prep$grouped_data$endpoint_preputial_separation__sex_male$study_ID %in%
excluded_studies_male),
]
table(hcd_male_step4$index_study)##
## no yes
## 15 1nrow(hcd_male_step4)## [1] 163.4 Step 4 - Modelling HCD Variability
In Step 4, the restricted HCD dataset selected in Step 3b is used to model the
between-study variability in the historical control means. For continuous
summary-level data, model_continuous_hcd_variability() uses the study-level
mean, standard deviation, standard error and sample size.
For this case study, the model is fitted on the original scale. In addition to an unadjusted model, an adjusted meta-regression model is fitted using offspring body weight at the endpoint as a continuous covariate. This reflects the case study objective of exploring body weight as a potential confounder for age at preputial separation.
step4_male <- model_continuous_hcd_variability(
data = hcd_male_step4,
mean_col = "mean_value",
sd_col = "sd_value",
n_col = "n_value",
se_col = "se_value",
study_id_col = "study_ID",
hcd_indicator_col = "index_study",
hcd_indicator_value = "no",
covariate_cols = "bw_offspring_endpoint",
covariate_types = c(bw_offspring_endpoint = "continuous"),
scale = "orig",
endpoint_label = "Age at preputial separation [days]",
make_forest = TRUE,
show_forest = TRUE,
save_forest = FALSE
)
step4_male$formatted_summary_table## model fitting_package adjusted covariates
## 1 meta_unadjusted meta FALSE none
## 2 metafor_unadjusted metafor FALSE none
## 3 metafor_adjusted metafor TRUE bw_offspring_endpoint
## prediction_at estimate ci_lower ci_upper pi_lower pi_upper tau2
## 1 overall HCD mean 43.8479 42.9522 44.7435 40.4543 47.2415 2.3292
## 2 overall HCD mean 43.8479 42.9731 44.7226 40.4597 47.2360 2.3292
## 3 concurrent control row 41.8971 39.7290 44.0652 38.2442 45.5500 1.8519
## tau Q Q_df Q_pvalue
## 1 1.5262 200.3994 14 0
## 2 1.5262 200.3994 14 0
## 3 1.3608 158.4161 13 0step4_male$formatted_manual_pi_table## model adjusted covariates
## 1 meta_unadjusted_manual_HK_style FALSE none
## 2 metafor_adjusted_manual_HK_style TRUE bw_offspring_endpoint
## prediction_at estimate pi_lower pi_upper var_pooled_HK tau2
## 1 overall HCD mean 43.8479 40.4296 47.2661 0.1744 2.3292
## 2 concurrent control row 41.8971 38.8234 44.9709 0.1724 1.8519step4_male$formatted_coefficient_table## method term estimate se statistic df p_value
## 1 metafor_rma.mv intrcpt 32.5611 5.4464 5.9784 13 0.0000
## 2 metafor_rma.mv bw_offspring_endpoint 0.0491 0.0237 2.0747 13 0.0584
## ci_lower ci_upper
## 1 20.7949 44.3274
## 2 -0.0020 0.1002The summary table reports the unadjusted and adjusted models, including the estimated mean, confidence interval, prediction interval and heterogeneity statistics. The manual prediction interval table provides the Step 4 prediction intervals used in the later comparison with the concurrent control. The displayed model labels reflect the confidence-interval method used in the Step 4 model. The coefficient table summarises the adjusted meta-regression effect of offspring body weight at the endpoint.
3.5 Step 5 - Comparison of Concurrent Control to HCD
In Step 5, the concurrent control from the index study is compared with the HCD prediction interval obtained in Step 4. For this continuous case study, the BfR report selected Approach 1, where the confidence interval for the concurrent control mean is compared with the HCD prediction interval for the population mean.
This approach is stricter than comparing only the observed concurrent-control mean. It is appropriate when the concurrent control is estimated reasonably precisely and one wants stronger evidence that the concurrent control population mean is compatible with the estimated between-study variability of the HCD.
The selected approach should normally be specified in the analysis protocol before inspecting the Step 5 results. For illustration, we also run Approach 2, where the concurrent-control sample mean is compared with the prediction interval for a new group mean of the same size.
step5_male_approach1 <- compare_continuous_cc_to_hcd(
step4_result = step4_male,
approach = "population_mean",
use_manual_pi = TRUE,
make_plots = TRUE
)
step5_male_approach2 <- compare_continuous_cc_to_hcd(
step4_result = step4_male,
approach = "new_group_mean",
use_manual_pi = TRUE,
make_plots = TRUE
)
step5_male_approach1$formatted_summary_table## Model Adjusted Covariates
## 1 meta_unadjusted_manual_HK_style FALSE none
## 2 metafor_adjusted_manual_HK_style TRUE bw_offspring_endpoint
## Covariate values used CC mean CC 90% CI lower CC 90% CI upper
## 1 None 41.2 40.7381 41.6619
## 2 Concurrent control 41.2 40.7381 41.6619
## HCD 95% PI lower HCD 95% PI upper Compatible
## 1 40.4296 47.2661 TRUE
## 2 38.8234 44.9709 TRUEstep5_male_approach2$formatted_summary_table## Model Adjusted Covariates
## 1 meta_unadjusted_manual_HK_style FALSE none
## 2 metafor_adjusted_manual_HK_style TRUE bw_offspring_endpoint
## Covariate values used CC mean HCD 95% PI lower HCD 95% PI upper Compatible
## 1 None 41.2 40.3806 47.3152 TRUE
## 2 Concurrent control 41.2 38.7690 45.0253 TRUEThe formatted Step 5 table is intentionally compact. For Approach 1, it shows the concurrent-control mean and 90% confidence interval, the HCD 95% prediction interval, and the compatibility conclusion. For Approach 2, it shows the concurrent-control mean, the HCD prediction interval, and the compatibility conclusion.
The plots below show the selected Approach 1 comparison. The dashed vertical lines show the HCD 95% prediction interval. The thick horizontal line shows the concurrent-control 90% confidence interval, and the square point shows the concurrent-control mean. Compatibility is assessed by comparing the concurrent control interval with the HCD prediction interval.
3.5.1 Step 5b - Sensitivity Around the Step 5 Comparison
Step 5b is used when the Step 5 comparison needs further investigation. Here, it is used to examine whether the conclusion from the population-mean approach is sensitive to the selected interval level or to the prediction-interval source.
This is an exploratory diagnostic step. It does not replace the Step 5 approach selected in the protocol, but helps assess how stable the compatibility conclusion is under alternative settings.
step5b_male <- investigate_continuous_cc_hcd_mismatch(
step4_result = step4_male,
approach = "population_mean",
alpha_values = c(0.10, 0.05, 0.025, 0.01),
pi_source = "both"
)
step5b_male$formatted_sensitivity_table## Approach CC interval level HCD PI level Tau2 scenario Adjusted
## 1 population_mean 0.80 0.95 central FALSE
## 2 population_mean 0.80 0.95 central TRUE
## 3 population_mean 0.90 0.95 central FALSE
## 4 population_mean 0.90 0.95 central TRUE
## 5 population_mean 0.95 0.95 central FALSE
## 6 population_mean 0.95 0.95 central TRUE
## 7 population_mean 0.98 0.95 central FALSE
## 8 population_mean 0.98 0.95 central TRUE
## 9 population_mean 0.80 0.95 as_returned FALSE
## 10 population_mean 0.80 0.95 as_returned FALSE
## 11 population_mean 0.80 0.95 as_returned TRUE
## 12 population_mean 0.90 0.95 as_returned FALSE
## 13 population_mean 0.90 0.95 as_returned FALSE
## 14 population_mean 0.90 0.95 as_returned TRUE
## 15 population_mean 0.95 0.95 as_returned FALSE
## 16 population_mean 0.95 0.95 as_returned FALSE
## 17 population_mean 0.95 0.95 as_returned TRUE
## 18 population_mean 0.98 0.95 as_returned FALSE
## 19 population_mean 0.98 0.95 as_returned FALSE
## 20 population_mean 0.98 0.95 as_returned TRUE
## Covariates CC mean CC CI lower CC CI upper HCD PI lower
## 1 none 41.2 40.8446 41.5554 40.4296
## 2 bw_offspring_endpoint 41.2 40.8446 41.5554 38.8234
## 3 none 41.2 40.7381 41.6619 40.4296
## 4 bw_offspring_endpoint 41.2 40.7381 41.6619 38.8234
## 5 none 41.2 40.6422 41.7578 40.4296
## 6 bw_offspring_endpoint 41.2 40.6422 41.7578 38.8234
## 7 none 41.2 40.5253 41.8747 40.4296
## 8 bw_offspring_endpoint 41.2 40.5253 41.8747 38.8234
## 9 none 41.2 40.8446 41.5554 40.4543
## 10 none 41.2 40.8446 41.5554 40.4597
## 11 bw_offspring_endpoint 41.2 40.8446 41.5554 38.2442
## 12 none 41.2 40.7381 41.6619 40.4543
## 13 none 41.2 40.7381 41.6619 40.4597
## 14 bw_offspring_endpoint 41.2 40.7381 41.6619 38.2442
## 15 none 41.2 40.6422 41.7578 40.4543
## 16 none 41.2 40.6422 41.7578 40.4597
## 17 bw_offspring_endpoint 41.2 40.6422 41.7578 38.2442
## 18 none 41.2 40.5253 41.8747 40.4543
## 19 none 41.2 40.5253 41.8747 40.4597
## 20 bw_offspring_endpoint 41.2 40.5253 41.8747 38.2442
## HCD PI upper Compatible
## 1 47.2661 TRUE
## 2 44.9709 TRUE
## 3 47.2661 TRUE
## 4 44.9709 TRUE
## 5 47.2661 TRUE
## 6 44.9709 TRUE
## 7 47.2661 TRUE
## 8 44.9709 TRUE
## 9 47.2415 TRUE
## 10 47.2360 TRUE
## 11 45.5500 TRUE
## 12 47.2415 TRUE
## 13 47.2360 TRUE
## 14 45.5500 TRUE
## 15 47.2415 TRUE
## 16 47.2360 TRUE
## 17 45.5500 TRUE
## 18 47.2415 TRUE
## 19 47.2360 TRUE
## 20 45.5500 TRUEThe formatted sensitivity table keeps only the main interpretation columns: approach, concurrent-control interval level, HCD prediction-interval level, tau-squared scenario, adjustment status, covariates, concurrent-control mean and interval, HCD prediction interval, and compatibility.
3.6 Step 6 - Analysis of the Index Study Including HCD
In Step 6, the index study is analysed with and without incorporation of HCD. For the continuous-data workflow, the HCD contribution is first converted into a pseudo-control row using the Step 4 model. This pseudo-control represents the HCD-derived information about the expected control response and its uncertainty. The raw pseudo-control object keeps the traceability columns, while the formatted table displays only the original-scale mean and SD, the log-scale mean and SD, and the effective sample size.
Important: The analyses shown with and without HCD are included only for illustrative purposes to demonstrate the corresponding Step 6 workflow and package outputs. In a real application of the EFSA workflow, the decision to include HCD in the analysis should be driven by the conclusions from the previous workflow steps.
The function prepare_continuous_index_study() prepares the dose-response data
from the index study. The function create_continuous_hcd_pseudo_control()
creates the HCD pseudo-control from the Step 4 results. Finally,
analyze_continuous_index_with_hcd() fits the index-study models with and
without the HCD pseudo-control.
# Read index-study dose-response data from package extdata
file_index_cont <- system.file(
"extdata",
"input_index study_case study 2.csv",
package = "efsaHCD"
)
doseresp <- read.csv(
file_index_cont,
header = TRUE,
stringsAsFactors = FALSE
)
# Prepare the male preputial-separation index-study data
index_male <- prepare_continuous_index_study(
data = doseresp,
endpoint_mean_col = "age_offspring_endpoint",
endpoint_sd_col = "sd_age_offspring_endpoint",
n_col = "No..of.animals_pups",
dose_col = "dose",
study_id_col = "study_ID",
endpoint_col = "endpoint",
endpoint_value = "preputial separation",
subgroup_col = "sex",
subgroup_value = "male",
control_dose = 0,
covariate_cols = "bw_offspring_endpoint",
date_col = "date",
date_format = "%d.%m.%Y",
endpoint_name = "Age at preputial separation [days]"
)##
## Prepared continuous index-study dataset
##
## Rows: 4
## Control rows: 1
## Treated rows: 3
## Endpoint mean column: age_offspring_endpoint
## Sample size column: No..of.animals_pupsindex_male$data## study_ID dose mean_value n_value sd_value se_value variance_value v_value
## 1 TOX2002-248 0 41.2 23 1.29 0.2689836 1.6641 1.6641
## 2 TOX2002-248 150 41.9 29 1.35 0.2506887 1.8225 1.8225
## 3 TOX2002-248 500 42.5 27 1.56 0.3002221 2.4336 2.4336
## 4 TOX2002-248 2500 47.9 27 2.75 0.5292377 7.5625 7.5625
## group_type index_code index_study endpoint sex
## 1 concurrent_control 2 yes preputial separation male
## 2 treated 3 yes preputial separation male
## 3 treated 3 yes preputial separation male
## 4 treated 3 yes preputial separation male
## bw_offspring_endpoint study_date study_year log_mean_value log_sd_value
## 1 190.1 1998-06-04 1998 3.717948 0.03130301
## 2 198.7 1998-06-04 1998 3.734767 0.03221121
## 3 191.9 1998-06-04 1998 3.748831 0.03669353
## 4 189.7 1998-06-04 1998 3.867470 0.05736405
## log_variance_value log_se_value original_study_ID original_endpoint
## 1 0.0009798784 0.006527129 TOX2002-248 preputial separation
## 2 0.0010375623 0.005981472 TOX2002-248 preputial separation
## 3 0.0013464150 0.007061673 TOX2002-248 preputial separation
## 4 0.0032906342 0.011039717 TOX2002-248 preputial separation
## original_sex original_No..of.animals_pups original_age_offspring_endpoint
## 1 male 23 41.2
## 2 male 29 41.9
## 3 male 27 42.5
## 4 male 27 47.9
## original_sd_age_offspring_endpoint original_bw_offspring_endpoint
## 1 1.29 190.1
## 2 1.35 198.7
## 3 1.56 191.9
## 4 2.75 189.7
## original_sd_bw_offspring_endpoint original_bw_offspring_3wPP
## 1 17.6 49.8
## 2 12.8 49.6
## 3 19.4 47.2
## 4 15.4 37.0
## original_sd_bw_offspring_3wPP original_No_parents original_bw_parent_3wPP
## 1 4.46 30 338.4
## 2 3.55 30 333.0
## 3 5.77 30 330.4
## 4 3.74 30 302.6
## original_sd_bw_parent_3wPP original_index_study original_study
## 1 20.60 yes TG416
## 2 22.67 yes TG416
## 3 18.20 yes TG416
## 4 16.52 yes TG416
## original_study_period_start_date original_date original_facility
## 1 04-Jun-98 04.06.1998 BCAD
## 2 04-Jun-98 04.06.1998 BCAD
## 3 04-Jun-98 04.06.1998 BCAD
## 4 04-Jun-98 04.06.1998 BCAD
## original_diet original_supplier_abb original_supplier_binned original_dose
## 1 PMRLC Crl_US Crl 0
## 2 PMRLC Crl_US Crl 150
## 3 PMRLC Crl_US Crl 500
## 4 PMRLC Crl_US Crl 2500# Create the HCD-derived pseudo-control rows from Step 4
hcd_pseudo_male <- create_continuous_hcd_pseudo_control(
step4_result = step4_male,
index_data = index_male,
pi_source = "manual"
)
hcd_pseudo_male$formatted_pseudo_control## Original mean Original SD Log mean Log SD Effective n
## 1 43.8479 8.0680 3.7807 0.1840 26
## 2 41.8971 7.2548 3.7352 0.1732 26# Analyse the index study with and without HCD
step6_male <- analyze_continuous_index_with_hcd(
index_data = index_male,
hcd_pseudo_control = hcd_pseudo_male,
dose_model = "both",
covariate_strategy = c("none", "pseudo_control_adjusted"),
include_cc_only = TRUE,
include_hcd = TRUE,
covariate_cols = "bw_offspring_endpoint",
make_plots = TRUE
)
step6_male$formatted_coefficient_table## Analysis HCD model Dose model Model term Estimate
## 1 cc_only none continuous intrcpt 41.3196
## 2 cc_only none continuous .dose_cont 0.0026
## 3 cc_only none categorical intrcpt 41.2000
## 4 cc_only none categorical .dose_fac150 0.7000
## 5 cc_only none categorical .dose_fac500 1.3000
## 6 cc_only none categorical .dose_fac2500 6.7000
## 7 with_hcd meta_unadjusted_manual_HK_style continuous intrcpt 42.0946
## 8 with_hcd meta_unadjusted_manual_HK_style continuous .dose_cont 0.0026
## 9 with_hcd meta_unadjusted_manual_HK_style categorical intrcpt 42.0649
## 10 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac150 0.6727
## 11 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac500 1.2727
## 12 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac2500 6.6727
## 13 with_hcd metafor_adjusted_manual_HK_style continuous intrcpt 41.3281
## 14 with_hcd metafor_adjusted_manual_HK_style continuous .dose_cont 0.0026
## 15 with_hcd metafor_adjusted_manual_HK_style categorical intrcpt 41.2241
## 16 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac150 0.6759
## 17 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac500 1.2759
## 18 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac2500 6.6759
## SE 95% CI lower 95% CI upper Term p-value Overall dose-effect p-value
## 1 0.1735 40.9796 41.6596 < 0.001 <0.001
## 2 0.0002 0.0022 0.0031 < 0.001 <0.001
## 3 0.2690 40.6728 41.7272 < 0.001 <0.001
## 4 0.3677 -0.0207 1.4207 0.05694 <0.001
## 5 0.4031 0.5099 2.0901 0.00126 <0.001
## 6 0.5937 5.5364 7.8636 < 0.001 <0.001
## 7 1.1657 39.8099 44.3793 < 0.001 <0.001
## 8 0.0002 0.0022 0.0031 < 0.001 <0.001
## 9 1.2418 39.6310 44.4987 < 0.001 <0.001
## 10 0.3667 -0.0460 1.3913 0.06658 <0.001
## 11 0.4022 0.4844 2.0609 0.00155 <0.001
## 12 0.5930 5.5103 7.8350 < 0.001 <0.001
## 13 0.1722 40.9905 41.6656 < 0.001 <0.001
## 14 0.0002 0.0022 0.0031 < 0.001 <0.001
## 15 0.2643 40.7060 41.7421 < 0.001 <0.001
## 16 0.3643 -0.0380 1.3899 0.06352 <0.001
## 17 0.4000 0.4920 2.0599 0.00142 <0.001
## 18 0.5916 5.5165 7.8354 < 0.001 <0.001step6_male$formatted_comparison_table## Analysis HCD model Dose model Model term Estimate
## 2 cc_only none continuous .dose_cont 0.0026
## 4 cc_only none categorical .dose_fac150 0.7000
## 5 cc_only none categorical .dose_fac500 1.3000
## 6 cc_only none categorical .dose_fac2500 6.7000
## 8 with_hcd meta_unadjusted_manual_HK_style continuous .dose_cont 0.0026
## 10 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac150 0.6727
## 11 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac500 1.2727
## 12 with_hcd meta_unadjusted_manual_HK_style categorical .dose_fac2500 6.6727
## 14 with_hcd metafor_adjusted_manual_HK_style continuous .dose_cont 0.0026
## 16 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac150 0.6759
## 17 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac500 1.2759
## 18 with_hcd metafor_adjusted_manual_HK_style categorical .dose_fac2500 6.6759
## SE 95% CI lower 95% CI upper Term p-value Overall dose-effect p-value
## 2 0.0002 0.0022 0.0031 < 0.001 <0.001
## 4 0.3677 -0.0207 1.4207 0.05694 <0.001
## 5 0.4031 0.5099 2.0901 0.00126 <0.001
## 6 0.5937 5.5364 7.8636 < 0.001 <0.001
## 8 0.0002 0.0022 0.0031 < 0.001 <0.001
## 10 0.3667 -0.0460 1.3913 0.06658 <0.001
## 11 0.4022 0.4844 2.0609 0.00155 <0.001
## 12 0.5930 5.5103 7.8350 < 0.001 <0.001
## 14 0.0002 0.0022 0.0031 < 0.001 <0.001
## 16 0.3643 -0.0380 1.3899 0.06352 <0.001
## 17 0.4000 0.4920 2.0599 0.00142 <0.001
## 18 0.5916 5.5165 7.8354 < 0.001 <0.001step6_male$formatted_descriptive_table## analysis_type hcd_model covariate_strategy
## 1 cc_only none none
## 2 cc_only none none
## 3 cc_only none none
## 4 cc_only none none
## 9 with_hcd meta_unadjusted_manual_HK_style none
## 10 with_hcd meta_unadjusted_manual_HK_style none
## 11 with_hcd meta_unadjusted_manual_HK_style none
## 12 with_hcd meta_unadjusted_manual_HK_style none
## 13 with_hcd meta_unadjusted_manual_HK_style none
## 19 with_hcd metafor_adjusted_manual_HK_style pseudo_control_adjusted
## 20 with_hcd metafor_adjusted_manual_HK_style pseudo_control_adjusted
## 21 with_hcd metafor_adjusted_manual_HK_style pseudo_control_adjusted
## 22 with_hcd metafor_adjusted_manual_HK_style pseudo_control_adjusted
## 23 with_hcd metafor_adjusted_manual_HK_style pseudo_control_adjusted
## study_ID dose mean_value se_value ci_lower ci_upper n_value
## 1 TOX2002-248 0 41.2000 0.2690 40.6728 41.7272 23
## 2 TOX2002-248 150 41.9000 0.2507 41.4087 42.3913 29
## 3 TOX2002-248 500 42.5000 0.3002 41.9116 43.0884 27
## 4 TOX2002-248 2500 47.9000 0.5292 46.8627 48.9373 27
## 9 HCD_pooled_1 0 43.8479 1.5823 40.7467 46.9491 26
## 10 TOX2002-248 0 41.2000 0.2690 40.6728 41.7272 23
## 11 TOX2002-248 150 41.9000 0.2507 41.4087 42.3913 29
## 12 TOX2002-248 500 42.5000 0.3002 41.9116 43.0884 27
## 13 TOX2002-248 2500 47.9000 0.5292 46.8627 48.9373 27
## 19 HCD_pooled_2 0 41.8971 1.4228 39.1085 44.6857 26
## 20 TOX2002-248 0 41.2000 0.2690 40.6728 41.7272 23
## 21 TOX2002-248 150 41.9000 0.2507 41.4087 42.3913 29
## 22 TOX2002-248 500 42.5000 0.3002 41.9116 43.0884 27
## 23 TOX2002-248 2500 47.9000 0.5292 46.8627 48.9373 27
## group_type hcd_adjusted hcd_covariates
## 1 concurrent_control NA <NA>
## 2 treated NA <NA>
## 3 treated NA <NA>
## 4 treated NA <NA>
## 9 historical_control_pooled FALSE none
## 10 concurrent_control NA <NA>
## 11 treated NA <NA>
## 12 treated NA <NA>
## 13 treated NA <NA>
## 19 historical_control_pooled TRUE bw_offspring_endpoint
## 20 concurrent_control NA <NA>
## 21 treated NA <NA>
## 22 treated NA <NA>
## 23 treated NA <NA>names(step6_male$plots)## [1] "forest__cc_only__continuous"
## [2] "dose_response__cc_only__continuous"
## [3] "forest__cc_only__categorical"
## [4] "dose_response__cc_only__categorical"
## [5] "forest__with_hcd__none__meta_unadjusted_manual_HK_style__continuous"
## [6] "dose_response__with_hcd__none__meta_unadjusted_manual_HK_style__continuous"
## [7] "forest__with_hcd__none__meta_unadjusted_manual_HK_style__categorical"
## [8] "dose_response__with_hcd__none__meta_unadjusted_manual_HK_style__categorical"
## [9] "forest__with_hcd__pseudo_control_adjusted__metafor_adjusted_manual_HK_style__continuous"
## [10] "dose_response__with_hcd__pseudo_control_adjusted__metafor_adjusted_manual_HK_style__continuous"
## [11] "forest__with_hcd__pseudo_control_adjusted__metafor_adjusted_manual_HK_style__categorical"
## [12] "dose_response__with_hcd__pseudo_control_adjusted__metafor_adjusted_manual_HK_style__categorical"plot_name <- grep(
"forest__with_hcd.*categorical",
names(step6_male$plots),
value = TRUE
)[1]
step6_male$plots[[plot_name]]
3.7 Step 7: Sensitivity Analysis
Step 7 of the EFSA HCD workflow evaluates the robustness of the conclusions under alternative assumptions, modelling strategies, or HCD selections. Typical sensitivity analyses may include repeating the analysis with different HCD subsets, alternative covariate adjustments, or different statistical models.
In the present vignette, a formal Step 7 sensitivity-analysis workflow was not implemented separately. This is because the package workflow already illustrated several alternative modelling approaches throughout Steps 4 to 6, including:
- unadjusted versus covariate-adjusted HCD models;
- different prediction-interval approaches;
- and analyses with and without incorporation of HCD pseudo-controls.
Therefore, the focus of this vignette was on demonstrating the main EFSA HCD workflow rather than performing an exhaustive sensitivity-analysis exercise. If needed, users can readily repeat the analysis using the available package functions under alternative assumptions, HCD selections, covariate structures, or modelling strategies.
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.
4 Count Data Case Study
This case study uses individual dam-level data from reproductive toxicity
studies in Wistar rats. The response variable is litter_size, defined as the
number of viable pups per dam at delivery.
Before applying the generic count-data HCD workflow, the raw dataset is prepared so that each row corresponds to one dam. This preparation step is case-specific, because it depends on the variable names, date format, generation coding, manual data checks and one duplicated/problematic study identified in the original analysis.
After this step, the cleaned dataset d can be passed to the generic count-data
functions in efsaHCD.
file_count <- system.file(
"extdata",
"input_case studies 5 and 6.csv",
package = "efsaHCD"
)
dat <- read.csv(
file_count,
header = TRUE,
na.strings = c("nd", "X", "")
)
# These checks/corrections were already addressed in the data file used here:
# - bw_birth: value 56 corrected to 5.6
# - bw_PND4: value 8..3 corrected to 8.3
# - extra spaces removed in selected fields
# - study 2019-12202 F1: index_study corrected from "No" to "Yes"
dat$date <- as.Date(dat$date, format = "%d.%m.%Y")
dat$n_date <- as.numeric(dat$date) / 365.25
dat$n_date <- dat$n_date + 1970
dat$year <- as.factor(as.numeric(format(dat$date, "%Y")))
dat$generation2 <- ifelse(dat$generation == "P", 1, 0)
dat$dam_ID <- paste(
dat$study_ID,
dat$generation,
dat$dam_ID_notunique
)
d <- dat[!duplicated(dat$dam_ID), ]
d <- d[d$study_ID != "ASB2010-12094", ]
table(d$study_ID, d$generation)##
## F1 P
## ASB2007-4840 27 27
## ASB2008-2016 20 19
## ASB2009-3519 26 21
## ASB2009-3719 20 22
## ASB2009-5909 28 28
## ASB2019-12202 23 89
## ASB2022-6248 26 27
## TOX2002-1278 28 26
## TOX2002-1333 0 10
## TOX2002-1334 27 24
## TOX2002-1679 20 23
## TOX2002-413 20 19
## TOX2006-1911 26 26
## TOX2006-839 19 28length(unique(d$dam_ID))## [1] 699length(unique(d$dam_ID[d$dose == 0]))## [1] 6334.1 Steps 1 and 2: Planning and Data Acceptability
In this case study, the endpoint is litter size, defined as the number of viable pups per dam at delivery. In Step 1, the analysis plan specifies that individual dam-level count data are required, together with relevant covariates such as generation, facility, diet, animal supplier, dam body weight and food consumption.
In Step 2, the reference studies are checked against the planned inclusion criteria. These include species, strain, study type, laboratory, animal supplier and study period. One duplicate/problematic study is excluded before proceeding with the HCD workflow.
The function prepare_count_hcd_data() standardises the count-data structure
for the package workflow. It creates a unique dam identifier, separates
historical control data from the index study and concurrent control, applies the
specified study exclusions, and prepares the selected covariates for later
analysis.
prep_count <- prepare_count_hcd_data(
data = dat,
response_col = "litter_size",
study_id_col = "study_ID",
dose_col = "dose",
control_dose = 0,
index_study_col = "index_study",
index_study_value = "Yes",
subject_id_cols = c("study_ID", "generation", "dam_ID_notunique"),
subject_id_new = "dam_ID",
deduplicate = TRUE,
date_col = "date",
date_format = "%d.%m.%Y",
factor_cols = c("supplier", "diet", "facility", "generation"),
relevel_factor_cols = c("supplier", "diet", "facility"),
covariate_cols = c(
"generation",
"food_consumption",
"bw_dam",
"bw_birth",
"bw_PND4"
),
exclude_study_ids = "ASB2010-12094",
restrict_expr = NULL,
keep_missing_response = FALSE
)
d <- prep_count$data
hcd <- prep_count$hcd
cc <- prep_count$cc
hcd_cc <- prep_count$hcd_cc
index_study <- prep_count$index_study
prep_count$settings## $response_col
## [1] "litter_size"
##
## $study_id_col
## [1] "study_ID"
##
## $dose_col
## [1] "dose"
##
## $control_dose
## [1] 0
##
## $index_study_col
## [1] "index_study"
##
## $index_study_value
## [1] "Yes"
##
## $subject_id_cols
## [1] "study_ID" "generation" "dam_ID_notunique"
##
## $subject_id_new
## [1] "dam_ID"
##
## $deduplicate
## [1] TRUE
##
## $date_col
## [1] "date"
##
## $date_format
## [1] "%d.%m.%Y"
##
## $factor_cols
## [1] "supplier" "diet" "facility" "generation"
##
## $relevel_factor_cols
## [1] "supplier" "diet" "facility"
##
## $covariate_cols
## [1] "generation" "food_consumption" "bw_dam" "bw_birth"
## [5] "bw_PND4"
##
## $exclude_study_ids
## [1] "ASB2010-12094"
##
## $restrict_expr
## NULL
##
## $keep_missing_response
## [1] FALSE
##
## $n_rows
## [1] 699
##
## $n_studies
## [1] 14
##
## $n_hcd_rows
## [1] 587
##
## $n_cc_rows
## [1] 46
##
## $n_index_rows
## [1] 112
##
## $n_hcd_cc_rows
## [1] 633prep_count$excluded_studies## [1] "ASB2010-12094"table(dat$litter_size, useNA = "ifany")##
## 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 12 24 50 66 91 400 640 1260 1397 1968 1209 546 270 80range(dat$litter_size, na.rm = TRUE)## [1] 3 16table(dat$litter_size > 6, useNA = "ifany")##
## FALSE TRUE
## 152 78614.2 Step 3 - Preliminary Investigation of HCD
4.2.1 Step 3a - Outlier Exploration
In Step 3, the historical control data are inspected before formal modelling. For this count-data case study, the first question is whether there are unusually small or large litter sizes, either at the individual-dam level or at the study-summary level.
The original report noted that very small litters occurred in several studies
and treated them as a special feature of the data-generating process. Therefore,
the function explore_count_hcd_outliers() is used here to inspect the
distribution of litter size and to evaluate the threshold used in the report,
namely keeping litters with more than six viable pups.
The exclusion of very small litters in this illustrative analysis was motivated by modelling considerations and by the possibility that these litters may represent a subgroup with different unmeasured characteristics. This should not be interpreted as implying that such litters are biologically implausible or invalid, since small litters can naturally occur in reproductive studies.
step3a_out <- explore_count_hcd_outliers(
data = hcd_cc,
hcd_data = hcd,
response_col = "litter_size",
study_id_col = "study_ID",
index_study_col = "index_study",
index_study_value = "Yes",
aggregate_by = c(
"study_ID",
"facility",
"supplier",
"diet",
"n_date",
"index_study"
),
aggregate_cols = c("litter_size", "bw_dam"),
threshold = 6,
threshold_rule = "greater",
fit_threshold_model = TRUE,
make_plots = TRUE,
save_plots = FALSE,
output_dir = "figures",
file_prefix = "step3a_count"
)
step3a_out$threshold_table##
## 0 1
## 28 605step3a_out$threshold_model_summary## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: threshold_indicator ~ 1 + (1 | study_ID)
## Data: dat
##
## AIC BIC logLik -2*log(L) df.resid
## 233.3 242.2 -114.7 229.3 631
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7591 0.2101 0.2133 0.2161 0.2222
##
## Random effects:
## Groups Name Variance Std.Dev.
## study_ID (Intercept) 0.03365 0.1834
## Number of obs: 633, groups: study_ID, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0883 0.2209 13.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1step3a_out$study_mean_shapiro_test##
## Shapiro-Wilk normality test
##
## data: response_mean
## W = 0.94839, p-value = 0.574head(step3a_out$aggregated_hcd)## study_ID facility supplier diet n_date index_study litter_size.mean
## 1 TOX2002-413 RVD SHR ALM 1997.822 No 10.589744
## 2 TOX2002-1333 ISG CRLUSA PMCD 1998.684 No 10.700000
## 3 TOX2002-1334 ISG CRLUSA PMCD 1999.090 No 10.470588
## 4 TOX2002-1278 ISG CRLUSA PMCD 1999.166 No 10.574074
## 5 TOX2006-839 ISG CRLUSA PMCD 2001.786 No 10.659574
## 6 TOX2002-1679 RVD SHR NAFAG 2002.906 No 10.418605
## litter_size.sd litter_size.length bw_dam.mean bw_dam.sd bw_dam.length
## 1 2.035448 39.000000 230.94872 17.00456 39.00000
## 2 2.213594 10.000000 215.52000 12.53580 10.00000
## 3 2.436004 51.000000 229.89020 19.52576 51.00000
## 4 1.754798 54.000000 228.56852 18.04464 54.00000
## 5 2.238755 47.000000 224.48936 16.57944 47.00000
## 6 2.342491 43.000000 221.95349 18.01977 43.00000The output summarises the number of observations below and above the selected threshold, fits the threshold model used to investigate whether small litters are concentrated in specific studies, and creates study-level summaries for further inspection. The Shapiro-Wilk test is applied to the study means as a simple check of their approximate distribution after aggregation.
4.2.2 Step 3a - Comparison of Count Model Families
After the initial outlier exploration, different model families are compared for the HCD litter-size data. This follows the logic of the original report, where a standard Poisson model was found to be unsuitable because the data were strongly under-dispersed.
The function compare_count_model_families() fits several candidate models with
study as a random effect and generation as an adjustment variable. The fitted
models are then compared using AIC and optional diagnostic checks.
model_comp <- compare_count_model_families(
data = hcd,
response_col = "litter_size",
study_id_col = "study_ID",
adjust_cols = "generation",
fit_poisson = TRUE,
fit_gaussian = TRUE,
fit_genpois_identity = TRUE,
fit_genpois_log = TRUE,
use_reml = TRUE,
run_performance_checks = TRUE,
make_plots = TRUE,
save_plots = FALSE,
output_dir = "figures",
file_prefix = "step3a2_model_compare"
)
model_comp$aic_table## model AIC fit_basis
## 2 gaussian_lmm 2602.629 REML
## 3 genpois_identity 2604.305 REML
## 4 genpois_log 2613.821 REML
## 1 poisson_glmm 2762.547 MLThe AIC table is used to identify which distributional assumption gives the best fit to the HCD. The output also stores the fitted model objects, so individual models can be inspected directly when needed. For example, the generalized Poisson model with identity link reports its estimated dispersion parameter. A value below 1 indicates under-dispersion relative to the standard Poisson model.
model_comp$models$genpois_identity## Formula: litter_size ~ 1 + generation + (1 | study_ID)
## Data: dat
## AIC BIC logLik -2*log(L) df.resid
## 2604.305 2621.805 -1298.153 2596.305 583
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## study_ID (Intercept) 0.0001031
##
## Number of obs: 587 / Conditional model: study_ID, 13
##
## Dispersion parameter for genpois family (): 0.415
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) generationP
## 10.5354 0.52414.2.3 Step 3a - Screening for Associations with Covariates
After selecting the generalized Poisson model family, we screen selected
covariates for association with litter size. The function
screen_count_hcd_covariate() fits one covariate model at a time, with study ID
as a random intercept and generation included as an adjustment variable.
Both categorical and continuous covariates can be assessed. For categorical covariates, the function compares levels of the covariate. For continuous covariates, the covariate can optionally be standardized before fitting. In this case study, facility and study date are fitted with an identity link, while dam body weight and food consumption are fitted with a log link, following the model choices used in the original analysis.
# Facility: categorical covariate, identity link
res_facility <- screen_count_hcd_covariate(
data = hcd,
plot_data = hcd_cc,
response_col = "litter_size",
covariate_col = "facility",
covariate_type = "categorical",
study_id_col = "study_ID",
adjust_cols = "generation2",
index_study_col = "index_study",
index_study_value = "Yes",
model_type = "glmmTMB_genpois",
link = "identity",
use_reml = FALSE,
standardize = FALSE,
make_plot = TRUE,
save_plot = FALSE
)
res_facility$coefficient_table## Estimate Std. Error z value Pr(>|z|) term
## 1 10.56546265 0.1832657 57.6510686 0.000000000 (Intercept)
## 2 -0.04156359 0.1929673 -0.2153919 0.829461785 facilityISG
## 3 0.52565462 0.1725657 3.0461122 0.002318213 generation2res_facility$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 10.2062685 10.9246568 10.5654626462 (Intercept)
## 2 -0.4197725 0.3366453 -0.0415635853 facilityISG
## 3 0.1874320 0.8638773 0.5256546226 generation2
## 4 0.0000000 Inf 0.0000891259 Std.Dev.(Intercept)|study_IDres_facility$aic_table## model AIC
## 1 reduced 2599.568
## 2 full 2601.521res_facility$plot
# Study date: continuous covariate, standardized, identity link
res_date <- screen_count_hcd_covariate(
data = hcd,
plot_data = hcd_cc,
response_col = "litter_size",
covariate_col = "n_date",
covariate_type = "continuous",
study_id_col = "study_ID",
adjust_cols = "generation2",
index_study_col = "index_study",
index_study_value = "Yes",
model_type = "glmmTMB_genpois",
link = "identity",
use_reml = FALSE,
standardize = TRUE,
make_plot = TRUE,
save_plot = FALSE
)
res_date$coefficient_table## Estimate Std. Error z value Pr(>|z|) term
## 1 10.5291855 0.12177950 86.461065 0.000000000 (Intercept)
## 2 0.1927035 0.08565525 2.249756 0.024464409 z_n_date
## 3 0.5378936 0.17156969 3.135132 0.001717767 generation2res_date$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 10.29050209 10.7678690 1.052919e+01 (Intercept)
## 2 0.02482225 0.3605847 1.927035e-01 z_n_date
## 3 0.20162323 0.8741641 5.378936e-01 generation2
## 4 0.00000000 Inf 5.794645e-05 Std.Dev.(Intercept)|study_IDres_date$aic_table ## model AIC
## 1 reduced 2599.568
## 2 full 2596.580res_date$plot
# Dam body weight: continuous covariate, standardized, log link
res_bw <- screen_count_hcd_covariate(
data = hcd,
plot_data = hcd_cc,
response_col = "litter_size",
covariate_col = "bw_dam",
covariate_type = "continuous",
study_id_col = "study_ID",
adjust_cols = "generation2",
index_study_col = "index_study",
index_study_value = "Yes",
model_type = "glmmTMB_genpois",
link = "log",
use_reml = FALSE,
standardize = TRUE,
make_plot = TRUE,
save_plot = FALSE
)
res_bw$coefficient_table## Estimate Std. Error z value Pr(>|z|) term
## 1 2.35925487 0.01138755 207.178463 0.000000e+00 (Intercept)
## 2 0.03681658 0.00788724 4.667866 3.043446e-06 z_bw_dam
## 3 0.03857518 0.01576291 2.447212 1.439660e-02 generation2res_bw$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 2.336935685 2.38157406 2.359255e+00 (Intercept)
## 2 0.021357872 0.05227529 3.681658e-02 z_bw_dam
## 3 0.007680449 0.06946991 3.857518e-02 generation2
## 4 0.000000000 Inf 7.858730e-06 Std.Dev.(Intercept)|study_IDres_bw$aic_table## model AIC
## 1 reduced 2599.568
## 2 full 2580.291res_bw$plot
# Food consumption: continuous covariate, standardized, log link
# Food consumption is not available for all dams, so complete cases are used.
hcd_food <- hcd[!is.na(hcd$food_consumption), ]
hcd_cc_food <- hcd_cc[!is.na(hcd_cc$food_consumption), ]
res_food <- screen_count_hcd_covariate(
data = hcd_food,
plot_data = hcd_cc_food,
response_col = "litter_size",
covariate_col = "food_consumption",
covariate_type = "continuous",
study_id_col = "study_ID",
adjust_cols = "generation2",
index_study_col = "index_study",
index_study_value = "Yes",
model_type = "glmmTMB_genpois",
link = "log",
use_reml = FALSE,
standardize = TRUE,
make_plot = TRUE,
save_plot = FALSE
)
res_food$coefficient_table## Estimate Std. Error z value Pr(>|z|) term
## 1 2.35328837 0.012576743 187.114293 0.0000000000 (Intercept)
## 2 0.01077932 0.008545557 1.261394 0.2071668146 z_food_consumption
## 3 0.05629761 0.016985682 3.314416 0.0009183483 generation2res_food$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 2.328638410 2.37793834 2.353288e+00 (Intercept)
## 2 -0.005969667 0.02752830 1.077932e-02 z_food_consumption
## 3 0.023006287 0.08958894 5.629761e-02 generation2
## 4 0.000000000 Inf 1.212949e-05 Std.Dev.(Intercept)|study_IDres_food$aic_table## model AIC
## 1 reduced 2322.592
## 2 full 2323.001res_food$plot
4.3 Step 3b - Selection of the Restricted HCD Dataset
Based on the exploratory analyses from Step 3a, the HCD dataset is restricted before proceeding to Step 4. In the original analysis, litters with six or fewer viable pups were treated separately because they may represent a subgroup with different unmeasured characteristics and created modelling difficulties.
The exploratory analyses also indicated:
- substantial under-dispersion relative to a standard Poisson model,
- improved fit for generalized Poisson models,
- associations with some continuous covariates, especially study date and dam body weight,
- weaker or negligible effects for several categorical covariates.
For consistency with the original modelling strategy, standardized versions of
continuous covariates are recreated here, and the binary generation indicator
(generation2) is retained for the later generalized Poisson models.
The resulting restricted dataset hcd_step4 is used in Step 4 for modelling
historical control variability.
# Restrict HCD to litters with more than six viable pups
hcd_step4 <- hcd[hcd$litter_size > 6, ]
# Standardized dam body weight
hcd_step4$z_bw_dam <-
(hcd_step4$bw_dam -
mean(hcd_step4$bw_dam, na.rm = TRUE)) /
sd(hcd_step4$bw_dam, na.rm = TRUE)
# Standardized study date
hcd_step4$z_n_date <-
(hcd_step4$n_date -
mean(hcd_step4$n_date, na.rm = TRUE)) /
sd(hcd_step4$n_date, na.rm = TRUE)
# Binary generation indicator used in the original scripts
hcd_step4$generation2 <-
ifelse(hcd_step4$generation == "P", 1, 0)
# Optional subset for analyses involving food consumption
hcd_step4_food <-
hcd_step4[!is.na(hcd_step4$food_consumption), ]
hcd_step4_food$z_food_consumption <-
(hcd_step4_food$food_consumption -
mean(hcd_step4_food$food_consumption, na.rm = TRUE)) /
sd(hcd_step4_food$food_consumption, na.rm = TRUE)
table(hcd_step4$litter_size <= 6)##
## FALSE
## 560length(unique(hcd_step4$study_ID))## [1] 134.4 Step 4: Modelling HCD Variability
In Step 4, the restricted HCD dataset from Step 3b is used to model historical control variability in litter size. The aim is to estimate the expected control response and the between-study variability, which will later be used to compare the concurrent control with the HCD.
For this count-data case study, the main model follows the original analysis and uses a generalized Poisson mixed model with identity link. Study ID is included as a random intercept and generation is included as a fixed adjustment variable.
An additional food-consumption and dam-body-weight-adjusted model is also fitted as a sensitivity analysis. This model uses the subset of HCD with available food-consumption data, includes standardized food consumption and dam body weight as continuous covariates, and uses a generalized Poisson model with a log link, matching the model structure used for these covariates in the original BfR analysis.
The main model estimates the overall HCD distribution after adjusting for generation only, whereas the alternative model additionally adjusts for food consumption and dam body weight to assess the sensitivity of the estimated HCD distribution to these covariates.
step4_main <- model_count_hcd_variability(
data = hcd_step4_food,
response_col = "litter_size",
study_id_col = "study_ID",
random_effects = "study_ID",
between_study_re = "study_ID",
covariate_cols = NULL,
adjust_cols = "generation2",
continuous_covariates = NULL,
standardize_continuous = FALSE,
model_type = "glmmTMB_genpois",
link = "identity",
use_reml = FALSE,
run_diagnostics = TRUE,
make_plots = TRUE,
save_plots = FALSE
)
step4_main$hcd_distribution## component mean sd lower upper link model_type
## 1 between_study 10.81202 0.154815 10.50859 11.11545 identity glmmTMB_genpois
## between_study_re
## 1 study_IDstep4_main$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 10.5701512 11.0538914 10.8120213 (Intercept)
## 2 0.3091271 0.9286613 0.6188942 generation2
## 3 0.0336754 0.7117266 0.1548150 Std.Dev.(Intercept)|study_ID# Alternative Step 4 model:
step4_main_food_subset <- model_count_hcd_variability(
data = hcd_step4_food,
response_col = "litter_size",
study_id_col = "study_ID",
random_effects = "study_ID",
between_study_re = "study_ID",
covariate_cols = c("food_consumption", "bw_dam"),
adjust_cols = "generation2",
continuous_covariates = c("food_consumption", "bw_dam"),
standardize_continuous = TRUE,
model_type = "glmmTMB_genpois",
link = "log",
use_reml = FALSE,
run_diagnostics = TRUE,
make_plots = TRUE,
save_plots = FALSE
)
step4_main_food_subset$hcd_distribution## component mean sd lower upper link model_type
## 1 between_study 10.89018 0.0204745 10.85005 10.93031 log glmmTMB_genpois
## between_study_re
## 1 study_IDstep4_main_food_subset$confidence_intervals## 2.5 % 97.5 % Estimate term
## 1 2.3643322007 2.41139126 2.38786173 (Intercept)
## 2 -0.0007147217 0.03042426 0.01485477 z_food_consumption
## 3 0.0234516309 0.05216685 0.03780924 z_bw_dam
## 4 0.0149686600 0.07026083 0.04261474 generation2
## 5 0.0075908855 0.05522480 0.02047450 Std.Dev.(Intercept)|study_ID4.5 Step 5 - Comparison of Concurrent Control to HCD
In Step 5, the concurrent control from the index study is compared with the HCD prediction distribution obtained from Step 4.
For this case study, the report used Approach 2. This approach compares the observed total number of pups in the concurrent-control group with the predicted distribution for a new control group of the same size. The concurrent-control dataset is restricted in the same way as the HCD used in Step 4, namely to litters with more than six viable pups.
For illustration, we also run Approach 1, where the concurrent-control mean and its confidence interval are compared with the HCD prediction interval for the population mean. However, the selected comparison for this case study is Approach 2.
# Concurrent-control data for Step 5:
# same endpoint restriction as used for the HCD in Step 4
cc <- d[
d$dose == 0 &
d$index_study == "Yes" &
d$litter_size > 6,
]
cc$generation2 <- ifelse(cc$generation == "P", 1, 0)
# Main Step 5 analysis: Approach 2
# Compare the observed total number of pups in the concurrent control
# with the predicted distribution for a new control group.
step5_main_a2 <- compare_count_cc_to_hcd(
cc_data = cc,
step4_result = step4_main,
response_col = "litter_size",
approach = "new_control_group_total",
sigma_source = "estimate",
random_effect_sampling = "observation",
N_MC = 10000,
x_range = 300:600,
pi_level = 0.95,
seed = 123,
make_plot = TRUE,
save_plot = FALSE
)
step5_main_a2$result$prediction_interval## lower upper
## 476 524step5_main_a2$result$observed_total## [1] 534step5_main_a2$result$inside_interval## [1] FALSEstep5_main_a2$plot
# Secondary illustration: Approach 1
# Compare the concurrent-control mean interval with the HCD
# prediction interval for the population mean.
step5_main_a1 <- compare_count_cc_to_hcd(
cc_data = cc,
step4_result = step4_main,
response_col = "litter_size",
approach = "population_mean",
sigma_source = "estimate",
pi_level = 0.95,
cc_ci_level = 0.95,
make_plot = TRUE,
save_plot = FALSE
)
step5_main_a1$result$hcd_prediction_interval## [1] 10.81116 11.41802step5_main_a1$result$cc_confidence_interval## [1] 11.28806 12.44527step5_main_a1$result$cc_mean## [1] 11.86667step5_main_a1$plot
4.6 Step 5b - Sensitivity Analyses
Because the selected Step 5 comparison indicated that the concurrent control was not compatible with the HCD prediction interval, the robustness of this result was explored further using sensitivity settings.
The function investigate_count_cc_hcd_mismatch() repeats the pre-specified
Step 5 approach under different prediction-interval levels. It does not compare
different Step 5 approaches and does not make regulatory or biological
conclusions. It is used only as a descriptive helper to assess how stable the
Step 5 result is under alternative assumptions.
Here, we apply the sensitivity analysis to the pre-specified Approach 2 using
the main Step 4 model. Instead of using the point estimate of the between-study
standard deviation, we use its upper confidence limit through
sigma_sources = "upper_ci". This is a more conservative setting because it
allows for greater between-study variability and therefore produces wider
prediction intervals.
The Step 4 model for this case study contains only the study-level random effect. To reproduce the original BfR calculation, this variability is included in the prediction but sampled independently for each simulated litter. This is the default behaviour of the updated functions, so no additional random-effect arguments are required.
step5b_main_a2 <- investigate_count_cc_hcd_mismatch(
cc_data = cc,
step4_result = step4_main,
response_col = "litter_size",
approach = "new_control_group_total",
pi_levels = c(0.90, 0.95, 0.99),
sigma_sources = "upper_ci",
N_MC = 10000,
seed = 123,
make_plot = TRUE
)
step5b_main_a2$formatted_sensitivity_table## pi_level sigma_source sigma_study observed_total observed_mean
## 1 0.90 upper_ci 0.7117266 534 11.86667
## 2 0.95 upper_ci 0.7117266 534 11.86667
## 3 0.99 upper_ci 0.7117266 534 11.86667
## prediction_lower prediction_upper prediction_lower_raw prediction_upper_raw
## 1 479 521 478.4082 520.9282
## 2 475 525 474.3296 524.9939
## 3 467 533 466.3486 532.9528
## observed_inside_prediction_interval mass_inside_interval
## 1 FALSE 0.9038931
## 2 FALSE 0.9515540
## 3 FALSE 0.9904574
## lower_tail_probability upper_tail_probability prediction_random_effects
## 1 0.046712345 0.049394548 study_ID
## 2 0.023475554 0.024970424 study_ID
## 3 0.004600342 0.004942232 study_ID
## group_level_random_effects
## 1 none
## 2 none
## 3 nonestep5b_main_a2$plot
Using the upper confidence limit of the between-study standard deviation produced slightly wider prediction intervals than the main analysis. Nevertheless, the observed concurrent-control total remained outside the HCD prediction interval at the 90%, 95%, and 99% levels. This indicates that the incompatibility between the concurrent control and the HCD was robust to both a more conservative estimate of between-study variability and the use of wider prediction intervals.
4.7 Step 6 - Analysis of the Index Study
Because Step 5 indicated that the concurrent control was not compatible with the HCD, the HCD are not incorporated into the Step 6 dose-response analysis. The index study is therefore analysed on its own.
The same restriction used in Step 5 is applied here: litters with six or fewer viable pups are excluded. We also restrict to records with available food consumption, following the model structure used in the report.
The function analyze_count_index_with_hcd() is used here with
hcd_data = NULL. This means that only the index-study data are analysed.
Several model specifications are fitted to inspect the dose-response
relationship and the role of dam body weight.
# Index-study data used in Step 6
index <- d[
d$litter_size > 6 &
!is.na(d$food_consumption) &
d$index_study == "Yes",
]
# Analyse the index study without HCD incorporation
step6_res <- analyze_count_index_with_hcd(
index_data = index,
hcd_data = NULL,
response_col = "litter_size",
dose_col = "dose",
study_id_col = "study_ID",
model_specs = list(
list(
name = "index_dose_only",
data_source = "index",
response_col = "litter_size",
predictor_cols = "dose",
continuous_covariates = "dose",
random_effect_cols = "study_ID",
model_type = "glmmTMB_genpois",
link = "log"
),
list(
name = "index_dose_bw",
data_source = "index",
response_col = "litter_size",
predictor_cols = c("bw_dam", "dose"),
continuous_covariates = "dose",
random_effect_cols = "study_ID",
model_type = "glmmTMB_genpois",
link = "log"
),
list(
name = "bw_dam_dose",
data_source = "index",
response_col = "bw_dam",
predictor_cols = "dose",
continuous_covariates = "dose",
random_effect_cols = "study_ID",
model_type = "glmmTMB_gaussian",
link = "identity"
)
),
make_plots = TRUE,
save_plots = FALSE
)
step6_res$model_summary_table## model data_source
## index_dose_only index_dose_only index
## index_dose_bw index_dose_bw index
## bw_dam_dose bw_dam_dose index
## formula model_type
## index_dose_only litter_size ~ z_dose + (1 | study_ID) glmmTMB_genpois
## index_dose_bw litter_size ~ bw_dam + z_dose + (1 | study_ID) glmmTMB_genpois
## bw_dam_dose bw_dam ~ z_dose + (1 | study_ID) glmmTMB_gaussian
## link n n_studies AIC BIC logLik
## index_dose_only log 109 1 441.1393 451.9047 -216.5696
## index_dose_bw log 109 1 435.0952 448.5519 -212.5476
## bw_dam_dose identity 109 1 934.7393 945.5047 -463.3696step6_res$formatted_coefficient_tables$index_dose_only## Term Estimate Standard error p-value
## 1 (Intercept) 2.4406 0.0148 <0.001
## 2 z_dose -0.0334 0.0147 0.024step6_res$formatted_coefficient_tables$index_dose_bw## Term Estimate Standard error p-value
## 1 (Intercept) 1.9166 0.1818 <0.001
## 2 bw_dam 0.0023 0.0008 0.004
## 3 z_dose -0.0183 0.0151 0.226step6_res$formatted_coefficient_tables$bw_dam_dose## Term Estimate Standard error p-value
## 1 (Intercept) 223.7982 1.6266 <0.001
## 2 z_dose -6.6492 1.6341 <0.001step6_res$plots$index_dose_only$response_by_dose
step6_res$plots$index_dose_bw$response_by_dose
step6_res$plots$bw_dam_dose$response_by_dose
4.8 Step 7 - Influence and Sensitivity Analysis
Step 7 is used to assess whether important decisions made earlier in the workflow influence the final interpretation. In the original report, several sensitivity analyses were performed for this count-data case study.
First, the Step 5 comparison was repeated using the food-consumption and dam-body-weight adjusted Step 4 model. The concurrent-control total still lay just outside the 95% prediction interval, so the conclusion of non-compatibility was unchanged.
Second, the upper confidence limit of the between-study variability was used instead of the point estimate. This widened the prediction interval, but the concurrent-control total still remained outside the 95% prediction interval.
Finally, the two generations were analysed separately. This showed that the parental-generation concurrent control was compatible with the HCD, whereas the F1-generation concurrent control was not. This indicated that the decision to pool both generations and adjust for generation influenced the Step 5 result for one generation, but not the other.
In this vignette, we do not repeat all Step 7 analyses separately. The same package functions shown in Steps 4 and 5 can be reused under alternative model specifications, prediction-interval assumptions, or HCD subsets when such sensitivity analyses are needed.
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.
5 Session information
sessionInfo()## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] readxl_1.5.0 patchwork_1.3.2 openxlsx_4.2.8.1 ggplot2_4.0.3
## [5] efsaHCD_0.0.7
##
## loaded via a namespace (and not attached):
## [1] mathjaxr_2.0-0 RColorBrewer_1.1-3 rstudioapi_0.19.0
## [4] jsonlite_2.0.0 magrittr_2.0.5 farver_2.1.2
## [7] nloptr_2.2.1 rmarkdown_2.31 fs_2.1.0
## [10] vctrs_0.7.3 minqa_1.2.8 mrfDepth_1.0.17
## [13] CompQuadForm_1.4.4 base64enc_0.1-6 htmltools_0.5.9
## [16] forcats_1.0.1 haven_2.5.5 cellranger_1.1.0
## [19] Formula_1.2-5 sass_0.4.10 bslib_0.11.0
## [22] htmlwidgets_1.6.4 plyr_1.8.9 sandwich_3.1-1
## [25] lubridate_1.9.5 rootSolve_1.8.2.4 zoo_1.8-15
## [28] cachem_1.1.0 TMB_1.9.21 lifecycle_1.0.5
## [31] pkgconfig_2.0.3 Matrix_1.7-5 R6_2.6.1
## [34] fastmap_1.2.0 rbibutils_2.4.1 magic_1.6-1
## [37] clue_0.3-68 digest_0.6.39 Exact_3.3
## [40] numDeriv_2016.8-1.1 colorspace_2.1-2 spatial_7.3-18
## [43] Hmisc_5.2-6 labeling_0.4.3 metadat_1.6-0
## [46] timechange_0.4.0 httr_1.4.8 abind_1.4-8
## [49] mgcv_1.9-4 compiler_4.5.1 proxy_0.4-29
## [52] withr_3.0.3 htmlTable_2.5.0 S7_0.2.2
## [55] backports_1.5.1 metafor_5.0-1 performance_0.17.1
## [58] meta_8.5-0 metabook_0.2-0 MASS_7.3-65
## [61] fBasics_4052.98 gld_2.6.8 tools_4.5.1
## [64] foreign_0.8-91 zip_3.0.0 nnet_7.3-20
## [67] glue_1.8.1 stabledist_0.7-2 nlme_3.1-169
## [70] grid_4.5.1 checkmate_2.3.4 cluster_2.1.8.2
## [73] reshape2_1.4.5 generics_0.1.4 glmmTMB_1.1.14
## [76] gtable_0.3.6 tzdb_0.5.0 class_7.3-23
## [79] tidyr_1.3.2 data.table_1.18.4 lmom_3.3
## [82] hms_1.1.4 univOutl_0.5.0 xml2_1.6.0
## [85] rmutil_1.1.10 pillar_1.11.1 stringr_1.6.0
## [88] robustbase_0.99-7 splines_4.5.1 dplyr_1.2.1
## [91] lattice_0.22-9 tidyselect_1.2.1 knitr_1.51
## [94] reformulas_0.4.4 gridExtra_2.3.1 bookdown_0.47
## [97] stats4_4.5.1 xfun_0.59 expm_1.0-0
## [100] timeDate_4052.112 matrixStats_1.5.0 LaplacesDemon_16.1.8
## [103] DEoptimR_1.2-0 stringi_1.8.7 VGAM_1.1-14
## [106] yaml_2.3.12 boot_1.3-32 evaluate_1.0.5
## [109] timeSeries_4052.112 tibble_3.3.1 cli_3.6.6
## [112] rpart_4.1.27 geometry_0.5.2 DescTools_0.99.60
## [115] Rdpack_2.6.6 jquerylib_0.1.4 Rcpp_1.1.1-1.1
## [118] stable_1.1.7 EnvStats_3.1.0 oaStyle_0.4.6
## [121] outliers_0.15 parallel_4.5.1 modeest_2.4.0
## [124] readr_2.2.0 CATTexact_0.1.1 lme4_2.0-1
## [127] mvtnorm_1.4-1 scales_1.4.0 e1071_1.7-17
## [130] statip_0.2.3 insight_1.5.2 purrr_1.2.2
## [133] rlang_1.2.0