Nothing
## ----setup, results = "hide", echo = FALSE, message = FALSE, warnings = FALSE----
set.seed(290875)
load(system.file("rda", "bloodloss.rda", package = "TH.data"))
### some packages
library("trtf")
library("tram")
library("rms")
library("coin")
library("survival")
library("ATR")
library("multcomp")
library("gridExtra")
library("vcd")
library("colorspace")
library("lattice")
tcols <- diverge_hcl(50, h = c(246, 40), c = 96, l = c(65, 90), alpha = .5)
cols <- qualitative_hcl(3, palette = "Harmonic")
vR <- paste(R.Version()$major, R.Version()$minor, sep = ".")
vtram <- packageDescription("tram")$Version
vtrtf <- packageDescription("trtf")$Version
### plot setup
trellis.par.set(list(plot.symbol = list(col="black", cex = .75),
box.rectangle = list(col=1),
box.umbrella = list(lty=1, col=1),
strip.background = list(col = "white")))
ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
lattice.options(default.theme = ltheme)
frmt1 <- function(x) formatC(round(x, 1), digits = 1, format = "f")
frmt3 <- function(x) {
if (!is.numeric(x)) return(x)
formatC(round(x, 3), digits = 3, format = "f")
}
### tree plots
ctrl <- ctree_control(alpha = 0.05, minbucket = 50,
teststat = "max", splitstat = "max", maxsurrogate = 3)
en <- function(obj, col = "black", bg = "white", fill = "transparent",
ylines = 2, id = TRUE, mainlab = NULL, gp = gpar(), K = 20,
type = c("trafo", "distribution", "survivor",
"density", "logdensity", "hazard",
"loghazard", "cumhazard", "quantile"),
flip = FALSE, axes = TRUE, xaxis = NULL, ...)
{
### panel function for ecdf in nodes
rval <- function(node) {
nid <- id_node(node)
dat <- data_party(obj, nid)
wn <- dat[["(weights)"]]
cf <- obj$coef[as.character(nid),c("Hb.prae", "F1.prae",
"F2.prae", "F13.Akt.prae")]
cf <- matrix(round(exp(cf), 3), nrow = 1)
# cf <- rbind(cf, obj$ci[as.character(nid),])
rownames(cf) <- c("OR")#, "CI")
colnames(cf) <- c("Hb.prae", "F1.prae", "F2.prae", "F13.Akt.prae")
colnames(cf) <- c("hemoglobin", "F. I", "F. II", "F. XIII")
top_vp <- viewport(layout = grid.layout(nrow = 1, ncol = 2,
widths = unit(c(1, ylines), c("null", "lines")),
heights = unit(1, "null")),
width = unit(1, "npc"),
height = unit(1, "npc"), # - unit(2, "lines"),
name = paste("node_mlt", nid, sep = ""), gp = gp)
pushViewport(top_vp)
grid.rect(gp = gpar(fill = bg, col = 0))
grid.draw(tableGrob(cf))
upViewport(1)
}
return(rval)
}
class(en) <- "grapcon_generator"
### format confidence intervals
ci <- function(m) {
cf <- coef(m)
idx <- 1:length(cf)
i <- grep("(Intercept)", names(cf))
if (length(i) > 0)
idx <- idx[-i]
cbind(exp(cf)[idx], exp(confint(m)[idx,]))
}
vlab <- function(x) {
lab <- code$desc_EN[code$varname == x]
lab <- paste0(toupper(substr(lab, 1, 1)), substr(lab, 2, nchar(lab)))
paste(lab, " (in ", code$unit[code$varname == x], ")", sep = "")
}
pvar <- function(x)
paste(paste(x[-length(x)], collapse = ", "), ", and ", x[length(x)], sep = "")
## ----vignette, eval = FALSE---------------------------------------------------
# library("knitr")
# knit("blood_loss_report.Rnw")
# library("tools")
# texi2pdf("blood_loss_report.tex")
## ----preproc, echo = FALSE----------------------------------------------------
x <- c("GA", "AGE", "MULTIPAR", "BMI", "TWIN", "FET.GEW", "IOL", "AIS")
z <- subset(code, type %in% c("confounder", "reason"))$varname
prae <- c("F1.prae", "F2.prae", "Hb.prae", ### "F13.Ag.prae",
"F13.Akt.prae")
### "DD.prae")
vars <- c("MBL", prae, z)
blood$mode <- with(blood, SECTIO.prim == "yes" | SECTIO.sek == "yes" |
SECTIO.not == "yes") + 1
blood$mode[blood$SECTIO.sek == "yes" | blood$SECTIO.not == "yes"] <- 3
blood$mode <- factor(blood$mode, levels = 1:3,
labels = c("Vaginal delivery", "Planned Cesarean", "Unplanned Cesarean"))
blood$VCmode <- blood$mode
levels(blood$VCmode) <- c("Vaginal delivery", "Cesarean Sectio", "Cesarean Sectio")
ct <- table(blood$mode)
## ----table-1, echo = FALSE, results = "asis"----------------------------------
frmtnum <- function(x)
paste(frmt1(x)[2], " (", paste(frmt1(x)[-2], collapse = "-"), ")", sep = "")
frmtrow <- function(x) {
if (is.factor(blood[[x]])) return(table(blood[[x]], blood$mode))
prb <- c(0.25, .5, .75)
ret <- c(vaginal = frmtnum(quantile(subset(blood, mode == "Vaginal delivery")[[x]],
prob = prb, na.rm = TRUE)),
planned = frmtnum(quantile(subset(blood, mode == "Planned Cesarean")[[x]],
prob = prb, na.rm = TRUE)),
unplanned = frmtnum(quantile(subset(blood, mode == "Unplanned Cesarean")[[x]],
prob = prb, na.rm = TRUE)))
ret <- matrix(ret, nrow = 1)
rownames(ret) <- "Med (IQR)"
if (all(!is.na(blood[[x]]))) return(ret)
ret <- rbind(table(is.na(blood[[x]]), blood$mode)["TRUE",], ret)
rownames(ret) <- c("missing", "Med (IQR)")
ret
}
tab <- lapply(vars, frmtrow)
names(tab) <- vars
desc <- code$desc_EN[match(vars, code$varname)]
desc <- rep(desc, sapply(tab, nrow))
desc[dd <- duplicated(desc)] <- ""
unit <- code$unit[match(vars, code$varname)]
unit <- gsub("\\%", "\\\\%", unit)
unit <- gsub("kg/m\\^2", "$\\\\text{kg} / \\\\text{m}^2$", unit)
unit <- rep(unit, sapply(tab, nrow))
unit[dd] <- ""
unit[is.na(unit)] <- ""
xtab <- do.call("rbind", tab)
xtab <- cbind(var = desc, unit = unit, item = rownames(xtab), xtab)
writeLines(paste(xtab[,1], " & ", xtab[,2], "&", xtab[, 3], " & ", xtab[,4],
" & ", xtab[,5], " & ", xtab[,6], " \\\\"))
write.csv(xtab, file = "table_1.csv")
## ----MBL-plot, echo = FALSE, fig.width = 6, fig.height = 5, cache = TRUE, warning = FALSE, dev = c("tiff", "pdf", "png"), dpi = 300----
### unconditional MBL
qy <- 0:max(blood$MBL)
### MBL outcome as interval censored
### interval length: 50 for MBL < 1000; 100 for MBL > 1000
off <- 25
tm1 <- with(blood, ifelse(MBL < 1000, MBL - off, MBL - 2 * off))
tm2 <- with(blood, ifelse(MBL >= 1000, MBL + 2 * off, MBL + off))
### some measurements are more precise, use length 10 here
ex <- !blood$MLB %in% seq(from = 100, to = 6000, by = 50)
stopifnot(sum(ex) == 0)
tm1[ex] <- blood$MBL[ex] - off / 5
tm2[ex] <- blood$MBL[ex] + off / 5
blood$MBLsurv <- Surv(time = tm1, time2 = tm2, type = "interval2")
MBLlim <- c(0, 2700)
nd <- data.frame(mode = sort(unique(blood$mode)))
### takes too long on Windows
if (FALSE) {
plot(m <- as.mlt(Colr(MBLsurv | mode ~ 1, data = blood, order = 15,
bounds = c(0, Inf), support = c(250, 2000),
extrapolate = TRUE)), newdata = nd,
q = qy, type = "distribution", col = cols, lwd = 2, xlim = MBLlim,
xlab = vlab("MBL"), ylab = "Probability", ylim = c(-.05, 1.05),
inset = 10)
rug(blood$MBL[blood$mode == "Vaginal delivery"], lwd = 2, col = cols[1])
rug(blood$MBL[blood$mode == "Planned Cesarean"], side = 3, lwd = 2, col = cols[2])
rug(blood$MBL[blood$mode == "Unplanned Cesarean"], side = 3, lwd = 2, col = cols[3])
legend("bottomright", lwd = 2, col = cols, legend = levels(blood$mode), bty = "n")
}
## ----MBL-plot-check, eval = FALSE, echo = FALSE, fig.width = 6, fig.height = 5, warning = FALSE, dev = c("tiff", "pdf", "png"), dpi = 300----
# plot(survfit(MBLsurv ~ mode, data = blood), col = cols, xlim = c(0, 2700),
# lty = 2, xlab = vlab("MBL"), ylab = "1 - Probability", ylim = c(-.05, 1.05))
# plot(m, newdata= nd, type = "survivor", add = TRUE, col = cols, lty = 1)
# rug(blood$MBL[blood$mode == "Vaginal delivery"], lwd = 2, col = cols[1])
# rug(blood$MBL[blood$mode == "Planned Cesarean"], side = 3, lwd = 2, col = cols[2])
# rug(blood$MBL[blood$mode == "Unplanned Cesarean"], side = 3, lwd = 2, col = cols[3])
# legend("topright", lwd = 2, col = cols, legend = levels(blood$mode), bty = "n")
## ----MBL-Colr, echo = TRUE, cache = TRUE--------------------------------------
mvars <- c("Hb.prae", "F1.prae", "F2.prae", "F13.Akt.prae")
fm <- paste(mvars, collapse = "+")
### continuous outcome logistic regression
m_MBL <- Colr(as.formula(paste("MBLsurv ~ ", fm)), data = blood,
bounds = c(0, Inf), support = c(250, 2000))
### number of observations
sum(complete.cases(model.frame(m_MBL)))
summary(m_MBL)
logLik(m_MBL)
## ----MBL-Colr-ci, echo = TRUE, cache = TRUE-----------------------------------
(ci_all <- ci(m_MBL))
## ----MBL-orm, echo = TRUE-----------------------------------------------------
(m_MBL_orm <- orm(as.formula(paste("MBL ~ ", fm)), data = blood))
### OR and confidence interval for F. XIII
### (sign of the coefficient is different in rms::orm and tram::Colr)
exp(-coef(m_MBL_orm)["F13.Akt.prae"])
exp(-rev(confint(m_MBL_orm)["F13.Akt.prae",]))
## ----MBL-Colr-Cesar, echo = TRUE, cache = TRUE--------------------------------
m_MBL_C <- Colr(as.formula(paste("MBLsurv | VCmode ~ VCmode:(", fm, ")")),
data = blood, bounds = c(0, Inf), support = c(250, 2000))
summary(m_MBL_C)
logLik(m_MBL_C)
## ----MBL-Colr-pre, echo = FALSE, fig.width = 6, fig.height = 5, dev = c("tiff", "pdf", "png"), dpi = 300----
nd <- blood[1,mvars, drop = FALSE]
nd[1,] <- cm <- round(apply(blood[, mvars], 2, median, na.rm = TRUE), 1)
F13 <- seq(from = min(blood$F13.Akt.prae, na.rm = TRUE),
to = max(blood$F13.Akt.prae, na.rm = TRUE),
length = 25)
nd <- nd[rep(1, length(F13)),]
nd[, "F13.Akt.prae"] <- F13
nd$MBLsurv <- 500
X <- model.matrix(as.mlt(m_MBL)$model, data = nd)
cf <- confint(glht(as.mlt(m_MBL), linfct = X))
plot(F13, plogis(cf$confint[, "Estimate"], lower.tail = FALSE), type = "l",
ylim = c(0, 1), xlab = "Prepartum F. XIII (%)",
ylab = expression(paste("Probability PPH (", MBL >= 500, "ml)")))
lwr <- plogis(cf$confint[, "lwr"], lower.tail = FALSE)
upr <- plogis(cf$confint[, "upr"], lower.tail = FALSE)
polygon(c(F13, rev(F13)), c(lwr, rev(upr)),
border = NA, col = "lightblue")
lines(F13, plogis(cf$confint[, "Estimate"], lower.tail = FALSE))
rug(blood$F13.Akt.prae, col = rgb(.1, .1, .1, .1))
## ----sample-size, echo = TRUE-------------------------------------------------
blood$PPH <- factor(blood$MBL >= 500, levels = c(FALSE, TRUE),
labels = c("no", "yes"))
summary(m_PPH <- lm(F13.Akt.prae ~ PPH, data = blood))
confint(m_PPH)["PPHyes",]
## ----sample-size-W, echo = TRUE, cache = TRUE---------------------------------
wilcox_test(F13.Akt.prae ~ PPH, data = blood,
distribution = approximate(10000), conf.int = TRUE)
## ----MBL-check-1, echo = TRUE, cache = TRUE, warning = FALSE, message = FALSE, results = "hide"----
### Tobit model
tll <- logLik(t_MBL <- Lm(as.formula(paste("MBLsurv ~ ", fm)),
data = blood))
### distribution regression
drll <- logLik(dr_MBL <- Colr(as.formula(paste("MBLsurv | ", fm, "~ 1")),
data = blood, bounds = c(0, Inf),
support = c(250, 2000)))
## ----dr-plot, echo = FALSE, warning = FALSE, dev = c("tiff", "pdf", "png"), dpi = 300----
nd0 <- data.frame("Hb.prae" = 0, "F1.prae" = 0, "F2.prae" = 0, "F13.Akt.prae" = 0)
q <- seq(from = MBLlim[1], to = MBLlim[2], length.out = 200)
p0 <- predict(as.mlt(dr_MBL), newdata = nd0, q = q, type = "trafo")
layout(matrix(1:4, ncol = 2))
for (i in 1:4) {
nd <- nd0
nd[[i]] <- 1
cim <- confint(m_MBL)[i,]
p <- c(predict(as.mlt(dr_MBL), newdata = nd, q = q, type = "trafo") - p0)
ylim <- exp(max(abs(c(p[is.finite(p)], cim))) * c(-1, 1))
cn <- code$desc_EN[code$varname == colnames(nd)[i]]
plot(q, exp(p), type = "l", main = cn,
ylim = ylim, xlab = vlab("MBL"), ylab = expression(exp(beta)))
polygon(c(-100, 3100, 3100, -100), rep(exp(cim), each = 2),
border = NA, col = rgb(.1, .1, .1, .1))
abline(h = exp(coef(m_MBL)[i]), lty = 2)
abline(h = 1, lty = 1, col = "darkred", lwd = 3)
rug(blood$MBL, col = rgb(.1, .1, .1, .1))
}
## ----MBL-check-2, echo = FALSE, results = "hide", cache = TRUE, warning = FALSE, message = FALSE----
### Binary logistic regression
lr500 <- glm(as.formula(paste("I(MBL < 500) ~ ", fm)),
data = blood, family = binomial())
(ci_500 <- ci(lr500))
lr750 <- glm(as.formula(paste("I(MBL < 750) ~ ", fm)),
data = blood, family = binomial())
(ci_750 <- ci(lr750))
lr1000 <- glm(as.formula(paste("I(MBL < 1000) ~ ", fm)),
data = blood, family = binomial())
(ci_1000 <- ci(lr1000))
## ----table-2, echo = FALSE, results = "asis"----------------------------------
p_all <- paste("$", format.pval(summary(m_MBL)$test$test$pvalues, eps = 0.001), "$")
p_500 <- paste("$", format.pval(summary(lr500)$coef[-1,4], eps = 0.001), "$")
p_750 <- paste("$", format.pval(summary(lr750)$coef[-1,4], eps = 0.001), "$")
p_1000 <- paste("$", format.pval(summary(lr1000)$coef[-1,4], eps = 0.001), "$")
ci_tab <- rbind(cbind(frmt3(ci_all), p_all),
cbind(frmt3(ci_500), p_500),
cbind(frmt3(ci_750), p_750),
cbind(frmt3(ci_1000), p_1000))
ci_tab <- cbind(c("All", rep("", 3),
"500", rep("", 3),
"750", rep("", 3),
"1000", rep("", 3)), rownames(ci_tab), ci_tab)
writeLines(paste(ci_tab[,1], " & ", ci_tab[,2], " & ",
ci_tab[,3], " & ", ci_tab[,4], " & ", ci_tab[,5], " & ", ci_tab[,6], " \\\\"))
write.csv(ci_tab, file = "table_2.csv")
## ----table-3, echo = FALSE, results = "asis"----------------------------------
ci_tab <- cbind(frmt3(ci(m_MBL_C)), format.pval(summary(m_MBL_C)$test$test$pvalues, eps = 0.001))
ci_tab <- cbind(names(coef(m_MBL_C)), ci_tab)
writeLines(paste(ci_tab[,1], " & ", ci_tab[,2], " & ",
ci_tab[,3], " & ", ci_tab[,4], " & ", ci_tab[,5], " \\\\"))
write.csv(ci_tab, file = "table_3.csv")
## ----MBL-xtree, echo = FALSE, cache = TRUE, fig.width = 12.25, fig.height = 4.5, dev = c("tiff", "pdf", "png"), dpi = 300----
### partitioning
xfm <- paste(x, collapse = "+")
zfm <- paste(z, collapse = "+")
xfm_MBL <- as.formula(paste("MBLsurv ~ ", fm, "|", xfm))
zfm_MBL <- as.formula(paste("MBLsurv ~ ", fm, "|", zfm))
blood$DAUER.ap[is.na(blood$DAUER.ap)] <- 0
yfm <- as.formula(paste("MBLsurv ~ ", fm))
xMBL_cc <- complete.cases(blood[, all.vars(yfm)])
zMBL_cc <- complete.cases(blood[, all.vars(yfm)])
xitr_MBL <- trafotree(as.mlt(m_MBL), formula = xfm_MBL, data = blood[xMBL_cc,],
parm = mvars, control = ctrl)
nodeid <- predict(xitr_MBL, newdata = blood, type = "node")
blood$nd <- factor(nodeid, levels = sort(unique(nodeid)),
labels = sort(unique(nodeid)))
m <- Colr(as.formula(paste("MBLsurv | nd ~ nd:(", fm, ")")),
data = blood[xMBL_cc,], bounds = c(0, Inf), support = c(250, 2000))
cf <- ci(m)
cf <- formatC(round(cf, 3), digits = 3, format = "f")
xitr_MBL$ci <- matrix(paste(matrix(cf[,2], ncol = 4),
matrix(cf[,3], ncol = 4), sep = "-"), ncol = 4)
rownames(xitr_MBL$ci) <- levels(blood$nd)
plot(rotate(xitr_MBL), terminal_panel = en)
## ----MBL-extreme, echo = FALSE------------------------------------------------
xn <- tapply(1:nrow(blood), blood$nd, function(i) colMeans(blood[i, mvars], na.rm = TRUE))
nd <- as.data.frame(do.call("rbind", xn))
nd$nd <- sort(unique(blood$nd))
cb <- confband(as.mlt(m), newdata = nd, type = "distribution",
K = 500, calpha = univariate_calpha())
### 500
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 500)^2),-1]))))[, c(1, 2, 4, 3)]
## ----MBL-extreme-, echo = FALSE-----------------------------------------------
### 750
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 750)^2),-1]))))[, c(1, 2, 4, 3)]
## ----MBL-extreme-1000, echo = FALSE-------------------------------------------
### 1000
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 1000)^2),-1]))))[, c(1, 2, 4, 3)]
## ----MBL-extreme-50, echo = FALSE---------------------------------------------
nd50 <- nd
nd50$F13.Akt.prae <- nd50$F13.Akt.prae + 50
cb <- confband(as.mlt(m), newdata = nd50, type = "distribution",
K = 500, calpha = univariate_calpha())
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 500)^2),-1]))))[, c(1, 2, 4, 3)]
## ----MBL-extreme-50-750, echo = FALSE-----------------------------------------
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 750)^2),-1]))))[, c(1, 2, 4, 3)]
## ----MBL-extreme-50-1000, echo = FALSE----------------------------------------
data.frame(Subgroup = nd$nd, do.call("rbind", lapply(cb, function(x) 100 * (1 - x[which.min((x[, "q"] - 1000)^2),-1]))))[, c(1, 2, 4, 3)]
## ----F13-trt, echo = FALSE, results = "hide", dev = c("tiff", "pdf", "png"), dpi = 300----
nd <- nd[rep(1:nrow(nd), 2),]
nd$type <- gl(2, nrow(nd) / 2, labels = c("prepartum F. XIII", "prepartum F. XIII + 50 units"))
nd$F13.Akt.prae <- nd$F13.Akt.prae + c(0, 50)[nd$type]
nd$med <- c(predict(as.mlt(m), newdata = nd, type = "quantile", prob = .5))
nd$MBL1000 <- (1 - c(predict(as.mlt(m), newdata = nd, type = "distribution", q = 1000))) * 100
print(nd)
qy <- 2:2000
p <- predict(as.mlt(m), newdata = nd, type = "distribution", q = qy)
nd <- nd[rep(1:nrow(nd), each = length(qy)),]
nd$MBL <- rep(qy, nrow(nd) / length(qy))
nd$p <- c(p)
key <- simpleKey(levels(nd$type), points = FALSE, lines = TRUE,
space = "top", lty = 1)
key$lines$col <- cols2 <- tcols[c(1, length(tcols))]
plt <- vector(mode = "list", length = 2)
# plt[[1]] <- plot(rotate(xitr_MBL), terminal_panel = en, pop = FALSE)
pfun <- function(x, y, ...) {
panel.abline(v = c(500, 750, 1000), col = "lightgrey")
panel.xyplot(x = x, y = y, ...)
}
levels(nd$nd) <- paste("Subgroup", levels(nd$nd))
plt[[2]] <- xyplot(p ~ MBL | nd, data = nd, groups = type, type = "l",
panel = pfun,
key = key, col = cols2, xlab = vlab("MBL"), ylab = "Probability")
# layout = matrix(1:nlevels(blood$nd), ncol = 1))
plot(plt[[2]])
# grid.arrange(plt[[1]], plt[[2]], ncol = 2)
## ----MBL-ztree, echo = FALSE, cache = TRUE, fig.width = 12.5, fig.height = 4.5, dev = c("tiff", "pdf", "png"), dpi = 300----
zitr_MBL <- trafotree(as.mlt(m_MBL), formula = zfm_MBL, data = blood[zMBL_cc,],
parm = mvars, control = ctrl)
blood$nd <- factor(predict(zitr_MBL, newdata = blood, type = "node"))
m <- Colr(as.formula(paste("MBLsurv | nd ~ nd:(", fm, ")")),
data = blood[zMBL_cc,], bounds = c(0, Inf), support = c(250, 2000))
cf <- ci(m)
cf <- formatC(round(cf, 3), digits = 3, format = "f")
zitr_MBL$ci <- matrix(paste(matrix(cf[,2], ncol = 4),
matrix(cf[,3], ncol = 4), sep = "-"), ncol = 4)
rownames(zitr_MBL$ci) <- levels(blood$nd)
plot(rotate(zitr_MBL), terminal_panel = en)
## ----session, echo = FALSE----------------------------------------------------
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.