Combined residual unexplained variability (RUV) implementation

[marianklose]

Hi everyone, thanks for developing PKPDposterior and PKPDsim. Both are great tools!

I have a minor concern regarding the implementation of the combined RUV error model. I am not entirely sure if my concern is justified but I'll just put it up for discussion.

The following will be based on the vancomycin example page from the docs. It relies on the Thomson et al. model (pub). They report a standard deviation for the additive RUV part of 1.6 mg/mL and 15% for the proportional RUV part. Thus, in the example the error terms are defined as list(prop = 0.15, add = 1.6) within

data <- new_stan_data(
  regimen,
  covariates, 
  tdm_data,
  dose_cmt = 2,
  parameters = parameters,
  fixed = c("TH_CRCL", "TDM_INIT"),
  iiv = iiv,
  ruv = list(
    prop = 0.15,
    add = 1.6
  ),
  ltbs = FALSE
)

and the resulting list object data has stored:

$ruv
$ruv$prop
[1] 0.15

$ruv$add
[1] 1.6

(...)

$stan_data$ruv_prop_pk
[1] 0.15

$stan_data$ruv_add_pk
[1] 1.6

The resulting stan model file (generated with write_stan_model()) directly takes the ruv_prop_pk and ruv_add_pk as input in the data{} block:

// error_parameters
int<lower = 0, upper = 1> ltbs_pk;
real<lower = 0> ruv_prop_pk;
real<lower = 0> ruv_add_pk;

There are no further manipulations or calculations involved and we directly jump to the likelihood statement and the simulated posterior distributions:

// likelihood_observed_data
if(ltbs_pk) {
  log_dv_pk ~ normal(log(ipred_obs_pk), ruv_add_pk); 
} else {
  dv_pk ~ normal(ipred_obs_pk, (ruv_prop_pk * ipred_obs_pk + ruv_add_pk));
}

(...)

// simulate_posterior_ruv
real ipred_ruv_pk[n_obs_pk];
for(i in 1:n_obs_pk){
  ipred_ruv_pk[i] = normal_rng(ipred_obs_pk[i], (ruv_prop_pk * ipred_obs_pk[i] + ruv_add_pk));
}

Since ruv_prop_pk is a CV in decimal units we should get a standard deviation for ruv_prop_pk * ipred_obs_pk. And ruv_add_pk is already a standard deviation. For the likelihood and sampling statements stan expects the standard deviation of the respective distribution (ref). So we need to find the overall standard deviation of the normal distribution at the respective concentration of ipred_obs_pk. However, variance, not standard deviation, is additive. Thus, to my understanding simply adding the proportional standard deviation (= ruv_prop_pk * ipred_obs_pk) and the additive standard deviation ruv_add_pk to get the joint standard deviation is invalid. We need to first transform the SD's back to variances, add them together and then transform the resulting variance back to the SD domain. Instead of

$$\sigma = cv \cdot ipred + add$$

we would need to calculate the standard deviation by

$$\sigma = \sqrt{(cv \cdot ipred)^2 + (add)^2}$$

I was initially not sure how big the difference between both calculations are so I perfomed a small simulation with varying magnitudes of additional and proportional errors:

# define additive implementation
add_fun <- function(prop, add, conc) {
  res <- prop * conc + add
  return(res)
}

# define square root implementation
sqrt_fun <- function(prop, add, conc) {
  res <- sqrt((prop*conc)^2+add^2)
  return(res)
}

# generate grid
prop <- seq(0, 1, 0.05)
add <- seq(0, 10, 0.5)
conc <- 100
grid <- expand.grid(prop=prop, add=add, conc=conc) |> as.tibble()

# delete 0s
grid <- grid |> 
  dplyr::filter(prop != 0 & add != 0)

# apply functions to grid
grid$add_fun <- apply(grid, 1, function(row) add_fun(row["prop"], row["add"], row["conc"]))
grid$sqrt_fun <- apply(grid, 1, function(row) sqrt_fun(row["prop"], row["add"], row["conc"]))

# calculate differences
grid$diff <- grid$add_fun - grid$sqrt_fun
grid$rel_diff <- (grid$diff / grid$sqrt_fun) * 100

# plot
plot_ly(x=grid$prop, y=grid$add, z=grid$rel_diff, type="scatter3d", mode="markers", color=grid$rel_diff) |> 
  layout(scene = list(xaxis = list(title = 'Proportional RUV [decimal]'),
                      yaxis = list(title = 'Additive RUV [mg/L]'),
                      zaxis = list(title = 'Difference [%]')))

The relative difference between both calculation methods ranges somewhere between 0.5 % and 41 % depending on the magnitude of prop and add:

The current implementation leads to higher standard deviations than the sqrt method, thus the MCMC sampler would be more prior driven and less influenced by the likelihood/data.

Again, I don't want to rule out that I just have a misunderstanding somewhere or that I may not have fully grasped the implementation in the code. Happy to hear your opinion on this.

[marianklose]

I have also did a small comparison (just out of interest) between both calculation methods for the vancomycin example using 20000 MCMC draws:

You can see some smaller shifts for CL and V1. The resulting predictions of the pk model are also slightly affected:

For this example the difference seems to be not that huge, but surely it depends on all those parameters and observations.