%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Applications of Nonparanormal Adjusted Marginal Inference} %\VignetteDepends{tram, TH.data, multcomp, survival, Stat2Data, tramME} \documentclass[article,nojss,shortnames]{jss} %% packages \usepackage{thumbpdf} \usepackage{amsfonts,amstext,amsmath,amssymb,amsthm} \usepackage{accents} \usepackage{color} \usepackage{rotating} \usepackage{verbatim} \usepackage[utf8]{inputenc} \usepackage{xspace} %% need no \usepackage{Sweave.sty} %%\usepackage[nolists]{endfloat} \newcommand{\cmd}[1]{\texttt{#1()}} \usepackage{tikz} \usetikzlibrary{shapes,arrows,chains} \usepackage{verbatim} <>= pkgs <- c("tram", "TH.data", "multcomp", "survival", "Stat2Data", "tramME") pkgs <- sapply(pkgs, require, character.only = TRUE) @ \newcommand{\TODO}[1]{{\color{red} #1}} \newcommand\Torsten[1]{{\color{blue}Torsten: ``#1''}} \newcommand{\THcite}[2]{\citeauthor{#2} (\citeyear{#2})} \newcommand\norm[1]{\left\lVert#1\right\rVert} \newcommand{\etc}{\textit{etc.}} \usepackage{booktabs} \newcommand{\NAMI}{nonparanormal adjusted marginal inference\xspace} \newcommand{\expit}{\text{expit}} \input{defs} \renewcommand{\thefootnote}{} %% code commands \newcommand{\Rclass}[1]{`\code{#1}'} %% JSS \author{Susanne Dandl \\ Universit\"at Z\"urich \And Torsten Hothorn \\ Universit\"at Z\"urich} \Plainauthor{Dandl and Hothorn} \title{Some Applications of Nonparanormal Adjusted Marginal Inference} \Plaintitle{Nonparanormal Adjusted Marginal Inference} \Shorttitle{Nonparanormal Adjusted Marginal Inference} \Abstract{ Covariate adjustment is recommended to improve precision of estimated treatment effects, even in randomised trials. For non-collapsible effect measures (such as Cohen's $d$, odds and hazard ratios), conditioning on covariates can alter the effect interpretation such that effects are not comparable given different adjustment sets. A novel nonparanormal adjusted marginal inference approach (NAMI) can be formulated for randomised controlled trials to derive marginal effects from the estimated joint distribution of covariates and outcomes. NAMI preserves the interpretation of treatment effects regardless of covariates used, while enhancing precision. We present applications of NAMI with continuous, binary and right-censored survival outcomes. Data is analysed by a proof-of-concept implementation of NAMI available in the \pkg{tram} add-on package to the \proglang{R} system for statistical computing. } \Keywords{marginal effect, noncollapsibility, covariate adjustment, randomised trial, transformation model} \Plainkeywords{marginal effect, noncollapsibility, covariate adjustment, randomised trial, transformation model} \Address{ Susanne Dandl \& Torsten Hothorn\\ Institut f\"ur Epidemiologie, Biostatistik und Pr\"avention \\ Universit\"at Z\"urich \\ Hirschengraben 84, CH-8001 Z\"urich, Switzerland \\ \texttt{Torsten.Hothorn@R-project.org} } \begin{document} <>= year <- substr(packageDescription("tram")$Date, 1, 4) version <- packageDescription("tram")$Version @ \footnote{Please cite this document as: Susanne Dandl and Torsten Hothorn (\Sexpr{year}) Some Applications of Nonparanormal Adjusted Marginal Inference. \textsf{R} package vignette version \Sexpr{version}, URL \url{https://doi.org/10.32614/CRAN.package.tram}.} <>= 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() } @ <>= 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) @ \section{Introduction} \label{sec:introduction} We present three applications from different domains illustrating the potential to estimate treatment effects more precisely with nonparanormal adjusted marginal inference \citep[NAMI,][]{Dandl_2025_nami}: A toxicity study demonstrates equivalence of a continuous outcome by Cohen's $d$, an efficiency trial compares a binary outcome between two arms by means of an odds ratio, and a survival study measures differences in terms of a hazard ratio. We apply the proof-of-concept implementation of \NAMI available in package \pkg{tram}. All results of this document can be reproduced from the \code{NAMI} demo available in the package: <>= install.packages("tram") demo("NAMI", package = "tram") @ \section{Continuous outcome: Immunotoxicity study on Chloramine} \label{subsec:appl-cont} To assess the toxicity of Chloramine, a study on mice was conducted as part of the \cite{toxicology_2000}. Female mice were randomly assigned to two groups of which one group received Chlor\-amine-dosed water and the other not. Repeated measurements of weights were conducted on days $1$, $8$, $15$, $22$, and $29$ in five dose groups ($0$, $2$, $10$, $20$, $100$ mg/kg). <>= 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") @ The following data set is a subset of the data obtained from the \cite{toxicology_2000}. It focuses on the comparison of the highest dose group ($W = 1, N_1 = 40$) with the control group ($W = 0, N_0 = 40$) with respect to the weight at day $29$ (denoted as outcome $\rY$). Information on the weight on day $1$ is also available and is used for adjustment, in the following. The research question of ``\emph{no} effect of Chloramine on weight'' is formulated in terms of the hypotheses $H_0: |\tau| \ge \delta$ and $H_1: |\tau| < \delta$ where Cohen's $d$ is the marginal treatment effect $\tau$. $H_0$ is rejected at $5\%$ when the $95\%$ confidence interval for $\tau$ is completely contained in the equivalence interval $(-\delta, \delta)$. We set $\delta = 0.36$ as recommended by \cite{wellek_equiv_2010} (Table 1.1.), acknowledging this as a potential oversimplification. \subsection{Unadjusted marginal inference} The following fits an unadjusted marginal outcome model featuring Cohen's $d$: <>= m0 <- Lm(y ~ w, data = immun) @ <>= coef(m0) ### marginal Cohen's d sqrt(vcov(m0)) ### observed sqrt(2/nrow(immun) * (coef(m0)^2/4 + 2)) ### expected confint(m0) ### Wald @ The unadjusted estimate of Cohen's $d$ is $\hat{\tau} = \Sexpr{frmt3(coef(m0))}$. The standard error based on the observed Fisher information is $\text{SE}(\hat{\tau}) = \Sexpr{frmt3(sqrt(vcov(m0)))}$. It is identical to the expected Fisher information \citep[see Lemma~1 in][]{Dandl_2025_nami} evaluated at $\hat{\tau}$. The resulting $95\%$ Wald interval obtained is $\Sexpr{frmtci(confint(m0))}$. Because the confidence interval is not fully contained in $(-0.36, 0.36)$, the unadjusted analysis cannot reject $H_0$. \subsection{Nonparanormal adjusted marginal inference} Adjusting for weight at baseline (variable \code{x}) in \NAMI (with transformation function in Bernstein form of order six) is conducted by fitting a marginal model \code{m1} for baseline weight in addition to the marginal model for the outcome \code{m0} and estimating all marginal and Copula parameters simultaneously with \code{Mmlt()}: <>= m1 <- BoxCox(x ~ 1, data = immun) m <- Mmlt(m0, m1, formula = ~ 1, data = immun) @ <>= (cf1 <- coef(m)["y.w100"]) ### marginal adjusted Cohen's d sqrt(diag(vcov(m))["y.w100"]) ### observed @ We obtain $\hat{\tau} = \Sexpr{frmt3(coef(m)["y.w100"])}$ as marginal adjusted Cohen's $d$. The standard error based on the observed Fisher information is $\Sexpr{frmt3(sqrt(diag(vcov(m))["y.w100"]))}$. Given an estimate of the copula parameters in $\mLambda$ \citep[see Sec.~2.2 in][]{Dandl_2025_nami}, we can confirm that the observed Fisher information is equal to the expected Fisher information evaluated at the maximum-likelihood estimates $\hat{\tau}$ and $\hat{\lambda}$: <>= lambda <- c(unclass(coef(m, type = "Lambdapar"))) sqrt(2/nrow(immun) * ((1 + lambda^2)*cf1^2 + 8)/(4*(1 + lambda^2))) ### expected @ The Wald interval $\Sexpr{frmtci(confint(m)["y.w100",])}$ is completely contained in the equivalence interval and thus the absence of an effect of Chloramine on weight can be inferred: <>= confint(m)["y.w100",] @ This finding is in accordance with previous research by \cite{guo_immuno_2011}. The reduction in standard error comes from the high correlation between the outcome $\rY$ and the covariate which can be estimated via <>= coef(m, type = "Corr") @ The coefficient of determination $R^2$ \citep[see Section~2.2 in][]{Dandl_2025_nami} obtained by <>= Omega <- as.array(coef(m, type = "Lambda"))[,,1] 1 - Omega[nrow(Omega), ncol(Omega)]^(-2) @ suggests an improvement of the conditional over the marginal model. \subsection{Model diagnosis} Model diagnosis can be conducted by fitting an additive transformation model for $\rY$ using the \pkg{tramME} package: <>= mad <- BoxCoxME(y ~ w + s(x), data = immun) @ Inspecting the monotonicity of the estimated smooth functions of the baseline covariate, allows to assess the copula fit. The following plot shows the smooth function for baseline weight estimated from the additive transformation model: <>= plot(smooth_terms(mad)) @ Because the functions are monotonically increasing, the copula fit seems to be appropriate. The marginal model \code{m0} for $\rY$ relies on the normality assumption and assumes that a linear transformation is sufficient. Via \code{plot()}, linearity of the transformation function of $\rY$ can be assessed using the additive transformation model: <>= plot(mad) @ The plot suggests that a more flexible model allowing a nonlinear transformation $h$ for $\rY$ might be better suited. This can be done in the \pkg{tram} package using \code{BoxCox()} for fitting a marginal model. The function estimates a nonlinear transformation function for $\rY$ in the form of polynomials in Bernstein form of order $M$ ($M = 6$ is the default): <>= m0 <- BoxCox(y ~ w, data = immun) m <- Mmlt(m0, m1, formula = ~ 1, data = immun) @ The coefficient of determination $R^2$ can be obtained via <>= Omega <- as.array(coef(m, type = "Lambda"))[,,1] 1 - Omega[nrow(Omega), ncol(Omega)]^(-2) @ suggesting an improved model fit over \NAMI with a linear transformation. This generalised version of Cohen's $d$ is then <>= (cf1 <- coef(m)["y.w100"]) confint(m)["y.w100",] @ The obtained treatment effect $\hat{\tau} = \Sexpr{frmt3(coef(m)["y.w100"])}$ with Wald interval $\Sexpr{frmtci(confint(m)["y.w100",])}$ is interpreted as the mean difference on the latent normal scale in this model. Transforming this estimate to the probabilistic index via $\Phi(\hat{\tau}/\sqrt{2})$, provides a more intuitive interpretation in terms of the probability of obtaining a lower outcome in the control than in the treatment group: <>= dnorm(cf1 / sqrt(2)) @ \section{Binary outcome: Efficacy study of new chemotherapy} \cite{roedel_oxaliplatin_2012} conducted an efficacy study of combining standard care with a new therapy for patients with rectal cancer. Data from the completed trial \citep{Roedel_Graeven_Fietkau_2015} are available in the \pkg{TH.data} package. <>= load(system.file("rda", "Primary_endpoint_data.rda", package = "TH.data")) @ <>= rt <- table(CAOsurv$randarm) @ $\Sexpr{rt["5-FU"]}$ patients received the Fluorouracil-based standard of care ($W = 0$), and $\Sexpr{rt["5-FU + Oxaliplatin"]}$ patients received a combination therapy adding Oxaliplatin ($W = 1$). Early effects were assessed based on a binary outcome \code{ypT0ypN0} -- the absence of viable tumour cells in the primary tumour and lymph nodes after surgery: <>= CAOsurv$ypT0ypN0 <- factor(CAOsurv$path_stad == "ypT0ypN0") @ \subsection{Unadjusted marginal inference} \cite{roedel_oxaliplatin_2012} reported an unadjusted odds ratio of $1.4$ with $95\%$ confidence interval $(1.02, 1.92)$ based on a Cochran-Mantel-Haenszel $\chi^2$ test (without continuity correction) stratified for lymph node involvement (positive vs.~negative) and clinical T category (1–3 vs.~4). Similar results can be obtained by fitting a binary logistic regression model: <>= 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]) @ An equivalent model can be obtained with \code{Polr()} using the \pkg{tram} package. The outcome variable \code{ypT0ypN0} contains $\Sexpr{sum(is.na(CAOsurv$ypT0ypN0))}$ missings. For \code{Polr()}, the missing values in the outcome are retained in the data set, but are ignored in the log-likelihood by choosing \code{na.action = na.pass}: <>= mpCR <- Polr(ypT0ypN0 ~ randarm, data = CAOsurv, na.action = na.pass, method = "logistic") exp(coef(mpCR)["randarm5-FU + Oxaliplatin"]) @ <>= exp(confint(glht(mpCR, coef. = function(...) coef(..., fixed = FALSE)), calpha = univariate_calpha())$confint) @ \subsection{Nonparanormal adjusted marginal inference} Information on six potentially prognostic covariates is available: \code{age}, sex (\code{geschlecht}), ECOG performance status (\code{ecog\_o}), distance to the anal verge of the tumour (\code{bentf}) and the two stratum variables lymph node involvement (\code{strat\_n}) and clinical T category (\code{strat\_t}). For \NAMI, we define one marginal model for each covariate: a Box-Cox-type model for continuous variable age, and binary or ordinal probit models for the remaining binary or ordinal covariates. The ECOG performance status and the distance to the anal verge have 14 missing values, respectively. As for the outcome model, all missing values are retained in the data set but ignored in the marginal log-likelihoods (\code{na.action = na.pass}): <>= 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") @ The joint model is fitted with \code{Mmlt}. It is important to note that, unlike the marginal models, this is not a complete cases analysis. If, for example, the outcome is missing for one observation, the corresponding datum is still taken into account. The log-likelihood contribution is obtained by integrating out the missing dimension from the model. Details on the estimation of the mixed continuous-discrete likelihood are given in \cite{hothorn_2024}. Here, we fit the model <>= 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]) @ and compute an adjusted confidence interval for the marginal log-odds ratio <>= ci <- confint(glht(m, coef. = function(...) coef(..., fixed = FALSE)), calpha = univariate_calpha())$confint exp(ci[prm,-1]) @ The corresponding odds ratio $\Sexpr{frmt3(exp(coef(m)[prm]))}$ with Wald interval $\Sexpr{frmtci(exp(ci[prm,-1]))}$ is very close to the initial results obtained by \cite{roedel_oxaliplatin_2012}. The reason is that none of the six variables (including the stratum variables) carries strong prognostic information. The covariate with the largest association with the binary outcome can be identified by ranking the covariates according to their linear correlation after transformation to normality: <>= mr <- as.array(coef(m, type = "Cor"))["ypT0ypN0",,1] i <- which.max(abs(mr[-length(mr)])) (ni <- names(mr)[i]) (mr <- mr[i]) @ \code{ecog\_o} is the variable with the largest association with the outcome but its linear correlation under transformation to normality only has a value of $\hat{\rho}_{J,\text{\Sexpr{toupper(substr(ni, 1, 4))}}} = \Sexpr{frmt3(mr)}$. Consequently, adjusting for covariates did not improve fit nor precision, reflected in a low $R^2$: <>= Omega <- as.array(coef(m, type = "Lambda"))[,,1] 1 - Omega[nrow(Omega), ncol(Omega)]^(-2) @ However, the standard error did not increase despite the adjustment for multiple covariates, as the following comparison to the standard error of the marginal model shows: <>= sqrt(vcov(mpCR)) sqrt(vcov(m)[prm, prm]) @ \section{Survival outcome: Longevity study of male fruit flies} \label{subsec:appl-surv} \cite{partridge_flies_1981} conducted a study on the sexual behavior of fruit flies. The aim was to investigate whether increased sexual activity leads to shorter life spans for male fruit flies. They randomly assigned a total of $125$ male fruit flies into five groups of $25$ flies: males that live alone, males that live with one or eight receptive females, and males that live with one or eight nonreceptive females. The data set \code{FruitFlies} is available in the \pkg{Stat2Data} package: <>= data("FruitFlies", package = "Stat2Data") @ We follow \cite{negassa_flies_2007} and compare only the two groups with eight female flies added that were either all nonreceptive ($W = 0$) or receptive ($W = 1$): <>= 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) @ Of interest is the survival time $\rY$ of male flies in days. Because thorax length of the male flies is strongly associated with longevity, it is used as covariate in this analysis. \subsection{Unadjusted marginal inference} The following code fits a marginal Cox proportional hazards model for time to death, with baseline log-cumulative hazard function in Bernstein form of order six: <>= coxph_w <- Coxph(survival ~ Treatment, data = flies) coef(coxph_w) ### log-hazard ratio confint(coxph_w) ### Wald @ The marginal treatment effect is defined as a log-hazard ratio $\hat{\tau} = \Sexpr{frmt3(coef(coxph_w))}$ with $95\%$ Wald interval $\Sexpr{frmtci(confint(coxph_w))}$. Consequently, the hazard of dying in the sexually active group is around $\exp(\hat{\tau}) = \Sexpr{frmt3(exp(coef(coxph_w)))}$ times higher than for the nonactive group. \subsection{Nonparanormal adjusted marginal inference} For \NAMI, we parameterise the transformation function $\h_1$ in the marginal model of thorax length in Bernstein form with order six. The joint distribution of both variables was expressed by a Gaussian copula model: <>= 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",]) @ The log-hazard ratio of \NAMI is $\Sexpr{frmt3(cf1)}$, with a shorter $95\%$ Wald interval of $\Sexpr{frmtci(ci1)}$. The coefficient of determination indicates that thorax length is highly prognostic: <>= Omega <- as.array(coef(m, type = "Lambda"))[,,1] 1 - Omega[nrow(Omega), ncol(Omega)]^(-2) @ \subsection{Model diagnosis} Model diagnosis is conducted on the level of the marginal model for the outcome and on the level of the Gaussian copula parameterising the joint model. On the level of the marginal model, we compare model-based and nonparametric estimators of the distribution functions of time-to-death on the complementary log-log scale. Plots of these functions are parallel under proportional hazards: <>= 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) @ From the obtained plot, we see that the smooth model-based transformation functions closely follow the transformed ECDFs (step function). This reflects that the transformation function fits well. The ECDF-based step functions are parallel, revealing that the proportional hazards function of the Cox model is appropriate. To assess the fit of the Gaussian copula, we fit a more flexible additive transformation model for the outcome given the covariate thorax length as splines for each treatment group: <>= 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) @ Because the smooth functions reflecting the effect of thorax length in both treatment group are monotone, the Gaussian copula seem to be appropriate. \section{Conclusion} The present vignette complements the theoretical work in \cite{Dandl_2025_nami} by discussing three application cases of \NAMI with diverse outcomes. It was showed that \NAMI can identify the marginal effect of the treatment of interest with higher precision (narrower confidence intervals) as compared to unadjusted analysis if at least one relevant prognostic covariate exists. In addition, model diagnosis tools were discussed that equip users with helpful tools to assess model fits and to examine underlying model assumptions. <>= 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) @ \bibliography{mlt} \end{document}