SD and geoSD of ratio of PK parameters for population simulations seem to be wrongly calculated

[AnnikaRPS]

The ratio of PK parameters for population workflows seems to be calculated as the ratio of the respective summary statistics of the population of interest and the control population, e.g. if the case example compares a pediatric population to an adult population the ratio summary parameters are calculated as:

mean of ratio = mean_pediatric/mean_adult
geomean of ratio = geomean_pediatric/geomean_adult
SD of ratio = sd_pediatric/sd_adult
geoSD of ratio = geosd_pediatric/geosd_adult

This calculation of ratios is an appropriate approximation in case of mean and geomean but wrong for SD and geoSD. Here, the delta method is needed.

For the arithmetic SD this would be:
z = x/y
dzdx = 1/meanY
dzdy = -meanX/meanY^2
meanZ = meanX/meanY
sdZ = sqrt(dzdx^2 times sdX^2+dzdy^2 times sdY^2)

For the geometric SD this would be:
z = x/y
geomeanZ = exp(log(geomeanX)-log(geomeanY))
geosdZ = exp(sqrt(log(geosdX)^2+log(geosdY)^2))

[msevestre]

@Yuri05 Can you confirm?

[Yuri05]

Can you confirm?

Wrong person to ask 😄

[AnnikaRPS]

@sfrechen

[sfrechen]

@AnnikaRPS
Where are we interested in these ratios? -> only in case the reference population is not the same as the comparator.
Otherwise we are interested in the statistics over the individual ratios!

[AnnikaRPS]

@AnnikaRPS
Where are we interested in these ratios? -> only in case the reference population is not the same as the comparator.
Otherwise we are interested in the statistics over the individual ratios!

That's true, in case the reference and the comparator population are the same (e.g. in DDI scenarios), statistics on the individual ratios are more precise

[sfrechen]

Can you confirm?

Yes. I do confirm. The equations given above calculate an estimate of the respective statistic of the distribution of the ratio. This is only applicable in case that indviduals differ between the two populations to be compared.

Please note that we need in all other cases (same population is compared under different scenarios, e.g. DDI), we DEFINITELY need the statitiscis over the individual ratios, i.e.

Let's say we are interested in the ratio of the AUC (AUCR):

Let AUC_x be n AUC values of scenario x, and AUC_y be n AUC values of the reference scenario y.

geoMean_AUCR = geomean(AUC_x/AUC_y) = exp(mean(log(AUC_x/AUC_y)))
geoSD_AUCR = geoSD(AUC_x/AUC_y) = exp(sd(log(AUC_x/AUC_y)))
geoCV_AUCR = geoCV(AUC_x/AUC_y) = sqrt(exp((sd(log(AUC_x/AUC_y)))^2-1))

arithMean_AUCR = mean(AUC_x/AUC_y)
arithSD_AUCR = sd(AUC_x/AUC_y)
arithCV_AUCR = sd(AUC_x/AUC_y)/mean(AUC_x/AUC_y)

quantile_q_AUCR = quantile(AUC_x/AUC_y, q) with q being the q th quantile between 0 and 1

[Yuri05]

OK, that means that for every calculation of PK-Ratios between 2 SimulationSets: RE must check if the used population is the same or not and depending on that use the appropriate calculation.

[sfrechen]

Derivations of equations provided by @AnnikaRPS

Being z = x/y introduces a new random variable in a nonlinear relationship of the primary random variables x and y.
Thus, we indeed need the so called Delta Method or better known as first order Taylor expansion to get rid of the ratio (with a cut after the first derivative, i.e. we linearize the issue).

Arith. Mean

Applying the first order taylor expansion to z = x/y, we obtain
arithMeanZ is indeed best estimated via arithMeanX/arithMeanY
with arithMeanX = mean(x) and arithMeanY = mean(y), i.e. we obtain:
arithMeanZ ≈ arithMeanX/arithMeanY

Arith. SD

Applying the first order taylor expansion to z = x/y, we really need to account for teh derivatives as @AnnikaRPS showsa and we obtain
with arithSdX = sd(x) and arithSdY = sd(y):
arithSdZ ≈ sqrt(dzdx^2*arithSdX^2+dzdy^2*arithSdY^2)
The derivatives are evaluated at arithMeanX and arithMeanY.

For the geometric statistics, things are easier since the logarithm already linearizes the issue.

GeoMean

Define new random variable z as the ratio x over y:
z = x/y

log transform both sides:
log(z) = log(x/y) = log(x)-log(y)

take the mean:
mean(log(z)) = mean(log(x)-log(y)) = mean(log(x))-mean(log(y))

exponentiate both sides:
exp(mean(log(z))) = exp(mean(log(x))-mean(log(y)))

Left hand side is geoMean of z: geoMeanZ := exp(mean(log(z)))
geoMeanZ = exp(mean(log(x))-mean(log(y)))

With geoMeanX := exp(mean(log(x))) and geoMeanY := exp(mean(log(y))):
geoMeanZ = exp(log(geoMeanX)-log(geoMeanY))

GeoSD

Define new random variable z as the ratio x over y:
z = x/y

log transform both sides:
log(z) = log(x/y) = log(x)-log(y)

take the variance:
var(log(z)) = var(log(x)-log(y)) = var(log(x))+var(log(y))

transform to standard deviation:
sd(log(z)) = sqrt(var(log(x))+var(log(x))) = sqrt((sd(log(x)))^2+(sd(log(x)))^2)

exponentiate both sides:
exp(sd(log(z))) = exp(sqrt((sd(log(x)))^2+(sd(log(y)))^2))

Left hand side is geoSD of z: geoSdZ := exp(sd(log(z)))
geoSdZ = exp(sqrt((sd(log(x)))^2+(sd(log(y)))^2))

With geoSdX := exp(sd(log(x))) and geoSdY := exp(mean(sd(y))):
geoSdZ = exp(sqrt((log(geoSdX))^2+(log(geoSdY))^2 ))

[sfrechen]

OK, that means that for every calculation of PK-Ratios between 2 SimulationSets: RE must check if the used population is the same or not and depending on that use the appropriate calculation.

YES!

[sfrechen]

CASE 1: same population is compared under different scenarios

Let's say we are interested in the ratio of the PK Parameter PKParam and the population is of size n:
Let PKParam_x be n PKParam values of scenario x, and PKParam_y be n PKParam values of the reference scenario y.

geoMean_PKParamR = geomean(PKParam_x/PKParam_y)
    = exp(mean(log(PKParam_x/PKParam_y)))

geoSd_PKParamR = geoSd(PKParam_x/PKParam_y)
    = exp(sd(log(PKParam_x/PKParam_y)))

geoCV_PKParamR = geoCV(PKParam_x/PKParam_y)
    = sqrt(exp((sd(log(PKParam_x/PKParam_y)))^2)-1)

arithMean_PKParamR = mean(PKParam_x/PKParam_y)

arithSD_PKParamR = sd(PKParam_x/PKParam_y)

arithCV_PKParamR = arithSD_PKParamR/arithMean_PKParamR 
    = sd(PKParam_x/PKParam_y)/mean(PKParam_x/PKParam_y)

quantile_q_PKParamR = quantile(PKParam_x/PKParam_y, q)

with q being the q th quantile between 0 and 1, i.e. 0.05, 0.25, 0.5, 0.75, 0.95

CASE 2: different populations (population 1 and 2) are compared

Let's say we are interested in the ratio of the PK Parameter PKParam and population 1 is of size n and population 2 is of size m:
Let PKParam_1 be n PKParam values of population 1, and PKParam_2 be m PKParam values of population 2.

Analytically (via First order Taylor expansion)

geoMean_PKParamR = exp(log(geoMean_PKParam_1)-log(geoMean_PKParam_2))
    = exp(log(exp(mean(log(PKParam_1))))-log(exp(mean(log(PKParam_2))))) 
    = exp(mean(log(PKParam_1))-mean(log(PKParam_2)))

geoSd_PKParamR = exp(sqrt((log(geoSd_PKParam_1))^2+(log(geoSd_PKParam_2))^2)) 
    = exp(sqrt((log(exp(sd(log(PKParam_1)))))^2+(log(exp(sd(log(PKParam_2)))))^2)) 
    = exp(sqrt((sd(log(PKParam_1)))^2+(sd(log(PKParam_2)))^2))

geoCV_PKParamR = sqrt(exp((log(geoSd_PKParamR))^2)-1)
    = sqrt(exp((log(exp(sqrt((sd(log(PKParam_1)))^2+(sd(log(PKParam_2)))^2))))^2)-1)
    = sqrt(exp((sd(log(PKParam_1)))^2+(sd(log(PKParam_2)))^2)-1)

arithMean_PKParamR = arithMean_PKParam1/arithMean_PKParam2
    = mean(PKParam_1)/mean(PKParam_2)

arithSD_PKParamR = sqrt((1/arithMean_PKParam_2)^2*(arithSd_PKParam_1)^2+(-arithMean_PKParam_1/(arithMean_PKParam_2)^2)^2*(arithSd_PKParam_2)^2) 
    = sqrt((1/mean(PKParam_2)*sd(PKParam_1))^2+(mean(PKParam_1)/(mean(PKParam_2))^2*sd(PKParam_2))^2)

arithCV_PKParamR = arithSD_PKParamR/arithMean_PKParamR

For quantiles, this is however not straight-forward, thus @abdullahhamadeh had in a brief chat about it the idea to solve the issue via MC sampling.

MC sampling

This method is needed for quantiles and might be even superior to first order Taylor expansion, in particular for arithMean and arithSD.

Let B the number of replicate samples (e.g. 10000).
Then, we create B samples of PKParam_1_ii and PKParam_1_ii (ii = 1,...,B) each of size min(n,m) and calculate all desired statistics using the equations from CASE 1 (see above).
Finally, we take the median over the distributions of the statistics of PKParamR.

Exemplary R code:

set.seed(123456) # for reproducibility
B <- 10000
PKParamR <- data.frame(geoMean=numeric(),geoSd=numeric(),geoCV=numeric(),
                       arithMean=numeric(),arithSD=numeric(),arithCV=numeric(),
                       q05=numeric(),q25=numeric(),q50=numeric(),q75=numeric(),q95=numeric())
sampleSize <- min(c(n,m))
for(ii in 1:B){
  s_1 <- sample(x = PKParam_1,size = sampleSize,replace = FALSE) # (sub-)sample of population 1
  s_2 <- sample(x = PKParam_2,size = sampleSize,replace = FALSE) # (sub-)sample of population 2
  PKParamR[ii,"geoMean"] <- exp(mean(log(s_1/s_2)))
  PKParamR[ii,"geoSd"] <- exp(sd(log(s_1/s_2)))
  PKParamR[ii,"geoCV"] <- sqrt(exp((sd(log(s_1/s_2)))^2)-1)
  PKParamR[ii,"arithMean"] <- mean(s_1/s_2)
  PKParamR[ii,"arithSD"] <- sd(s_1/s_2)
  PKParamR[ii,"arithCV"] <- sd(s_1/s_2)/mean(s_1/s_2)
  PKParamR[ii,c("q05","q25","q50","q75","q95")] <- quantile(s_1/s_2,c(0.05,0.25,0.5,0.75,0.95))
}
apply(PKParamR,MARGIN = 2,median)

Are we ready to implement? Any objections against the MC sampling method?

@Yuri05 @abdullahhamadeh @pchelle @msevestre @PavelBal @AnnikaRPS @Christoph27

[sfrechen]

SUMMARY

After some discussione with @AnnikaRPS, we here now present the desired implementation.
Sorry for any confusions.

Population workflow: PARALLEL COMPARISON (with 1 reference population)

Let's say we are interested in PK Parameter PKParam and the populations are of size n1, n2, n3.
Let Population 1 be our reference population.

The following tables should be reported:

Table 1: PK_Param - Population Statisctics

Population	N	5th percentile	25th percentile	50th percentile	75th percentile	95th percentile	arith mean	arith standard deviation	arith CV	geo mean	geo standard deviation	geo CV
Population 1	n1
Population 2	n2
Population 3	n3

Table 2: PK_Param - Population Statistics compared to the reference: Population 1 (of size n1)

Population	N	5th percentile	25th percentile	50th percentile	75th percentile	95th percentile	arith mean	arith standard deviation	arith CV	geo mean	geo standard deviation	geo CV
Population 2	n2
Population 3	n3

In Table 1 the statistics reflect simply the given statistic over the respective PKParam values in each population,
i.e. for Population 1 :

• let PKParam_1 be n1 PKParam values of population 1

geoMean_PKParam = exp(mean(log(PKParam_1)))

geoSd_PKParam = exp(sd(log(PKParam_1)))

geoCV_PKParam = sqrt(exp((log(sd(PKParam_1)))^2)-1)

arithMean_PKParam = mean(PKParam_1)

arithSD_PKParam = sd(PKParam_1)

arithCV_PKParam = sd(PKParam_1)/mean(PKParam_1)

percentile_PKParam = quantile(PKParam_1,q)

with q being the q th quantile between 0 and 1, i.e. 0.05, 0.25, 0.5, 0.75, 0.95

In Table 2 the statistics reflect the relative change in the respective statistic,
i.e. for Population 2 in comparison to Population 1:

• let PKParam_1 be n1 PKParam values of population 1
• let PKParam_2 be n2 PKParam values of population 2

geoMean_PKParam_relChange = geoMean_PKParam_Population2/geoMean_PKParam_Population1
    = exp(mean(log(PKParam_2)))/exp(mean(log(PKParam_1)))

geoSd_PKParam_relChange = geoSd_PKParam_Population2/geoSd_PKParam_Population1
    = exp(sd(log(PKParam_2)))/exp(sd(log(PKParam_1)))

geoCV_PKParam_relChange = geoCV_PKParam_Population2/geoCV_PKParam_Population1
    = sqrt(exp((log(sd(PKParam_2)))^2)-1)/sqrt(exp((log(sd(PKParam_1)))^2)-1)

arithMean_PKParam_relChange = arithMean_PKParam_Population2/arithMean_PKParam_Population1
    = mean(PKParam_2)/mean(PKParam_1)

arithSD_PKParam_relChange = arithSD_PKParam_Population2/arithSD_PKParam_Population1
    = sd(PKParam_2)/sd(PKParam_1)

arithCV_PKParam_relChange = arithCV_PKParam_Population2/arithCV_PKParam_Population1
    = (sd(PKParam_2)/mean(PKParam_2))/(sd(PKParam_1)/mean(PKParam_1))

percentile_PKParam_relChange = percentile_PKParam_Population2/percentile_PKParam_Population1
    = quantile(PKParam_2, q)/quantile(PKParam_1,q)

with q being the q th quantile between 0 and 1, i.e. 0.05, 0.25, 0.5, 0.75, 0.95

Thus, the comparisons of the statistics itself becomes a standard output (not only for ratio comparison workflows).

Population workflow: RATIO COMPARISON

The tables Table 1 and Table 2 from the "Population workflow: PARALLEL COMPARISON (with 1 reference population)" as descibed above should also be part of the output for this workflow.
Additionally, we now provide statistics on the RATIO PKParamR of the PK parameter PKParam in an additional Table 3.

CASE 1: same population is compared under different scenarios

Let's say we are interested in the ratio PKParamR of the PK Parameter PKParam and the population is of size n.
Let Population Scenario x be our reference population.

Table 3: PK_Param Ratio - Population Statistics on the ratio of PK_Param with the reference: Scenario x (of size n).

Population	N	5th percentile	25th percentile	50th percentile	75th percentile	95th percentile	arith mean	arith standard deviation	arith CV	geo mean	geo standard deviation	geo CV
Scenario y	n
Scenario z	n

In Table 3 the statistics are calculated as follows on the "individual" ratios,
e.g. for scenario y in comparison to scenario x:

• let PKParam_x be n PKParam values of scenario x,
• let PKParam_y be n PKParam values of scenario y,

geoMean_PKParamR = geomean(PKParam_y/PKParam_x)
    = exp(mean(log(PKParam_y/PKParam_x)))

geoSd_PKParamR = geoSd(PKParam_y/PKParam_x)
    = exp(sd(log(PKParam_y/PKParam_x)))

geoCV_PKParamR = geoCV(PKParam_y/PKParam_x)
    = sqrt(exp((sd(log(PKParam_y/PKParam_x)))^2)-1)

arithMean_PKParamR = mean(PKParam_y/PKParam_x)

arithSD_PKParamR = sd(PKParam_y/PKParam_x)

arithCV_PKParamR = arithSD_PKParamR/arithMean_PKParamR 
    = sd(PKParam_y/PKParam_x)/mean(PKParam_y/PKParam_x)

quantile_q_PKParamR = quantile(PKParam_y/PKParam_x, q)

with q being the q th quantile between 0 and 1, i.e. 0.05, 0.25, 0.5, 0.75, 0.95

CASE 2: different populations are compared

Let's say we are interested in the ratio of the PK Parameter PKParam and the populations are of size n1, n2, n3.
Let Population 1 be our reference population.

Table 3: PK_Param Ratio - Population Statistics on the ratio of PK_Param with the reference: Population 1 (of size n1).

Population	N	5th percentile	25th percentile	50th percentile	75th percentile	95th percentile	arith mean	arith standard deviation	arith CV	geo mean	geo standard deviation	geo CV
Population 2	n2
Population 3	n3

Now, in Table 3 the statistics are calculated as follows,
e.g. for Population 2 in comparison to Population 1:

• let PKParam_1 be n1 PKParam values of population 1
• let PKParam_2 be n2 PKParam values of population 2

For the geometric statistics, well-defined analytical solutions exist:

geoMean_PKParamR = exp(log(geoMean_PKParam_2)-log(geoMean_PKParam_1))
    = exp(log(exp(mean(log(PKParam_2))))-log(exp(mean(log(PKParam_1))))) 
    = exp(mean(log(PKParam_2))-mean(log(PKParam_1)))

geoSd_PKParamR = exp(sqrt((log(geoSd_PKParam_2))^2+(log(geoSd_PKParam_1))^2)) 
    = exp(sqrt((log(exp(sd(log(PKParam_2)))))^2+(log(exp(sd(log(PKParam_1)))))^2)) 
    = exp(sqrt((sd(log(PKParam_2)))^2+(sd(log(PKParam_1)))^2))

geoCV_PKParamR = sqrt(exp((log(geoSd_PKParamR))^2)-1)
    = sqrt(exp((log(exp(sqrt((sd(log(PKParam_2)))^2+(sd(log(PKParam_1)))^2))))^2)-1)
    = sqrt(exp((sd(log(PKParam_2)))^2+(sd(log(PKParam_1)))^2)-1)

For the arithmetic and percentile statistics, we approximate the ratio statistics via MC Sampling:

Let B the number of replicate samples (e.g. 10000).
Then, we create B samples of PKParam_1_ii and PKParam_2_ii (ii = 1,...,B) each of size min(n1,n2) and calculate all desired statistics using the equations from CASE 1 (see above).
Finally, we take the median over the distributions of the statistics of PKParamR.

Exemplary R code:

set.seed(123456) # for reproducibility
B <- 10000
PKParamR_temp <- data.frame(q05=numeric(),q25=numeric(),q50=numeric(),q75=numeric(),q95=numeric(),
                            arithMean=numeric(),arithSD=numeric(),arithCV=numeric())
sampleSize <- min(c(n1,n2))
for(ii in 1:B){
  s_1 <- sample(x = PKParam_1,size = sampleSize,replace = FALSE) # (sub-)sample of population 1
  s_2 <- sample(x = PKParam_2,size = sampleSize,replace = FALSE) # (sub-)sample of population 2
  PKParamR_temp[ii,"arithMean"] <- mean(s_2/s_1)
  PKParamR_temp[ii,"arithSd"] <- sd(s_2/s_1)
  PKParamR_temp[ii,"arithCV"] <- sd(s_2/s_1)/mean(s_2/s_1)
  PKParamR_temp[ii,c("q05","q25","q50","q75","q95")] <- quantile(s_2/s_1,c(0.05,0.25,0.5,0.75,0.95))
}
PKParamR1 <- apply(PKParamR_temp,MARGIN = 2,median)
PKParamR2 <- data.frame(geoMean=exp(mean(log(PKParam_2))-mean(log(PKParam_1))),
                        geoSd=exp(sqrt((sd(log(PKParam_2)))^2+(sd(log(PKParam_1)))^2)),
                        geoCV=sqrt(exp((sd(log(PKParam_2)))^2+(sd(log(PKParam_1)))^2)-1))
PKParamR <- data.frame(c(PKParamR1,PKParamR2))
print(PKParamR)

@Yuri05 @abdullahhamadeh @pchelle @msevestre @PavelBal @AnnikaRPS @Christoph27

[AnnikaRPS]

Looks good, no objections from my side.

For pediatric workflows, parallel comparison (with 1 reference population) would be the default

[sfrechen]

@Yuri05 : not only a feature. For workflow "ratioComparison", it is a bug. The current results are not at all what would be expected by any user.

[pchelle]

@sfrechen
I am currently implementing the feature and have few questions:

1- How do you assert that 2 populations are "the same" ?
So far, I am assessing whether the population and study design files are identical. Is this okay ?

2- In the report, do you want the analytical solution for geomean, geosd, etc. or would Monte Carlo solution be better ?
Regardless, I intend to calculate both and include them in the debug logs (this way users can check if both are consistent)

3- I will add 2 settings for the Monte Carlo sampling in the task
- The number of repetitions: plotPKParameters$settings$mcRepetitions
- The Random Seed for reproducibility: plotPKParameters$settings$mcRandomSeed
Would this be okay ?

[sfrechen]

@pchelle @Yuri05

1- How do you assert that 2 populations are "the same" ?
So far, I am assessing whether the population and study design files are identical. Is this okay ?

Ideally, we would check some kind of checksum / hash of the population file. Sometimes (due to automization) we have different file names with the same content...
(Study design check should not be implemented. Or is it possible to modify population properties via study design?)

2- In the report, do you want the analytical solution for geomean, geosd, etc. or would Monte Carlo solution be better ?
Regardless, I intend to calculate both and include them in the debug logs (this way users can check if both are consistent)

Great. Monte Carlo solution should be default.

3- I will add 2 settings for the Monte Carlo sampling in the task

• The number of repetitions: plotPKParameters$settings$mcRepetitions
• The Random Seed for reproducibility: plotPKParameters$settings$mcRandomSeed
  Would this be okay ?

Sounds very reasonable. Great!

[Yuri05]

@sfrechen

Ideally, we would check some kind of checksum / hash of the population file. Sometimes (due to automization) we have different file names with the same content...

Can you give an example of a use case, where this could happen? 2 different population files with de-facto the same population?

Or is it possible to modify population properties via study design?

Yes, sure :)

[sfrechen]

Ideally, we would check some kind of checksum / hash of the population file. Sometimes (due to automization) we have different file names with the same content...

We have plenty of them (I would say)

[Yuri05]

We have plenty of them (I would say)

Can you give some concrete examples? For me it sounds more like a user error.
But if we really need such a comparison, I would not compare any files (which could be different in the 20th digit somewhere, e.g. if the same population was exported on different computers etc.).

Instead I would e.g. load populations into dataframes and compare them numerically.
And if it turns out, that 2 populations (from 2 different files) are in fact the same - I would produce a warning.

[sfrechen]

Can you give some concrete examples? For me it sounds more like a user error.

No, no use error. Concretely: We have pool of virtual individuals. Then a script generates the different work packages with certain instructions which individuals to pull from this pool, creates the pkml files etc. Finally, the script saves everything in folders and gives matching names to poulation and pkml file...