This document accompanies the “A method to estimate probability of disease and vaccine efficacy from clinical trial immunogenicity data.” publication. It describes the application of PoDBAY package in the PoDBAY general simulation example.
The goal of PoDBAY simulation analysis is to simulate the clinical trial and validate the PoDBAY method by comparing “True” efficacy, Case-count efficacy and PoDBAY efficacy. By running the simulations we obtain following information:
PoDBAY simulation is summarized by three subsequent steps as described in the publication, section Methods.
method
parameter in the simulation inputs part below for further information.vignette("efficacyestimation", package = "PoDBAY")
.Steps 2 and 3 are repeated 1,000 times (1,000 simulations). Here we describe how to setup one simulation for the example based on Scenario A from the publication, see section Results.
To run a simulation, following data needs to be provided:
method
parameters for creating immunogenicity sample - for details see ?ImmunogenicitySubset()
and ?BlindSampling()
repeatCount
parameter - represents size of the PoDBAY efficacy set \(e^{``}\) we want to estimate. In our simulations we use repeatCount
= 500 (in vignettes the number is set to 50 because of time constrains for vignettes building). For details see PoDBAY efficacy estimation
section in vignette("efficacyestimation", package = "PoDBAY")
. adjustTiters
, adjustFrom
, adjustTo
. For details see ?PoD()
.This parameter decides how the Immunogenicity sample in the clinical trial is created (if any).
It is a named list with two parameters “name” and “value” which are discussed below.
method
= “Full”No immunogenicity sample is created. We have full titer information about diseased and non-diseased populations. See example of inactivated influenza vaccines for further details, section Results in the publication.
method
= “Ratio”Immunogenicity sample is created. We have full titer information about diseased but only partial titer information about nondiseased in the immunogenicity sample.
method
= “Fixed”Immunogenicity sample is created. We have full titer information about diseased but only partial titer information about nondiseased in the immunogenicity sample. See example of zoster vaccine and dengue vaccine for further details, section Results in the publication.
In some cases we might not be able to measure titer values below certain level (detection limit). In this case we might want to adjust level of titers below this detection limit to certain value.
There three parameters which needs to be set:
adjustTiters
= TRUE or FALSE - TRUE if titer values should be adjustedadjustFrom
- detection limit. All values below detection limit will be adjusted to adjustTo
parameter valueadjustTo
- value to which titers below detection limit will be adjustedDengue vaccine case might serve as an example where the detection limit for \(log_2\) titers is set to \(log_{2}10\). All \(log_2\) titers below detection limit are assigned to \(log_{2}5\) value. In this specific example we would set the parameters in the following way:
adjustTiters
= TRUEadjustFrom
= \(log_{2}10\)adjustTo
= \(log_{2}5\)Required libraries
vaccinated <- list()
vaccinated$N = 2000
vaccinated$mean = 8
vaccinated$stdDev = 2
control <- list()
control$N = 1000
control$mean = 5
control$stdDev = 2
PoDParams$pmax = 0.03
PoDParams$et50 = 7
PoDParams$slope = 7
#methods: "Full", "Ratio", "Fixed"
method <- list(name = "Full",
value = NA)
# method <- list(name = "Ratio",
# value = 4)
# method <- list(name = "Fixed",
# value = 300)
parameters <- list(vaccinated = vaccinated,
control = control,
PoDParams = PoDParams,
method = method,
repeatCount = 50,
adjustTiters = FALSE,
adjustFrom = NA,
adjustTo = NA)
True efficacy is calculated based on the provided parameters for vaccinated and control populations and PoD curve parameters.
Vaccinated and control populations are generated using generatePopulation()
function based on the simulation input parameters.
vaccinated <- generatePopulation(parameters$vaccinated$N,
parameters$vaccinated$mean,
parameters$vaccinated$stdDev)
control <- generatePopulation(parameters$control$N,
parameters$control$mean,
parameters$control$stdDev)
str(vaccinated)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 2000
#> $ mean : num 8
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:2000] 5.59 8.55 10.17 3.31 8.86 ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi(0)
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
str(control)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 1000
#> $ mean : num 5
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:1000] 3.05 4.8 4.78 7.38 1.69 ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi(0)
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
Probability of Disease (PoD) is assigned to each patient based on individual titer values and PoD curve parameters. Function PoD()
and population class method asignPoD()
are used.
vaccinated$assignPoD(
PoD(titer = vaccinated$popX(),
pmax = PoDParams$pmax,
et50 = PoDParams$et50,
slope = PoDParams$slope,
adjustTiters = parameters$adjustTiters,
adjustFrom = parameters$adjustFrom,
adjustTo = parameters$adjustTo)
)
control$assignPoD(
PoD(titer = control$popX(),
pmax = PoDParams$pmax,
et50 = PoDParams$et50,
slope = PoDParams$slope,
adjustTiters = parameters$adjustTiters,
adjustFrom = parameters$adjustFrom,
adjustTo = parameters$adjustTo)
)
str(vaccinated)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 2000
#> $ mean : num 8
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:2000] 5.59 8.55 10.17 3.31 8.86 ...
#> $ PoDs : num [1:2000] 0.02487 0.00591 0.00205 0.02984 0.00484 ...
#> $ diseaseStatus : logi(0)
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
str(control)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 1000
#> $ mean : num 5
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:1000] 3.05 4.8 4.78 7.38 1.69 ...
#> $ PoDs : num [1:1000] 0.0299 0.028 0.0281 0.0122 0.03 ...
#> $ diseaseStatus : logi(0)
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
Disease status (DS) is assigned based on the PoD of each patient. Case-count efficacy with its confidence intervals (80%, 90% and 95% level of significance) is estimated based on the DS. ClinicalTrialCoverage()
or ClinicalTrial()
function is used depending which CIs we want to estimate.
CaseCountEfficacy <- ClinicalTrialCoverage(vaccinated,
control)
list(CaseCountEfficacy = CaseCountEfficacy$efficacy,
confidenceInterval95 = unlist(CaseCountEfficacy$confidenceInterval95),
confidenceInterval90 = unlist(CaseCountEfficacy$confidenceInterval90),
confidenceInterval80 = unlist(CaseCountEfficacy$confidenceInterval80))
#> $CaseCountEfficacy
#> [1] 0.4565217
#>
#> $confidenceInterval95
#> lowerBound upperBound
#> 0.04743002 0.68992449
#>
#> $confidenceInterval90
#> lowerBound upperBound
#> 0.1296098 0.6606481
#>
#> $confidenceInterval80
#> lowerBound upperBound
#> 0.2155981 0.6234474
str(vaccinated)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 2000
#> $ mean : num 8
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:2000] 5.59 8.55 10.17 3.31 8.86 ...
#> $ PoDs : num [1:2000] 0.02487 0.00591 0.00205 0.02984 0.00484 ...
#> $ diseaseStatus : logi [1:2000] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
str(control)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : num 1000
#> $ mean : num 5
#> $ stdDev : num 2
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : num [1:1000] 3.05 4.8 4.78 7.38 1.69 ...
#> $ PoDs : num [1:1000] 0.0299 0.028 0.0281 0.0122 0.03 ...
#> $ diseaseStatus : logi [1:1000] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
Diseased (diseased) and non-diseased (nondiseased) populations are created based on the DS of each patient. ExtractDiseased()
and ExtractNondiseased
functions with population class methods getDiseasedTiters()
, getNondiseasedTiters()
, getDiseasedCount()
and getNondiseasedCount()
are used.
diseasedAll <- ExtractDiseased(CaseCountEfficacy$vaccinated,
CaseCountEfficacy$control)
nondiseasedAll <- ExtractNondiseased(CaseCountEfficacy$vaccinated,
CaseCountEfficacy$control)
str(diseasedAll)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : int 48
#> $ mean : num 5.53
#> $ stdDev : num 1.94
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : Named num [1:48] 6.13 5.73 6.7 6.38 6.84 ...
#> ..- attr(*, "names")= chr [1:48] "vacc" "vacc" "vacc" "vacc" ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi [1:48] TRUE TRUE TRUE TRUE TRUE TRUE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
str(nondiseasedAll)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : int 2952
#> $ mean : num 7.04
#> $ stdDev : num 2.43
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : Named num [1:2952] 5.59 8.55 10.17 3.31 8.86 ...
#> ..- attr(*, "names")= chr [1:2952] "vacc" "vacc" "vacc" "vacc" ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi [1:2952] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
Depending on the chosen method, a required immunogenicity sample (IS) is created out of the full population in the clinical trial.
method
= “Full” - Immunogenicity sample is NOT created. Immunogenicity sample = whole clinical trialmethod
= “Ratio” - Immunogenicity sample is created.method
= “Fixed” - Immunogenicity sample is created.Then \(Nondiseased_{IS}\), \(Vaccinated_{IS}\) and \(Control_{IS}\) populations are identified within this immunogenicity sample.
In reality, in the case of “Fixed”" method patients in the Immunogenicity sample are picked before enrolling to the clinical trial thus both disease statuses (diseased and non-diseased) are possible to appear in the IS.
In the case of “Ratio”" method patients in the Immunogenicity sample are picked based on the number of diseased in the clinical trial and the whole Immunogenicity sample has non-diseased disease status.
\(Nondiseased_{IS}\)
Subjects from immunogenicity sample with “non-diseased” disease status.
str(ImmunogenicitySample$ImmunogenicityNondiseased)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : int 2952
#> $ mean : num 7.04
#> $ stdDev : num 2.43
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : Named num [1:2952] 7.92 6.33 10.95 9.81 8.67 ...
#> ..- attr(*, "names")= chr [1:2952] "vacc" "vacc" "vacc" "vacc" ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi [1:2952] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
\(Vaccinated_{IS}\)
Subjects from immunogenicity sample with “vaccinated” vaccination status.
str(ImmunogenicitySample$ImmunogenicityVaccinated)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : int 2000
#> $ mean : num 7.99
#> $ stdDev : num 1.98
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : Named num [1:2000] 7.92 6.33 10.95 9.81 8.67 ...
#> ..- attr(*, "names")= chr [1:2000] "vacc_FALSE" "vacc_FALSE" "vacc_FALSE" "vacc_FALSE" ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi [1:2000] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
\(Control_{IS}\)
Subjects from immunogenicity sample with “control” vaccination status.
str(ImmunogenicitySample$ImmunogenicityControl)
#> Reference class 'Population' [package "PoDBAY"] with 8 fields
#> $ N : int 1000
#> $ mean : num 5.06
#> $ stdDev : num 2.02
#> $ unknownDistribution: logi FALSE
#> $ UDFunction :function ()
#> $ titers : Named num [1:1000] 3.07 5.11 9.23 5.95 4.58 ...
#> ..- attr(*, "names")= chr [1:1000] "control_FALSE" "control_FALSE" "control_FALSE" "control_FALSE" ...
#> $ PoDs : num(0)
#> $ diseaseStatus : logi [1:1000] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> and 24 methods, of which 10 are possibly relevant:
#> assignPoD, getDiseasedCount, getDiseasedTiters, getNondiseasedCount,
#> getNondiseasedTiters, getTiters, getUnknown, initialize, popFun, popX
PoDBAY efficacy follows the structure of EfficacyEstimation vignette - see vignette("efficacyestimation", package = "PoDBAY")
for further details.
PoD curve is estimated (Point estimate together with confidence intervals) in three steps - further details can be found in the publication, section Methods.
Titers of all diseased and all non-diseased subjects are used for estimation of PoD curve parameters. Parameter estimates \(p_{max}^`\), \(et_{50}^`\) and \(\gamma^`\) are obtained.
Titers of all diseased and all non-diseased subjects are put together and bootstrapped. For each individual titer a probability of disease is calculated using the PoD curve with parameter values \(p_{max}^`\), \(et_{50}^`\) and \(\gamma^`\). New disease status is assigned to each titer based on the probability of disease.
Titers of all new diseased and all new non-diseased subjects are used for re-estimation of PoD curve parameters. Parameter estimates \(p_{max}^{``}\), \(et_{50}^{``}\) and \(\gamma^{``}\) are obtained.
estimatedParameters <- PoDParamEstimation(diseasedAll$titers,
ImmunogenicitySample$ImmunogenicityNondiseased$titers,
nondiseasedAll$N,
parameters$repeatCount,
adjustTiters = parameters$adjustTiters,
adjustFrom = parameters$adjustFrom,
adjustTo = parameters$adjustTo)
# Confidence intervals
PoDParamsCI <- PoDParamsCICoverage(estimatedParameters$results)
unlist(PoDParamsCI)
#> PmaxCILow95 PmaxCIHigh95 PmaxCILow90 PmaxCIHigh90 PmaxCILow80
#> 0.02066134 0.05089425 0.02135638 0.04454633 0.02313595
#> PmaxCIHigh80 Et50CILow95 Et50CIHigh95 Et50CILow90 Et50CIHigh90
#> 0.03740940 5.25016770 8.52981399 6.17865695 8.41597546
#> Et50CILow80 Et50CIHigh80 SlopeCILow95 SlopeCIHigh95 SlopeCILow90
#> 6.46806621 8.28328580 4.26560267 31.41464940 5.85036433
#> SlopeCIHigh90 SlopeCILow80 SlopeCIHigh80
#> 30.57352127 6.06343116 29.55021314
# Point estimate
PoDParamsPointEst <- PoDParamPointEstimation(estimatedParameters$resultsPriorReset)
unlist(PoDParamsPointEst)
#> pmax slope et50
#> 0.02781535 11.26514623 7.53419999
PoDBAY Efficacy (Point estimate together with confidence intervals) is estimated - further details can be found in the publication, section Methods.
PoDParamsPointEst
from step 1 PoD-titer relationship estimation
estimatedParameters$results
from step 1 PoD-titer relationship estimation
step 2
) and standard deviationsstep 3
# Point estimate
meansBlind <- list("vaccinated" = ImmunogenicitySample$ImmunogenicityVaccinated$mean,
"control" = ImmunogenicitySample$ImmunogenicityControl$mean)
standardDeviationsBlind <- list("vaccinated" = ImmunogenicitySample$ImmunogenicityVaccinated$stdDev,
"control" = ImmunogenicitySample$ImmunogenicityControl$stdDev)
EfficacyPointEst <- efficacyComputation(PoDParamsPointEst, meansBlind, standardDeviationsBlind)
# PoDBAY efficacy set
efficacySet <- PoDBAYEfficacy(estimatedParameters$results,
ImmunogenicitySample$ImmunogenicityVaccinated,
ImmunogenicitySample$ImmunogenicityControl,
adjustTiters = parameters$adjustTiters,
adjustFrom = parameters$adjustFrom,
adjustTo = parameters$adjustTo)
# PoDBAY efficacy confidence intervals
CI <- EfficacyCICoverage(efficacySet)
EfficacyPointEst
#> [1] 0.4983887
unlist(CI)
#> mean median CILow95 CIHigh95 CILow90 CIHigh90 CILow80 CIHigh80
#> 0.4966027 0.4954539 0.3630324 0.6286025 0.3868161 0.6183307 0.3989669 0.6128655
Analysis provides following results:
result <- list(
TrueEfficacy = TrueEfficacy,
CaseCountEfficacy = CaseCountEfficacy$efficacy,
confidenceInterval95 = unlist(CaseCountEfficacy$confidenceInterval95),
confidenceInterval90 = unlist(CaseCountEfficacy$confidenceInterval90),
confidenceInterval80 = unlist(CaseCountEfficacy$confidenceInterval80),
EfficacyPointEst = EfficacyPointEst,
efficacyCI = unlist(CI),
PoDParamsPointEst = unlist(PoDParamsPointEst),
PoDParamsCI = unlist(PoDParamsCI))
result
#> $TrueEfficacy
#> [1] 0.5336551
#>
#> $CaseCountEfficacy
#> [1] 0.4565217
#>
#> $confidenceInterval95
#> lowerBound upperBound
#> 0.04743002 0.68992449
#>
#> $confidenceInterval90
#> lowerBound upperBound
#> 0.1296098 0.6606481
#>
#> $confidenceInterval80
#> lowerBound upperBound
#> 0.2155981 0.6234474
#>
#> $EfficacyPointEst
#> [1] 0.4983887
#>
#> $efficacyCI
#> mean median CILow95 CIHigh95 CILow90 CIHigh90 CILow80 CIHigh80
#> 0.4966027 0.4954539 0.3630324 0.6286025 0.3868161 0.6183307 0.3989669 0.6128655
#>
#> $PoDParamsPointEst
#> pmax slope et50
#> 0.02781535 11.26514623 7.53419999
#>
#> $PoDParamsCI
#> PmaxCILow95 PmaxCIHigh95 PmaxCILow90 PmaxCIHigh90 PmaxCILow80
#> 0.02066134 0.05089425 0.02135638 0.04454633 0.02313595
#> PmaxCIHigh80 Et50CILow95 Et50CIHigh95 Et50CILow90 Et50CIHigh90
#> 0.03740940 5.25016770 8.52981399 6.17865695 8.41597546
#> Et50CILow80 Et50CIHigh80 SlopeCILow95 SlopeCIHigh95 SlopeCILow90
#> 6.46806621 8.28328580 4.26560267 31.41464940 5.85036433
#> SlopeCIHigh90 SlopeCILow80 SlopeCIHigh80
#> 30.57352127 6.06343116 29.55021314
In a special case when serum samples at baseline and after vaccination are collected and assayed only in a subset of subjects (“immunogenicity sample/ subset”) and the assay value of titer is obtained also for all disease cases at the same time points, the general method for PoD curve estimation described above can be extended. Further details can be found in the publication Appendix A. Details about the PoD curve estimation can be found in the vignette("efficacyestimation", package = "PoDBAY")
Appendix.
If you are interested in running the simulation with this setup, change the method
parameter in the simulations parameter input and run the simulation. The only change would be that immunogenicity sample is created and PoD curve is estimated accordingly inside the simulations. The analysis steps remain the same.
vaccinated <- list()
vaccinated$N = 2000
vaccinated$mean = 8
vaccinated$stdDev = 2
control <- list()
control$N = 1000
control$mean = 5
control$stdDev = 2
PoDParams$pmax = 0.03
PoDParams$et50 = 7
PoDParams$slope = 7
#methods: "Full", "Ratio", "Fixed"
# method <- list(name = "Full",
# value = NA)
# method <- list(name = "Ratio",
# value = 4)
method <- list(name = "Fixed",
value = 300)
parameters <- list(vaccinated = vaccinated,
control = control,
PoDParams = PoDParams,
method = method,
repeatCount = 50,
adjustTiters = FALSE,
adjustFrom = NA,
adjustTo = NA)