Nothing
## ----nami-pkgs, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
pkgs <- c("tram", "TH.data", "multcomp", "survival", "Stat2Data", "tramME")
pkgs <- sapply(pkgs, require, character.only = TRUE)
## ----nami-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()
}
if (!interactive() && .Platform$OS.type != "unix")
{
cat("Vignette only compiled under Unix alikes.")
knitr::knit_exit()
}
## ----nami-setup, echo = FALSE, results = "hide", message = FALSE, warning = FALSE----
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)
frmt1 <- function(x, digits = 1) {
if (!is.numeric(x)) return(x)
formatC(round(x, digits = digits), digits = digits, format = "f")
}
frmt3 <- function(x)
frmt1(x, digits = 3)
frmtci <- function(x, digits = 3) {
if (!is.numeric(x)) return(x)
if (length(x) != 2) stop("not a confidence interval")
return(paste("(", frmt1(x[1], digits = digits),
",", frmt1(x[2], digits = digits), ")"))
}
## discrete nonparanormal likelihood relies on Monte Carlo, set seed
set.seed(221224L)
## ----nami-demo, eval = FALSE---------------------------------------------
# install.packages("tram")
# demo("NAMI", package = "tram")
## ----immun, echo = FALSE, message = FALSE, results = "hide"--------------
immun <- structure(list(y = c(21.699999999999999, 23.899999999999999,
22.699999999999999, 23.399999999999999, 26.800000000000001, 24.800000000000001,
23.399999999999999, 25.100000000000001, 24.100000000000001, 23.300000000000001,
25.399999999999999, 25.100000000000001, 23.800000000000001, 23.100000000000001,
24, 24.199999999999999, 27.399999999999999, 23.300000000000001,
22.600000000000001, 23.300000000000001, 24.800000000000001, 23.899999999999999,
22.199999999999999, 21.399999999999999, 22.800000000000001, 22.300000000000001,
22.399999999999999, 30.399999999999999, 30.600000000000001, 21.699999999999999,
24.800000000000001, 25.600000000000001, 21.399999999999999, 24.300000000000001,
23.5, 25.800000000000001, 21.600000000000001, 22.899999999999999,
23.800000000000001, 22.600000000000001, 24.199999999999999, 24.300000000000001,
25.699999999999999, 23.199999999999999, 24.600000000000001, 24.5,
22.699999999999999, 26.300000000000001, 27.199999999999999, 27.100000000000001,
22.699999999999999, 24.600000000000001, 23, 23.199999999999999,
23.899999999999999, 23, 20.800000000000001, 23.399999999999999,
24.300000000000001, 24.399999999999999, 22.600000000000001, 22.100000000000001,
22.199999999999999, 24.100000000000001, 28.100000000000001, 23.399999999999999,
26.800000000000001, 24, 25.899999999999999, 24.699999999999999,
24.100000000000001, 26.899999999999999, 23.899999999999999, 24.399999999999999,
25.199999999999999, 22.899999999999999, 25.699999999999999, 24.300000000000001,
25.199999999999999, 24.100000000000001), w = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("0",
"100"), class = "factor"), x = c(19, 22, 19.899999999999999,
20.800000000000001, 22.899999999999999, 21, 21.800000000000001,
22.600000000000001, 21.5, 22.600000000000001, 21.5, 22, 20.699999999999999,
21.699999999999999, 20.800000000000001, 22.300000000000001, 21.800000000000001,
19.800000000000001, 20.699999999999999, 20.600000000000001, 19.899999999999999,
20.199999999999999, 20.199999999999999, 19, 19.600000000000001,
18.699999999999999, 18.699999999999999, 21.699999999999999, 21.300000000000001,
19.199999999999999, 21, 21.100000000000001, 18.699999999999999,
20.300000000000001, 21.600000000000001, 21.899999999999999, 18.100000000000001,
19.199999999999999, 20, 18.899999999999999, 20.399999999999999,
20.5, 21.300000000000001, 19.199999999999999, 19.399999999999999,
20.199999999999999, 18.199999999999999, 21, 21.5, 22.100000000000001,
19.600000000000001, 21.399999999999999, 19.199999999999999, 21,
20.899999999999999, 19.600000000000001, 19, 21.300000000000001,
22.199999999999999, 21.699999999999999, 20.800000000000001, 20.600000000000001,
19.199999999999999, 20.800000000000001, 23.600000000000001, 21.199999999999999,
22.600000000000001, 21.100000000000001, 23.300000000000001, 22.300000000000001,
21.600000000000001, 22.600000000000001, 21.199999999999999, 23.899999999999999,
23.800000000000001, 20.699999999999999, 22.199999999999999, 22.5,
22.199999999999999, 22.699999999999999)), row.names = c("1.5",
"2.10", "3.15", "4.20", "5.25", "6.30", "7.35", "8.40", "33.165",
"34.170", "35.175", "36.180", "37.185", "38.190", "39.195", "40.200",
"57.285", "58.290", "59.295", "60.300", "61.305", "62.310", "63.315",
"64.320", "89.445", "90.450", "91.455", "92.460", "93.465", "94.470",
"95.475", "96.480", "121.605", "122.610", "123.615", "124.620",
"125.625", "126.630", "127.635", "128.640", "153.765", "154.770",
"155.775", "156.780", "157.785", "158.790", "159.795", "160.800",
"177.885", "178.890", "179.895", "180.900", "181.905", "182.910",
"183.915", "184.920", "209.1045", "210.1050", "211.1055", "212.1060",
"213.1065", "214.1070", "215.1075", "216.1080", "233.1165", "234.1170",
"235.1175", "236.1180", "237.1185", "238.1190", "239.1195", "240.1200",
"265.1325", "266.1330", "267.1335", "268.1340", "269.1345", "270.1350",
"271.1355", "272.1360"), class = "data.frame")
## ----marginal_outcome, echo = TRUE, message = FALSE, results = "hide"----
m0 <- Lm(y ~ w, data = immun)
## ----cf0_cont, echo = TRUE, message = FALSE------------------------------
coef(m0) ### marginal Cohen's d
sqrt(vcov(m0)) ### observed
sqrt(2/nrow(immun) * (coef(m0)^2/4 + 2)) ### expected
confint(m0) ### Wald
## ----ami_cont, echo = TRUE, message = FALSE, results = "hide"------------
m1 <- BoxCox(x ~ 1, data = immun)
m <- Mmlt(m0, m1, formula = ~ 1, data = immun)
## ----cfseci_cont, echo = TRUE, message = FALSE---------------------------
(cf1 <- coef(m)["y.w100"]) ### marginal adjusted Cohen's d
sqrt(diag(vcov(m))["y.w100"]) ### observed
## ----setheo_cont, echo = TRUE, message = FALSE---------------------------
lambda <- c(unclass(coef(m, type = "Lambdapar")))
sqrt(2/nrow(immun) * ((1 + lambda^2)*cf1^2 + 8)/(4*(1 + lambda^2))) ### expected
## ----ci_cont, echo = TRUE, message = FALSE-------------------------------
confint(m)["y.w100",]
## ----corr_cont, echo = TRUE, message = FALSE-----------------------------
coef(m, type = "Corr")
## ----r2_cont, echo = TRUE, message = FALSE-------------------------------
Omega <- as.array(coef(m, type = "Lambda"))[,,1]
1 - Omega[nrow(Omega), ncol(Omega)]^(-2)
## ----cfsec_cont, echo = TRUE, message = FALSE----------------------------
mad <- BoxCoxME(y ~ w + s(x), data = immun)
## ----hx_smooth, echo = TRUE, fig.width = 4, fig.height = 3.5, out.width='.6\\linewidth'----
plot(smooth_terms(mad))
## ----hy_cont, echo = TRUE, fig.width = 4, fig.height = 3.5, out.width='.6\\linewidth'----
plot(mad)
## ----boxcox_cont, echo = TRUE, message = FALSE---------------------------
m0 <- BoxCox(y ~ w, data = immun)
m <- Mmlt(m0, m1, formula = ~ 1, data = immun)
## ----r2_boxcox, echo = TRUE, message = FALSE-----------------------------
Omega <- as.array(coef(m, type = "Lambda"))[,,1]
1 - Omega[nrow(Omega), ncol(Omega)]^(-2)
## ----coef_boxcox_cont, echo = TRUE, message = FALSE----------------------
(cf1 <- coef(m)["y.w100"])
confint(m)["y.w100",]
## ----trt, echo = TRUE, message = FALSE-----------------------------------
dnorm(cf1 / sqrt(2))
## ----CAOdata, echo = FALSE, message = FALSE------------------------------
load(system.file("rda", "Primary_endpoint_data.rda", package = "TH.data"))
## ----CAOno, echo = FALSE, message = FALSE--------------------------------
rt <- table(CAOsurv$randarm)
## ----CAOoutcome, echo = TRUE, message = FALSE----------------------------
CAOsurv$ypT0ypN0 <- factor(CAOsurv$path_stad == "ypT0ypN0")
## ----CAO-glm, echo = TRUE, message = FALSE-------------------------------
mg_w <- glm(ypT0ypN0 ~ randarm,
data = CAOsurv, family = binomial())
exp(coef(mg_w)["randarm5-FU + Oxaliplatin"])
exp(confint(glht(mg_w), calpha = univariate_calpha())$confint[2,-1])
## ----CAO-marg, echo = TRUE, message = FALSE------------------------------
mpCR <- Polr(ypT0ypN0 ~ randarm, data = CAOsurv, na.action = na.pass,
method = "logistic")
exp(coef(mpCR)["randarm5-FU + Oxaliplatin"])
## ----CAO-margci, echo = TRUE, message = FALSE----------------------------
exp(confint(glht(mpCR, coef. = function(...) coef(..., fixed = FALSE)),
calpha = univariate_calpha())$confint)
## ----CAO-covariates, echo = TRUE-----------------------------------------
mage <- BoxCox(age ~ 1, data = CAOsurv)
msex <- Polr(geschlecht ~ 1, data = CAOsurv, method = "probit")
CAOsurv$ecog_o <- as.ordered(CAOsurv$ecog_b)
mecog <- Polr(ecog_o ~ 1, data = CAOsurv, na.action = na.pass,
method = "probit")
mentf <- Polr(bentf ~ 1, data = CAOsurv, na.action = na.pass,
method = "probit")
mT <- Polr(strat_t ~ 1, data = CAOsurv, method = "probit")
mN <- Polr(strat_n ~ 1, data = CAOsurv, method = "probit")
## ----CAO-Mmlt, echo = TRUE-----------------------------------------------
m <- Mmlt(mT, mN, mentf, mage, msex, mecog, mpCR,
data = CAOsurv, args = list(type = "ghalton", M = 250),
optim = mltoptim(hessian = TRUE)["constrOptim"])
prm <- "ypT0ypN0.randarm5-FU + Oxaliplatin"
exp(coef(m)[prm])
## ----CAO-mmltci, echo = TRUE---------------------------------------------
ci <- confint(glht(m, coef. = function(...) coef(..., fixed = FALSE)),
calpha = univariate_calpha())$confint
exp(ci[prm,-1])
## ----CAO-corrhide, echo = TRUE-------------------------------------------
mr <- as.array(coef(m, type = "Cor"))["ypT0ypN0",,1]
i <- which.max(abs(mr[-length(mr)]))
(ni <- names(mr)[i])
(mr <- mr[i])
## ----CAO-r2, echo = TRUE, message = FALSE--------------------------------
Omega <- as.array(coef(m, type = "Lambda"))[,,1]
1 - Omega[nrow(Omega), ncol(Omega)]^(-2)
## ----CAO-se, echo = TRUE, message = FALSE--------------------------------
sqrt(vcov(mpCR))
sqrt(vcov(m)[prm, prm])
## ----flies, echo = TRUE, results = "hide"--------------------------------
data("FruitFlies", package = "Stat2Data")
## ----flies_subset, echo = TRUE, message = FALSE--------------------------
flies <- FruitFlies
flies <- flies[flies$Treatment %in% c("8 virgin", "8 pregnant"),]
flies$Treatment <- flies$Treatment[, drop = TRUE]
flies$Longevity <- as.double(flies$Longevity)
flies$survival <- Surv(flies$Longevity)
## ----flies_marg, echo = TRUE---------------------------------------------
coxph_w <- Coxph(survival ~ Treatment, data = flies)
coef(coxph_w) ### log-hazard ratio
confint(coxph_w) ### Wald
## ----flies_mmlt, echo = TRUE, message = FALSE----------------------------
xmod <- BoxCox(Thorax ~ 1, data = flies)
m <- Mmlt(xmod, coxph_w, data = flies, formula = ~ 1)
(cf1 <- coef(m)["survival.Treatment8 virgin"])
(ci1 <- confint(m)["survival.Treatment8 virgin",])
## ----flies_r2, echo = TRUE, message = FALSE------------------------------
Omega <- as.array(coef(m, type = "Lambda"))[,,1]
1 - Omega[nrow(Omega), ncol(Omega)]^(-2)
## ----trans_y, echo = TRUE, fig.width = 4, fig.height = 3.5, out.width='.6\\linewidth'----
q <- 0:100
cols <- c("grey20", "grey70")
### nonparametric
plot(q, log(-log(1 - ecdf(subset(flies, Treatment == "8 pregnant")$survival[,1])(q))),
main = "", xlab = "Time", type = "S", lwd = 1,
ylab = "cloglog(Probability)")
lines(q, log(-log(1 - ecdf(subset(flies, Treatment == "8 virgin")$survival[,1])(q))),
type = "S", lty = 2, lwd = 1)
legend("bottomright", lty = c(1, 2),
legend = levels(flies$Treatment), bty = "n")
### model-based
nd <- expand.grid(survival = q, Treatment = sort(unique(flies$Treatment)))
nd$h <- predict(as.mlt(coxph_w), newdata = nd, type = "trafo")
fm <- nd$Treatment == "8 virgin"
lines(nd$survival[fm], nd$h[fm], lty = 2)
fm <- nd$Treatment == "8 pregnant"
lines(nd$survival[fm], nd$h[fm], lty = 1)
## ----copula, echo = TRUE, fig.with = 4, fig.height = 4, out.width='.6\\linewidth'----
m <- CoxphME(survival ~ s(Thorax, k = 5), data = flies,
subset = Treatment == levels(Treatment)[1])
plot(smooth_terms(m))
m <- CoxphME(survival ~ s(Thorax, k = 5), data = flies,
subset = Treatment == levels(Treatment)[2])
plot(smooth_terms(m), add = TRUE, lty = 2)
## ----nami-pkgs-funs, echo = FALSE, results = "hide"----------------------
if (file.exists("packages.bib")) file.remove("packages.bib")
## sentence style for titles
toLower <- function(text) {
parts <- strsplit(unname(text), split = ":")[[1]]
w1 <- paste0(parts[1], ":")
p2 <- strsplit(parts[2], split = "}")[[1]]
p3 <- strsplit(p2[1], " ")
p4 <- strsplit(p2[1], " ")[[1]]
w2 <- p4[2]
w3 <- paste(p4[3:length(p4)], collapse = " ")
w3 <- tolower(w3)
paste(w1, w2, w3, "},")
}
sentence_style <- FALSE
## R
x <- citation()[[1]]
b <- toBibtex(x)
b <- gsub("R:", paste0("\\\\proglang{R}:"), b)
b <- gsub("R ", paste0("\\\\proglang{R} "), b)
if (sentence_style) b["title"] <- toLower(b["title"])
b[1] <- "@Manual{R,"
cat(b, sep = "\n", file = "packages.bib", append = TRUE)
pkgv <- function(pkg) packageVersion(pkg)
pkgbib <- function(pkg) {
x <- citation(package = pkg, auto = TRUE)[[1]]
b <- toBibtex(x)
b <- gsub("packaging by", "", b)
b <- gsub("with contributions from", "", b)
b <- gsub("Gruen", "Gr{\\\\\"u}n", 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 <- b[names(b) != "url"]
if (is.na(b["doi"])) {
b[length(b)] <- paste(" DOI = {10.32614/CRAN.package.",
pkg, "}", sep = "")
b <- c(b, "}")
}
b["note"] <- gsub("R package", "\\\\proglang{R} package", b["note"])
cat(b, sep = "\n", file = "packages.bib", append = TRUE)
}
pkg <- function(pkg)
paste("\\pkg{", pkg, "} \\citep[version~",
pkgv(pkg), ",][]{pkg:", pkg, "}", sep = "")
pkgs <- c("tram")
sapply(pkgs, pkgbib)
out <- sapply(pkgs, pkg)
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.