Nothing
demo/mtram.R
In tram: Transformation Models
## ----mtram-pkgs, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
set.seed(290875)
pkgs <- c("colorspace", "survival", "lme4", "tram", "gridExtra",
"lattice", "latticeExtra", "mvtnorm", "ordinalCont", "tramME")
pkgs <- pkgs %in% installed.packages()
## ----mtram-citation, echo = FALSE---------------------------------------------
year <- substr(packageDescription("tram")$Date, 1, 4)
version <- packageDescription("tram")$Version
## ----fail, results = "asis", echo = FALSE-------------------------------------
if (any(!pkgs))
{
cat(paste("Package(s)", paste(names(pkgs)[!pkgs], collapse = ", "),
"not available, stop processing.",
"\\end{document}\n"))
knitr::knit_exit()
}
## ----mtram-setup, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
library("lattice")
library("latticeExtra")
library("gridExtra")
trellis.par.set(list(plot.symbol = list(col=1,pch=20, cex=0.7),
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)
knitr::opts_chunk$set(echo = TRUE, results = 'markup', error = FALSE,
warning = FALSE, message = FALSE,
tidy = FALSE, cache = FALSE, size = "small",
fig.width = 6, fig.height = 4, fig.align = "center",
out.width = NULL, ###'.6\\linewidth',
out.height = NULL,
fig.scap = NA)
knitr::render_sweave() # use Sweave environments
knitr::set_header(highlight = '') # do not \usepackage{Sweave}
## R settings
options(prompt = "R> ", continue = "+ ", useFancyQuotes = FALSE) # JSS style
options(width = 75)
### ecdf plots
myprepanel <- function (x, y, f.value = NULL, ...)
{
ans <- prepanel.default.qqmath(x, f.value = f.value, distribution = qunif)
with(ans, list(xlim = ylim, ylim = c(0, 1), dx = dy, dy = dx))
}
mypanel <- function (x, y, f.value = NULL, type = "s", groups = NULL, qtype = 7,
ref = TRUE, ...)
{
if (ref) {
reference.line <- trellis.par.get("reference.line")
do.call(panel.abline, c(list(h = c(0, 1)), reference.line))
}
x <- as.numeric(x)
distribution <- qunif
nobs <- sum(!is.na(x))
if (!is.null(groups)) {
panel.superpose(x, y = NULL, f.value = f.value, type = type,
distribution = distribution, qtype = qtype, groups = groups,
panel.groups = panel.ecdfplot, ...)
}
else if (nobs) {
if (is.null(f.value)) {
panel.xyplot(x = sort(x), y = cumsum(y[order(x)]) / sum(y),
type = type, ...)
}
else {
p <- if (is.numeric(f.value))
f.value
else f.value(nobs)
panel.xyplot(x = quantile(x, p, names = FALSE, type = qtype,
na.rm = TRUE), y = distribution(p), type = type,
...)
}
}
}
library("colorspace")
col <- diverge_hcl(2, h = c(246, 40), c = 120, l = c(65, 90), alpha = .75)
## ----mtram-vignette, eval = FALSE---------------------------------------------
## install.packages("tram")
## demo("mtram", package = "tram")
## ----mtram-sleep-plot, echo = FALSE-------------------------------------------
library("lme4")
xyplot(Reaction ~ Days | Subject, data = sleepstudy,
xlab = "Days of sleep deprivation", ylab = "Average reaction time (in ms)")
## ----mtram-sleep_lmer, cache = FALSE------------------------------------------
sleep_lmer <- lmer(Reaction ~ Days + (Days | Subject),
data = sleepstudy, REML = FALSE)
## ----mtram-tram---------------------------------------------------------------
library("tram")
## ----mtram-sleep_mtram, cache = FALSE-----------------------------------------
sleep_LM <- Lm(Reaction ~ Days, data = sleepstudy)
sleep_LMmer <- mtram(sleep_LM, ~ (Days | Subject), data = sleepstudy)
## ----mtram-sleep_cmp----------------------------------------------------------
logLik(sleep_lmer)
logLik(sleep_LMmer)
## ----mtram-sleep_sd-----------------------------------------------------------
(sdinv <- 1 / summary(sleep_lmer)$sigma)
coef(sleep_LMmer)["Reaction"]
## ----mtram-sleep_beta---------------------------------------------------------
fixef(sleep_lmer) * c(-1, 1) * sdinv
coef(sleep_LMmer)[c("(Intercept)", "Days")]
## ----mtram-sleep_varparm------------------------------------------------------
sleep_lmer@theta
coef(sleep_LMmer)[-(1:3)]
## ----mtram-sleep_varcorr------------------------------------------------------
sleep_LMmer$G * (1 / sdinv)^2
cov2cor(sleep_LMmer$G * (1 / sdinv)^2)
unclass(VarCorr(sleep_lmer))$Subject
## ----mtram-sleep-Surv---------------------------------------------------------
library("survival")
sleepstudy$Reaction_I <- with(sleepstudy, Surv(Reaction - 20, Reaction + 20,
type = "interval2"))
sleepstudy$Reaction_I[1:5]
## ----mtram-sleep-interval, cache = FALSE--------------------------------------
sleep_LM_I <- Lm(Reaction_I ~ Days, data = sleepstudy)
sleep_LMmer_I <- mtram(sleep_LM_I, ~ (Days | Subject), data = sleepstudy)
## ----mtram-sleep-interval-results---------------------------------------------
logLik(sleep_LMmer_I)
coef(sleep_LMmer_I)
coef(sleep_LMmer)
## ----mtram-sleep_BoxCox, cache = FALSE----------------------------------------
sleep_BC <- BoxCox(Reaction ~ Days, data = sleepstudy)
sleep_BCmer <- mtram(sleep_BC, ~ (Days | Subject), data = sleepstudy,
Hessian = TRUE)
logLik(sleep_BCmer)
## ----mtram-sleep_BoxCoxPlot, echo = FALSE, fig.height = 5---------------------
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
plot(tmp, newdata = data.frame(Days = 0), type = "trafo", col = "black",
xlab = "Average reaction time (in ms)", ylab = expression(h(y)))
rug(sleepstudy$Reaction, col = rgb(.1, .1, .1, .1))
## ----mtram-sleep_marginal, fig.height = 5, fig.width = 7, echo = FALSE--------
days <- 0:9
q <- seq(from = min(sleepstudy$Reaction), to = max(sleepstudy$Reaction),
length.out = 100)
meth <- c("Normal linear mixed-effects model", "Non-normal linear transformation model")
ex <- expand.grid(Reaction = q, Days = days,
Method = factor(meth, levels = meth, labels = meth))
U <- cbind(1, days)
### Linear
tmp <- as.mlt(sleep_LM)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_LMmer)[names(cf)]
SLM <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SLM))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_LM <- pnorm(t(t(h) / sd ))
### BoxCox
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_BC <- pnorm(t(t(h) / sd ))
ex$prob <- c(prob_LM, prob_BC)
plotfun <- function(prob, data, ...) {
fm <- as.formula(paste(prob, "~ Reaction | Days"))
xyplot(fm, data = data, type = "l",
panel = function(x, y, subscripts, ...) {
tmp <- subset(sleepstudy, Days == unique(nd[subscripts, "Days"]))
mypanel(tmp$Reaction, rep(1, nrow(tmp)), lwd = 3, col = grey)
panel.xyplot(x, y, subscripts = subscripts, ...)
}, col = col, xlab = "Average reaction time (in ms)",
ylab = "Marginal distribution function", lwd = 2, groups = Method, ...)
}
grey <- rgb(.75, .75, .75)
nd <- ex
plotfun("prob", ex, layout = c(5, 2), par.settings = simpleTheme(col=c(grey, col), lwd = 3),
auto.key = list(text = c("Empirical cumulative distribution function", levels(nd$Method)),
points = FALSE, lines = TRUE))
## ----mtram-sleep_corr---------------------------------------------------------
cov2cor(sleep_BCmer$G)
## ----mtram-sleep_vcov---------------------------------------------------------
library("mvtnorm")
VC <- solve(sleep_BCmer$Hessian)
idx <- (nrow(VC) - 2):nrow(VC)
Rcoef <- rmvnorm(1000, mean = coef(sleep_BCmer), sigma = VC)[,idx]
ret <- apply(Rcoef, 1, function(gamma) {
L <- matrix(c(gamma[1:2], 0, gamma[3]), nrow = 2)
V <- tcrossprod(L)
c(diag(V), cov2cor(V)[1,2])
})
## ----mtram-sleep_ci-----------------------------------------------------------
### variance random intercept
quantile(ret[1,], c(.025, .5, .975))
### variance random slope
quantile(ret[2,], c(.025, .5, .975))
### correlation random intercept / random slope
quantile(ret[3,], c(.025, .5, .975))
## ----mtram-sleep_Colr, cache = FALSE------------------------------------------
sleep_C <- Colr(Reaction ~ Days, data = sleepstudy)
sleep_Cmer <- mtram(sleep_C, ~ (Days | Subject), data = sleepstudy)
logLik(sleep_Cmer)
## ----mtram-sleep_marginal-2, fig.height = 5, fig.width = 7, echo = FALSE------
days <- 0:9
q <- seq(from = min(sleepstudy$Reaction), to = max(sleepstudy$Reaction),
length.out = 100)
meth <- c("Normal linear mixed-effects model", "Probit transformation model",
"Marginal logit transformation model")
ex <- expand.grid(Reaction = q, Days = days,
Method = factor(meth, levels = meth, labels = meth))
U <- cbind(1, days)
### Linear
tmp <- as.mlt(sleep_LM)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_LMmer)[names(cf)]
SLM <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SLM))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_LM <- pnorm(t(t(h) / sd ))
### BoxCox
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_BC <- pnorm(t(t(h) / sd ))
### Colr
tmp <- as.mlt(sleep_C)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_Cmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_Cmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_C <- plogis(t(t(h) / sd ))
ex$prob <- c(prob_LM, prob_BC, prob_C)
plotfun <- function(prob, data, ...) {
fm <- as.formula(paste(prob, "~ Reaction | Days"))
xyplot(fm, data = data, type = "l",
panel = function(x, y, subscripts, ...) {
tmp <- subset(sleepstudy, Days == unique(nd[subscripts, "Days"]))
mypanel(tmp$Reaction, rep(1, nrow(tmp)), lwd = 3, col = grey)
panel.xyplot(x, y, subscripts = subscripts, ...)
}, col = c(col, col[2]), xlab = "Average reaction time (in ms)",
ylab = "Marginal distribution function", lwd = 2, groups = Method, lty =
c(1, 1, 3), ...)
}
grey <- rgb(.75, .75, .75)
nd <- ex
plotfun("prob", ex, layout = c(5, 2), par.settings = simpleTheme(col=c(grey, col, col[2]), lwd =
3, lty = c(1, 1, 1, 3)),
auto.key = list(text = c("Empirical cumulative distribution function", levels(nd$Method)),
points = FALSE, lines = TRUE))
## ----mtram-toenail-plot, echo = FALSE, cache = FALSE--------------------------
data("toenail", package = "HSAUR3")
layout(matrix(1:2, ncol = 2))
trt <- levels(toenail$treatment)
cdplot(outcome ~ time, data = subset(toenail, treatment == trt[1]),
main = trt[1], xlab = "Time", ylab = "Toe nail infection")
cdplot(outcome ~ time, data = subset(toenail, treatment == trt[2]),
main = trt[2], xlab = "Time", ylab = "")
## ----mtram-toenail_glmer_RI, cache = FALSE------------------------------------
### Laplace
toenail_glmer_RI_1 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 1)
summary(toenail_glmer_RI_1)
toenail_glmer_RI_1@theta
### Adaptive Gaussian Quadrature
toenail_glmer_RI_2 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 20)
summary(toenail_glmer_RI_2)
toenail_glmer_RI_2@theta
## ----mtram-toenail_glmmTMB_RI, cache = FALSE----------------------------------
library("glmmTMB")
toenail_glmmTMB_RI_3 <-
glmmTMB(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"))
summary(toenail_glmmTMB_RI_3)
## ----mtram-toenail_mtram_RI, cache = FALSE------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_RI <-
mtram(toenail_probit, ~ (1 | patientID),
data = toenail, Hessian = TRUE)
coef(toenail_mtram_RI)
## ----mtram-toenail-hessian----------------------------------------------------
vcov(toenail_glmer_RI_2)
solve(toenail_mtram_RI$Hessian)[1:4, 1:4]
## ----mtram-toenail-coef-------------------------------------------------------
cf <- coef(toenail_mtram_RI)
cf[2:4] / sqrt(1 + cf["gamma1"]^2)
## ----mtram-toenail-gee-probit-------------------------------------------------
library("geepack")
gin <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, corstr = "independence",
family = binomial(link = "probit"))
gex <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "exchangeable",
family = binomial(link = "probit"))
gun <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "unstructured",
family = binomial(link = "probit"))
## ----mtram-toenail-gee-probit-coef--------------------------------------------
cbind(mtram = cf[2:4] / sqrt(1 + cf["gamma1"]^2),
indep = coef(gin)[-1],
excha = coef(gex)[-1],
unstr = coef(gun)[-1])
## ----mtram-toenail-gee-logit--------------------------------------------------
gin <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, corstr = "independence",
family = binomial())
gex <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "exchangeable",
family = binomial())
gun <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "unstructured",
family = binomial())
## ----mtram-toenail-gee-logit-coef---------------------------------------------
coef(gin)
coef(gex)
coef(gun)
## ----mtram-toenail_logit, cache = FALSE---------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Logistic", negative = TRUE)
toenail_logit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_logit <- mtram(toenail_logit, ~ (1 | patientID),
data = toenail, Hessian = TRUE)
## ----mtram-toenail_logit-coef-------------------------------------------------
cf <- coef(toenail_mtram_logit)
cf[2:4] / sqrt(1 + cf["gamma1"]^2)
## ----mtram-toenail-trt--------------------------------------------------------
S <- rmvnorm(10000, mean = coef(toenail_mtram_logit),
sigma = solve(toenail_mtram_logit$Hessian))
(ci <- quantile(S[,"treatmentterbinafine:time"] / sqrt(1 + S[, "gamma1"]^2),
prob = c(.025, .975)))
## ----mtram-toenail-gee-logit-mcoef--------------------------------------------
cbind(mtram = cf[2:4] / sqrt(1 + cf["gamma1"]^2),
indep = coef(gin)[-1],
excha = coef(gex)[-1],
unstr = coef(gun)[-1])
## ----mtram-GEE-CI-------------------------------------------------------------
exp(coef(gun)["treatmentterbinafine:time"] +
c(-1, 1) * qnorm(.975) * sqrt(diag(vcov(gun)))["treatmentterbinafine:time"])
## ----mtram-toenail_marginal_logit_s-------------------------------------------
tmp <- toenail_logit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_logit))[1]
cf <- coef(toenail_mtram_logit)[names(cf)] / sqrt(sdrf^2 + 1)
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
time <- 0:180/10
treatment <- sort(unique(toenail$treatment))
nd <- expand.grid(time = time, treatment = treatment)
nd$prob_logit <- predict(tmp, newdata = nd, type = "distribution")[1,]
nd$odds <- exp(predict(tmp, newdata = nd, type = "trafo")[1,])
## ----mtram-toenail_OR_2, dev = "png", cache = FALSE, echo = FALSE, dpi = 300----
X <- model.matrix(~ treatment * time, data = nd)
rbeta <- rmvnorm(10000, mean = coef(toenail_mtram_logit),
sigma = solve(toenail_mtram_logit$Hessian))
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
odds <- exp(-X %*% t(rbeta))
OR <- odds[1:length(time),] / odds[-(1:length(time)),]
plot(time, rep(0, length(time)), ylim = range(OR), type = "n",
xlab = "Time", ylab = "Odds ratio")
colgrey <- rgb(.1, .1, .1, .01)
out <- apply(OR, 2, function(x) lines(time, x, col = colgrey))
ORest <- nd$odds[1:length(time)] / nd$odds[-(1:length(time))]
lines(time, ORest, col = col[1], lwd = 2)
## ----mtram-toenail_marginal_logit---------------------------------------------
tmp <- toenail_logit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_logit))[1]
cf <- coef(toenail_mtram_logit)[names(cf)]
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
pr <- predict(tmp, newdata = nd, type = "distribution")[1,]
nd$prob_logit <- pnorm(qnorm(pr) / sdrf)
## ----mtram-toenail_marginal_probit--------------------------------------------
tmp <- toenail_probit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_RI))[1]
cf <- coef(toenail_mtram_RI)[names(cf)] / sqrt(sdrf^2 + 1)
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
nd$prob_probit <- predict(tmp, newdata = nd, type = "distribution")[1,]
## ----mtram-toenail_probplot, echo = FALSE-------------------------------------
nd2 <- nd[rep(1:nrow(nd), 2),]
nd2$prob <- c(nd$prob_probit, nd$prob_logit)
lev <- c("Probit", "Logit")
nd2$model <- rep(factor(lev, labels = lev, levels = lev), each = nrow(nd))
xyplot(prob ~ time | model, data = nd2, group = treatment, ylim = c(0, 1),
xlab = "Time",
par.settings = simpleTheme(col = col),
auto.key = list(text = levels(nd2$treatment),
points = FALSE, lines = TRUE),
col = col, type = "l", ylab = "Probability (none or mild)")
## ----mtram-toenail-subset-----------------------------------------------------
(rlev <- levels(toenail$patientID)[xtabs(~ patientID,
data = toenail) == 1])
toenail_gr1 <- subset(toenail, !patientID %in% rlev)
toenail_gr1$patientID <- toenail_gr1$patientID[, drop = TRUE]
## ----mtram-toenail_glmer_RS, cache = FALSE------------------------------------
toenail_glmer_RS <-
glmer(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit"))
summary(toenail_glmer_RS)
toenail_glmer_RS@theta
## ----mtram-toenail_glmmTMB_RS, cache = FALSE----------------------------------
toenail_glmmTMB_RS_1 <-
glmmTMB(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit"))
summary(toenail_glmmTMB_RS_1)
## ----mtram-toenail_mtram_RS, cache = FALSE------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail_gr1),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail_gr1,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_RS <-
mtram(toenail_probit, ~ (1 + time | patientID),
data = toenail_gr1)
logLik(toenail_mtram_RS)
coef(toenail_mtram_RS)
## ----toenail-comparisons, cache = FALSE, echo = FALSE, results = "hide"-------
t1 <- system.time(toenail_glmer_RI_1 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 1))
t2 <- system.time(toenail_glmer_RI_2 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 20))
t3 <- system.time(toenail_glmmTMB_RI_3 <-
glmmTMB(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit")))
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
t4 <- system.time(toenail_mtram_RI <-
mtram(toenail_probit, ~ (1 | patientID),
data = toenail, Hessian = TRUE))
t5 <- system.time(toenail_glmer_RS <-
glmer(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit")))
t6 <- system.time(toenail_glmmTMB_RS_1 <-
glmmTMB(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit")))
m <- ctm(as.basis(~ outcome, data = toenail_gr1),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail_gr1,
fixed = c("outcomemoderate or severe" = 0))
t7 <- system.time(toenail_mtram_RS <-
mtram(toenail_probit, ~ (1 + time | patientID),
data = toenail_gr1))
## ----output, echo = FALSE------------------------------------------------
tn_RI_glmer_L <- c(fixef(toenail_glmer_RI_1), toenail_glmer_RI_1@theta, 0, 0)
tn_RI_glmer_A <- c(fixef(toenail_glmer_RI_2), toenail_glmer_RI_2@theta, 0, 0)
tn_RI_glmmTMB <- c(fixef(toenail_glmmTMB_RI_3)$cond, sqrt(VarCorr(toenail_glmmTMB_RI_3)$cond$patientID), 0, 0)
tn_RI_mlt <- c(coef(toenail_mtram_RI), 0, 0)
tn_RS_glmer <- c(fixef(toenail_glmer_RS), toenail_glmer_RS@theta)
tn_RS_glmmTMB <- c(fixef(toenail_glmer_RS), chol(VarCorr(toenail_glmmTMB_RS_1)$cond$patientID)[c(1,3, 4)])
tn_RS_mlt <- coef(toenail_mtram_RS)
tn <- cbind(tn_RI_glmer_L, tn_RI_glmer_A , tn_RI_glmmTMB, tn_RI_mlt ,
tn_RS_glmer, tn_RS_glmmTMB, tn_RS_mlt)
logLik(toenail_glmer_RI_1)
logLik(toenail_glmer_RI_2)
logLik(toenail_glmmTMB_RI_3)
logLik(toenail_mtram_RI)
logLik(toenail_glmer_RS)
logLik(toenail_glmmTMB_RS_1)
logLik(toenail_mtram_RS)
ll <- c(
### logLik of transformation model for glmer (Laplace) parameters
logLik(toenail_mtram_RI, tn_RI_glmer_L[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model for glmer (AGQ) parameters
logLik(toenail_mtram_RI, tn_RI_glmer_A[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model for glmmTMB (Laplace) parameters
logLik(toenail_mtram_RI, tn_RI_glmmTMB[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model
logLik(toenail_mtram_RI),
### logLik of transformation model for glmer (Laplace) parameters
logLik(toenail_mtram_RS, tn_RS_glmer * c(-1, rep(1, 6))),
### logLik of transformation model for glmmTMB (Laplace) parameters
logLik(toenail_mtram_RS, tn_RS_glmmTMB * c(-1, rep(1, 6))),
### logLik of transformation model
logLik(toenail_mtram_RS))
tm <- c(t1["user.self"],
t2["user.self"],
t3["user.self"],
t4["user.self"],
t5["user.self"],
t6["user.self"],
t7["user.self"])
tm <- formatC(tm, format = "f", digits = 2, width = 5)
tn <- formatC(tn, format = "f", digits = 2, width = 5)
ll <- formatC(ll, format = "f", digits = 2, width = 6)
tn <- cbind(c("$\\alpha$", "$\\eshiftparm_1$", "$\\eshiftparm_2$", "$\\eshiftparm_3$", "$\\gamma_1$", "$\\gamma_2$", "$\\gamma_3$"), tn)
ret <- c("
\\begin{tabular}{lrrrr|rrr} \\\\ \\hline
& \\multicolumn{4}{c|}{RI} & \\multicolumn{3}{c}{RI + RS} \\\\
& \\texttt{glmer} & \\texttt{glmer} & \\texttt{glmmTMB} & & \\texttt{glmer} & \\texttt{glmmTMB} & \\\\
& L & AGQ & L & (7) & L & L & (7) \\\\ \\hline")
ret <- c(ret, apply(tn, 1, function(x) c(paste(x, collapse = " & "), "\\\\")))
ret <- c(ret, "\\hline")
ret <- c(ret,
paste("LogLik &", paste(ll, collapse = "&"), "\\\\ "),
paste("Time (sec) &", paste(tm, collapse = "&"), "\\\\ \\hline"),
"\\end{tabular}")
## ----table, echo = FALSE, results = "asis"------------------------------------
cat(ret, sep = "\n")
## ----mtram-neck_plot, echo = FALSE, fig.height = 4, fig.width = 7-------------
data("neck_pain", package = "ordinalCont")
pain_df <- neck_pain
idt <- xtabs(~ id, data = pain_df)
miss <- names(idt)[idt < 3]
pain_df <- subset(pain_df, !id %in% miss)
pain_df$id <- factor(pain_df$id)
levels(pain_df$laser) <- c("Placebo", "Active")
levels(pain_df$time) <- c("Baseline", "7 weeks", "12 weeks")
pain <- rbind(subset(pain_df, laser == levels(pain_df$laser)[1]),
subset(pain_df, laser == levels(pain_df$laser)[2]))
p1 <- xyplot(vas ~ time | laser, data = pain,
groups = id, type = "l", col = rgb(.1, .1, .1, .1),
lwd = 2, layout = c(2, 1),
ylab = "Neck pain (on visual analog scale)", xlab = "Examinations")
plot(p1)
## ----mtram-ordinalCont--------------------------------------------------------
library("ordinalCont")
## ----mtram-neck_ocm, cache = FALSE, results = "hide"--------------------------
neck_ocm <- ocm(vas ~ laser * time + (1 | id), data = pain_df,
scale = c(0, 1))
## ----mtram-neck_ocm_summary---------------------------------------------------
summary(neck_ocm)
## ----mtram-neck_ocm_ci--------------------------------------------------------
exp(cbind(coef(neck_ocm)[2:6], confint(neck_ocm)[2:6,]))
## ----tramME-neck--------------------------------------------------------------
library("tramME")
neck_ColrME <- ColrME(vas ~ laser * time + (1 | id), data = pain_df,
bounds = c(0, 1), support = c(0, 1))
## ----tramME-neck_ci-----------------------------------------------------------
exp(coef(neck_ColrME))
## ----mtram-neck_Colr, cache = FALSE-------------------------------------------
neck_Colr <- Colr(vas ~ laser * time, data = pain_df,
bounds = c(0, 1), support = c(0, 1),
extrapolate = TRUE)
neck_Colrmer <- mtram(neck_Colr, ~ (1 | id), data = pain_df,
Hessian = TRUE)
## ----mtram-neck_Colr_distr, echo = FALSE, fig.height = 4, fig.width = 7-------
nd <- expand.grid(laser = unique(pain_df$laser),
time = unique(pain_df$time))
tmp <- as.mlt(neck_Colr)
coef(tmp)[] <- coef(neck_Colrmer)[1:12]
q <- 1:99/100
nd2 <- expand.grid(vas = q, laser = unique(pain_df$laser),
time = unique(pain_df$time))
# tmp <- as.mlt(neck_Colr)
sd <- sqrt(coef(neck_Colrmer)[13]^2 + 1)
prb <- predict(tmp, newdata = nd, type = "distribution", q = q)
nd2$prob <- c(pnorm(qnorm(prb) / sd))
p2 <- xyplot(prob ~ vas | time, data = nd2, groups = laser, type = "l",
col = col,
layout = c(3, 1),
xlab = "Neck pain (on visual analog scale)",
ylab = "Marginal distribution",
par.settings = simpleTheme(col=col),
auto.key = list(text = levels(nd2$laser),
points = FALSE, lines = TRUE))
plot(p2)
## M1
# neck_Colrmer <- mtram(neck_Colr, ~ (1 | id), data = pain_df,
# Hessian = TRUE, standardise = FALSE)
# logLik(neck_Colrmer)
#
# nd <- expand.grid(laser = unique(pain_df$laser),
# time = unique(pain_df$time))
# q <- 1:99/100
# nd2 <- expand.grid(vas = q, laser = unique(pain_df$laser),
# time = unique(pain_df$time))
# tmp <- as.mlt(neck_Colr)
# coef(tmp)[] <- coef(neck_Colrmer)[1:12]
# sd <- sqrt(coef(neck_Colrmer)[13]^2 + 1)
# prb <- predict(tmp, newdata = nd, type = "distribution", q = q)
# nd2$prob <- c(pnorm(qnorm(prb) / sd))
# p2 <- xyplot(prob ~ vas | time, data = nd2, groups = laser, type = "l",
# col = col, ylim = c(-0.05, 1.05),
# layout = c(3, 1),
# xlab = "Neck pain (on visual analog scale)",
# ylab = "Marginal distribution",
# par.settings = simpleTheme(col=col),
# auto.key = list(text = levels(nd2$laser),
# points = FALSE, lines = TRUE))
# plot(p2)
## ----mtram-neck_Colr-CI, echo = TRUE, eval=TRUE-------------------------------
S <- solve(neck_Colrmer$Hessian)
rbeta <- rmvnorm(10000, mean = coef(neck_Colrmer), sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
t(apply(rbeta[, 8:12], 2, function(x) {
quantile(exp(x),prob = c(.025, .5, .975))}))
## ----mtram-neck_Colr-std_beta-------------------------------------------------
beta <- coef(neck_Colrmer)[8:12]
alpha <- coef(neck_Colrmer)[13]
(std_beta <- cbind(beta / sqrt(1 + alpha^2)))
## ----mtram-neck_Colr-ctr_mat--------------------------------------------------
ctr_mat <- matrix(c(1, 0, 0, 0, 0,
1, 0, 0, 1, 0,
1, 0, 0, 0, 1), nrow = 3, byrow = TRUE)
ctr_mat %*% std_beta
## ----mtram-neck_PImanual, eval=TRUE-------------------------------------------
(ci_emp <- t(apply(ctr_mat %*% t(rbeta[, 8:12]), 1, function(x) {
quantile(x, prob = c(.025, .5, .975))})))
PI(-ci_emp, link = "logistic")
## ----neck_Colr-PI, echo = TRUE------------------------------------------------
nd <- expand.grid(time = unique(pain_df$time),
laser = unique(pain_df$laser))
neck_Colr_marg <- neck_Colr
neck_Colr_marg$coef <- coef(neck_Colrmer)[1:12] /
sqrt(coef(neck_Colrmer)[13]^2 + 1)
(neck_Colr_PI <- PI(neck_Colr_marg, newdata = nd[1:3, ],
reference = nd[4:6, ],
one2one = TRUE, conf.level = .95))[1:3, 1:3]
## ----mtram-CAO, echo = FALSE--------------------------------------------------
dir <- system.file("rda", package = "TH.data")
load(file.path(dir, "Primary_endpoint_data.rda"))
## ----mtram-CAO-plot, cache = FALSE, echo = FALSE------------------------------
ra <- sort(unique(CAOsurv$randarm))
st <- sort(unique(CAOsurv$strat_t))
sn <- sort(unique(CAOsurv$strat_n))
su <- c(1, 1700)
add <- c(0, max(CAOsurv$iDFS[, "time2"]) - su[2])
ylim <- c(-.05, 1.05)
tmp <- as.mlt(Coxph(iDFS | 0 + strat_n:strat_t:randarm ~ 1, data = CAOsurv,
support = su, add = add, log_first = TRUE))
nd <- expand.grid(strat_n = sn, strat_t = st, randarm = ra)
q <- mkgrid(tmp, 100)[[1]]
surv <- predict(tmp, newdata = nd, type = "survivor", q = q)
nd <- nd[rep(1:nrow(nd), each = nrow(surv)),]
nd$time <- q
nd$surv <- c(surv)
xyplot(surv ~ time | strat_t + strat_n, data = nd, groups = randarm,
type = "l", ylim = c(0, 1), col = col, ylab = "Probability",
xlab = "Time (in days)",
par.settings = simpleTheme(col=col),
auto.key = list(text = levels(nd$randarm),
points = FALSE, lines = TRUE))
## ----mtram-CAO_DFS------------------------------------------------------------
### convert "exact" event dates to interval-censoring (+/- one day)
tmp <- CAOsurv$iDFS
exact <- tmp[,3] == 1
tmp[exact,2] <- tmp[exact,1] + 2
tmp[exact,1] <- pmax(tmp[exact,1] - 2, 0)
tmp[exact,3] <- 3
CAOsurv$iDFS2 <- tmp
## ----mtram-CAO_SR, cache = TRUE-----------------------------------------------
CAO_SR <- Survreg(iDFS2 ~ randarm, data = CAOsurv)
CAO_SR_mtram <- mtram(CAO_SR, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_SR_mtram)
(cf <- coef(CAO_SR_mtram))
(OR <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
## ----mtram-CAO-CI-------------------------------------------------------------
S <- solve(CAO_SR_mtram$Hessian)
# sqrt(diag(S))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_mtram),
sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[, -ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_SR_2, cache = TRUE---------------------------------------------
CAO_SR_2 <- Survreg(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv)
CAO_SR_2_mtram <- mtram(CAO_SR_2, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_SR_2_mtram)
(cf <- coef(CAO_SR_2_mtram))
(OR_2 <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
## ----mtram-CAO-CI-2, echo = FALSE---------------------------------------------
H <- CAO_SR_2_mtram$Hessian
S <- solve(H + .01 * diag(nrow(H)))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_2_mtram),
sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[, -ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_Cox_2, cache = TRUE--------------------------------------------
CAO_Cox_2 <- Coxph(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv,
support = c(1, 1700), log_first = TRUE, order = 4)
logLik(CAO_Cox_2)
CAO_Cox_2_mtram <- mtram(CAO_Cox_2, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_Cox_2_mtram)
coef(CAO_Cox_2_mtram)
## ----mtram-CAO-CI-3, echo = FALSE---------------------------------------------
H <- CAO_Cox_2_mtram$Hessian
S <- solve(H + .1 * diag(nrow(H)))
rbeta <- rmvnorm(10000, mean = coef(CAO_Cox_2_mtram),
sigma = S)
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_PI-------------------------------------------------------------
nd <- CAOsurv[1:2, ]
tmp <- CAO_Cox_2
tmp$coef <- coef(CAO_Cox_2_mtram)[-22] / sqrt(coef(CAO_Cox_2_mtram)[22]^2 + 1)
(CAO_Cox_PI <- PI(tmp, newdata = nd[2, ], reference = nd[1, ],
one2one = TRUE, conf.level = .95))[1, ]
## ----mtram-CAO_PI_man---------------------------------------------------------
ci_man <- quantile(-rbeta[, ncol(rbeta)], prob = c(.025, .5, .975))
(CAO_Cox_PIm <- PI(ci_man, link = "minimum extreme value"))
## ----tramME-CAO_SR, cache = TRUE----------------------------------------------
CAO_Cox_2_tramME <- CoxphME(iDFS2 | 0 + strat_n:strat_t ~ randarm + (1 | Block),
data = CAOsurv, log_first = TRUE)
## ----tramME-CAO_SR-hr, cache = TRUE-------------------------------------------
exp(coef(CAO_Cox_2_tramME))
exp(confint(CAO_Cox_2_tramME, parm = "randarm5-FU + Oxaliplatin",
estimate = TRUE))
## ----echo=FALSE, eval=FALSE---------------------------------------------------
## sqrt(VarCorr(CAO_Cox_2_tramME)$Block$var)
## coef(CAO_Cox_2_mtram)["gamma1"]
## ----mtram-CAO-coxme, echo = FALSE, eval = FALSE------------------------------
## library("coxme")
## m <- coxme(DFS ~ randarm + (1 | Block), data = CAOsurv)
## summary(m)
## sd <- sqrt(diag(vcov(m)))
## exp(coef(m) + c(-1, 0, 1) * qnorm(.975) * sd)
## ----sim, eval = FALSE--------------------------------------------------------
## source(system.file("simulations", "mtram_sim.R", package = "tram"), echo = TRUE)
## ----mtram-sessionInfo, echo = FALSE, results = "hide"------------------------
sessionInfo()
## ----mtram-funs, echo = FALSE, results = "hide"-------------------------------
if (file.exists("packages.bib")) file.remove("packages.bib")
pkgversion <- function(pkg) {
pkgbib(pkg)
packageDescription(pkg)$Version
}
pkgbib <- function(pkg) {
x <- citation(package = pkg, auto = TRUE)[[1]]
b <- toBibtex(x)
b <- gsub("Buehlmann", "B{\\\\\"u}hlmann", b)
b[1] <- paste("@Manual{pkg:", pkg, ",", sep = "")
if (is.na(b["url"])) {
b[length(b)] <- paste(" URL = {http://CRAN.R-project.org/package=",
pkg, "}", sep = "")
b <- c(b, "}")
}
cat(b, sep = "\n", file = "packages.bib", append = TRUE)
}
pkg <- function(pkg) {
vrs <- try(pkgversion(pkg))
if (inherits(vrs, "try-error")) return(NA)
paste("\\\\pkg{", pkg, "} \\\\citep[version~",
vrs, ",][]{pkg:", pkg, "}", sep = "")
}
pkg("mlt")
pkg("tram")
pkg("SparseGrid")
cat(c("@Manual{vign:mlt.docreg,",
" title = {Most Likely Transformations: The mlt Package},",
" author = {Torsten Hothorn},",
paste(" year = ", substr(packageDescription("mlt.docreg")$Date, 1, 4), ",", sep = ""),
paste(" note = {R package vignette version ", packageDescription("mlt.docreg")$Version, "},", sep = ""),
" url = {https://CRAN.R-project.org/package=mlt.docreg},",
"}"), file = "packages.bib", append = TRUE, sep = "\n")
Try the tram package in your browser
Any scripts or data that you put into this service are public.
tram documentation built on Aug. 25, 2023, 5:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----mtram-pkgs, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
set.seed(290875)
pkgs <- c("colorspace", "survival", "lme4", "tram", "gridExtra",
"lattice", "latticeExtra", "mvtnorm", "ordinalCont", "tramME")
pkgs <- pkgs %in% installed.packages()
## ----mtram-citation, echo = FALSE---------------------------------------------
year <- substr(packageDescription("tram")$Date, 1, 4)
version <- packageDescription("tram")$Version
## ----fail, results = "asis", echo = FALSE-------------------------------------
if (any(!pkgs))
{
cat(paste("Package(s)", paste(names(pkgs)[!pkgs], collapse = ", "),
"not available, stop processing.",
"\\end{document}\n"))
knitr::knit_exit()
}
## ----mtram-setup, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
library("lattice")
library("latticeExtra")
library("gridExtra")
trellis.par.set(list(plot.symbol = list(col=1,pch=20, cex=0.7),
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)
knitr::opts_chunk$set(echo = TRUE, results = 'markup', error = FALSE,
warning = FALSE, message = FALSE,
tidy = FALSE, cache = FALSE, size = "small",
fig.width = 6, fig.height = 4, fig.align = "center",
out.width = NULL, ###'.6\\linewidth',
out.height = NULL,
fig.scap = NA)
knitr::render_sweave() # use Sweave environments
knitr::set_header(highlight = '') # do not \usepackage{Sweave}
## R settings
options(prompt = "R> ", continue = "+ ", useFancyQuotes = FALSE) # JSS style
options(width = 75)
### ecdf plots
myprepanel <- function (x, y, f.value = NULL, ...)
{
ans <- prepanel.default.qqmath(x, f.value = f.value, distribution = qunif)
with(ans, list(xlim = ylim, ylim = c(0, 1), dx = dy, dy = dx))
}
mypanel <- function (x, y, f.value = NULL, type = "s", groups = NULL, qtype = 7,
ref = TRUE, ...)
{
if (ref) {
reference.line <- trellis.par.get("reference.line")
do.call(panel.abline, c(list(h = c(0, 1)), reference.line))
}
x <- as.numeric(x)
distribution <- qunif
nobs <- sum(!is.na(x))
if (!is.null(groups)) {
panel.superpose(x, y = NULL, f.value = f.value, type = type,
distribution = distribution, qtype = qtype, groups = groups,
panel.groups = panel.ecdfplot, ...)
}
else if (nobs) {
if (is.null(f.value)) {
panel.xyplot(x = sort(x), y = cumsum(y[order(x)]) / sum(y),
type = type, ...)
}
else {
p <- if (is.numeric(f.value))
f.value
else f.value(nobs)
panel.xyplot(x = quantile(x, p, names = FALSE, type = qtype,
na.rm = TRUE), y = distribution(p), type = type,
...)
}
}
}
library("colorspace")
col <- diverge_hcl(2, h = c(246, 40), c = 120, l = c(65, 90), alpha = .75)
## ----mtram-vignette, eval = FALSE---------------------------------------------
## install.packages("tram")
## demo("mtram", package = "tram")
## ----mtram-sleep-plot, echo = FALSE-------------------------------------------
library("lme4")
xyplot(Reaction ~ Days | Subject, data = sleepstudy,
xlab = "Days of sleep deprivation", ylab = "Average reaction time (in ms)")
## ----mtram-sleep_lmer, cache = FALSE------------------------------------------
sleep_lmer <- lmer(Reaction ~ Days + (Days | Subject),
data = sleepstudy, REML = FALSE)
## ----mtram-tram---------------------------------------------------------------
library("tram")
## ----mtram-sleep_mtram, cache = FALSE-----------------------------------------
sleep_LM <- Lm(Reaction ~ Days, data = sleepstudy)
sleep_LMmer <- mtram(sleep_LM, ~ (Days | Subject), data = sleepstudy)
## ----mtram-sleep_cmp----------------------------------------------------------
logLik(sleep_lmer)
logLik(sleep_LMmer)
## ----mtram-sleep_sd-----------------------------------------------------------
(sdinv <- 1 / summary(sleep_lmer)$sigma)
coef(sleep_LMmer)["Reaction"]
## ----mtram-sleep_beta---------------------------------------------------------
fixef(sleep_lmer) * c(-1, 1) * sdinv
coef(sleep_LMmer)[c("(Intercept)", "Days")]
## ----mtram-sleep_varparm------------------------------------------------------
sleep_lmer@theta
coef(sleep_LMmer)[-(1:3)]
## ----mtram-sleep_varcorr------------------------------------------------------
sleep_LMmer$G * (1 / sdinv)^2
cov2cor(sleep_LMmer$G * (1 / sdinv)^2)
unclass(VarCorr(sleep_lmer))$Subject
## ----mtram-sleep-Surv---------------------------------------------------------
library("survival")
sleepstudy$Reaction_I <- with(sleepstudy, Surv(Reaction - 20, Reaction + 20,
type = "interval2"))
sleepstudy$Reaction_I[1:5]
## ----mtram-sleep-interval, cache = FALSE--------------------------------------
sleep_LM_I <- Lm(Reaction_I ~ Days, data = sleepstudy)
sleep_LMmer_I <- mtram(sleep_LM_I, ~ (Days | Subject), data = sleepstudy)
## ----mtram-sleep-interval-results---------------------------------------------
logLik(sleep_LMmer_I)
coef(sleep_LMmer_I)
coef(sleep_LMmer)
## ----mtram-sleep_BoxCox, cache = FALSE----------------------------------------
sleep_BC <- BoxCox(Reaction ~ Days, data = sleepstudy)
sleep_BCmer <- mtram(sleep_BC, ~ (Days | Subject), data = sleepstudy,
Hessian = TRUE)
logLik(sleep_BCmer)
## ----mtram-sleep_BoxCoxPlot, echo = FALSE, fig.height = 5---------------------
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
plot(tmp, newdata = data.frame(Days = 0), type = "trafo", col = "black",
xlab = "Average reaction time (in ms)", ylab = expression(h(y)))
rug(sleepstudy$Reaction, col = rgb(.1, .1, .1, .1))
## ----mtram-sleep_marginal, fig.height = 5, fig.width = 7, echo = FALSE--------
days <- 0:9
q <- seq(from = min(sleepstudy$Reaction), to = max(sleepstudy$Reaction),
length.out = 100)
meth <- c("Normal linear mixed-effects model", "Non-normal linear transformation model")
ex <- expand.grid(Reaction = q, Days = days,
Method = factor(meth, levels = meth, labels = meth))
U <- cbind(1, days)
### Linear
tmp <- as.mlt(sleep_LM)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_LMmer)[names(cf)]
SLM <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SLM))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_LM <- pnorm(t(t(h) / sd ))
### BoxCox
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_BC <- pnorm(t(t(h) / sd ))
ex$prob <- c(prob_LM, prob_BC)
plotfun <- function(prob, data, ...) {
fm <- as.formula(paste(prob, "~ Reaction | Days"))
xyplot(fm, data = data, type = "l",
panel = function(x, y, subscripts, ...) {
tmp <- subset(sleepstudy, Days == unique(nd[subscripts, "Days"]))
mypanel(tmp$Reaction, rep(1, nrow(tmp)), lwd = 3, col = grey)
panel.xyplot(x, y, subscripts = subscripts, ...)
}, col = col, xlab = "Average reaction time (in ms)",
ylab = "Marginal distribution function", lwd = 2, groups = Method, ...)
}
grey <- rgb(.75, .75, .75)
nd <- ex
plotfun("prob", ex, layout = c(5, 2), par.settings = simpleTheme(col=c(grey, col), lwd = 3),
auto.key = list(text = c("Empirical cumulative distribution function", levels(nd$Method)),
points = FALSE, lines = TRUE))
## ----mtram-sleep_corr---------------------------------------------------------
cov2cor(sleep_BCmer$G)
## ----mtram-sleep_vcov---------------------------------------------------------
library("mvtnorm")
VC <- solve(sleep_BCmer$Hessian)
idx <- (nrow(VC) - 2):nrow(VC)
Rcoef <- rmvnorm(1000, mean = coef(sleep_BCmer), sigma = VC)[,idx]
ret <- apply(Rcoef, 1, function(gamma) {
L <- matrix(c(gamma[1:2], 0, gamma[3]), nrow = 2)
V <- tcrossprod(L)
c(diag(V), cov2cor(V)[1,2])
})
## ----mtram-sleep_ci-----------------------------------------------------------
### variance random intercept
quantile(ret[1,], c(.025, .5, .975))
### variance random slope
quantile(ret[2,], c(.025, .5, .975))
### correlation random intercept / random slope
quantile(ret[3,], c(.025, .5, .975))
## ----mtram-sleep_Colr, cache = FALSE------------------------------------------
sleep_C <- Colr(Reaction ~ Days, data = sleepstudy)
sleep_Cmer <- mtram(sleep_C, ~ (Days | Subject), data = sleepstudy)
logLik(sleep_Cmer)
## ----mtram-sleep_marginal-2, fig.height = 5, fig.width = 7, echo = FALSE------
days <- 0:9
q <- seq(from = min(sleepstudy$Reaction), to = max(sleepstudy$Reaction),
length.out = 100)
meth <- c("Normal linear mixed-effects model", "Probit transformation model",
"Marginal logit transformation model")
ex <- expand.grid(Reaction = q, Days = days,
Method = factor(meth, levels = meth, labels = meth))
U <- cbind(1, days)
### Linear
tmp <- as.mlt(sleep_LM)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_LMmer)[names(cf)]
SLM <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SLM))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_LM <- pnorm(t(t(h) / sd ))
### BoxCox
tmp <- as.mlt(sleep_BC)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_BCmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_BCmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_BC <- pnorm(t(t(h) / sd ))
### Colr
tmp <- as.mlt(sleep_C)
cf <- coef(tmp)
coef(tmp) <- coef(sleep_Cmer)[names(cf)]
SBC <- tcrossprod(U %*% as(sleep_Cmer$G, "matrix"), U) + diag(nrow(U))
sd <- sqrt(diag(SBC))
h <- predict(tmp, newdata = data.frame(Days = days), q = q, type = "trafo")
prob_C <- plogis(t(t(h) / sd ))
ex$prob <- c(prob_LM, prob_BC, prob_C)
plotfun <- function(prob, data, ...) {
fm <- as.formula(paste(prob, "~ Reaction | Days"))
xyplot(fm, data = data, type = "l",
panel = function(x, y, subscripts, ...) {
tmp <- subset(sleepstudy, Days == unique(nd[subscripts, "Days"]))
mypanel(tmp$Reaction, rep(1, nrow(tmp)), lwd = 3, col = grey)
panel.xyplot(x, y, subscripts = subscripts, ...)
}, col = c(col, col[2]), xlab = "Average reaction time (in ms)",
ylab = "Marginal distribution function", lwd = 2, groups = Method, lty =
c(1, 1, 3), ...)
}
grey <- rgb(.75, .75, .75)
nd <- ex
plotfun("prob", ex, layout = c(5, 2), par.settings = simpleTheme(col=c(grey, col, col[2]), lwd =
3, lty = c(1, 1, 1, 3)),
auto.key = list(text = c("Empirical cumulative distribution function", levels(nd$Method)),
points = FALSE, lines = TRUE))
## ----mtram-toenail-plot, echo = FALSE, cache = FALSE--------------------------
data("toenail", package = "HSAUR3")
layout(matrix(1:2, ncol = 2))
trt <- levels(toenail$treatment)
cdplot(outcome ~ time, data = subset(toenail, treatment == trt[1]),
main = trt[1], xlab = "Time", ylab = "Toe nail infection")
cdplot(outcome ~ time, data = subset(toenail, treatment == trt[2]),
main = trt[2], xlab = "Time", ylab = "")
## ----mtram-toenail_glmer_RI, cache = FALSE------------------------------------
### Laplace
toenail_glmer_RI_1 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 1)
summary(toenail_glmer_RI_1)
toenail_glmer_RI_1@theta
### Adaptive Gaussian Quadrature
toenail_glmer_RI_2 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 20)
summary(toenail_glmer_RI_2)
toenail_glmer_RI_2@theta
## ----mtram-toenail_glmmTMB_RI, cache = FALSE----------------------------------
library("glmmTMB")
toenail_glmmTMB_RI_3 <-
glmmTMB(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"))
summary(toenail_glmmTMB_RI_3)
## ----mtram-toenail_mtram_RI, cache = FALSE------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_RI <-
mtram(toenail_probit, ~ (1 | patientID),
data = toenail, Hessian = TRUE)
coef(toenail_mtram_RI)
## ----mtram-toenail-hessian----------------------------------------------------
vcov(toenail_glmer_RI_2)
solve(toenail_mtram_RI$Hessian)[1:4, 1:4]
## ----mtram-toenail-coef-------------------------------------------------------
cf <- coef(toenail_mtram_RI)
cf[2:4] / sqrt(1 + cf["gamma1"]^2)
## ----mtram-toenail-gee-probit-------------------------------------------------
library("geepack")
gin <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, corstr = "independence",
family = binomial(link = "probit"))
gex <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "exchangeable",
family = binomial(link = "probit"))
gun <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "unstructured",
family = binomial(link = "probit"))
## ----mtram-toenail-gee-probit-coef--------------------------------------------
cbind(mtram = cf[2:4] / sqrt(1 + cf["gamma1"]^2),
indep = coef(gin)[-1],
excha = coef(gex)[-1],
unstr = coef(gun)[-1])
## ----mtram-toenail-gee-logit--------------------------------------------------
gin <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, corstr = "independence",
family = binomial())
gex <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "exchangeable",
family = binomial())
gun <- geeglm(I((0:1)[outcome]) ~ treatment * time,
id = patientID, data = toenail, cor = "unstructured",
family = binomial())
## ----mtram-toenail-gee-logit-coef---------------------------------------------
coef(gin)
coef(gex)
coef(gun)
## ----mtram-toenail_logit, cache = FALSE---------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Logistic", negative = TRUE)
toenail_logit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_logit <- mtram(toenail_logit, ~ (1 | patientID),
data = toenail, Hessian = TRUE)
## ----mtram-toenail_logit-coef-------------------------------------------------
cf <- coef(toenail_mtram_logit)
cf[2:4] / sqrt(1 + cf["gamma1"]^2)
## ----mtram-toenail-trt--------------------------------------------------------
S <- rmvnorm(10000, mean = coef(toenail_mtram_logit),
sigma = solve(toenail_mtram_logit$Hessian))
(ci <- quantile(S[,"treatmentterbinafine:time"] / sqrt(1 + S[, "gamma1"]^2),
prob = c(.025, .975)))
## ----mtram-toenail-gee-logit-mcoef--------------------------------------------
cbind(mtram = cf[2:4] / sqrt(1 + cf["gamma1"]^2),
indep = coef(gin)[-1],
excha = coef(gex)[-1],
unstr = coef(gun)[-1])
## ----mtram-GEE-CI-------------------------------------------------------------
exp(coef(gun)["treatmentterbinafine:time"] +
c(-1, 1) * qnorm(.975) * sqrt(diag(vcov(gun)))["treatmentterbinafine:time"])
## ----mtram-toenail_marginal_logit_s-------------------------------------------
tmp <- toenail_logit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_logit))[1]
cf <- coef(toenail_mtram_logit)[names(cf)] / sqrt(sdrf^2 + 1)
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
time <- 0:180/10
treatment <- sort(unique(toenail$treatment))
nd <- expand.grid(time = time, treatment = treatment)
nd$prob_logit <- predict(tmp, newdata = nd, type = "distribution")[1,]
nd$odds <- exp(predict(tmp, newdata = nd, type = "trafo")[1,])
## ----mtram-toenail_OR_2, dev = "png", cache = FALSE, echo = FALSE, dpi = 300----
X <- model.matrix(~ treatment * time, data = nd)
rbeta <- rmvnorm(10000, mean = coef(toenail_mtram_logit),
sigma = solve(toenail_mtram_logit$Hessian))
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
odds <- exp(-X %*% t(rbeta))
OR <- odds[1:length(time),] / odds[-(1:length(time)),]
plot(time, rep(0, length(time)), ylim = range(OR), type = "n",
xlab = "Time", ylab = "Odds ratio")
colgrey <- rgb(.1, .1, .1, .01)
out <- apply(OR, 2, function(x) lines(time, x, col = colgrey))
ORest <- nd$odds[1:length(time)] / nd$odds[-(1:length(time))]
lines(time, ORest, col = col[1], lwd = 2)
## ----mtram-toenail_marginal_logit---------------------------------------------
tmp <- toenail_logit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_logit))[1]
cf <- coef(toenail_mtram_logit)[names(cf)]
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
pr <- predict(tmp, newdata = nd, type = "distribution")[1,]
nd$prob_logit <- pnorm(qnorm(pr) / sdrf)
## ----mtram-toenail_marginal_probit--------------------------------------------
tmp <- toenail_probit
cf <- coef(tmp)
cf <- cf[names(cf) != "outcomemoderate or severe"]
sdrf <- rev(coef(toenail_mtram_RI))[1]
cf <- coef(toenail_mtram_RI)[names(cf)] / sqrt(sdrf^2 + 1)
cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
coef(tmp) <- cf
nd$prob_probit <- predict(tmp, newdata = nd, type = "distribution")[1,]
## ----mtram-toenail_probplot, echo = FALSE-------------------------------------
nd2 <- nd[rep(1:nrow(nd), 2),]
nd2$prob <- c(nd$prob_probit, nd$prob_logit)
lev <- c("Probit", "Logit")
nd2$model <- rep(factor(lev, labels = lev, levels = lev), each = nrow(nd))
xyplot(prob ~ time | model, data = nd2, group = treatment, ylim = c(0, 1),
xlab = "Time",
par.settings = simpleTheme(col = col),
auto.key = list(text = levels(nd2$treatment),
points = FALSE, lines = TRUE),
col = col, type = "l", ylab = "Probability (none or mild)")
## ----mtram-toenail-subset-----------------------------------------------------
(rlev <- levels(toenail$patientID)[xtabs(~ patientID,
data = toenail) == 1])
toenail_gr1 <- subset(toenail, !patientID %in% rlev)
toenail_gr1$patientID <- toenail_gr1$patientID[, drop = TRUE]
## ----mtram-toenail_glmer_RS, cache = FALSE------------------------------------
toenail_glmer_RS <-
glmer(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit"))
summary(toenail_glmer_RS)
toenail_glmer_RS@theta
## ----mtram-toenail_glmmTMB_RS, cache = FALSE----------------------------------
toenail_glmmTMB_RS_1 <-
glmmTMB(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit"))
summary(toenail_glmmTMB_RS_1)
## ----mtram-toenail_mtram_RS, cache = FALSE------------------------------------
m <- ctm(as.basis(~ outcome, data = toenail_gr1),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail_gr1,
fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_RS <-
mtram(toenail_probit, ~ (1 + time | patientID),
data = toenail_gr1)
logLik(toenail_mtram_RS)
coef(toenail_mtram_RS)
## ----toenail-comparisons, cache = FALSE, echo = FALSE, results = "hide"-------
t1 <- system.time(toenail_glmer_RI_1 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 1))
t2 <- system.time(toenail_glmer_RI_2 <-
glmer(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit"),
nAGQ = 20))
t3 <- system.time(toenail_glmmTMB_RI_3 <-
glmmTMB(outcome ~ treatment * time + (1 | patientID),
data = toenail, family = binomial(link = "probit")))
m <- ctm(as.basis(~ outcome, data = toenail),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail,
fixed = c("outcomemoderate or severe" = 0))
t4 <- system.time(toenail_mtram_RI <-
mtram(toenail_probit, ~ (1 | patientID),
data = toenail, Hessian = TRUE))
t5 <- system.time(toenail_glmer_RS <-
glmer(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit")))
t6 <- system.time(toenail_glmmTMB_RS_1 <-
glmmTMB(outcome ~ treatment * time + (1 + time | patientID),
data = toenail_gr1, family = binomial(link = "probit")))
m <- ctm(as.basis(~ outcome, data = toenail_gr1),
shifting = ~ treatment * time,
data = toenail, todistr = "Normal", negative = TRUE)
toenail_probit <- mlt(m, data = toenail_gr1,
fixed = c("outcomemoderate or severe" = 0))
t7 <- system.time(toenail_mtram_RS <-
mtram(toenail_probit, ~ (1 + time | patientID),
data = toenail_gr1))
## ----output, echo = FALSE------------------------------------------------
tn_RI_glmer_L <- c(fixef(toenail_glmer_RI_1), toenail_glmer_RI_1@theta, 0, 0)
tn_RI_glmer_A <- c(fixef(toenail_glmer_RI_2), toenail_glmer_RI_2@theta, 0, 0)
tn_RI_glmmTMB <- c(fixef(toenail_glmmTMB_RI_3)$cond, sqrt(VarCorr(toenail_glmmTMB_RI_3)$cond$patientID), 0, 0)
tn_RI_mlt <- c(coef(toenail_mtram_RI), 0, 0)
tn_RS_glmer <- c(fixef(toenail_glmer_RS), toenail_glmer_RS@theta)
tn_RS_glmmTMB <- c(fixef(toenail_glmer_RS), chol(VarCorr(toenail_glmmTMB_RS_1)$cond$patientID)[c(1,3, 4)])
tn_RS_mlt <- coef(toenail_mtram_RS)
tn <- cbind(tn_RI_glmer_L, tn_RI_glmer_A , tn_RI_glmmTMB, tn_RI_mlt ,
tn_RS_glmer, tn_RS_glmmTMB, tn_RS_mlt)
logLik(toenail_glmer_RI_1)
logLik(toenail_glmer_RI_2)
logLik(toenail_glmmTMB_RI_3)
logLik(toenail_mtram_RI)
logLik(toenail_glmer_RS)
logLik(toenail_glmmTMB_RS_1)
logLik(toenail_mtram_RS)
ll <- c(
### logLik of transformation model for glmer (Laplace) parameters
logLik(toenail_mtram_RI, tn_RI_glmer_L[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model for glmer (AGQ) parameters
logLik(toenail_mtram_RI, tn_RI_glmer_A[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model for glmmTMB (Laplace) parameters
logLik(toenail_mtram_RI, tn_RI_glmmTMB[1:5] * c(-1, 1, 1, 1, 1)),
### logLik of transformation model
logLik(toenail_mtram_RI),
### logLik of transformation model for glmer (Laplace) parameters
logLik(toenail_mtram_RS, tn_RS_glmer * c(-1, rep(1, 6))),
### logLik of transformation model for glmmTMB (Laplace) parameters
logLik(toenail_mtram_RS, tn_RS_glmmTMB * c(-1, rep(1, 6))),
### logLik of transformation model
logLik(toenail_mtram_RS))
tm <- c(t1["user.self"],
t2["user.self"],
t3["user.self"],
t4["user.self"],
t5["user.self"],
t6["user.self"],
t7["user.self"])
tm <- formatC(tm, format = "f", digits = 2, width = 5)
tn <- formatC(tn, format = "f", digits = 2, width = 5)
ll <- formatC(ll, format = "f", digits = 2, width = 6)
tn <- cbind(c("$\\alpha$", "$\\eshiftparm_1$", "$\\eshiftparm_2$", "$\\eshiftparm_3$", "$\\gamma_1$", "$\\gamma_2$", "$\\gamma_3$"), tn)
ret <- c("
\\begin{tabular}{lrrrr|rrr} \\\\ \\hline
& \\multicolumn{4}{c|}{RI} & \\multicolumn{3}{c}{RI + RS} \\\\
& \\texttt{glmer} & \\texttt{glmer} & \\texttt{glmmTMB} & & \\texttt{glmer} & \\texttt{glmmTMB} & \\\\
& L & AGQ & L & (7) & L & L & (7) \\\\ \\hline")
ret <- c(ret, apply(tn, 1, function(x) c(paste(x, collapse = " & "), "\\\\")))
ret <- c(ret, "\\hline")
ret <- c(ret,
paste("LogLik &", paste(ll, collapse = "&"), "\\\\ "),
paste("Time (sec) &", paste(tm, collapse = "&"), "\\\\ \\hline"),
"\\end{tabular}")
## ----table, echo = FALSE, results = "asis"------------------------------------
cat(ret, sep = "\n")
## ----mtram-neck_plot, echo = FALSE, fig.height = 4, fig.width = 7-------------
data("neck_pain", package = "ordinalCont")
pain_df <- neck_pain
idt <- xtabs(~ id, data = pain_df)
miss <- names(idt)[idt < 3]
pain_df <- subset(pain_df, !id %in% miss)
pain_df$id <- factor(pain_df$id)
levels(pain_df$laser) <- c("Placebo", "Active")
levels(pain_df$time) <- c("Baseline", "7 weeks", "12 weeks")
pain <- rbind(subset(pain_df, laser == levels(pain_df$laser)[1]),
subset(pain_df, laser == levels(pain_df$laser)[2]))
p1 <- xyplot(vas ~ time | laser, data = pain,
groups = id, type = "l", col = rgb(.1, .1, .1, .1),
lwd = 2, layout = c(2, 1),
ylab = "Neck pain (on visual analog scale)", xlab = "Examinations")
plot(p1)
## ----mtram-ordinalCont--------------------------------------------------------
library("ordinalCont")
## ----mtram-neck_ocm, cache = FALSE, results = "hide"--------------------------
neck_ocm <- ocm(vas ~ laser * time + (1 | id), data = pain_df,
scale = c(0, 1))
## ----mtram-neck_ocm_summary---------------------------------------------------
summary(neck_ocm)
## ----mtram-neck_ocm_ci--------------------------------------------------------
exp(cbind(coef(neck_ocm)[2:6], confint(neck_ocm)[2:6,]))
## ----tramME-neck--------------------------------------------------------------
library("tramME")
neck_ColrME <- ColrME(vas ~ laser * time + (1 | id), data = pain_df,
bounds = c(0, 1), support = c(0, 1))
## ----tramME-neck_ci-----------------------------------------------------------
exp(coef(neck_ColrME))
## ----mtram-neck_Colr, cache = FALSE-------------------------------------------
neck_Colr <- Colr(vas ~ laser * time, data = pain_df,
bounds = c(0, 1), support = c(0, 1),
extrapolate = TRUE)
neck_Colrmer <- mtram(neck_Colr, ~ (1 | id), data = pain_df,
Hessian = TRUE)
## ----mtram-neck_Colr_distr, echo = FALSE, fig.height = 4, fig.width = 7-------
nd <- expand.grid(laser = unique(pain_df$laser),
time = unique(pain_df$time))
tmp <- as.mlt(neck_Colr)
coef(tmp)[] <- coef(neck_Colrmer)[1:12]
q <- 1:99/100
nd2 <- expand.grid(vas = q, laser = unique(pain_df$laser),
time = unique(pain_df$time))
# tmp <- as.mlt(neck_Colr)
sd <- sqrt(coef(neck_Colrmer)[13]^2 + 1)
prb <- predict(tmp, newdata = nd, type = "distribution", q = q)
nd2$prob <- c(pnorm(qnorm(prb) / sd))
p2 <- xyplot(prob ~ vas | time, data = nd2, groups = laser, type = "l",
col = col,
layout = c(3, 1),
xlab = "Neck pain (on visual analog scale)",
ylab = "Marginal distribution",
par.settings = simpleTheme(col=col),
auto.key = list(text = levels(nd2$laser),
points = FALSE, lines = TRUE))
plot(p2)
## M1
# neck_Colrmer <- mtram(neck_Colr, ~ (1 | id), data = pain_df,
# Hessian = TRUE, standardise = FALSE)
# logLik(neck_Colrmer)
#
# nd <- expand.grid(laser = unique(pain_df$laser),
# time = unique(pain_df$time))
# q <- 1:99/100
# nd2 <- expand.grid(vas = q, laser = unique(pain_df$laser),
# time = unique(pain_df$time))
# tmp <- as.mlt(neck_Colr)
# coef(tmp)[] <- coef(neck_Colrmer)[1:12]
# sd <- sqrt(coef(neck_Colrmer)[13]^2 + 1)
# prb <- predict(tmp, newdata = nd, type = "distribution", q = q)
# nd2$prob <- c(pnorm(qnorm(prb) / sd))
# p2 <- xyplot(prob ~ vas | time, data = nd2, groups = laser, type = "l",
# col = col, ylim = c(-0.05, 1.05),
# layout = c(3, 1),
# xlab = "Neck pain (on visual analog scale)",
# ylab = "Marginal distribution",
# par.settings = simpleTheme(col=col),
# auto.key = list(text = levels(nd2$laser),
# points = FALSE, lines = TRUE))
# plot(p2)
## ----mtram-neck_Colr-CI, echo = TRUE, eval=TRUE-------------------------------
S <- solve(neck_Colrmer$Hessian)
rbeta <- rmvnorm(10000, mean = coef(neck_Colrmer), sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
t(apply(rbeta[, 8:12], 2, function(x) {
quantile(exp(x),prob = c(.025, .5, .975))}))
## ----mtram-neck_Colr-std_beta-------------------------------------------------
beta <- coef(neck_Colrmer)[8:12]
alpha <- coef(neck_Colrmer)[13]
(std_beta <- cbind(beta / sqrt(1 + alpha^2)))
## ----mtram-neck_Colr-ctr_mat--------------------------------------------------
ctr_mat <- matrix(c(1, 0, 0, 0, 0,
1, 0, 0, 1, 0,
1, 0, 0, 0, 1), nrow = 3, byrow = TRUE)
ctr_mat %*% std_beta
## ----mtram-neck_PImanual, eval=TRUE-------------------------------------------
(ci_emp <- t(apply(ctr_mat %*% t(rbeta[, 8:12]), 1, function(x) {
quantile(x, prob = c(.025, .5, .975))})))
PI(-ci_emp, link = "logistic")
## ----neck_Colr-PI, echo = TRUE------------------------------------------------
nd <- expand.grid(time = unique(pain_df$time),
laser = unique(pain_df$laser))
neck_Colr_marg <- neck_Colr
neck_Colr_marg$coef <- coef(neck_Colrmer)[1:12] /
sqrt(coef(neck_Colrmer)[13]^2 + 1)
(neck_Colr_PI <- PI(neck_Colr_marg, newdata = nd[1:3, ],
reference = nd[4:6, ],
one2one = TRUE, conf.level = .95))[1:3, 1:3]
## ----mtram-CAO, echo = FALSE--------------------------------------------------
dir <- system.file("rda", package = "TH.data")
load(file.path(dir, "Primary_endpoint_data.rda"))
## ----mtram-CAO-plot, cache = FALSE, echo = FALSE------------------------------
ra <- sort(unique(CAOsurv$randarm))
st <- sort(unique(CAOsurv$strat_t))
sn <- sort(unique(CAOsurv$strat_n))
su <- c(1, 1700)
add <- c(0, max(CAOsurv$iDFS[, "time2"]) - su[2])
ylim <- c(-.05, 1.05)
tmp <- as.mlt(Coxph(iDFS | 0 + strat_n:strat_t:randarm ~ 1, data = CAOsurv,
support = su, add = add, log_first = TRUE))
nd <- expand.grid(strat_n = sn, strat_t = st, randarm = ra)
q <- mkgrid(tmp, 100)[[1]]
surv <- predict(tmp, newdata = nd, type = "survivor", q = q)
nd <- nd[rep(1:nrow(nd), each = nrow(surv)),]
nd$time <- q
nd$surv <- c(surv)
xyplot(surv ~ time | strat_t + strat_n, data = nd, groups = randarm,
type = "l", ylim = c(0, 1), col = col, ylab = "Probability",
xlab = "Time (in days)",
par.settings = simpleTheme(col=col),
auto.key = list(text = levels(nd$randarm),
points = FALSE, lines = TRUE))
## ----mtram-CAO_DFS------------------------------------------------------------
### convert "exact" event dates to interval-censoring (+/- one day)
tmp <- CAOsurv$iDFS
exact <- tmp[,3] == 1
tmp[exact,2] <- tmp[exact,1] + 2
tmp[exact,1] <- pmax(tmp[exact,1] - 2, 0)
tmp[exact,3] <- 3
CAOsurv$iDFS2 <- tmp
## ----mtram-CAO_SR, cache = TRUE-----------------------------------------------
CAO_SR <- Survreg(iDFS2 ~ randarm, data = CAOsurv)
CAO_SR_mtram <- mtram(CAO_SR, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_SR_mtram)
(cf <- coef(CAO_SR_mtram))
(OR <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
## ----mtram-CAO-CI-------------------------------------------------------------
S <- solve(CAO_SR_mtram$Hessian)
# sqrt(diag(S))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_mtram),
sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[, -ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_SR_2, cache = TRUE---------------------------------------------
CAO_SR_2 <- Survreg(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv)
CAO_SR_2_mtram <- mtram(CAO_SR_2, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_SR_2_mtram)
(cf <- coef(CAO_SR_2_mtram))
(OR_2 <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
## ----mtram-CAO-CI-2, echo = FALSE---------------------------------------------
H <- CAO_SR_2_mtram$Hessian
S <- solve(H + .01 * diag(nrow(H)))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_2_mtram),
sigma = S)
s <- rbeta[, ncol(rbeta)]
rbeta <- rbeta[, -ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_Cox_2, cache = TRUE--------------------------------------------
CAO_Cox_2 <- Coxph(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv,
support = c(1, 1700), log_first = TRUE, order = 4)
logLik(CAO_Cox_2)
CAO_Cox_2_mtram <- mtram(CAO_Cox_2, ~ (1 | Block), data = CAOsurv,
Hessian = TRUE)
logLik(CAO_Cox_2_mtram)
coef(CAO_Cox_2_mtram)
## ----mtram-CAO-CI-3, echo = FALSE---------------------------------------------
H <- CAO_Cox_2_mtram$Hessian
S <- solve(H + .1 * diag(nrow(H)))
rbeta <- rmvnorm(10000, mean = coef(CAO_Cox_2_mtram),
sigma = S)
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
## ----mtram-CAO_PI-------------------------------------------------------------
nd <- CAOsurv[1:2, ]
tmp <- CAO_Cox_2
tmp$coef <- coef(CAO_Cox_2_mtram)[-22] / sqrt(coef(CAO_Cox_2_mtram)[22]^2 + 1)
(CAO_Cox_PI <- PI(tmp, newdata = nd[2, ], reference = nd[1, ],
one2one = TRUE, conf.level = .95))[1, ]
## ----mtram-CAO_PI_man---------------------------------------------------------
ci_man <- quantile(-rbeta[, ncol(rbeta)], prob = c(.025, .5, .975))
(CAO_Cox_PIm <- PI(ci_man, link = "minimum extreme value"))
## ----tramME-CAO_SR, cache = TRUE----------------------------------------------
CAO_Cox_2_tramME <- CoxphME(iDFS2 | 0 + strat_n:strat_t ~ randarm + (1 | Block),
data = CAOsurv, log_first = TRUE)
## ----tramME-CAO_SR-hr, cache = TRUE-------------------------------------------
exp(coef(CAO_Cox_2_tramME))
exp(confint(CAO_Cox_2_tramME, parm = "randarm5-FU + Oxaliplatin",
estimate = TRUE))
## ----echo=FALSE, eval=FALSE---------------------------------------------------
## sqrt(VarCorr(CAO_Cox_2_tramME)$Block$var)
## coef(CAO_Cox_2_mtram)["gamma1"]
## ----mtram-CAO-coxme, echo = FALSE, eval = FALSE------------------------------
## library("coxme")
## m <- coxme(DFS ~ randarm + (1 | Block), data = CAOsurv)
## summary(m)
## sd <- sqrt(diag(vcov(m)))
## exp(coef(m) + c(-1, 0, 1) * qnorm(.975) * sd)
## ----sim, eval = FALSE--------------------------------------------------------
## source(system.file("simulations", "mtram_sim.R", package = "tram"), echo = TRUE)
## ----mtram-sessionInfo, echo = FALSE, results = "hide"------------------------
sessionInfo()
## ----mtram-funs, echo = FALSE, results = "hide"-------------------------------
if (file.exists("packages.bib")) file.remove("packages.bib")
pkgversion <- function(pkg) {
pkgbib(pkg)
packageDescription(pkg)$Version
}
pkgbib <- function(pkg) {
x <- citation(package = pkg, auto = TRUE)[[1]]
b <- toBibtex(x)
b <- gsub("Buehlmann", "B{\\\\\"u}hlmann", b)
b[1] <- paste("@Manual{pkg:", pkg, ",", sep = "")
if (is.na(b["url"])) {
b[length(b)] <- paste(" URL = {http://CRAN.R-project.org/package=",
pkg, "}", sep = "")
b <- c(b, "}")
}
cat(b, sep = "\n", file = "packages.bib", append = TRUE)
}
pkg <- function(pkg) {
vrs <- try(pkgversion(pkg))
if (inherits(vrs, "try-error")) return(NA)
paste("\\\\pkg{", pkg, "} \\\\citep[version~",
vrs, ",][]{pkg:", pkg, "}", sep = "")
}
pkg("mlt")
pkg("tram")
pkg("SparseGrid")
cat(c("@Manual{vign:mlt.docreg,",
" title = {Most Likely Transformations: The mlt Package},",
" author = {Torsten Hothorn},",
paste(" year = ", substr(packageDescription("mlt.docreg")$Date, 1, 4), ",", sep = ""),
paste(" note = {R package vignette version ", packageDescription("mlt.docreg")$Version, "},", sep = ""),
" url = {https://CRAN.R-project.org/package=mlt.docreg},",
"}"), file = "packages.bib", append = TRUE, sep = "\n")
Try the tram package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.