Low-density lipoprotein analysis
In the following paragraphs we show a real life healthcare example for LDL monitoring. Two approaches will be presented: one with and one without taking the standard deviation into account. The aim is to minimise the healthcare burden generated by patients with high cardiovascular (CV) event risk. The model is built upon the relationship between the low-density lipoprotein level and the risk of CV events, thus the LDL level is the process of interest (Boekholdt et al. 2012).
Parameters were estimated from several sources. The table below gives information about the meaning of the parameter values in the healthcare setting and shows the source of the parameter estimations.| Parameter value | Meaning | Parameter source |
|---|---|---|
| 3 mmol/l | Target value | Set according to the European guideline for patients at risk |
| 0.1 mmol/l | Process standard deviation | Estimated using real life data from Hungary, namely registry data through Healthware Consulting Ltd. |
| 0.8/3 | Expected shift size, 0.8 increase in LDL per year on average | Estimated with the help of a health professional |
| 1/120 | Expected number of shifts in a day, 3 shifts per year on average | Estimated with the help of a health professional |
| 0.027, 1.15 | Parameters of the repair size beta distribution | Estimated using an international study which included Hungary |
| 0.1, 30 | Parameters of the sampling probability logistic function | Patient non-compliance in LDL controlling medicine is quite high, and this is represented through the parametrisation of the logistic function |
| 5.01 EUR | Sampling cost | Estimated using the LDL testing cost and visit cost in Hungary |
| 4.60 EUR | Shift-proportional daily out-of-control (OOC) cost | Estimated using real world data of cardiovascular event costs from Hungary |
| 9.97 EUR | Base repair cost | Estimated using the simvastatin therapy costs in Hungary |
| 7.48 EUR | Shift-proportional repair cost | Estimated using the simvastatin therapy costs in Hungary |
It is very difficult to give a good estimate for the type and parameters of a distribution that properly models the non-compliance, thus the results here can at best be regarded as close approximations to a real life situation. This is not a limiting factor, as patients themselves can have vast differences in their behaviour, so evaluation of different scenarios are often required, and will also be presented here. Since high LDL levels rarely produce noticeable symptoms, the sampling probability only depends on the time between samplings (Institute for Quality and Efficiency in Health Care 2017). Thus, the sampling probability was modelled by the logistic function and not by the beta distribution (see the help page of the Markovstat function for more details). It is important to note that the proportional costs increase according to a Taguchi-type (squared) loss function, thus huge, expenses can be generated if the patient’s health is highly out-of-control.
Optimisation using only the cost expectation
The LDL aggregated dataset is included in the package with all of the above-mentioned parameter estimations. For cost estimation and optimisation purposes we need to use two function: Markovstat for stationary distribution estimation and Markovchart for cost calculation. The code below feeds the proper parameters to these functions. Notice that exponential shift size distribution, random repair and sampling is used and Markovchart is set for optimisation.
data("LDL")
stat_LDL_cost <- Markovstat(
shiftfun = 'exp', h = 50, k = 0.15,
sigma = LDL$standard_deviation, s = LDL$num_exp_daily_shifts,
delta = LDL$expected_shift_size,
RanRep = TRUE, alpha = LDL$rep_size_first, beta = LDL$rep_size_second,
RanSam = TRUE, q = LDL$samp_prob_first, z = LDL$samp_prob_second,
Vd = 50, V = 3)
#Defining parallel_opt parallel settings.
#parallel_opt can also be left empty to be defined automatically by the function.
num_workers <- min(c(detectCores(),2))
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res_LDL_cost <- Markovchart(
statdist = stat_LDL_cost,
OPTIM=TRUE, p = 1,
cs = LDL$sampling_cost,
cf = LDL$base_rep_cost,
coparams = c(0,LDL$OOC_cost),
crparams = c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
res_LDL_cost
#> $Results
#> G-value Expected cost Total cost sd Cost sd due to process var.
#> 1 0.4258544 0.4258544 0.5041891 0.5041891
#> 2nd process moment 3rd process moment 4th process moment
#> 1 0.4355586 3.703977 77.98069
#>
#> $Subcosts
#> Sampling cost Repair cost OOC cost
#> 1 0.08998974 0.08023143 0.262028
#>
#> $Parameters
#> Time between samplings (h) Critical value (k)
#> 1 55.70492 0.1583099The optimal parameters for the case when the cost standard deviation was not taken into account were 55.7 days and 0.158 mmol/l for the time between samplings and the critical increase in the LDL level from the guideline value respectively. These parameters entailed an average daily cost of 0.43€ and standard deviation of 0.5€. This result is interesting, because according to the optimisation we should use a somewhat higher critical value than the one in the guideline – 0 mmol/l critical increase would be the 3 mmol/l value in the guideline. At the same time patients should be monitored more often than usual as times of the LDL measurements are usually several months or years apart. However, this is a strictly cost effective viewpoint which could be overwritten by a health professional. Nonetheless, the results provide a possible new approach to the therapeutic regime for controlling LDL level. Often, it is good to look at the interaction between the parameters and the resulting average cost, especially in situations where the optimal parameters cannot be used because of real life reasons.
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res_LDL_cost_grid <- Markovchart(
statdist = stat_LDL_cost,
h=seq(50,75,2.5),
k=seq(0.05,0.25,0.02),
p = 1,
cs = LDL$sampling_cost,
cf = LDL$base_rep_cost,
coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
plot(res_LDL_cost_grid,
y = 'Expected \ndaily cost',
mid = '#ff9494',
high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical LDL increase') +
geom_point(aes(x = res_LDL_cost$Parameters[[1]],
y = res_LDL_cost$Parameters[[2]]))
The figure shows the expected cost as function of the time between samplings and the critical value.
The heat map shows the average cost as the function of the different parameter values. The dot in the lightest area of Figure \(\ref{expectedcost}\) corresponds to the optimal cost. Any other point entails a higher average cost. It can be seen that too low or high critical values will both increase the average daily cost. Note that the change in time between samplings entails relatively low change in the critical LDL increase: even if the time between control visits is changed the critical value should stay around 0.12 - 0.18 mmol/l.
Optimisation using cost expectation and cost standard deviation
The code below is largely the same as before, but notice that p=0.9 instead of p=1, thus prompting the 10% use of the standard deviation during the optimisation process.
stat_LDL_cost <- Markovstat(
shiftfun = 'exp', h = 50, k = 0.15,
sigma = LDL$standard_deviation, s = LDL$num_exp_daily_shifts,
delta = LDL$expected_shift_size,
RanRep = TRUE, alpha = LDL$rep_size_first, beta = LDL$rep_size_second,
RanSam=TRUE, q = LDL$samp_prob_first, z = LDL$samp_prob_second,
Vd = 50, V = 3)
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res_LDL_G <- Markovchart(
statdist = stat_LDL_cost,
OPTIM=TRUE, p = 0.9,
cs = LDL$sampling_cost,
cf = LDL$base_rep_cost,
coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
res_LDL_G
#> $Results
#> G-value Expected cost Total cost sd Cost sd due to process var.
#> 1 0.4266033 0.4323714 0.3746903 0.3746903
#> 2nd process moment 3rd process moment 4th process moment
#> 1 0.3273378 2.038361 42.26514
#>
#> $Subcosts
#> Sampling cost Repair cost OOC cost
#> 1 0.07838653 0.08099193 0.2755366
#>
#> $Parameters
#> Time between samplings (h) Critical value (k)
#> 1 63.95067 0.145161In this part, the cost standard deviation is also taken into account as p=0.9, thus the weight of the standard deviation in the calculation of \(G\) is 0.1. The optimal parameters found by our approach were 63.95 days and 0.145 mmol/l for the time between samplings and critical increase in the LDL level respectively. These parameters entailed an average daily cost of 0.43€ and standard deviation of 0.37€. The inclusion of the cost standard deviation into the model has somewhat increased the time between control visits and decreased the critical value. The expected cost increased, while the cost standard deviation has moderately decreased. The figure below shows the previous heat map with non-compliance included in the model.
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res_LDL_G_grid <- Markovchart(
statdist = stat_LDL_cost,
h=seq(50,75,2.5),
k=seq(0.05,0.25,0.02),
p = 0.9,
cs = LDL$sampling_cost,
cf = LDL$base_rep_cost,
coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
plot(res_LDL_G_grid,
y = 'Expected \ndaily cost',
mid = '#ff9494',
high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical LDL increase',
nbreaks = 14) +
geom_point(aes(x = res_LDL_G$Parameters[[1]],
y = res_LDL_G$Parameters[[2]]))
The figure shows \(G\) value (with p=0.9) as function of the time between samplings and the critical value. It can be seen that the elliptical shape of the heat map has not changed: the change in the time between control visits still does not entail great change in the critical value.
Sensitivity Analysis
As there were uncertainties about the estimation of several parameters, it is important to assess the effect of different parameter setups. The results for different out-of-control costs are plotted for both approaches. The results can be seen in the figure below. The code runs the optimisation repeatedly to collect the necessary data and then gathers the results into a ggplot-compatible data format.
results <- NULL
statds <- NULL
for (i in 3:10)
{
for (j in c(0.9,1))
{
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res_LDL_optim <- Markovchart(
statdist = stat_LDL_cost,
OPTIM = TRUE,
p = j,
cs = LDL$sampling_cost,
cf = LDL$base_rep_cost,
coparams = c(0,i),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
results <- rbind(results,c(j,i,unlist(c(res_LDL_optim$Results[c(2,3)],res_LDL_optim$Parameters))))
statds <- rbind(statds,res_LDL_optim$Stationary_distribution)
}
}
results <- resultsbck <- as.data.frame(results)
colnames(results) <- c("p","co","EC","SDC","h","k")
results$p <- as.character(results$p)
results <- cbind(results[,1:2],melt(results[,3:6]))
results$variable <- as.character(results$variable)
results$variable[results$variable=="h"] <- "Days between samplings"
results$variable[results$variable=="k"] <- "Critical LDL increase"
results$variable[results$variable=="EC"] <- "Expected cost"
results$variable[results$variable=="SDC"] <- "Cost standard deviation"
results$variable = factor(results$variable, levels=c("Critical LDL increase","Days between samplings","Expected cost","Cost standard deviation"))
results$min <- NA
results$min[results$variable=="Days between samplings"] <- min(results$value[results$variable=="Days between samplings"])
results$max[results$variable=="Days between samplings"] <- max(results$value[results$variable=="Days between samplings"])
results$min[results$variable=="Critical LDL increase"] <- min(results$value[results$variable=="Critical LDL increase"])
results$max[results$variable=="Critical LDL increase"] <- max(results$value[results$variable=="Critical LDL increase"])
results$min[results$variable=="Expected cost"] <- min(results$value[results$variable=="Expected cost"])
results$max[results$variable=="Expected cost"] <- max(results$value[results$variable=="Expected cost"])
results$min[results$variable=="Cost standard deviation"] <- min(results$value[results$variable=="Cost standard deviation"])
results$max[results$variable=="Cost standard deviation"] <- max(results$value[results$variable=="Cost standard deviation"])
app_LDL_params_plot <- ggplot(results, aes(x=co, y=value, group=p)) +
geom_line(aes(colour=p),size = 1.1) +
facet_wrap(~variable, labeller = label_bquote(rows=.(variable)), scales="free_y", nrow=2) +
scale_colour_manual(expression(italic(p)),values=c("black","#800000")) +
geom_blank(aes(y = min)) +
geom_blank(aes(y = max)) +
xlab("Out-of-control cost") +
ylab("Value") +
theme_bw() +
theme(strip.text.x = element_text(size = 11),
strip.text.y = element_text(size = 11,angle = 0), legend.position="top")
app_LDL_params_plot
The figure shows parameters, average total cost and cost standard deviation as function of the out-of-control cost. The critical value and time between samplings decrease and the average cost and cost standard deviation increase with higher out-of-control costs. Just as on the heat maps, one can observe here, that if the cost standard deviation is taken into account in the optimisation, then the critical value should be lowered and the time between samplings increased. It can be observed that a substantial decrease can be achieved in the cost standard deviation while the cost expectation barely changes.
