# NOT RUN {
data("diabetes")
str(diabetes)
# }
# NOT RUN {
##### Example of data assessment for control chart use #####
### Packages
require(zoo)
require(reshape2)
RANDOMISED_DATA <- diabetes
### Functions
weighted.var <- function(x, w, na.rm = FALSE) {
if (na.rm) {
w <- w[i <- !is.na(x)]
x <- x[i]
}
sum.w <- sum(w)
sum.w2 <- sum(w^2)
mean.w <- sum(x * w) / sum(w)
(sum.w / (sum.w^2 - sum.w2)) * sum(w * (x - mean.w)^2, na.rm =
na.rm)
}
estBetaParams <- function(mu, var) {
alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
beta <- alpha * (1 / mu - 1)
return(params = list(alpha = alpha, beta = beta))
}
### Setting up data
# Way too high HbA1C levels are discarded as outliers
RANDOMISED_DATA$HBA1C_AVG[RANDOMISED_DATA$HBA1C_AVG>20 & !is.na(RANDOMISED_DATA$HBA1C_AVG)] <- NA
# Lowest HbA1c level taken into account
lowest <- 4
### Gathering data for several estimates
RANDOMISED_DATA <- RANDOMISED_DATA[RANDOMISED_DATA$ID %in%
RANDOMISED_DATA$ID[grepl("E11",RANDOMISED_DATA$ICD)],]
# Process standard deviation
sigma_param <- sigma <- sqrt(weighted.mean((RANDOMISED_DATA$HBA1C_SD[
RANDOMISED_DATA$SAMPLING_IN_MONTH>=2 & !is.na(RANDOMISED_DATA$SAMPLING_IN_MONTH)])^2,
RANDOMISED_DATA$SAMPLING_IN_MONTH[RANDOMISED_DATA$SAMPLING_IN_MONTH>=2 &
!is.na(RANDOMISED_DATA$SAMPLING_IN_MONTH)]))
IDLIST <- unique(RANDOMISED_DATA$ID[!is.na(RANDOMISED_DATA$HBA1C_AVG)][
duplicated(RANDOMISED_DATA$ID[!is.na(RANDOMISED_DATA$HBA1C_AVG)])])
IDLIST <- unique(RANDOMISED_DATA$ID[(RANDOMISED_DATA$ID %in% IDLIST) & RANDOMISED_DATA$AGE>39])
shiftdat <- NULL
stimedat <- NULL
repaidat <- NULL
deltats <- NULL
deltaATC <- NULL
for(i in IDLIST)
{
smalldat <- RANDOMISED_DATA[RANDOMISED_DATA$ID==i,c("DATE","HBA1C_AVG","THERAPY_VECTOR")]
smalldat <- smalldat[!is.na(smalldat$DATE) & !is.na(smalldat$HBA1C_AVG),]
patshiftdat <- as.data.frame(cbind(i,smalldat$DATE[2:dim(smalldat)[1]],diff(smalldat$DATE),
diff(smalldat$HBA1C_AVG))[diff(smalldat$HBA1C_AVG)>2*sigma,,drop=FALSE])
if(dim(patshiftdat)[1]>1) stimedat <- rbind(stimedat,cbind(i,diff(as.Date(patshiftdat$V2))))
patrepaidat <- as.data.frame(cbind(i,diff(smalldat$DATE),(smalldat$HBA1C_AVG-lowest)[2:
length(smalldat$HBA1C_AVG)]/(smalldat$HBA1C_AVG-lowest)[1:(length(smalldat$HBA1C_AVG)-1)],
as.character(smalldat$THERAPY_VECTOR[1:(length(smalldat$THERAPY_VECTOR)-1)]))[
(which(diff(smalldat$HBA1C_AVG)<(-2*sigma) &
smalldat$HBA1C_AVG[1:(length(smalldat$HBA1C_AVG)-1)]>6 &
smalldat$HBA1C_AVG[1:(length(smalldat$HBA1C_AVG)-1)]<=20)),,drop=FALSE])
shiftdat <- rbind(shiftdat,patshiftdat)
repaidat <- rbind(repaidat,patrepaidat)
deltats <- rbind(deltats,cbind(i,diff(as.Date(RANDOMISED_DATA$DATE[
!is.na(RANDOMISED_DATA$HBA1C_AVG) & RANDOMISED_DATA$ID==i]))))
try(deltaATC <- rbind(deltaATC,cbind(i,diff(as.Date(RANDOMISED_DATA$DATE[
!is.na(RANDOMISED_DATA$THERAPY) & RANDOMISED_DATA$ID==i])))), silent=TRUE)
}
colnames(shiftdat) <- c("ID","TIME","TIMEDIFF","SHIFTSIZE")
colnames(deltats) <- c("ID","DeltaT")
colnames(deltaATC) <- c("ID","deltaATC")
# delta: expected shift size
delta_param <- mean(shiftdat$SHIFTSIZE[shiftdat$TIMEDIFF>=90 & shiftdat$TIMEDIFF<184])
# Empirical distribution of elapsed times (between samplings)
summary(deltats[,2])
mean(deltats[,2])
median(deltats[,2])
sd(deltats[,2])
# s: expected number of shifts per unit time
stimedat <- as.data.frame(stimedat)
colnames(stimedat) <- c("ID","TIMEDIFF")
s_param <- 1/mean(stimedat$TIMEDIFF[stimedat$TIMEDIFF<367])
# alphas, betas: therapy effectiveness parameters
colnames(repaidat) <- c("ID","TIMEDIFF","REMAIN","THERAP")
repaidat$REMAIN <- as.numeric(as.character(repaidat$REMAIN))
repaidat$TIMEDIFF <- as.numeric(as.character(repaidat$TIMEDIFF))
repaidat$THERAP <- as.character(repaidat$THERAP)
repaidat <- repaidat[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184,]
repaidat$REMAIN[repaidat$REMAIN<0] <- 0
ther_eff <- as.data.frame(rbind(
cbind("ANALOGUE",repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("ANALOGUE",repaidat$THERAP) & !grepl("-H",repaidat$THERAP) & !grepl("GLP",repaidat$THERAP)]),
cbind("GLP",repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("GLP",repaidat$THERAP) & !grepl("-H",repaidat$THERAP)])))
ther_eff[,1] <- factor(ther_eff[,1], levels = c("ANALOGUE", "GLP"))
ther_eff[,2] <- as.numeric(as.character(ther_eff[,2]))
colnames(ther_eff) <- c("Type","Effectiveness")
ANALOGUE <- estBetaParams(mean(repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("ANALOGUE",repaidat$THERAP) & !grepl("-H",repaidat$THERAP) & !grepl("GLP",repaidat$THERAP)]),
var(repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("ANALOGUE",repaidat$THERAP) & !grepl("-H",repaidat$THERAP) &
!grepl("GLP",repaidat$THERAP)]))
GLP <- estBetaParams(mean(repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("GLP",repaidat$THERAP) & !grepl("-H",repaidat$THERAP)]),
var(repaidat$REMAIN[repaidat$TIMEDIFF>=90 & repaidat$TIMEDIFF<184 &
grepl("GLP",repaidat$THERAP) & !grepl("-H",repaidat$THERAP)]))
### Repair cost
HBA1C_fill <- NULL
for (i in unique(RANDOMISED_DATA$ID[!is.na(RANDOMISED_DATA$HBA1C_AVG)]))
{
vec <- RANDOMISED_DATA$HBA1C_AVG[RANDOMISED_DATA$ID==i]
vec[which(is.na(vec))[which(is.na(vec))<which(!is.na(vec))[1]]] <- vec[which(!is.na(vec))[1]]
vec[which(is.na(vec))[which(is.na(vec))>which(!is.na(vec))[length(which(!is.na(vec)))]]] <-
vec[which(!is.na(vec))[length(which(!is.na(vec)))]]
vec <- na.approx(vec)
HBA1C_fill <- rbind(HBA1C_fill,cbind(i,vec))
smaldat <- RANDOMISED_DATA[RANDOMISED_DATA$ID==i,]
smaldat$THERAPY_COST_EUR[smaldat$THERAPY_COST_EUR==0 & smaldat$THERAPY_VECTOR!=""] <- NA
if(is.na(smaldat$THERAPY_COST_EUR[1])) smaldat$THERAPY_COST_EUR[1] <- 0
new_burden <- na.locf(smaldat$THERAPY_COST_EUR)
seged <- cbind(rle(is.na(smaldat$THERAPY_COST_EUR))[[2]],
rle(is.na(smaldat$THERAPY_COST_EUR))[[1]])
seged[,2][seged[,1]==0] <- seged[,2][seged[,1]==0]-1
seged[,2][seged[,1]==1] <- seged[,2][seged[,1]==1]+1
if(seged[length(seged[,1]),1]==0) seged[length(seged[,2]),2] <- seged[length(seged[,2]),2]+1
seged2 <- cbind(rep(seged[,1], seged[,2]),rep(seged[,2], seged[,2]))
new_burden[seged2[,1]==1] <- new_burden[seged2[,1]==1]/seged2[,2][seged2[,1]==1]
RANDOMISED_DATA$THERAPY_COST_EUR[RANDOMISED_DATA$ID==i] <- new_burden
}
RANDOMISED_DATA$HBA1C_fill <- NA
RANDOMISED_DATA$HBA1C_fill[RANDOMISED_DATA$ID%in%HBA1C_fill[,1]] <- HBA1C_fill[,2]
RANDOMISED_DATA$HBA1C_fill_filter <- RANDOMISED_DATA$HBA1C_fill
RANDOMISED_DATA$HBA1C_fill_filter[RANDOMISED_DATA$HBA1C_fill_filter>=10] <- NA
discparam <- 150
cutHBA1C_AVG <- cut(na.omit(RANDOMISED_DATA$HBA1C_fill_filter),breaks=discparam)
newlvls <- seq(min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)),
max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)),
(max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter))-
min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)))/discparam)[1:discparam] +
(max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter))-
min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)))/discparam/2
levels(cutHBA1C_AVG) <- newlvls
costs <- cbind(as.numeric(as.character(cutHBA1C_AVG)),
RANDOMISED_DATA$THERAPY_COST_EUR[!is.na(
RANDOMISED_DATA$HBA1C_fill_filter)]/30,
as.character(RANDOMISED_DATA$THERAPY[
!is.na(RANDOMISED_DATA$HBA1C_fill_filter)]))
costs <- as.data.frame(costs)
colnames(costs) <- c("HBA1C","HC_BURDEN","THERAP")
costs$HBA1C <- as.numeric(as.character(costs$HBA1C))
costs$HC_BURDEN <- as.numeric(as.character(costs$HC_BURDEN))
costs$THERAP <- as.character(costs$THERAP)
costs$THERAP[grepl("ANALOGUE", costs$THERAP) & !grepl("GLP", costs$THERAP)] <- "ANALOGUE"
costs$THERAP[grepl("GLP",costs$THERAP)] <- "GLP"
cost.ANALOGUE <- as.data.frame(cbind(sort(unique(costs$HBA1C[costs$THERAP=="ANALOGUE"])),
as.numeric(tapply(costs$HC_BURDEN[costs$THERAP=="ANALOGUE"],
costs$HBA1C[costs$THERAP=="ANALOGUE"],mean))))
colnames(cost.ANALOGUE) <- c("HBA1C","HC_BURDEN")
cost.GLP <- as.data.frame(cbind(sort(unique(costs$HBA1C[costs$THERAP=="GLP"])),
as.numeric(tapply(costs$HC_BURDEN[costs$THERAP=="GLP"],
costs$HBA1C[costs$THERAP=="GLP"],mean))))
colnames(cost.GLP) <- c("HBA1C","HC_BURDEN")
## ANALOGUE therapy
# Mean
cost.ANALOGUE <- na.omit(as.data.frame(cbind(as.numeric(
costs$HBA1C[costs$THERAP=="ANALOGUE"]),
costs$HC_BURDEN[costs$THERAP=="ANALOGUE"])))
colnames(cost.ANALOGUE) <- c("HBA1C","HC_BURDEN")
cost.ANALOGUE <- cost.ANALOGUE[order(cost.ANALOGUE$HBA1C),]
cost.ANALOGUE <- cost.ANALOGUE[cost.ANALOGUE$HBA1C>lowest,]
cost.ANALOGUE$HBA1C <- cost.ANALOGUE$HBA1C-min(lowest)
# Fit non-linear model
mod.ANALOGUE <- nls(HC_BURDEN ~ a + b/(HBA1C + c),
start = list(a = 5, b = -5, c = 1), cost.ANALOGUE,
control = list(maxiter = 50000, minFactor = 0.000000000000001))
cost_ANALOGUE_plotdat <- cbind("ANALOGUE",as.data.frame(cbind(seq(0,6,6/99),
predict(mod.ANALOGUE,
newdata=data.frame(HBA1C = seq(0,6,6/99))))))
# Variance
cost_var.ANALOGUE <- na.omit(as.data.frame(cbind(sort(unique(
costs$HBA1C[costs$THERAP=="ANALOGUE"])),
as.numeric(tapply(costs$HC_BURDEN[costs$THERAP=="ANALOGUE"],
costs$HBA1C[costs$THERAP=="ANALOGUE"],var)))))
colnames(cost_var.ANALOGUE) <- c("HBA1C","HC_BURDEN")
cost_var.ANALOGUE <- cost_var.ANALOGUE[cost_var.ANALOGUE$HBA1C>lowest,]
cost_var.ANALOGUE$HBA1C <- cost_var.ANALOGUE$HBA1C-min(lowest)
# Fit non-linear model
mod_var.ANALOGUE <- nls(HC_BURDEN ~ a + b/(HBA1C + c),
start = list(a = 5, b = -3, c = 0.1),
cost_var.ANALOGUE[cost_var.ANALOGUE$HBA1C<10-lowest,],
control = list(maxiter = 50000, minFactor = 0.000000000000001))
cost_ANALOGUE_plotdat <- cbind(cost_ANALOGUE_plotdat,
cost_ANALOGUE_plotdat[,3] -
sqrt(predict(mod_var.ANALOGUE,
newdata=data.frame(HBA1C = seq(0,6,6/99)))),
cost_ANALOGUE_plotdat[,3] +
sqrt(predict(mod_var.ANALOGUE,
newdata=data.frame(HBA1C = seq(0,6,6/99)))))
colnames(cost_ANALOGUE_plotdat) <- c("Data","HbA1c","Therapy cost","low","high")
## GLP
# Mean
cost.GLP <- na.omit(as.data.frame(cbind(as.numeric(
costs$HBA1C[costs$THERAP=="GLP"]),
costs$HC_BURDEN[costs$THERAP=="GLP"])))
colnames(cost.GLP) <- c("HBA1C","HC_BURDEN")
cost.GLP <- cost.GLP[order(cost.GLP$HBA1C),]
cost.GLP <- cost.GLP[cost.GLP$HBA1C>lowest,]
cost.GLP$HBA1C <- cost.GLP$HBA1C-min(lowest)
# Simple linear
mod.GLP <- nls(HC_BURDEN ~ a + b * HBA1C,
start = list(a = 1, b = 1), cost.GLP,
control = list(maxiter = 50000, minFactor = 0.000000000000001))
cost_GLP_plotdat <- cbind("GLP",as.data.frame(cbind(seq(0,6,6/99),
predict(mod.GLP, newdata=data.frame(HBA1C = seq(0,6,6/99))))))
# Variance
cost_var.GLP <- na.omit(as.data.frame(cbind(sort(unique(
costs$HBA1C[costs$THERAP=="GLP"])),
as.numeric(tapply(costs$HC_BURDEN[costs$THERAP=="GLP"],
costs$HBA1C[costs$THERAP=="GLP"],var)))))
colnames(cost_var.GLP) <- c("HBA1C","HC_BURDEN")
cost_var.GLP <- cost_var.GLP[cost_var.GLP$HBA1C>lowest,]
cost_var.GLP$HBA1C <- cost_var.GLP$HBA1C-min(lowest)
# Simple linear
mod_var.GLP <- nls(HC_BURDEN ~ a + b*(HBA1C),
start = list(a = 5, b = -3), cost_var.GLP,
control = list(maxiter = 50000, minFactor = 0.000000000000001))
cost_GLP_plotdat <- cbind(cost_GLP_plotdat,
cost_GLP_plotdat[,3] -
sqrt(predict(mod_var.GLP, newdata=data.frame(HBA1C = seq(0,6,6/99)))),
cost_GLP_plotdat[,3] +
sqrt(predict(mod_var.GLP, newdata=data.frame(HBA1C = seq(0,6,6/99)))))
colnames(cost_GLP_plotdat) <- c("Data","HbA1c","Therapy cost","low","high")
### Out-of-control cost
COST_CUMU<-NULL
for (i in unique(RANDOMISED_DATA$ID[!is.na(RANDOMISED_DATA$HEALTHCARE_BURDEN_EUR)]))
{
vec <- RANDOMISED_DATA$HEALTHCARE_BURDEN_EUR[RANDOMISED_DATA$ID==i]
vec2 <- rollapply(vec,min(6,length(vec)),
sum,align="left",partial=TRUE)/
(pmin(length(vec)-(1:length(vec))+1,6)*30)
COST_CUMU <- rbind(COST_CUMU,cbind(i,vec2))
}
RANDOMISED_DATA$COST_CUMU <- NA
RANDOMISED_DATA$COST_CUMU[RANDOMISED_DATA$ID%in%COST_CUMU[,1]] <- COST_CUMU[,2]
discparam <- 150
cutHBA1C_AVG <- cut(na.omit(RANDOMISED_DATA$HBA1C_fill_filter),breaks=discparam)
newlvls <- seq(min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)),
max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)),
(max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter))-
min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)))/discparam)[1:discparam] +
(max(na.omit(RANDOMISED_DATA$HBA1C_fill_filter))-
min(na.omit(RANDOMISED_DATA$HBA1C_fill_filter)))/discparam/2
levels(cutHBA1C_AVG) <- newlvls
ooc_costs <- cbind(round(as.numeric(as.character(cutHBA1C_AVG)),1),
RANDOMISED_DATA$COST_CUMU[!is.na(RANDOMISED_DATA$HBA1C_fill_filter)])
ooc_costs <- as.data.frame(ooc_costs)
# Mean
disc_ooc_cost <- as.data.frame(cbind(as.numeric(ooc_costs[,1]),ooc_costs[,2]))
colnames(disc_ooc_cost) <- c("HBA1C","COST")
disc_ooc_cost <- disc_ooc_cost[order(disc_ooc_cost$HBA1C),]
disc_ooc_cost <- disc_ooc_cost[disc_ooc_cost$HBA1C >= lowest,]
disc_ooc_cost$HBA1C <- disc_ooc_cost$HBA1C - lowest
mod.COST <- nls(COST ~ a + c*HBA1C^2, start = list(a = 1, c = 1), disc_ooc_cost)
cost_COST_plotdat <- cbind("Complications",as.data.frame(cbind(seq(0, 6, 6/99),
predict(mod.COST, newdata=data.frame(HBA1C = seq(0, 6, 6/99))))))
# Variance
disc_ooc_cost_var <- as.data.frame(cbind(sort(unique(ooc_costs[,1])),
as.numeric(tapply(ooc_costs[,2],ooc_costs[,1],var))))
colnames(disc_ooc_cost_var) <- c("HBA1C","COST")
disc_ooc_cost_var <- disc_ooc_cost_var[disc_ooc_cost_var$HBA1C>=lowest,]
disc_ooc_cost_var$HBA1C <- disc_ooc_cost_var$HBA1C-lowest
mod_var.COST <- nls(COST ~ a + c*HBA1C^2,
start = list(a = 0.5, c = 0.5), disc_ooc_cost_var,
control = list(maxiter = 50000, minFactor = 0.000000000000001))
cost_COST_plotdat <- cbind(cost_COST_plotdat,
cost_COST_plotdat[,3] -
sqrt(predict(mod_var.COST,
newdata=data.frame(HBA1C = seq(0,6,6/99)))),
cost_COST_plotdat[,3] +
sqrt(predict(mod_var.COST,
newdata=data.frame(HBA1C = seq(0,6,6/99)))))
colnames(cost_COST_plotdat) <- c("Data","HbA1c","Therapy cost","low","high")
cost_plots <- rbind(cost_ANALOGUE_plotdat,cost_GLP_plotdat,cost_COST_plotdat)
cost_plots$HbA1c <- cost_plots$HbA1c + lowest
cost_plots[,"Therapy cost"] <- cost_plots[,"Therapy cost"]/1
cost_plots[,"low"] <- cost_plots[,"low"]/1
cost_plots[,"low"][cost_plots[,"low"]<0] <- 0
cost_plots[,"high"] <- cost_plots[,"high"]/1
cost_plots <- cost_plots
### Sampling cost: official, government-regulated prices related to sampling
### converted from HUF to EUR
sampling_cost=2875/369.05
### Empirical costs for comparison
# GLP
mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30 +
mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU) +
sampling_cost/mean(deltats[,2])
sd(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR/30 +
RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU +
sampling_cost/mean(deltats[,2]))
# ANALOGUE
mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30 +
mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU) +
sampling_cost/mean(deltats[,2])
sd(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR/30 +
RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU +
sampling_cost/mean(deltats[,2]))
### Empirical HbA1c distribution
# GLP
empi.GLP <- RANDOMISED_DATA$HBA1C_fill[grepl("GLP", RANDOMISED_DATA$THERAPY) &
RANDOMISED_DATA$HBA1C_fill>=4 & RANDOMISED_DATA$HBA1C_fill<=20]
cutHBA1C <- cut(na.omit(empi.GLP),breaks=100)
newlvls <- seq(min(na.omit(empi.GLP)),max(na.omit(empi.GLP)),
(max(na.omit(empi.GLP))-min(na.omit(empi.GLP)))/100)[1:100] +
(max(na.omit(empi.GLP))-min(na.omit(empi.GLP)))/100/2
levels(cutHBA1C) <- newlvls
empi.GLP <- as.data.frame(table(cutHBA1C)/sum(table(cutHBA1C)))
empi.GLP$cutHBA1C <- as.numeric(as.character(empi.GLP$cutHBA1C))
empi.GLP <- cbind("GLP", empi.GLP)
colnames(empi.GLP) <- c("Therapy", "HbA1c", "Probability")
# ANALOGUE
empi.ANALOGUE <- RANDOMISED_DATA$HBA1C_fill[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY) &
RANDOMISED_DATA$HBA1C_fill>=4 & RANDOMISED_DATA$HBA1C_fill<=20]
cutHBA1C <- cut(na.omit(empi.ANALOGUE),breaks=100)
newlvls <- seq(min(na.omit(empi.ANALOGUE)),
max(na.omit(empi.ANALOGUE)),
(max(na.omit(empi.ANALOGUE))-
min(na.omit(empi.ANALOGUE)))/100)[1:100] +
(max(na.omit(empi.ANALOGUE))-
min(na.omit(empi.ANALOGUE)))/100/2
levels(cutHBA1C) <- newlvls
empi.ANALOGUE <- as.data.frame(table(cutHBA1C)/sum(table(cutHBA1C)))
empi.ANALOGUE$cutHBA1C <- as.numeric(as.character(empi.ANALOGUE$cutHBA1C))
empi.ANALOGUE <- cbind("ANALOGUE", empi.ANALOGUE)
colnames(empi.ANALOGUE) <- c("Therapy", "HbA1c", "Probability")
# Merging datasets
hba1c_empir <- rbind(empi.ANALOGUE, empi.GLP)
# The list of gathered data and statistics:
# sigma_param, s_param, delta_param, ANALOGUE
# GLP, mod.ANALOGUE, mod_var.ANALOGUE
# mod.GLP, mod_var.GLP, mod.COST, mod_var.COST
# cost_plots, sampling_cost, hba1c_empir
# }
Run the code above in your browser using DataLab