expirest_wisle: What-if (approach for) shelf life estimation (wisle)

Description

Based on a linear regression model fitted to a stability data set the function expirest_wisle() estimates the shelf life, or retest period, for the specified release and specification limit following the ARGPM guidance “Stability testing for prescription medicines”. The abbreviation “wisle” stands for “what-if shelf life estimation” (because it estimates the shelf life (“what”) for a given release limit (“if”)).

Usage

expirest_wisle(
  data,
  response_vbl,
  time_vbl,
  batch_vbl,
  rl,
  rl_sf,
  sl,
  sl_sf,
  srch_range,
  alpha = 0.05,
  alpha_pool = 0.25,
  xform = c("no", "no"),
  shift = c(0, 0),
  sf_option = "tight",
  ivl = "confidence",
  ivl_type = "one.sided",
  ivl_side = "lower",
  ...
)

Value

An object of class ‘expirest_wisle’ is returned, containing the following elements:

Data: Data frame of the original data including new columns with transformed variables, if applicable.
Parameters: A list of the parameters with the elements alpha, alpha.pool, ivl, ivl.type and ivl.side.
Variables: A list of the variable names, i.e. the original names of batch_vbl, time_vbl and response_vbl and, if applicable, of the transformed variables.
Model.Type: A list of two elements that specifies which model, based on the ANCOVA analysis, suits best. The first element (type.spec) is a numeric vector of length 2 that specifies the best model accepted at the significance level specified by alpha.pool. The first number represents the decision on the intercept and the second on the slope, where 1 stands for “common” and 2 stands for “different”. The second element (type.acronym) is an acronym representing the first item.
Models: A list of four elements named cics, dics, dids.pmse and dids. The first three elements contain the ‘lm’ objects of the “common intercept / common slope” (cics), “different intercept / common slope” (dics) and “different intercept / different slope” (dids) models. The fourth element is a list of the ‘lm’ objects that is obtained from fitting a regression model to the data of each level of the categorical variable separately. The cics, dics and dids.pmse elements are NA if data of only a single batch is available.
AIC: A numeric named vector of the Akaike Information Criterion (AIC) values of the cics, dics and dids.pmse models.
BIC: A numeric named vector of the Bayesian Information Criterion (BIC) values of each of the cics, dics and dids.pmse models.
wc.icpt: A data frame of the worst case intercepts of each of the four fitted models.
wc.batch: A list of numeric value(s) of the worst case batch(es) per model type.
Limits: A list of all limits.
POI: A data frame of the intercepts, the differences between release and shelf life limits, the WCSLs, the expiry and release specification limits, the shelf lives and POI values.

Structure of the POI data frame:

Intercept.cics: The intercept of the worst case batch of the cics model.
Intercept.dics: The intercept of the worst case batch of the dics model.
Intercept.dids.pmse: The intercept of the worst case batch of the dids model with pooled mean square error (pmse).
Intercept.dids: The intercept of the worst case batch of the dids model obtained by fitting individual models to the data of each batch.
Delta.cics: Absolute difference between the release and and the shelf life specification of the cics model.
Delta.dics: Absolute difference between the release and and the shelf life specification of the dics model.
Delta.dids.pmse: Absolute difference between the release and and the shelf life specification of the dids model with pooled mean square error (pmse).
Delta.dids: Absolute difference between the release and and the shelf life specification of the dids model obtained by fitting individual models to the data of each batch.
WCSL.cics: WCSL of the cics model.
WCSL.dics: WCSL of the dics model.
WCSL.dids.pmse: WCSL of the dids model with pooled mean square error (pmse).
WCSL.dids: WCSL of the dids model obtained by fitting individual models to the data of each batch.
Exp.Spec: The (expiry) specification, i.e. the specification which is relevant for the determination of the expiry.
Rel.Spec: The calculated release specification.
Shelf.Life.cics: The estimated shelf life of the cics model.
Shelf.Life.dics: The estimated shelf life of the dics model.
Shelf.Life.dids.pmse: The estimated shelf life of the dids model with pooled mean square error (pmse).
Shelf.Life.dids: The estimated shelf life of the dids model obtained by fitting individual models to the data of each batch.
POI.Model.cics: The POI of the cics model.
POI.Model.dics: The POI of the dics model.
POI.Model.dids.pmse: The POI of the dids model with pooled mean square error (pmse).
POI.Model.dids: The POI of the dids model obtained by fitting individual models to the data of each batch.

Arguments

data: A data frame with the columns specified by response_vbl, time_vbl and batch_vbl.
response_vbl: A character string that specifies the response variable name that must be a column of data.
time_vbl: A character string that specifies the time variable name that must be a column of data.
batch_vbl: A character string that specifies the column in data with the grouping information (i.e. a factorial variable) for the differentiation of the observations of the various batches.
rl: A numeric value or a numeric vector that specifies the release specification limit(s) for which the corresponding expiry should be estimated.
rl_sf: A positive integer or a vector of positive integers that specifies the number of “significant figures” (sf) of rl. It must have the same length as rl.
sl: A numeric value or a numeric vector of length 2 that specifies the specification limit or limits. If a vector is provided it must be of the form c(lower limit, upper limit).
sl_sf: A positive integer or a vector of positive integers that specifies the number of “significant figures” (sf) of sl. It must have the same length as sl.
srch_range: A vector of length 2 that specifies the end-points of the (time) range that is supposed to contain the shelf life or retest period.
alpha: A numeric value between 0 and 1 that specifies the significance level for the calculation of confidence or prediction intervals. The default is 0.05.
alpha_pool: A numeric value between 0 and 1 that specifies the type I error rate for the test of the poolability of the batches. The default is 0.25, following ICH Q1E guideline. The value should not be changed unless supported by well-founded reasons.
xform: A vector of two character strings that specifies the transformation of the response and the time variable. The default is “no” transformation, i.e. c("no", "no"), where the first element specifies the transformation of the \(x\) variable and the second element the transformation of the \(y\) variable. Valid alternatives for \(x\) and/or \(y\) variable transformation are "log" (natural logarithm), "sqrt" (square root) and "sq" (square).
shift: A vector of two values which will be added to the variables \(x\) and/or \(y\) before they are transformed as specified by the xform parameter, where the first element will be added to the \(x\) variable and the second element to the \(y\) variable. The purpose is to prevent an undefined state which could arise when variables with values of \(\leq 0\) are log or square root transformed. The default is c(0, 0).
sf_option: A character string that specifies if the limits (rl or sl) should be regarded as “tight” or “loose”, i.e. either "tight" or "loose", respectively. The default is "tight". The option "tight" means that the limits are rounded to the number of significant digits specified by the parameters rl_sf and sl_sf. If sf_option = "loose", four on the order of the last significant digit + 1 is added if ivl_side = "upper" or five on the order of the last significant digit + 1 is subtracted if ivl_side = "upper".
ivl: A character string of either "confidence" or "prediction" that specifies the type of interval of interest. The default is "confidence".
ivl_type: A character string that specifies if a “one sided” or a “two sided” interval should be calculated, i.e. either "one.sided" or "two.sided", respectively. The default is "one.sided".
ivl_side: A character string that specifies if the “upper” or the “lower” limit is the relevant limit, i.e. either "upper" or "lower", respectively. The default is "lower". Since this parameter additionally specifies the relationship of rl with sl, i.e. which of the two sides of sl the rl is compared to, only either either "upper" or "lower" is possible. In this respect, the usage of ivl_side differs from its usage in the expirest_osle() function where ivl_side in addition can be "both".
...: Additional named or unnamed arguments passed on to uniroot().

Checking batch poolability

According to ICH Q1E guideline, construction of the 95% confidence interval on the basis of the combined data of all test batches is allowed only if it has been confirmed by aid of a statistical test whether the regression lines from the different batches have a common slope and a common intercept. A significance level of alpha_pool = 0.25 should to be used for both batch-related terms, and the test of the slopes has to precede the test of the intercepts. From these tests, three possible models may be appropriate, i.e.

a common intercept / common slope model (cics),
a different intercept / common slope model (dics) or
a different intercept / different slope model (dids).

The common intercept / different slope model is not of practical relevance because the corresponding model is missing an effect. If the slopes are significantly different, there is no point comparing intercepts. The dids model has individual intercepts and individual slopes, and the calculation of confidence intervals is based on the corresponding individual mean square errors. The different intercept / different slope model where the mean square error is pooled across batches is reported as dids.pmse.

These requirements can be checked by aid of an “ANalysis of COVAriance” (ANCOVA) including the batch variable as main effect and as \(batch \times time\) interaction term. The full ANCOVA model simultaneously tests all the effects, and non-significant effects can be identified and removed for fitting of the final regression model that is used for the estimation of the shelf life or retest period.

The significance level (alpha_pool = 0.25, Type I error) is used to increase the power of the test to detect cases where the data should not be pooled. Setting alpha_pool = 0.25 decreases the probability of incorrectly concluding that stability data from multiple batches can be pooled. On the other hand, though, it increases the probability of using a single batch to determine expiry when pooling batches would be more appropriate.

Details

For the shelf life estimation for submissions in Australia the Australian Regulatory Guidelines for Prescription Medicines (ARGPM), i.e. the guidance “Stability testing for prescription medicines”, is binding. In chapter 14.3.1, “Predicting shelf life from stability data”, it is described how the estimation should be done. It recommends to take the worst case situation at batch release into account. The following examples are listed:

1a): For medicine that has a lower Assay release limit of 95 per cent and a lower assay expiry limit of 90 per cent, the maximum decrease in assay allowed over the shelf-life is 5 per cent.
2a): Similarly, for a medicine that has a release limit for an individual impurity of 0.2 per cent and an expiry limit of 0.5 per cent, the maximum increase in the impurity allowed over the shelf life is 0.3 per cent.
1b): For the same medicine, if the actual release assay result is 101 per cent, then the shelf life should be determined at the time the medicine (or confidence interval of the regression line) decreases by 5 per cent and reaches 96 per cent, rather than the expiry limit (90 per cent). This takes into account the possibility of batches being released at the lower release limit (i.e. 95 per cent) and ensures they will comply with the expiry limit throughout the shelf life.
2b): Similarly, if the actual impurity content is 0.1 per cent at batch release, the shelf life should be determined as the time the 95 per cent confidence interval reaches 0.4 per cent (i.e. increases by 0.3 per cent).

Consequently, it is necessary to define

a release limit and
an expiry limit or shelf life limit (usually the specification limit).

If both of these limits were the same the shelf life would be zero, i.e. no change allowed.

For the estimation of the shelf life or expiry limit it is necessary to find the point where the upper or lower 95% confidence interval limit of the linear model fitted to the data (by aid of lm) intersects the “worst case scenario limit” (WCSL), which is defined by the ARGPM guidance “Stability testing for prescription medicines” as the intercept \(\pm\) difference between the specification limit and the release limit, where this differences is added (\(+\)) if the upper limit is relevant or it is subtracted (\(-\)) if the lower limit is relevant. In this package, this point is called the “point of intersection” or “point of interest”, abbreviated POI.

Before performing the retest period or shelf life estimation the most suitable model should be determined. It should particularly be verified if data of all test batches are poolable or not. Details on this are described in section “Checking batch poolability” below.

References

Therapeutic Goods Administration (TGA) of the Department of Health of the Australian Government, Australian Regulatory Guidelines for Prescription Medicines (ARGPM), Stability testing for prescription medicines, Version 1.1, March 2017

International Council for Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human (ICH), Harmonised Tripartite Guideline, Evaluation of Stability Data Q1E, step 4, February 2003 (CPMP/ICH/420/02).

Examples

Run this code

# Successful estimations
# A model with common intercepts / common slopes (cics)
res1 <-
  expirest_wisle(data = exp1[exp1$Batch %in% c("b2", "b5", "b7"), ],
                 response_vbl = "Potency", time_vbl = "Month",
                 batch_vbl = "Batch", rl = 98, rl_sf = 3, sl = 95,
                 sl_sf = 3, srch_range = c(0, 500), sf_option = "loose")
res1$Model.Type
res1$POI

# Expected results in res1$Model.Type
# $type.spec
# common.icpt  common.slp
#           1           1
#
# $type.acronym
# [1] "cics"

# (Expected) results in res1$POI
# Intercept.cics Intercept.dics Intercept.dids Intercept.dids.pmse Delta.cics
#       100.5669       100.3638       100.2491            100.2491          3
# Delta.dics Delta.dids Delta.dids.pmse WCSL.cics WCSL.dics WCSL.dids
#          3          3               3  97.56688  97.36375  97.24914
# WCSL.dids.pmse Exp.Spec.Report Exp.Spec Rel.Spec.Report Rel.Spec
#       97.24914              95    94.95              98    97.95
# Shelf.Life.cics Shelf.Life.dics Shelf.Life.dids Shelf.Life.dids.pmse
#        14.07398        13.23176        13.34561             13.80695
# POI.Model.cics POI.Model.dics POI.Model.dids POI.Model.dids.pmse
#        26.2241        24.8003       23.34184            23.66724

# The parameter settings sf_option = "loose" and ivl_side = "lower" (the
# default setting of ivl_side) cause the specification limit of 95.0
# (sl_sf = 3, i.e. 3 significant digits) to be reduced by 0.05, i.e. the
# actual specification limit is 94.95. Analogously, the release limit of 98.0
# is reduced to 97.95.

# A model with different intercepts / different slopes (dids)
res2 <-
  expirest_wisle(data = exp1[exp1$Batch %in% c("b4", "b5", "b8"), ],
                 response_vbl = "Potency", time_vbl = "Month",
                 batch_vbl = "Batch", rl = 98, rl_sf = 3, sl = 95,
                 sl_sf = 3, srch_range = c(0, 500), sf_option = "loose")
res2$Model.Type
res2$POI

# Expected results in res2$Model.Type
# $type.spec
# common.icpt  common.slp
#           0           0
#
# $type.acronym
# [1] "dids"

# (Expected) results in res2$POI
# Intercept.cics Intercept.dics Intercept.dids Intercept.dids.pmse Delta.cics
#       101.5498       100.4882       101.2594            101.2594          3
# Delta.dics Delta.dids Delta.dids.pmse WCSL.cics WCSL.dics WCSL.dids
#          3          3               3  98.54976  97.48822  98.25938
# WCSL.dids.pmse Exp.Spec.Report Exp.Spec Rel.Spec.Report Rel.Spec
#       98.25938              95    94.95              98    97.95
# Shelf.Life.cics Shelf.Life.dics Shelf.Life.dids Shelf.Life.dids.pmse
#        13.03332        11.42141        7.619661             7.483223
# POI.Model.cics POI.Model.dics POI.Model.dids POI.Model.dids.pmse
#       28.12518       22.47939       15.96453            15.72348

# Analysis with a single batch (i.e. the worst case batch of res2 above)
res3 <-
  expirest_wisle(data = exp1[exp1$Batch == "b8", ],
                 response_vbl = "Potency", time_vbl = "Month",
                 batch_vbl = "Batch", rl = 98, rl_sf = 3, sl = 95,
                 sl_sf = 3, srch_range = c(0, 500), sf_option = "loose")
res3$Model.Type
res3$POI

# Since only one batch is involved there is no model type. Nevertheless, the
# result is reported under the dids model name.

# Expected results in res3$Model.Type
# $type.spec
# common.icpt  common.slp
#         NA           NA
#
# $type.acronym
# [1] "n.a."

# (Expected) results in res3$POI
# Intercept.cics Intercept.dics Intercept.dids Intercept.dids.pmse Delta.cics
#             NA             NA       101.2594                  NA         NA
# Delta.dics Delta.dids Delta.dids.pmse WCSL.cics WCSL.dics WCSL.dids
#         NA          3              NA        NA        NA  98.25938
# WCSL.dids.pmse Exp.Spec.Report Exp.Spec Rel.Spec.Report Rel.Spec
#             NA              95    94.95              98    97.95
# Shelf.Life.cics Shelf.Life.dics Shelf.Life.dids Shelf.Life.dids.pmse
#              NA              NA        7.619661                   NA
# POI.Model.cics POI.Model.dics POI.Model.dids POI.Model.dids.pmse
#             NA             NA       15.96453                  NA

# Unsuccessful estimation
if (FALSE) {
  # The interval does not cross the limit (i.e. not within srch_range).
  res4 <-
    expirest_wisle(data = exp1[exp1$Batch %in% c("b2", "b5", "b7"), ],
                   response_vbl = "Potency", time_vbl = "Month",
                   batch_vbl = "Batch", rl = 98, rl_sf = 3, sl = 105,
                   sl_sf = 4, srch_range = c(0, 500), sf_option = "loose",
                   ivl_side = "upper")
  res4$POI

  # (Expected) results in res3$POI
  # Intercept.cics Intercept.dics Intercept.dids Intercept.dids.pmse Delta.cics
  #       100.5669             NA             NA                  NA          3
  # Delta.dics Delta.dids Delta.dids.pmse WCSL.cics WCSL.dics WCSL.dids
  #         NA         NA              NA  107.5669        NA        NA
  # WCSL.dids.pmse Exp.Spec.Report Exp.Spec Rel.Spec.Report Rel.Spec
  #             NA             105   105.04              98    98.04
  # Shelf.Life.cics Shelf.Life.dics Shelf.Life.dids Shelf.Life.dids.pmse
  #              NA              NA              NA                   NA
  # POI.Model.cics POI.Model.dics POI.Model.dids POI.Model.dids.pmse
  #             NA             NA             NA                  NA

  # Estimation may also fail because of an inappropriate 'srch_range' setting.
  res5 <-
    expirest_wisle(data = exp1[exp1$Batch %in% c("b2", "b5", "b7"), ],
                   response_vbl = "Potency", time_vbl = "Month",
                   batch_vbl = "Batch", rl = 98, rl_sf = 3, sl = 95,
                   sl_sf = 3, srch_range = c(0, 5), sf_option = "loose")
  res5$POI

  # (Expected) results in res4$POI
  # Intercept.cics Intercept.dics Intercept.dids Intercept.dids.pmse Delta.cics
  #       100.5669             NA             NA                  NA          3
  # Delta.dics Delta.dids Delta.dids.pmse WCSL.cics WCSL.dics WCSL.dids
  #         NA         NA              NA  97.56688        NA        NA
  # WCSL.dids.pmse Exp.Spec.Report Exp.Spec Rel.Spec.Report Rel.Spec
  #             NA              95    94.95              98    97.95
  # Shelf.Life.cics Shelf.Life.dics Shelf.Life.dids Shelf.Life.dids.pmse
  #              NA              NA              NA                   NA
  # POI.Model.cics POI.Model.dics POI.Model.dids POI.Model.dids.pmse
  #             NA             NA             NA                  NA
}

Run the code above in your browser using DataLab