Skip to content

Commit

Permalink
Merge pull request #50 from certara/support_binless_pctalq
Browse files Browse the repository at this point in the history 
generate pctalq if censoring uloq for binless()
  • Loading branch information
@certara-jcraig
certara-jcraig authored Sep 7, 2023
2 parents 816bb59 + 1c1bb38 commit 23bd54a
Showing 1 changed file with 91 additions and 38 deletions.
129 changes: 91 additions & 38 deletions R/binless.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)

}

.binlessvpcstats <- function(o, qpred=c(0.05, 0.5, 0.95), ..., conf.level=0.95, quantile.type=7, vpc.type){
.binlessvpcstats <- function(o, qpred=c(0.05, 0.5, 0.95), conf.level=0.95, quantile.type=7, vpc.type, ...){
y <- x <- blq <- fit <- . <- repl <- cprop <- rqssmed <- llam.med <- c.rqssmed <- NULL

obs.fits <- o$rqss.obs.fits
Expand All @@ -295,17 +295,34 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
qpred <- o$qpred
conf.level <- o$conf.level
qconf <- c(0, 0.5, 1) + c(1, 0, -1)*(1 - conf.level)/2

if(!is.null(obs$blq) && any(obs$blq)) {
if(!is.null(o$strat)) {
stratlloq <- c(names(o$strat), "lloq")
lloq <- obs[, stratlloq, with = FALSE]
lloq <- unique(lloq)
obs.fits <- obs.fits[lloq, on = names(o$strat)]
} else {
obs.fits[, lloq := rep(obs$lloq, len=.N)]

censor_vars <- c()
if (!is.null(obs$blq) && any(obs$blq)) {
censor_vars <- c(censor_vars, lloq = "blq")
}
if (!is.null(obs$alq) && any(obs$alq)) {
censor_vars <- c(censor_vars, uloq = "alq")
}


for (i in seq_along(censor_vars)) {
censor_var <- censor_vars[[i]]
censor_var_name <- names(censor_vars)[[i]]
if (!is.null(obs[[censor_var]]) && any(obs[[censor_var]])) {
if (!is.null(o$strat)) {
stratloq <- c(names(o$strat), censor_var_name)
loq <- obs[, stratloq, with = FALSE]
loq <- unique(loq)
obs.fits <- obs.fits[loq, on = names(o$strat)]
} else {
obs.fits[, (censor_var_name) := rep(obs[[censor_var_name]], len = .N)]
}
if(censor_var == "blq") {
obs.fits[, (censor_var) := ifelse(fit < lloq, TRUE, FALSE)]
} else {
obs.fits[, (censor_var) := ifelse(fit > uloq, TRUE, FALSE)]
}
}
obs.fits[, blq := ifelse(fit < lloq, TRUE, FALSE)]
}

obs.fits <- setnames(obs.fits[, lapply(.SD, median), by = x.binless], "fit", "y")
Expand All @@ -315,19 +332,33 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
obs.fits[, blq := ifelse(y < lloq, TRUE, FALSE)]
obs.fits[, y := ifelse(blq == TRUE, NA, y)]
}
if(!is.null(obs$alq) && any(obs$alq)) {
obs.fits[, alq := ifelse(y > uloq, TRUE, FALSE)]
obs.fits[, y := ifelse(alq == TRUE, NA, y)]
}

if (!is.null(o$strat)) {
stats <- obs.fits[sim.fits, on = c("x", "qname", names(o$strat))]
} else {
stats <- obs.fits[sim.fits, on = c("x", "qname")]
}

if (!is.null(obs$blq) && any(obs$blq)) {
if(is.null(o$strat)) {
sim[, lloq := rep(obs$lloq, len=.N)]
sim[, blq := (y < lloq)]
pctblq <- pctalq <- NULL

if (length(censor_vars) > 0) {
for (i in seq_along(censor_vars)) {
censor_var <- censor_vars[[i]]
censor_var_name <- names(censor_vars)[[i]]
if(is.null(o$strat)) {
sim[, (censor_var_name) := rep(obs[[censor_var_name]], len=.N)]

if (censor_var == "blq") {
sim[, (censor_var) := (y < lloq)]
} else {
sim[, (censor_var) := (y > uloq)]
}
setorder(obs, cols = x)
cprop.obs <- obs[, cprop := cumsum(blq) / 1:length(blq)]
cprop.obs <- obs[, cprop := cumsum(get(censor_var)) / 1:length(get(censor_var))]

sic.cprop <- function(llam) {
a <- AIC(
Expand All @@ -347,26 +378,41 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
)
cprop.obs$med <- fitted(med.obs.cprop)

setorder(sim, repl, x)[, cprop := cumsum(blq) / 1:length(blq), by=list(repl)]
setorder(sim, repl, x)[, cprop := cumsum(get(censor_var)) / 1:length(get(censor_var)), by=list(repl)]

suppressWarnings(sim[, rqssmed := fitted(rqss(cprop ~ qss(x, lambda = exp(llam.med.cprop)),
tau = 0.5, na.action = na.exclude, .SD)), by = .(repl)])

u.x <- unique(cprop.obs$x) #%#
u.x <- unique(cprop.obs$x)
med.obs.cprop <- tapply(cprop.obs$med, cprop.obs$x, median)
med.sims.blq <- tapply(sim$rqssmed, sim$x, median)
med.sims.blq.lb <- tapply(sim$rqssmed, sim$x, quantile, probs=c(qconf[[1]]))
med.sims.blq.ub <- tapply(sim$rqssmed, sim$x, quantile, probs=c(qconf[[3]]))
pctblq <- data.table(cbind(u.x,med.obs.cprop, med.sims.blq.lb, med.sims.blq, med.sims.blq.ub))

if (censor_var == "blq") {
pctblq <- data.table(cbind(u.x,med.obs.cprop, med.sims.blq.lb, med.sims.blq, med.sims.blq.ub))
setnames(pctblq, c("x", "y", "lo", "md", "hi"))
pctblq[, (setdiff(names(pctblq), "x")) := lapply(.SD, function(col) col * 100), .SDcols = -"x"]
}

if (censor_var == "alq") {
pctalq <- data.table(cbind(u.x,med.obs.cprop, med.sims.blq.lb, med.sims.blq, med.sims.blq.ub))
setnames(pctalq, c("x", "y", "lo", "md", "hi"))
pctalq[, (setdiff(names(pctalq), "x")) := lapply(.SD, function(col) col * 100), .SDcols = -"x"]
}

setnames(pctblq, c("x", "y", "lo", "md", "hi"))
} else {
strat <- o$strat
strat.split <- split(obs, strat)
strat.split <- strat.split[lapply(strat.split,NROW)>0]
x.strat <- c("x", names(strat))
sim[, lloq := rep(obs$lloq, len=.N), by = names(strat)]
sim[, blq := (y < lloq)]

sim[, (censor_var_name ) := rep(obs[[censor_var_name]], len=.N), by = names(strat)]
if (censor_var == "blq") {
sim[, blq := (y < lloq)]
} else {
sim[, alq := (y > uloq)]
}
stratx.binless <- obs[, list(x, o$strat)]
stratxrepl <- data.table(stratx.binless, sim[, .(repl)])
strat.split.sim <- split(sim, strat)
Expand All @@ -376,42 +422,49 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
rqss(
cprop ~
qss(x, lambda=exp(llam)),
tau=.5, na.action=na.exclude, data = strat.split[[i]]
tau=.5, na.action=na.exclude, data = strat.split[[j]]
),
k=-1
)
}
llam.strat.med.cprop <- vector("list", length(strat.split))
for (i in seq_along(strat.split)) {
setorder(strat.split[[i]], cols = x)
strat.split[[i]] <- strat.split[[i]][, cprop := cumsum(blq) / 1:length(blq)]
llam.strat.med.cprop[[i]] <- strat.split[[i]][, .(llam.med = optimize(sic.strat.cprop, interval=c(0, 7))$min)][,.(med = unlist(llam.med))]
strat.split[[i]][, c.rqssmed := fitted(rqss(cprop ~ qss(x, lambda = exp(llam.strat.med.cprop[[i]][[1]])),tau= .5, na.action = na.exclude, data = strat.split[[i]]))]
for (j in seq_along(strat.split)) {
setorder(strat.split[[j]], cols = x)
strat.split[[j]] <- strat.split[[j]][, cprop := cumsum(get(censor_var)) / 1:length(get(censor_var))]
llam.strat.med.cprop[[j]] <- strat.split[[j]][, .(llam.med = optimize(sic.strat.cprop, interval=c(0, 7))$min)][,.(med = unlist(llam.med))]
strat.split[[j]][, c.rqssmed := fitted(rqss(cprop ~ qss(x, lambda = exp(llam.strat.med.cprop[[j]][[1]])),tau= .5, na.action = na.exclude, data = strat.split[[j]]))]
}

obs.cprop <- rbindlist(strat.split)
obs.cprop <- setnames(obs.cprop[, lapply(.SD, median, na.rm = TRUE), by = x.strat, .SDcols = "c.rqssmed"], "c.rqssmed", "y")

for (i in seq_along(strat.split.sim)) {
setorder(strat.split.sim[[i]], cols = repl, x)
strat.split.sim[[i]] <- strat.split.sim[[i]][, cprop := cumsum(blq) / 1:length(blq), by = .(repl)]
strat.split.sim[[i]][, c.rqssmed := fitted(rqss(cprop ~ qss(x, lambda = exp(llam.strat.med.cprop[[i]][[1]])),tau= .5, na.action = na.exclude, .SD)), by = .(repl)]
for (j in seq_along(strat.split.sim)) {
setorder(strat.split.sim[[j]], cols = repl, x)
strat.split.sim[[j]] <- strat.split.sim[[j]][, cprop := cumsum(get(censor_var)) / 1:length(get(censor_var)), by = .(repl)]
strat.split.sim[[j]][, c.rqssmed := fitted(rqss(cprop ~ qss(x, lambda = exp(llam.strat.med.cprop[[j]][[1]])),tau= .5, na.action = na.exclude, .SD)), by = .(repl)]
}

sim.cprop <- rbindlist(strat.split.sim)
sim.med <- setnames(sim.cprop[, lapply(.SD, median, na.rm = TRUE), by = x.strat, .SDcols = "c.rqssmed"], "c.rqssmed", "md")
sim.lb <- setnames(sim.cprop[, lapply(.SD, quantile, probs = qconf[[1]]), by = x.strat, .SDcols = "c.rqssmed"], "c.rqssmed", "lo")
sim.ub <- setnames(sim.cprop[, lapply(.SD, quantile, probs = qconf[[3]]), by = x.strat, .SDcols = "c.rqssmed"], "c.rqssmed", "hi")

pctblq <- obs.cprop[sim.med, on = x.strat]
pctblq <- pctblq[sim.lb, on = x.strat]
pctblq <- pctblq[sim.ub, on = x.strat]
pctlq <- obs.cprop[sim.med, on = x.strat]
pctlq <- pctlq[sim.lb, on = x.strat]
pctlq <- pctlq[sim.ub, on = x.strat]
pctlq[, (setdiff(names(pctlq), x.strat)) := lapply(.SD, function(col) col * 100), .SDcols = -x.strat]

if (censor_var == "blq") {
pctblq <- pctlq
}
if (censor_var == "alq") {
pctalq <- pctlq
}
}
} else {
pctblq <- NULL
}
}

update(o, stats = stats, pctblq = pctblq, vpc.type = vpc.type)
update(o, stats = stats, pctblq = pctblq, pctalq = pctalq, vpc.type = vpc.type)
}

binlessaugment <- function(o, qpred = c(0.05, 0.50, 0.95), interval = c(0,7), loess.ypc = FALSE, ...) {
Expand Down

0 comments on commit 23bd54a

Please sign in to comment.