
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{mtram}
%\VignetteDepends{variables, basefun, mlt, tram, survival, lme4, gridExtra, lattice, latticeExtra, colorspace, HSAUR3, mvtnorm, ordinalCont}

\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}
%% need no \usepackage{Sweave.sty}
%%\usepackage[nolists]{endfloat}

\newcommand{\cmd}[1]{\texttt{#1()}}


\usepackage{tikz}
\usetikzlibrary{shapes,arrows,chains}
\usepackage{verbatim}

<<mtram-setup, echo = FALSE, results = "hide", message = FALSE, warning = FALSE>>=
set.seed(290875)

pkgs <- sapply(c("mlt", "survival", "tram", "lme4", "gridExtra", 
         "lattice", "latticeExtra", "mvtnorm", "ordinalCont"), require, char = TRUE)

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)
@

\newcommand{\TODO}[1]{{\color{red} #1}}

\newcommand\Torsten[1]{{\color{blue}Torsten: ``#1''}}

\newcommand\norm[1]{\left\lVert#1\right\rVert}

\newcommand{\CTM}{CTM Likelihood Boosting}
\newcommand{\STM}{STM Likelihood Boosting}

\newcommand{\etc}{\textit{etc.}}

\usepackage{booktabs}

\newcommand{\expit}{\text{expit}}

%%% mlt
%% rv
\newcommand{\rZ}{Z}
\newcommand{\rY}{Y}
\newcommand{\rX}{\mX}
\newcommand{\rz}{z}
\newcommand{\ry}{y}
\newcommand{\rx}{\xvec}
\newcommand{\ru}{\uvec}
\newcommand{\erx}{x}
%% sigma algebra
\newcommand{\sA}{\mathfrak{A}}
\newcommand{\sAZ}{\mathfrak{B}}
\newcommand{\sAY}{\mathfrak{C}}
\newcommand{\esA}{A}
\newcommand{\esAZ}{B}
\newcommand{\esAY}{C}
%% sample spaces
\newcommand{\sam}{\Omega}
\newcommand{\samZ}{\RR}
\newcommand{\samY}{\Xi}
\newcommand{\samX}{\chi}
%% measureable spaces
\newcommand{\ms}{(\sam, \sA)}
\newcommand{\msZ}{(\samZ, \sAZ)}
\newcommand{\msY}{(\samY, \sAY)}
%% probability spaces
\newcommand{\ps}{(\sam, \sA, \Prob)}
\newcommand{\psZ}{(\samZ, \sAZ, \Prob_\rZ)}
\newcommand{\psY}{(\samY, \sAY, \Prob_\rY)}
%% distributions
\newcommand{\pZ}{F}
\newcommand{\pY}{F_\rY}
\newcommand{\oY}{O_\rY}
\newcommand{\oYx}{O_{\rY \mid \rX = \rx}}
\newcommand{\hatpY}{\hat{F}_{\rY,N}}
\newcommand{\hatpYx}{\hat{F}_{\rY \mid \rX = \rx, N}}
\newcommand{\PYx}{\Prob_{\rY \mid \rX = \rx}}
\newcommand{\PYX}{\Prob_{\rY, \rX}}
\newcommand{\PY}{\Prob_{\rY}}
\newcommand{\pN}{\Phi}
\newcommand{\pSL}{F_{\SL}}
\newcommand{\pMEV}{F_{\MEV}}
\newcommand{\pExp}{F_{\ExpD}}
\newcommand{\pSW}{F_{\SW}}
\newcommand{\pYx}{F_{\rY \mid \rX = \rx}}
\newcommand{\pYA}{F_{\rY \mid \rX = A}}
\newcommand{\pYB}{F_{\rY \mid \rX = B}}
\newcommand{\qZ}{F^{-1}_\rZ}
\newcommand{\qY}{F^{-1}_\rY}
\newcommand{\dZ}{f}
\newcommand{\dY}{f_\rY}
\newcommand{\hatdY}{\hat{f}_{\rY, N}}
\newcommand{\dYx}{f_{\rY \mid \rX = \rx}}
\newcommand{\hazY}{\lambda_\rY}
\newcommand{\HazY}{\Lambda_\rY}
\newcommand{\HazYx}{\Lambda_{\rY \mid \rX = \rx}}

\newcommand{\hathazY}{\hat{\lambda}_{\rY, N}}
\newcommand{\hatHazY}{\hat{\Lambda}_{\rY, N}}
%% measures
\newcommand{\measureY}{\mu}
\newcommand{\lebesgue}{\mu_L}
\newcommand{\counting}{\mu_C}
%% trafo
\newcommand{\g}{g}
\newcommand{\h}{h}
\newcommand{\s}{\svec}
\newcommand{\hY}{h_\rY}
\newcommand{\hx}{h_\rx}
\newcommand{\hs}{\mathcal{H}}
\newcommand{\basisy}{\avec}
\newcommand{\bern}[1]{\avec_{\text{Bs},#1}}
\newcommand{\bernx}[1]{\bvec_{\text{Bs},#1}}
\newcommand{\basisx}{\bvec}
\newcommand{\basisyx}{\cvec}
\newcommand{\m}{m}
\newcommand{\lik}{\mathcal{L}}
\newcommand{\parm}{\varthetavec}
\newcommand{\eparm}{\vartheta}
\newcommand{\dimparm}{P}
\newcommand{\dimparmx}{Q}
\newcommand{\dimparmvar}{M}
\newcommand{\dimparmrand}{R}
\newcommand{\shiftparm}{\betavec}
\newcommand{\scaleparm}{\xivec}
\newcommand{\varparm}{\gammavec}
\newcommand{\eshiftparm}{\beta}
\newcommand{\escaleparm}{\xi}

\newcommand{\ie}{\textit{i.e.}~}
\newcommand{\eg}{\textit{e.g.}~}

\renewcommand{\Prob}{\mathbb{P}}
\newcommand{\Ex}{\mathbb{E}}
\newcommand{\RR}{\mathbb{R}}
\newcommand{\NN}{\mathbb{N}}
\newcommand{\eps}{\varepsilon}
\newcommand{\prodname}{tensor }
\newcommand{\Null}{\mathbf{0}}
\newcommand{\FI}{\mF}

\usepackage{dsfont}
\newcommand{\I}{\mathds{1}}



\def \dsP {\text{$\mathds{P}$}}
\def \dsE {\text{$\mathds{E}$}}
\def \dsR {\text{$\mathds{R}$}}
\def \dsN {\text{$\mathds{N}$}}


% Math Operators

 \DeclareMathOperator{\logit}{logit}
 \DeclareMathOperator{\LRT}{LRT}
 \DeclareMathOperator{\RLRT}{RLRT}
 \DeclareMathOperator{\Cov}{Cov}
 \DeclareMathOperator{\Cor}{Cor}
 \DeclareMathOperator{\Var}{Var}
 \DeclareMathOperator{\EW}{\dsE}
 \DeclareMathOperator{\D}{D}
 \DeclareMathOperator{\Bias}{Bias}
 \DeclareMathOperator{\MSE}{MSE}
 \DeclareMathOperator{\PLS}{PLS}
 \DeclareMathOperator{\rank}{rank}
 \DeclareMathOperator{\ncol}{ncol}
 \DeclareMathOperator{\pen}{pen}
 \DeclareMathOperator{\const}{const}
 \DeclareMathOperator{\diag}{diag}
 \DeclareMathOperator{\blockdiag}{blockdiag}
 \DeclareMathOperator{\df}{df}
 \DeclareMathOperator{\trace}{tr}
 \DeclareMathOperator{\iid}{i.i.d.}
 \DeclareMathOperator{\ind}{ind.}
 \DeclareMathOperator{\obs}{obs}
 \DeclareMathOperator{\acos}{acos}
 \DeclareMathOperator{\spat}{spat}
 \DeclareMathOperator{\fix}{{fix}}
 \DeclareMathOperator{\ran}{{ran}}
 \DeclareMathOperator*{\argmin}{{arg\,min}}
 \DeclareMathOperator*{\argmax}{{arg\,max}}
 \DeclareMathOperator{\BIC}{{BIC}}
 \DeclareMathOperator{\DIC}{{DIC}}
 \DeclareMathOperator{\AIC}{{AIC}}
 \DeclareMathOperator{\mAIC}{{mAIC}}
 \DeclareMathOperator{\cAIC}{{cAIC}}

% Distributions

 \DeclareMathOperator{\ND}{N}
 \DeclareMathOperator{\TND}{TN}
 \DeclareMathOperator{\UD}{U}
 \DeclareMathOperator{\GaD}{Ga}
 \DeclareMathOperator{\tD}{t}
 \DeclareMathOperator{\IGD}{IG}
 \DeclareMathOperator{\IWD}{IW}
 \DeclareMathOperator{\PoD}{Po}
 \DeclareMathOperator{\ExpD}{Exp}
 \DeclareMathOperator{\LapD}{Lap}
 \DeclareMathOperator{\MuD}{Mu}
 \DeclareMathOperator{\DirD}{Dir}
 \DeclareMathOperator{\PDD}{PD}
 \DeclareMathOperator{\BeD}{Be}
 \DeclareMathOperator{\BD}{B}
 \DeclareMathOperator{\DPD}{DP}
 \DeclareMathOperator{\KSD}{KS}
 \DeclareMathOperator{\SL}{SL}
 \DeclareMathOperator{\MEV}{MEV}
 \DeclareMathOperator{\SW}{SW}
 \DeclareMathOperator{\Chi1}{\chi^2_1}
 \DeclareMathOperator{\WD}{W}



% Boldface vectors and matrices

\def \avec {\text{\boldmath$a$}}    \def \mA {\text{\boldmath$A$}}
\def \bvec {\text{\boldmath$b$}}    \def \mB {\text{\boldmath$B$}}
\def \cvec {\text{\boldmath$c$}}    \def \mC {\text{\boldmath$C$}}
\def \dvec {\text{\boldmath$d$}}    \def \mD {\text{\boldmath$D$}}
\def \evec {\text{\boldmath$e$}}    \def \mE {\text{\boldmath$E$}}
\def \fvec {\text{\boldmath$f$}}    \def \mF {\text{\boldmath$F$}}
\def \gvec {\text{\boldmath$g$}}    \def \mG {\text{\boldmath$G$}}
\def \hvec {\text{\boldmath$h$}}    \def \mH {\text{\boldmath$H$}}
\def \ivec {\text{\boldmath$i$}}    \def \mI {\text{\boldmath$I$}}
\def \jvec {\text{\boldmath$j$}}    \def \mJ {\text{\boldmath$J$}}
\def \kvec {\text{\boldmath$k$}}    \def \mK {\text{\boldmath$K$}}
\def \lvec {\text{\boldmath$l$}}    \def \mL {\text{\boldmath$L$}}
\def \mvec {\text{\boldmath$m$}}    \def \mM {\text{\boldmath$M$}}
\def \nvec {\text{\boldmath$n$}}    \def \mN {\text{\boldmath$N$}}
\def \ovec {\text{\boldmath$o$}}    \def \mO {\text{\boldmath$O$}}
\def \pvec {\text{\boldmath$p$}}    \def \mP {\text{\boldmath$P$}}
\def \qvec {\text{\boldmath$q$}}    \def \mQ {\text{\boldmath$Q$}}
\def \rvec {\text{\boldmath$r$}}    \def \mR {\text{\boldmath$R$}}
\def \svec {\text{\boldmath$s$}}    \def \mS {\text{\boldmath$S$}}
\def \tvec {\text{\boldmath$t$}}    \def \mT {\text{\boldmath$T$}}
\def \uvec {\text{\boldmath$u$}}    \def \mU {\text{\boldmath$U$}}
\def \vvec {\text{\boldmath$v$}}    \def \mV {\text{\boldmath$V$}}
\def \wvec {\text{\boldmath$w$}}    \def \mW {\text{\boldmath$W$}}
\def \xvec {\text{\boldmath$x$}}    \def \mX {\text{\boldmath$X$}}
\def \yvec {\text{\boldmath$y$}}    \def \mY {\text{\boldmath$Y$}}
\def \zvec {\text{\boldmath$z$}}    \def \mZ {\text{\boldmath$Z$}}

 \def \calA {\mathcal A}
 \def \calB {\mathcal B}
 \def \calC {\mathcal C}
 \def \calD {\mathcal D}
 \def \calE {\mathcal E}
 \def \calF {\mathcal F}
 \def \calG {\mathcal G}
 \def \calH {\mathcal H}
 \def \calI {\mathcal I}
 \def \calJ {\mathcal J}
 \def \calK {\mathcal K}
 \def \calL {\mathcal L}
 \def \calM {\mathcal M}
 \def \calN {\mathcal N}
 \def \calO {\mathcal O}
 \def \calP {\mathcal P}
 \def \calQ {\mathcal Q}
 \def \calR {\mathcal R}
 \def \calS {\mathcal S}
 \def \calT {\mathcal T}
 \def \calU {\mathcal U}
 \def \calV {\mathcal V}
 \def \calW {\mathcal W}
 \def \calX {\mathcal X}
 \def \calY {\mathcal Y}
 \def \calZ {\mathcal Z}

\def \ahatvec {\text{\boldmath$\hat a$}}    \def \mhatA {\text{\boldmath$\hat A$}}
\def \bhatvec {\text{\boldmath$\hat b$}}    \def \mhatB {\text{\boldmath$\hat B$}}
\def \chatvec {\text{\boldmath$\hat c$}}    \def \mhatC {\text{\boldmath$\hat C$}}
\def \dhatvec {\text{\boldmath$\hat d$}}    \def \mhatD {\text{\boldmath$\hat D$}}
\def \ehatvec {\text{\boldmath$\hat e$}}    \def \mhatE {\text{\boldmath$\hat E$}}
\def \fhatvec {\text{\boldmath$\hat f$}}    \def \mhatF {\text{\boldmath$\hat F$}}
\def \ghatvec {\text{\boldmath$\hat g$}}    \def \mhatG {\text{\boldmath$\hat G$}}
\def \hhatvec {\text{\boldmath$\hat h$}}    \def \mhatH {\text{\boldmath$\hat H$}}
\def \ihatvec {\text{\boldmath$\hat i$}}    \def \mhatI {\text{\boldmath$\hat I$}}
\def \jhatvec {\text{\boldmath$\hat j$}}    \def \mhatJ {\text{\boldmath$\hat J$}}
\def \khatvec {\text{\boldmath$\hat k$}}    \def \mhatK {\text{\boldmath$\hat K$}}
\def \lhatvec {\text{\boldmath$\hat l$}}    \def \mhatL {\text{\boldmath$\hat L$}}
\def \mhatvec {\text{\boldmath$\hat m$}}    \def \mhatM {\text{\boldmath$\hat M$}}
\def \nhatvec {\text{\boldmath$\hat n$}}    \def \mhatN {\text{\boldmath$\hat N$}}
\def \ohatvec {\text{\boldmath$\hat o$}}    \def \mhatO {\text{\boldmath$\hat O$}}
\def \phatvec {\text{\boldmath$\hat p$}}    \def \mhatP {\text{\boldmath$\hat P$}}
\def \qhatvec {\text{\boldmath$\hat q$}}    \def \mhatQ {\text{\boldmath$\hat Q$}}
\def \rhatvec {\text{\boldmath$\hat r$}}    \def \mhatR {\text{\boldmath$\hat R$}}
\def \shatvec {\text{\boldmath$\hat s$}}    \def \mhatS {\text{\boldmath$\hat S$}}
\def \thatvec {\text{\boldmath$\hat t$}}    \def \mhatT {\text{\boldmath$\hat T$}}
\def \uhatvec {\text{\boldmath$\hat u$}}    \def \mhatU {\text{\boldmath$\hat U$}}
\def \vhatvec {\text{\boldmath$\hat v$}}    \def \mhatV {\text{\boldmath$\hat V$}}
\def \whatvec {\text{\boldmath$\hat w$}}    \def \mhatW {\text{\boldmath$\hat W$}}
\def \xhatvec {\text{\boldmath$\hat x$}}    \def \mhatX {\text{\boldmath$\hat X$}}
\def \yhatvec {\text{\boldmath$\hat y$}}    \def \mhatY {\text{\boldmath$\hat Y$}}
\def \zhatvec {\text{\boldmath$\hat z$}}    \def \mhatZ {\text{\boldmath$\hat Z$}}


\def \atildevec {\text{\boldmath$\tilde a$}}    \def \mtildeA {\text{\boldmath$\tilde A$}}
\def \btildevec {\text{\boldmath$\tilde b$}}    \def \mtildeB {\text{\boldmath$\tilde B$}}
\def \ctildevec {\text{\boldmath$\tilde c$}}    \def \mtildeC {\text{\boldmath$\tilde C$}}
\def \dtildevec {\text{\boldmath$\tilde d$}}    \def \mtildeD {\text{\boldmath$\tilde D$}}
\def \etildevec {\text{\boldmath$\tilde e$}}    \def \mtildeE {\text{\boldmath$\tilde E$}}
\def \ftildevec {\text{\boldmath$\tilde f$}}    \def \mtildeF {\text{\boldmath$\tilde F$}}
\def \gtildevec {\text{\boldmath$\tilde g$}}    \def \mtildeG {\text{\boldmath$\tilde G$}}
\def \htildevec {\text{\boldmath$\tilde h$}}    \def \mtildeH {\text{\boldmath$\tilde H$}}
\def \itildevec {\text{\boldmath$\tilde i$}}    \def \mtildeI {\text{\boldmath$\tilde I$}}
\def \jtildevec {\text{\boldmath$\tilde j$}}    \def \mtildeJ {\text{\boldmath$\tilde J$}}
\def \ktildevec {\text{\boldmath$\tilde k$}}    \def \mtildeK {\text{\boldmath$\tilde K$}}
\def \ltildevec {\text{\boldmath$\tilde l$}}    \def \mtildeL {\text{\boldmath$\tilde L$}}
\def \mtildevec {\text{\boldmath$\tilde m$}}    \def \mtildeM {\text{\boldmath$\tilde M$}}
\def \ntildevec {\text{\boldmath$\tilde n$}}    \def \mtildeN {\text{\boldmath$\tilde N$}}
\def \otildevec {\text{\boldmath$\tilde o$}}    \def \mtildeO {\text{\boldmath$\tilde O$}}
\def \ptildevec {\text{\boldmath$\tilde p$}}    \def \mtildeP {\text{\boldmath$\tilde P$}}
\def \qtildevec {\text{\boldmath$\tilde q$}}    \def \mtildeQ {\text{\boldmath$\tilde Q$}}
\def \rtildevec {\text{\boldmath$\tilde r$}}    \def \mtildeR {\text{\boldmath$\tilde R$}}
\def \stildevec {\text{\boldmath$\tilde s$}}    \def \mtildeS {\text{\boldmath$\tilde S$}}
\def \ttildevec {\text{\boldmath$\tilde t$}}    \def \mtildeT {\text{\boldmath$\tilde T$}}
\def \utildevec {\text{\boldmath$\tilde u$}}    \def \mtildeU {\text{\boldmath$\tilde U$}}
\def \vtildevec {\text{\boldmath$\tilde v$}}    \def \mtildeV {\text{\boldmath$\tilde V$}}
\def \wtildevec {\text{\boldmath$\tilde w$}}    \def \mtildeW {\text{\boldmath$\tilde W$}}
\def \xtildevec {\text{\boldmath$\tilde x$}}    \def \mtildeX {\text{\boldmath$\tilde X$}}
\def \ytildevec {\text{\boldmath$\tilde y$}}    \def \mtildeY {\text{\boldmath$\tilde Y$}}
\def \ztildevec {\text{\boldmath$\tilde z$}}    \def \mtildeZ {\text{\boldmath$\tilde Z$}}

\def \alphavec        {\text{\boldmath$\alpha$}}
\def \betavec         {\text{\boldmath$\beta$}}
\def \gammavec        {\text{\boldmath$\gamma$}}
\def \deltavec        {\text{\boldmath$\delta$}}
\def \epsilonvec      {\text{\boldmath$\epsilon$}}
\def \varepsilonvec   {\text{\boldmath$\varepsilon$}}
\def \zetavec         {\text{\boldmath$\zeta$}}
\def \etavec          {\text{\boldmath$\eta$}}
\def \thetavec        {\text{\boldmath$\theta$}}
\def \varthetavec     {\text{\boldmath$\vartheta$}}
\def \iotavec         {\text{\boldmath$\iota$}}
\def \kappavec        {\text{\boldmath$\kappa$}}
\def \lambdavec       {\text{\boldmath$\lambda$}}
\def \muvec           {\text{\boldmath$\mu$}}
\def \nuvec           {\text{\boldmath$\nu$}}
\def \xivec           {\text{\boldmath$\xi$}}
\def \pivec           {\text{\boldmath$\pi$}}
\def \varpivec        {\text{\boldmath$\varpi$}}
\def \rhovec          {\text{\boldmath$\rho$}}
\def \varrhovec       {\text{\boldmath$\varrho$}}
\def \sigmavec        {\text{\boldmath$\sigma$}}
\def \varsigmavec     {\text{\boldmath$\varsigma$}}
\def \tauvec          {\text{\boldmath$\tau$}}
\def \upsilonvec      {\text{\boldmath$\upsilon$}}
\def \phivec          {\text{\boldmath$\phi$}}
\def \varphivec       {\text{\boldmath$\varphi$}}
\def \psivec          {\text{\boldmath$\psi$}}
\def \chivec          {\text{\boldmath$\chi$}}
\def \omegavec        {\text{\boldmath$\omega$}}

\def \alphahatvec        {\text{\boldmath$\hat \alpha$}}
\def \betahatvec         {\text{\boldmath$\hat \beta$}}
\def \gammahatvec        {\text{\boldmath$\hat \gamma$}}
\def \deltahatvec        {\text{\boldmath$\hat \delta$}}
\def \epsilonhatvec      {\text{\boldmath$\hat \epsilon$}}
\def \varepsilonhatvec   {\text{\boldmath$\hat \varepsilon$}}
\def \zetahatvec         {\text{\boldmath$\hat \zeta$}}
\def \etahatvec          {\text{\boldmath$\hat \eta$}}
\def \thetahatvec        {\text{\boldmath$\hat \theta$}}
\def \varthetahatvec     {\text{\boldmath$\hat \vartheta$}}
\def \iotahatvec         {\text{\boldmath$\hat \iota$}}
\def \kappahatvec        {\text{\boldmath$\hat \kappa$}}
\def \lambdahatvec       {\text{\boldmath$\hat \lambda$}}
\def \muhatvec           {\text{\boldmath$\hat \mu$}}
\def \nuhatvec           {\text{\boldmath$\hat \nu$}}
\def \xihatvec           {\text{\boldmath$\hat \xi$}}
\def \pihatvec           {\text{\boldmath$\hat \pi$}}
\def \varpihatvec        {\text{\boldmath$\hat \varpi$}}
\def \rhohatvec          {\text{\boldmath$\hat \rho$}}
\def \varrhohatvec       {\text{\boldmath$\hat \varrho$}}
\def \sigmahatvec        {\text{\boldmath$\hat \sigma$}}
\def \varsigmahatvec     {\text{\boldmath$\hat \varsigma$}}
\def \tauhatvec          {\text{\boldmath$\hat \tau$}}
\def \upsilonhatvec      {\text{\boldmath$\hat \upsilon$}}
\def \phihatvec          {\text{\boldmath$\hat \phi$}}
\def \varphihatvec       {\text{\boldmath$\hat \varphi$}}
\def \psihatvec          {\text{\boldmath$\hat \psi$}}
\def \chihatvec          {\text{\boldmath$\hat \chi$}}
\def \omegahatvec        {\text{\boldmath$\hat \omega$}}

\def \alphatildevec        {\text{\boldmath$\tilde \alpha$}}
\def \betatildevec         {\text{\boldmath$\tilde \beta$}}
\def \gammatildevec        {\text{\boldmath$\tilde \gamma$}}
\def \deltatildevec        {\text{\boldmath$\tilde \delta$}}
\def \epsilontildevec      {\text{\boldmath$\tilde \epsilon$}}
\def \varepsilontildevec   {\text{\boldmath$\tilde \varepsilon$}}
\def \zetatildevec         {\text{\boldmath$\tilde \zeta$}}
\def \etatildevec          {\text{\boldmath$\tilde \eta$}}
\def \thetatildevec        {\text{\boldmath$\tilde \theta$}}
\def \varthetatildevec     {\text{\boldmath$\tilde \vartheta$}}
\def \iotatildevec         {\text{\boldmath$\tilde \iota$}}
\def \kappatildevec        {\text{\boldmath$\tilde \kappa$}}
\def \lambdatildevec       {\text{\boldmath$\tilde \lambda$}}
\def \mutildevec           {\text{\boldmath$\tilde \mu$}}
\def \nutildevec           {\text{\boldmath$\tilde \nu$}}
\def \xitildevec           {\text{\boldmath$\tilde \xi$}}
\def \pitildevec           {\text{\boldmath$\tilde \pi$}}
\def \varpitildevec        {\text{\boldmath$\tilde \varpi$}}
\def \rhotildevec          {\text{\boldmath$\tilde \rho$}}
\def \varrhotildevec       {\text{\boldmath$\tilde \varrho$}}
\def \sigmatildevec        {\text{\boldmath$\tilde \sigma$}}
\def \varsigmatildevec     {\text{\boldmath$\tilde \varsigma$}}
\def \tautildevec          {\text{\boldmath$\tilde \tau$}}
\def \upsilontildevec      {\text{\boldmath$\tilde \upsilon$}}
\def \phitildevec          {\text{\boldmath$\tilde \phi$}}
\def \varphitildevec       {\text{\boldmath$\tilde \varphi$}}
\def \psitildevec          {\text{\boldmath$\tilde \psi$}}
\def \chitildevec          {\text{\boldmath$\tilde \chi$}}
\def \omegatildevec        {\text{\boldmath$\tilde \omega$}}

\def \mGamma   {\mathbf{\Gamma}}
\def \mDelta   {\mathbf{\Delta}}
\def \mTheta   {\mathbf{\Theta}}
\def \mLambda  {\mathbf{\Lambda}}
\def \mXi      {\mathbf{\Xi}}
\def \mPi      {\mathbf{\Pi}}
\def \mSigma   {\mathbf{\Sigma}}
\def \mUpsilon {\mathbf{\Upsilon}}
\def \mPhi     {\mathbf{\Phi}}
\def \mPsi     {\mathbf{\Psi}}
\def \mOmega   {\mathbf{\Omega}}

\def \mhatGamma   {\mathbf{\hat \Gamma}}
\def \mhatDelta   {\mathbf{\hat \Delta}}
\def \mhatTheta   {\mathbf{\hat \Theta}}
\def \mhatLambda  {\mathbf{\hat \Lambda}}
\def \mhatXi      {\mathbf{\hat \Xi}}
\def \mhatPi      {\mathbf{\hat \Pi}}
\def \mhatSigma   {\mathbf{\hat \Sigma}}
\def \mhatUpsilon {\mathbf{\hat \Upsilon}}
\def \mhatPhi     {\mathbf{\hat \Phi}}
\def \mhatPsi     {\mathbf{\hat \Psi}}
\def \mhatOmega   {\mathbf{\hat \Omega}}

\def \nullvec {\mathbf{0}}
\def \onevec {\mathbf{1}}

%%% theorems
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{coro}{Corollary}
\newtheorem{defn}{Definition}
\newtheorem{remark}{Remark}

\newcommand{\ubar}[1]{\underaccent{\bar}{#1}}


\renewcommand{\thefootnote}{}

%% code commands
\newcommand{\Rclass}[1]{`\code{#1}'}
%% JSS
\author{Torsten Hothorn \\ Universit\"at Z\"urich}
\Plainauthor{Hothorn}

\title{Some Applications of Marginally Interpretable Linear 
       Transformation Models for Clustered Observations}
\Plaintitle{Marginally Interpretable Transformation Models}
\Shorttitle{Marginally Interpretabel Transformation Models}

\Abstract{
Owing to their generality, transformation models can be used to set-up and
compute many interesting regression models for discrete and continuous responses.  This
document focuses on the analysis of clustered observations.  Marginal
predictive distributions are defined by transformation models and their
joint normal distribution depends on a structured covariance matrix. 
Applications with skewed, bounded, and survival continuous outcomes as well
as binary and ordered categorical responses are presented. Data is analysed
by a proof-of-concept implementation of parametric linear transformation models for
clustered observations available in the \pkg{tram} add-on package to the
\proglang{R} system for statistical computing.
}

\Keywords{conditional mixed models, marginal models, marginal predictive
distributions, survival analysis, categorical data analysis}
\Plainkeywords{conditional mixed models, marginal models, marginal predictive
distributions, survival analysis, categorical data analysis}

\Address{
  Torsten Hothorn\\
  Institut f\"ur Epidemiologie, Biostatistik und Pr\"avention \\
  Universit\"at Z\"urich \\
  Hirschengraben 84, CH-8001 Z\"urich, Switzerland \\
  \texttt{Torsten.Hothorn@uzh.ch} \\
  \url{http://tiny.uzh.ch/bh}
}

\begin{document}

<<mtram-citation, echo = FALSE>>=
year <- substr(packageDescription("tram")$Date, 1, 4)
version <- packageDescription("tram")$Version
@
\footnote{Please cite this document as: Torsten Hothorn (\Sexpr{year})
Some Applications of Marginally Interpretable Linear 
Transformation Models for Clustered Observations.
\textsf{R} package vignette version \Sexpr{version}, 
URL \url{https://CRAN.R-project.org/package=tram}.}

% \input{todo}

<<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()
}
@

\section{Introduction}

The purpose of this document is to compare marginally interpretable linear
transformation models for clustered observations 
\citep{Hothorn_2019_mtram} to conventional conditional
formulations of mixed-effects models where such an overlap exists.  In
addition, novel transformation models going beyond the capabilities of
convential mixed-effects models are estimated and interpreted.  A
proof-of-concept implementation (meaning: the algorithms work but need
optimisation, the user interface is very rough) available in package
\pkg{tram} \citep{pkg:tram} is applied. The results presented in this
document can be reproduced from the \code{mtram} package vignette
<<mtram-vignette, eval = FALSE>>=
install.packages("tram")
vignette("mtram", package = "tram")
@

\section{Normal and Non-normal Mixed-effects Models}

First we consider mixed-effects models for reaction times in the sleep
deprivation study \citep{Belenky_Wesensten_Thorne_2003}.  The average
reaction times to a specific task over several days of sleep deprivation are
given for $18$ subjects in Figure~\ref{fig:sleepstudy}.

\begin{figure}[t]
<<mtram-sleep-plot, echo = FALSE>>=
xyplot(Reaction ~ Days | Subject, data = sleepstudy, 
       xlab = "Days of sleep deprivation", ylab = "Average reaction time (in ms)")
@
\caption{Sleep deprivation: Average reaction times to a specific task over several days 
         of sleep deprivation for $18$ subjects from
         \cite{Belenky_Wesensten_Thorne_2003}. \label{fig:sleepstudy}}
\end{figure}

The classical normal linear random-intercept/random-slope model, treating
the study participants as independent observations, is fitted by maximum
likelihood to the data using the \cmd{lmer} function from the \pkg{lme4}
add-on package \citep{pkg:lme4}:
%
<<mtram-sleep_lmer, cache = TRUE>>=
library("lme4")
sleep_lmer <- lmer(Reaction ~ Days + (Days | Subject), 
                   data = sleepstudy, REML = FALSE)
@
%
The corresponding conditional model for subject $i$ reads
%
\begin{eqnarray*}
\Prob(\text{Reaction} \le \ry \mid \text{day}, i) = \Phi\left(\frac{\ry -
\alpha - \beta \text{day} - \alpha_i - \beta_i \text{day}}{\sigma}\right),
\quad (\alpha_i, \beta_i) \sim \ND_2(\nullvec, \mG(\varparm))
\end{eqnarray*}
%
with $\sigma^{-2}\mG = \mLambda(\varparm) \mLambda(\varparm)^\top$ and
%
\begin{eqnarray*}
\mLambda(\varparm) = \left( \begin{array}{cc}
    \gamma_1 & 0  \\
    \gamma_2 & \gamma_3
\end{array} \right), \quad \varparm = (\gamma_1, \gamma_2, \gamma_3)^\top.
\end{eqnarray*}

The same model, however using the alternative parameterisation and an
independent (of \pkg{lme4}, only the \cmd{update} method for Cholesky
factors is reused) gradient-based maximisation of the log-likelihood, is estimated 
in a two-step approach as
<<mtram-sleep_mtram, cache = TRUE>>=
sleep_LM <- Lm(Reaction ~ Days, data = sleepstudy)
sleep_LMmer <- mtram(sleep_LM, ~ (Days | Subject), data = sleepstudy)
@
%
The first call to \cmd{Lm} computes the equivalent of a normal linear
regression model parameterised as a linear transformation model
\emph{ignoring} the longitudinal nature of the observations. The purpose if
to set-up the necessary model infrastructure (model matrices, inverse link
functions, etc.) and to compute reasonable starting values for the fixed
effects. The second call to \cmd{mtram} specifies the random effects
structure (here a correlated pair of random intercept for subject 
and random slope for days) and optimises the likelihood for all model
parameters $\eparm_1, \tilde{\alpha}, \tilde{\beta}$, and $\varparm$
in the model (here also looking at the conditional model for subject $i$)
%
\begin{eqnarray*}
\Prob(\text{Reaction} \le \ry \mid \text{day}, i) = \Phi\left(\eparm_1 \ry + \tilde{\alpha} - \tilde{\beta} \text{day} - \tilde{\alpha}_i - \tilde{\beta}_i \text{day}\right),
\quad (\tilde{\alpha}_i, \tilde{\beta}_i) \sim \ND_2(\nullvec, \mLambda(\varparm) \mLambda(\varparm))
\end{eqnarray*}
%
that is, all fixed and random effect parameters are 
divided by the residual standard deviation $\sigma$ (this is the
reparameterisation applied by \cmd{Lm}).
Of course, the parameter $\eparm_1$, the inverse residual standard
deviation, is ensured to be positive via an additional constraint in the
optimiser maximising the log-likelihood.
%

The log-likelihoods of the two models fitted by \cmd{lmer} and \cmd{mtram} 
are very close
<<mtram-sleep_cmp>>=
logLik(sleep_lmer)
logLik(sleep_LMmer)
@
Looking at the model coefficients, the two procedures lead to almost
identical inverse residual standard deviations
<<mtram-sleep_sd>>=
(sdinv <- 1 / summary(sleep_lmer)$sigma)
coef(sleep_LMmer)["Reaction"]
@
and fixed effects (the slope can be interpreted as inverse coefficient of
variation)
<<mtram-sleep_beta>>=
fixef(sleep_lmer) * c(-1, 1) * sdinv
coef(sleep_LMmer)[c("(Intercept)", "Days")]
@
The random-effect parameters $\varparm$ are also reasonably close
<<mtram-sleep_varparm>>=
sleep_lmer@theta
coef(sleep_LMmer)[-(1:3)]
@
Consequently, the variance-covariance and correlation matrices
<<mtram-sleep_varcorr>>=
sleep_LMmer$G * (1 / sdinv)^2
cov2cor(sleep_LMmer$G * (1 / sdinv)^2)
unclass(VarCorr(sleep_lmer))$Subject
@
are practically equivalent. This result indicates the correctness of the
alternative implementation of normal linear mixed-effects models in the
transformation model framework: \cmd{mtram} reuses some infrastructure from
\pkg{lme4} and \pkg{Matrix}, most importantly fast update methods for
Cholesky factors, but the likelihood and corresponding optimisation relies
on an independent implementation. So why are we doing this? Because
\cmd{mtram} is able to deal with models or likelihoods 
not available in \pkg{lme4}, for example the likelihood for
interval-censored observations.

Let's assume that the timing of the reaction times was less accurate than suggested by the
numerical representation of the results. The following code
<<mtram-sleep-Surv>>=
library("survival")
sleepstudy$Reaction_I <- with(sleepstudy, Surv(Reaction - 20, Reaction + 20, 
                                               type = "interval2"))
sleepstudy$Reaction_I[1:5]
@
converts the outcome to interval-censored values, where each interval has
length $40$. The above mixed model can now be estimated by maximising the
likelihood corresponding to interval-censored observations:
<<mtram-sleep-interval, cache = TRUE>>=
sleep_LM_I <- Lm(Reaction_I ~ Days, data = sleepstudy)
sleep_LMmer_I <- mtram(sleep_LM_I, ~ (Days | Subject), data = sleepstudy)
@
Of course, the log-likelihood changes (because this is a log-probability and
not a log-density of a continuous distribution) but the parameter estimates are reasonably close
<<mtram-sleep-interval-results>>=
logLik(sleep_LMmer_I)
coef(sleep_LMmer_I)
coef(sleep_LMmer)
@

The next question is if the normal assumption for reaction times is
appropriate. In the transformation world, this assumption is simple to
assess because we can easily (theoretically and in-silico) switch to 
the non-normal linear mixed-effects transformation model
%
\begin{eqnarray*}
\Prob(\text{Reaction} \le \ry \mid \text{day}, i) =
\Phi\left(\h(\ry) - \tilde{\beta} \text{day} - \tilde{\alpha}_i - \tilde{\beta}_i \text{day}\right),
\quad (\tilde{\alpha}_i, \tilde{\beta}_i) \sim \ND_2(\nullvec, \mLambda(\varparm) \mLambda(\varparm))
\end{eqnarray*}
%
where $\h(\ry) = \basisy(\ry)^\top \parm$ represents a monotone non-decreasing
transformation function. The function implementing such a more flexible
model in named in honor of the first paper on the analysis of
transformed responses by \cite{BoxCox_1964} but it \emph{does not} simply apply
what is known as a Box-Cox transformation. Bernstein polynomials
$\h(\ry) = \basisy(\ry)^\top \parm$ under suitable constraints are applied instead by
<<mtram-sleep_BoxCox, cache = TRUE>>=
sleep_BC <- BoxCox(Reaction ~ Days, data = sleepstudy)
sleep_BCmer <- mtram(sleep_BC, ~ (Days | Subject), data = sleepstudy, 
                     Hessian = TRUE)
logLik(sleep_BCmer)
@
%
The increase in the log-likelihood compared to the normal model is not a big
surprise.  Plotting the transformation function $\h(\ry) = \basisy(\ry)^\top \parm$ as
a function of reaction time can help to assess deviations from normality
because the latter assumption implies a linear transformation function. 
Figure~\ref{fig:sleepstudy_trafo} clearly indicates that models allowing a
certain skewness of reaction times will provide a better fit to the data.
This might also not come as a big surprise to experienced data analysts.

\begin{figure}[t]
<<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))
@
\caption{Sleep deprivation: Data-driven transformation $\hat{\h}$ 
         of average reaction times to sleep deprivation. The non-linearity
         induces a non-normal marginal distribution function of reaction times.
         \label{fig:sleepstudy_trafo}}
\end{figure}

However, what does this finding mean in terms of a direct comparison of the
model and the data?  Looking at the marginal cumulative distribution
functions of average reaction time conditional on days of sleep deprivation
in Figure~\ref{fig:sleepstudy_ecdf} one finds that the non-normal marginal
transformation models provided a better fit to the marginal empirical
cumulative distribution functions than the normal marginal models. 
Especially for short reaction times in the first week of sleep deprivation,
the yellowish marginal cumulative distribution is much closer to the
empirical cumulative distribution function representing the marginal
distribution of reaction times at each single day of study participation.

\begin{figure}[t]
<<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-norma 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))
@
\caption{Sleep deprivation: Marginal distribution of reaction times, separately for each day of
         study participation. The grey step-function corresponds to the
         empirical cumulative distribution function, the blue line to the
         marginal cumulative distribution of a normal linear mixed-effects model, and
         the yellowish line to a non-normal linear mixed-effects transformation
         model. \label{fig:sleepstudy_ecdf}}
\end{figure}

It should be noted that the small positive correlation between random intercept and random slope
observed in the normal linear mixed-effects model turned into a negative
correlation in this non-normal model 
<<mtram-sleep_corr>>=
cov2cor(sleep_BCmer$G)
@
What is the uncertainty associated with this parameter? The correlation is a
non-linear function of $\varparm$ and therefore the direct computation of
confidence intervals questionable. However, we
can extract an estimate of the covariance of the estimated model parameters
from the model and, relying on the asymptotic normality of the maximum likelihood
estimators, we can sample from the asymptotic distribution
of the variance of the random intercept $\tilde{\alpha}$, the random slope
$\tilde{\beta}$, and their correlation
<<mtram-sleep_vcov>>=
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])
})
@
The $95\%$ confidence intervals
<<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))
@
indicate rather strong unobserved heterogeneity affecting the intercept and
less pronouned variability in the slope. There is only weak information
about the correlation of the two random effects contained in the data.

The downside of this approach is that, although the model is nicely
interpretable on the scale of marginal or conditional distribution
functions, the direct interpretation of the fixed effect $\tilde{\beta}$ is
not very straightforward because it corresponds to the conditional mean
\emph{after} transforming the outcome.  This interpretability issue can be
addressed by exchanging the probit link to a logit link in
Section~\ref{sec:logit}.

\section{Binary Probit Mixed-effects Models}

Here we compare different implementations of binary probit mixed models for
the notoriously difficult toe nail data \citep{Backer_Vroey_1998}. The
outcome was categorised to two levels (this being probably the root of all
troubles) and a conditional density plot (Figure~\ref{fig:toenail}) suggests
an improvement in both treatment groups over time, however with a more rapid
advance in patients treated with terbinafine.

\begin{figure}[t]
<<mtram-toenail-plot, echo = FALSE>>=
data("toenail", package = "HSAUR3")
rlev <- levels(toenail$patientID)[xtabs(~ patientID, 
                                        data = toenail) == 1]
toenail <- subset(toenail, !patientID %in% rlev)
toenail$patientID <- toenail$patientID[, drop = TRUE]
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 = "")
@
\caption{Toe nail data: Conditional density plot of two outcome classes
(none or mild vs.~moderate or severe) under two treatments.
\label{fig:toenail}}
\end{figure}

The random intercept probit model fitted by Laplace and Adaptive
Gauss-Hermite Quadrature (AGQ) approximations to the
likelihood give quite different results:
<<mtram-toenail_glmer_RI, cache = TRUE>>=
### 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
@

%%%% this takes too long for CRAN
%The sequential reduction (SR) algorithm \citep{Ogden_2015} gives results close to
%AGQ
%<<mtram-toenail_glmmsr_RI, cache = TRUE>>=
%library("glmmsr")
%toenail_glmm_RI_3 <- 
%    glmm(outcome ~ treatment * time + (1 | patientID),
%         data = toenail, family = binomial(link = "probit"), 
%         method = "SR", control = list(nSL = 3))
%summary(toenail_glmm_RI_3)
%@

%%% start of static output
The sequential reduction (SR) algorithm \citep{Ogden_2015} gives results close to
AGQ
\begin{Schunk}
\begin{Sinput}
R> library("glmmsr")
R> toenail_glmm_RI_3 <- 
+      glmm(outcome ~ treatment * time + (1 | patientID),
+           data = toenail, family = binomial(link = "probit"), 
+           method = "SR", control = list(nSL = 3))
\end{Sinput}
\begin{Soutput}
Fitting the model......................... done.
\end{Soutput}
\begin{Sinput}
R> summary(toenail_glmm_RI_3)
\end{Sinput}
\begin{Soutput}
Generalized linear mixed model fit by maximum likelihood [glmmFit] 
Likelihood approximation: Sequential reduction at level 3 
 
Family: binomial ( probit ) 
Formula: outcome ~ treatment * time + (1 | patientID)

Random effects:
 Groups    Name        Estimate Std.Error
 patientID (Intercept) 2.119    0.1954   
Number of obs: 1903, groups: patientID, 289; 

Fixed effects: 
                          Estimate Std. Error z value  Pr(>|z|)
(Intercept)               -0.93105    0.23217  4.0102 6.066e-05
treatmentterbinafine      -0.07618    0.30945  0.2462 8.055e-01
time                      -0.19076    0.02060  9.2618 2.010e-20
treatmentterbinafine:time -0.06420    0.03099  2.0713 3.834e-02
\end{Soutput}
\end{Schunk}
%%% end of static output

Because of the probit link, this binary generalised linear model is
equivalent to a linear transformation model and we can thus use the exact
likelihood implemented for the latter model in \cmd{mtram} for parameter
estimation (it is still a bit nasty to set-up a constant transformation function
$\h(\ry) = \alpha$, I plan to add a more convenient interface later)
<<mtram-toenail_mtram_RI, cache = TRUE>>=
m <- ctm(as.basis(~ outcome, data = toenail), 
         shifting = ~ treatment * time, 
         data = toenail, todistr = "Normal")
toenail_probit <- mlt(m, data = toenail, 
                      fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_RI <- 
    mtram(toenail_probit, ~ (1 | patientID), 
          data = toenail, Hessian = TRUE)
logLik(toenail_mtram_RI)
coef(toenail_mtram_RI)
@
For this random intercept model, the exact likelihood is defined as a
one-dimensional integral over the unit interval.  We use sparse grids
\citep{Heiss_Winschel_2008, pkg:SparseGrid} to approximate this integral. 
The integrant is defined by products of normal probabilities, which are
approximated as described by \cite{Matic_Radoicic_2018}.  It is important to
note that this likelihood can be computed as accurately as necessary whereas
Laplace, AGQ, and SR are approximations of limited accuracy. 

The results are very close to SR and AGQ, indicating a very good quality of
the AGQ and SR approximations. We can also compare the corresponding covariances
<<mtram-toenail-hessian>>=
vcov(toenail_glmer_RI_2)
solve(toenail_mtram_RI$Hessian)[1:4, 1:4]
@
%We can compare the in-sample log-likelihood (as defined by the
%transformation model) for parameter estimates obtained from 
%the different procedures
%<<mtram-toenail_RI_mlt_glmmsr>>=
%### logLik of transformation model for glmer (Laplace) parameters
%logLik(toenail_mtram_RI, c(-fixef(toenail_glmer_RI_1), 
%                            toenail_glmer_RI_1@theta))
%### logLik of transformation model for glmer (AGQ) parameters
%logLik(toenail_mtram_RI, c(-fixef(toenail_glmer_RI_2), 
%                            toenail_glmer_RI_2@theta))
%### logLik of transformation model for glmmsr (SR) parameters
%logLik(toenail_mtram_RI, parm = c(-toenail_glmm_RI_3$estim[-1], 
%                                   toenail_glmm_RI_3$estim[1]))
%### logLik of transformation model
%logLik(toenail_mtram_RI)
%@
%Except for the Laplace approximation, the in-sample log-likelihoods
%coincide.

Things get a bit less straightforward when a random slope is added to the
model. The two implementations of the Laplace approximation in packages
\pkg{lme4}
<<mtram-toenail_glmer_RS, cache = TRUE>>=
toenail_glmer_RS <- 
    glmer(outcome ~ treatment * time + (1 + time | patientID),
          data = toenail, family = binomial(link = "probit"))
summary(toenail_glmer_RS)
toenail_glmer_RS@theta
@
and \pkg{glmmsr}
%%%% this takes too long for CRAN
%<<mtram-toenail_glmmsr_RS, cache = TRUE>>=
%toenail_glmm_RS_1 <- 
%    glmm(outcome ~ treatment * time + (1 + time | patientID),
%         data = toenail, family = binomial(link = "probit"), 
%         method = "Laplace")
%toenail_glmm_RS_1$estim[1:3]
%toenail_glmm_RS_1$estim[-(1:3)]
%@
%do not quite agree. Note that the standard deviation of the random intercept
%is twice as large in the \cmd{glmer} output.
%%% start of static output
\begin{Schunk}
\begin{Sinput}
R> toenail_glmm_RS_1 <- 
+      glmm(outcome ~ treatment * time + (1 + time | patientID),
+           data = toenail, family = binomial(link = "probit"), 
+           method = "Laplace")
\end{Sinput}
\begin{Soutput}
Fitting the model..... done.
\end{Soutput}
\begin{Sinput}
R> toenail_glmm_RS_1$estim[1:3]
\end{Sinput}
\begin{Soutput}
[1]  4.9992215 -0.5644628  0.4110986
\end{Soutput}
\begin{Sinput}
R> toenail_glmm_RS_1$estim[-(1:3)]
\end{Sinput}
\begin{Soutput}
[1] -3.49232631  0.04197115 -0.06769437 -0.12100940
\end{Soutput}
\end{Schunk}
do not quite agree. Note that the standard deviation of the random intercept
is twice as large in the \cmd{glmer} output.
%%% end of static output


The optimisation of the exact discrete likelihood in the transformation
framework gives
<<mtram-toenail_mtram_RS, cache = TRUE>>=
toenail_mtram_RS <- 
    mtram(toenail_probit, ~ (1 + time | patientID), 
          data = toenail)
logLik(toenail_mtram_RS)
coef(toenail_mtram_RS)
@
The variance parameters are not too far off the results reported by
\cmd{glmm}, but the fixed effects differ quite a bit. 
%Again, we can
%compare the in-sample log-likelihoods:
%<<mtram-toenail_RS_mlt_glmmsr>>=
%### logLik of transformation model for glmer (Laplace) parameters
%logLik(toenail_mtram_RS, c(-fixef(toenail_glmer_RS), 
%                            toenail_glmer_RS@theta))
%### logLik of transformation model for glmmsr (Laplace) parameters
%logLik(toenail_mtram_RS, parm = c(-toenail_glmm_RS_1$estim[-(1:3)], 
%                                   toenail_glmm_RS_1$estim[1:3]))
%### logLik of transformation model
%logLik(toenail_mtram_RS)
%@

At least in biostatistics, the probit model is less popular than the logit
model owing to the better interpretability of the fixed effects as
conditional log-odds ratios in the latter. Using a logit link, we can use
the transformation approach to compute marginally interpretable
time-dependent log-odds ratios from random intercept transformation logit
models \citep[\cmd{standardise = TRUE} computes model (M2) instead of the
default (M1), see][]{Hothorn_2019_mtram}:
<<mtram-toenail_logit, cache = TRUE>>=
m <- ctm(as.basis(~ outcome, data = toenail), 
         shifting = ~ treatment * time, 
         data = toenail, todistr = "Logistic")
toenail_logit <- mlt(m, data = toenail, 
                     fixed = c("outcomemoderate or severe" = 0))
toenail_mtram_logit <- mtram(toenail_logit, ~ (1 | patientID), 
                              data = toenail)
toenail_mtram_logit_s <- mtram(toenail_logit, ~ (1 | patientID), 
                               data = toenail, standardise = TRUE, 
                               Hessian = TRUE)
@
It is important to note that this model is \emph{not} a logistic 
mixed-effects model and thus we can't expect to obtain identical 
results from \cmd{glmer} as it was (partially) the case for the probit
model.

From the standardised model, we can compute marginally interpretable 
probabilities and odds ratios over 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_s))[1]
cf <- coef(toenail_mtram_logit_s)[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_s <- predict(tmp, newdata = nd, type = "distribution")[1,]
nd$odds <- exp(predict(tmp, newdata = nd, type = "trafo")[1,])
@

We can also sample from the distribution of the maximum likelihood
estimators to obtain an idea about the uncertainty
(Figure~\ref{fig:toenailOR}).

\begin{figure}[t]
<<mtram-toenail_OR_2, dev = "png", cache = TRUE, echo = FALSE, dpi = 300>>=
X <- model.matrix(~ treatment * time, data = nd)
rbeta <- rmvnorm(10000, mean = coef(toenail_mtram_logit_s), 
                 sigma = solve(toenail_mtram_logit_s$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)
@
\caption{Toe nail data: Marginal odds ratio over time (from a logistic
  random intercept model). The blue line represents the maximum likelihood
  estimator, the grey lines are samples from the corresponding distribution.
  \label{fig:toenailOR}}
\end{figure}

From the unstandardised logit and probit models, we can also obtain
marginally interpretable probabilities as (probit)
<<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)
@
and (logit)
<<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,]
@
The marginal time-dependent probabilities obtained from all three models are
very similar as shown in Figure~\ref{fig:toenailprob}.

\begin{figure}[t]
<<mtram-toenail_probplot, echo = FALSE>>=
nd2 <- nd[rep(1:nrow(nd), 3),]
nd2$prob <- c(nd$prob_probit, nd$prob_logit_s, nd$prob_logit)
lev <- c("Probit (M1) = (M2)", "Logit (M1)", "Logit (M2)")
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)")
@
\caption{Toe nail data: Comparison of marginal probabilities obtained from
         a probit linear mixed-effects model and two logistic transformation
         models (M2: with or M1: without marginal log-odds ratio treatment
         effect). \label{fig:toenailprob}}
\end{figure}


\section{Proportional Odds Models for Bounded Responses} \label{sec:logit}

\cite{Manuguerra_Heller_2010} proposed a mixed-effects model for bounded
responses  whose fixed effects can be interpreted as log-odds ratios. 
We fit a transformation model to data from
a randomised controlled trial on chronic neck pain treatment
\citep{Chow_Heller_2006}. The data are visualised in
Figure~\ref{fig:neck_pain}. Subjective neck pain levels were assessed on a 
visual analog scale, that is, on a bounded interval.

\begin{figure}[t]
<<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)
@
\caption{Neck pain: Trajectories of neck pain assessed on a visual analog
         scale with and without low-level laser therapy.
         \label{fig:neck_pain}}
\end{figure}

\cite{Manuguerra_Heller_2010} suggested the conditional model
\begin{eqnarray*}
& & \logit(\Prob(\text{pain} \le \ry \mid \text{treatment}, \text{time}, i)) =
\\
& &  \quad \h(\ry) + \eshiftparm_\text{Active} + \eshiftparm_\text{7 weeks} + 
  \eshiftparm_\text{12 weeks} + \eshiftparm_\text{7 weeks, Active} + 
  \eshiftparm_\text{12 weeks, Active} + \alpha_i
\end{eqnarray*}
with random intercepts $\tilde{\alpha}_i$ such that the odds at baseline, for example, are given by
\begin{eqnarray*}
\frac{\Prob(\text{pain} \le \ry \mid \text{Active}, \text{baseline}, i)}
     {\Prob(\text{pain} > \ry \mid \text{Active}, \text{baseline}, i)} = 
\exp(\eshiftparm_\text{Active}) 
\frac{\Prob(\text{pain} \le \ry \mid \text{Placebo}, \text{baseline}, i)}
     {\Prob(\text{pain} > \ry \mid \text{Placebo}, \text{baseline}, i)}
\end{eqnarray*}

<<mtram-neck_ocm, cache = TRUE, results = "hide">>=
library("ordinalCont")
neck_ocm <- ocm(vas ~ laser * time + (1 | id), data = pain_df, 
                scale = c(0, 1))
@
The results
<<mtram-neck_ocm_summary>>=
summary(neck_ocm)
@
suggest that there is a difference at baseline; the pain distribution of 
subjects in the placebo group on the odds scale is only 
$\Sexpr{round(exp(coef(neck_ocm)[2]) * 100)}\%$ of the odds in the active 
group for any cut-off $\ry$:
<<mtram-neck_ocm_ci>>=
exp(cbind(coef(neck_ocm)[2:6], confint(neck_ocm)[2:6,]))
@
In contrast, there seems to be a very large treatment effect (at week 7, the
odds in the placebo group is $\Sexpr{round(exp(coef(neck_ocm)[4]))}$ times
larger than in the active group. This levels off after 12 weeks, but the
effect is still significant at the $5\%$ level.

The corresponding transformation model with a different parameterisation of the
transformation function $\h$ and a completely different optimisation
procedure for maximising the log-likelihood, can be estimated by
<<mtram-neck_Colr, cache = TRUE>>=
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)
logLik(neck_Colrmer)
@
Based on this model of form (M1), it is possible to derive the marginal 
distribution functions in the two groups, see Figure~\ref{fig:distr_pain}.

\begin{figure}[t]
<<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)
@
\caption{Neck pain: Marginal distribution functions of chronic neck pain
         evaluated at three different time points under placebo or active
         low-level laser therapy. \label{fig:distr_pain}}
\end{figure}

\section{Marginally Interpretable Weibull and Cox Models}

The CAO/ARO/AIO-04 randomised clinical trial
\citep{Roedel_Graeven_Fietkau_2015} compared Oxaliplatin added to
fluorouracil-based preoperative chemoradiotherapy and postoperative
chemotherapy to the same therapy using fluorouracil only for rectal cancer
patients.  Patients were
randomised in the two treatment arms by block randomisation taking the study
center, the lymph node involvement (negative, positive), and tumour grading
(T1-3 vs.~T4) into account.  The primary endpoint was disease-free survival,
defined as the time between randomisation and non-radical surgery of the
primary tumour (R2 resection), locoregional recurrence after R0/1 resection,
metastatic disease or progression, or death from any cause, whichever
occurred first. The observed outcomes are a mix of exact dates (time to
death or incomplete removal of the primary tumour), right-censoring (end of
follow-up or drop-out), and interval-censoring (local or distant
metastases). We are interested in a clustered Cox or Weibull model for 
interval-censored survival times. The survivor functions, estimated
separately for each of the four strata defined by lymph node involvement and
tumour grading, are given in Figure~\ref{fig:CAO}.

<<mtram-CAO, echo = FALSE>>=
dir <- system.file("rda", package = "TH.data")
load(file.path(dir, "Primary_endpoint_data.rda"))
@

\begin{figure}
<<mtram-CAO-plot, cache = TRUE, 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))
@
\caption{Rectal cancer: Distribution of disease-free surival times for two
         treatments in the four strata defined by lymph node involvement
         (negative or positive) and tumour grading (T1-3 or T4). \label{fig:CAO}}
\end{figure}

The implementation of mixed transformation models is currently not able to
deal with mixed exact and censored outcomes in the same cluster. We
therefore recode exact event times as being interval-censored by adding a
4-day window to each exact event time (variable \code{iDFS2}).

<<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
@

We start with the random intercept model
\begin{eqnarray*}
\Prob(\rY > \ry \mid \text{treatment}) = 
\exp\left(-\exp\left(\frac{\eparm_1 + \eparm_2 \log(\ry) - 
                     \eshiftparm_\text{5-FU + Ox}}{\sqrt{\gamma_1^2 + 1}}\right)\right)
\end{eqnarray*}
assuming a marginal Weibull model whose effects are scaled depending on the
variance $\gamma_1^2$ of a block-specific (interaction of lymph node involvement,
tumour grading, and study center) random intercept:
<<mtram-CAO_SR, cache = TRUE>>=
CAO_SR <- Survreg(iDFS2 ~ randarm, data = CAOsurv)
CAO_SR_mtram_s <- mtram(CAO_SR, ~ (1 | Block), data = CAOsurv,
                        standardise = TRUE, Hessian = TRUE)
logLik(CAO_SR_mtram_s)
(cf <- coef(CAO_SR_mtram_s))
(OR <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
@
We are, of course, interested in the marginal treatment effect, that is, the
odds ratio 
%
\begin{eqnarray*}
\exp(-\eshiftparm_\text{5-FU + Ox} / \sqrt{\gamma_1^2 + 1}).
\end{eqnarray*}
%
We simply sample from the joint normal distribution of the maximum likelihood estimators
$\hat{\eparm}_1, \hat{\eparm}_2, \hat{\eshiftparm}_\text{5-FU + Ox},
\hat{\gamma}_1$ and compute confidence intervals for the marginal treatment
effect $\Sexpr{round(OR, 2)}$ as
<<mtram-CAO-CI>>=
S <- solve(CAO_SR_mtram_s$Hessian)
sqrt(diag(S))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_mtram_s), 
                 sigma = S)
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
@

In a next step, we stratify with respect to lymph node involvement and tumour
grading: For each of the four strata, the parameters $\eparm_1$ and
$\eparm_2$ are estimated separately:
<<mtram-CAO_SR_2, cache = TRUE>>=
CAO_SR_2 <- Survreg(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv)
CAO_SR_2_mtram_s <- mtram(CAO_SR_2, ~ (1 | Block), data = CAOsurv,
                          standardise = TRUE, Hessian  = TRUE)
logLik(CAO_SR_2_mtram_s)
(cf <- coef(CAO_SR_2_mtram_s))
(OR_2 <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
@
The corresponding confidence interval for the marginal treatment effect is
then
<<mtram-CAO-CI-2, echo = FALSE>>=
H <- CAO_SR_2_mtram_s$Hessian
S <- solve(H + .01 * diag(nrow(H)))
sqrt(diag(S))
rbeta <- rmvnorm(10000, mean = coef(CAO_SR_2_mtram_s), 
                 sigma = S)
s <- rbeta[,ncol(rbeta)]
rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
@
%%% this takes too long for CRAN, so just print the text
%We now relax the Weibull assumption in the Cox model
%\begin{eqnarray*}
%\Prob(\rY > \ry \mid \text{treatment}) = 
%\exp\left(-\exp\left(\frac{\basisy(\log(\ry))^\top \parm + 
%                     \eshiftparm_\text{5-FU + Ox}}{\sqrt{\gamma_1^2 + 1}}\right)\right)
%\end{eqnarray*}
%(note the positive sign of the treatment effect).
%<<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_s <- mtram(CAO_Cox_2, ~ (1 | Block), data = CAOsurv, 
%                           standardise = TRUE, Hessian = TRUE)
%logLik(CAO_Cox_2_mtram_s)
%coef(CAO_Cox_2_mtram_s)
%@
%with confidence interval
%<<mtram-CAO-CI-3, echo = FALSE>>=
%H <- CAO_Cox_2_mtram_s$Hessian
%S <- solve(H + .1 * diag(nrow(H)))
%rbeta <- rmvnorm(10000, mean = coef(CAO_Cox_2_mtram_s), 
%                 sigma = S)
%s <- rbeta[,ncol(rbeta)]
%rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
%quantile(exp(rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
%@

%%% start of static output
We now relax the Weibull assumption in the Cox model
\begin{eqnarray*}
\Prob(\rY > \ry \mid \text{treatment}, i) = 
\exp\left(-\exp\left(\frac{\basisy(\log(\ry))^\top \parm + 
                     \eshiftparm_\text{5-FU + Ox}}{\sqrt{\gamma_1^2 + 1}}\right)\right)
\end{eqnarray*}
(note the positive sign of the treatment effect).
\begin{Schunk}
\begin{Sinput}
R> CAO_Cox_2 <- Coxph(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv, 
+                     support = c(1, 1700), log_first = TRUE, order = 4)
R> logLik(CAO_Cox_2)
\end{Sinput}
\begin{Soutput}
'log Lik.' -2021.875 (df=21)
\end{Soutput}
\begin{Sinput}
R> CAO_Cox_2_mtram_s <- mtram(CAO_Cox_2, ~ (1 | Block), data = CAOsurv, 
+                               standardise = TRUE, Hessian = TRUE)
R> logLik(CAO_Cox_2_mtram_s)
\end{Sinput}
\begin{Soutput}
'log Lik.' -2031.051 (df=22)
\end{Soutput}
\begin{Sinput}
R> coef(CAO_Cox_2_mtram_s)
\end{Sinput}
\begin{Soutput}
Bs1(iDFS2):strat_ncN0:strat_tcT1-3 Bs2(iDFS2):strat_ncN0:strat_tcT1-3 
                     -5.832796e+01                      -3.163149e+00 
Bs3(iDFS2):strat_ncN0:strat_tcT1-3 Bs4(iDFS2):strat_ncN0:strat_tcT1-3 
                     -3.161654e+00                      -2.181838e+00 
Bs5(iDFS2):strat_ncN0:strat_tcT1-3 Bs1(iDFS2):strat_ncN+:strat_tcT1-3 
                     -7.656357e-01                      -1.653595e+01 
Bs2(iDFS2):strat_ncN+:strat_tcT1-3 Bs3(iDFS2):strat_ncN+:strat_tcT1-3 
                     -8.143542e+00                      -2.071762e+00 
Bs4(iDFS2):strat_ncN+:strat_tcT1-3 Bs5(iDFS2):strat_ncN+:strat_tcT1-3 
                     -1.792342e+00                      -7.614334e-01 
  Bs1(iDFS2):strat_ncN0:strat_tcT4   Bs2(iDFS2):strat_ncN0:strat_tcT4 
                     -2.525002e+00                      -2.519727e+00 
  Bs3(iDFS2):strat_ncN0:strat_tcT4   Bs4(iDFS2):strat_ncN0:strat_tcT4 
                     -2.330486e+00                      -3.587793e-01 
  Bs5(iDFS2):strat_ncN0:strat_tcT4   Bs1(iDFS2):strat_ncN+:strat_tcT4 
                     -1.456236e-01                      -4.096546e+01 
  Bs2(iDFS2):strat_ncN+:strat_tcT4   Bs3(iDFS2):strat_ncN+:strat_tcT4 
                     -2.000965e+00                      -1.987891e+00 
  Bs4(iDFS2):strat_ncN+:strat_tcT4   Bs5(iDFS2):strat_ncN+:strat_tcT4 
                     -3.570476e-01                       1.919546e-14 
         randarm5-FU + Oxaliplatin                             gamma1 
                     -1.806504e-01                      -1.215448e-06 
\end{Soutput}
\end{Schunk}
with confidence interval
\begin{Schunk}
\begin{Soutput}
     2.5%       50%     97.5% 
0.6911201 0.8325653 1.0084011 
\end{Soutput}
\end{Schunk}
%%% end of static output

Because the estimated variance parameter $\gamma_1$ is not very large, we
would expect to see similar results in a conditional Cox model with normal
frailty term
<<mtram-CAO-coxme>>=
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)
@


\section{Assessment of Unexplained Variability}

<<mtram-CHFLS, echo = FALSE>>=
library("TH.data")
load(file.path(path.package(package="TH.data"), "rda", "CHFLS.rda"))

### choose neccessary variables
org <- chfls1[, c("REGION6", "ZJ05", "ZJ06", "A35", "ZJ07", "ZJ16M", "INCRM",
                  "JK01", "JK02", "JK20", "HY04", "HY07", "A02", "AGEGAPM", 
                  "A07M", "A14", "A21", "A22M", "A23", "AX16", "INCAM", "SEXNOW", "ZW04")]

names(org) <- c("Region",
                "Rgender",               ### gender of respondent
                "Rage",                  ### age of respondent
                "RagestartA",            ### age of respondent at beginning of relationship with partner A
                "Redu",                  ### education of respondent
                "RincomeM",              ### rounded monthly income of respondent
                "RincomeComp",           ### inputed monthly income of respondent
                "Rhealth",               ### health condition respondent
                "Rheight",               ### respondent's height
                "Rhappy",                ### respondent's happiness
                "Rmartial",              ### respondent's marital status
                "RhasA",                 ### R has current A partner
                "Agender",               ### gender of partner A
                "RAagegap",              ### age gap
                "RAstartage",            ### age at marriage
                "Aheight",               ### height of partner A
                "Aedu",                  ### education of partner A
                "AincomeM",              ### rounded partner A income
                "AincomeEst",            ### estimated partner A income
                "orgasm",                ### orgasm frequency
                "AincomeComp",           ### imputed partner A income
                "Rsexnow",               ### has sex last year
                "Rhomosexual")           ### R is homosexual

### duration of partnership 
org$RAduration <- org$Rage - org$RagestartA

### code missing values
org$AincomeM[org$AincomeM < 0] <- NA
org$RincomeM[org$RincomeM < 0] <- NA
org$Aheight[org$Aheight < 0] <- NA

olevels <- c("never", "rarely", "sometimes", "often", "always")
orgA <- subset(org, Rgender == "female" & Rhomosexual != "yes" & orgasm %in% olevels)

orgA$orgasm <- ordered(as.character(orgA$orgasm),
        levels = c("never", "rarely", "sometimes", "often", "always"))

orgA$Redu <- factor(as.character(orgA$Redu),
        levels = c("univ/grad", "j col", "up mid", "low mid", "primary", "no school"))
levels(orgA$Redu) <-  c("univ", "jcol", "upmid", "lowmid", "primary", "noschool")

orgA$Aedu <- factor(as.character(orgA$Aedu),
        levels = c("univ/grad", "j col", "up mid", "low mid", "primary", "no school"))

orgA$Rhappy <- factor(as.character(orgA$Rhappy),
        levels = c("v unhappy", "not too", "relatively", "very"))

orgA$Rhealth <- factor(as.character(orgA$Rhealth),
        levels = c("poor", "not good", "fair", "good", "excellent"))

orgA$Region <- factor(as.character(orgA$Region),
        levels = c("CentralW", "Northeast", "North", "InlandS", "CoastalE", "CoastalS"))

orgA$AincomeSD <- orgA$AincomeComp/sd(orgA$AincomeComp)
orgA$AheightSD <- orgA$Aheight/sd(orgA$Aheight)
orgA$RageSD <- orgA$Rage/sd(orgA$Rage)
orgA$edudiff <- as.numeric(orgA$Aedu) - as.numeric(orgA$Redu)
orgA$edudiffSD <- orgA$edudiff/sd(orgA$edudiff, na.rm=TRUE)
orgA$wealthdiff <- orgA$RincomeComp - orgA$AincomeComp
orgA$wealthdiffSD <- orgA$wealthdiff/sd(orgA$wealthdiff, na.rm=TRUE)
orgA$RAdurationSD <- orgA$RAduration/sd(orgA$RAduration, na.rm=TRUE)
orgAtmp <- orgA[, c("orgasm", "AincomeSD", "AheightSD", "RAdurationSD",
                 "RageSD", "edudiffSD", "wealthdiffSD", "Redu", "Rhealth",
                 "Rhappy", "Region")]
cc <- complete.cases(orgAtmp)
orgAcc <- subset(orgA, cc)
@

\cite{Pollet2009} reported on an association between partner wealth and
female self-reported orgasm frequency.  It was later
\citep{Herberich+Hothorn+Nettle+Pollet2010} pointed out that the finding was
due to an incorrectly implemented variable selection procedure based on a
proportional odds (cumulative logit) model for the ordinal variable
corresponding to the question ``When having sex with your current partner,
how often did you have orgasm?'' with possible answer categories $\ry_1 =
\text{Always}$, $\ry_2 = \text{Often}$, $\ry_3 = \text{Sometimes}$, $\ry_4 =
\text{Rarely}$, or $\ry_5 = \text{Never}$.  The full model explains the
conditional distribution of orgasm frequency by $\rx = $ partner income,
partner height, the duration of the relationship, the respondents age, the
difference between both partners regarding education and wealth, the
respondents education, health, happiness, and place of living (regions in
China) of the form
%
\begin{eqnarray*}
\Prob(\text{orgasm} \le \ry_k \mid \rx) = \text{expit}(\eparm_k + 
                                                       \rx^\top \shiftparm)
\end{eqnarray*}
%
for $i = 1, \dots, N = \Sexpr{nrow(orgAcc)}$ independent heterosexual 
couples. In this model, the regression coefficients $\shiftparm$ can be
interpreted as log-odds ratios and we question the appropriateness of this
model here by including a subject-specific random intercept with standard
deviation $\gamma_1$. This changes the model to
\begin{eqnarray*}
\Prob(\text{orgasm} \le \ry_k \mid \rx) = \Phi\left(
  \frac{\Phi^{-1}(\text{expit}(\eparm_k + \rx^\top
  \shiftparm))}{\sqrt{\gamma_1^2 + 1}}\right)
\end{eqnarray*}
A value of $\gamma_1$ close to zero corresponds to marginal distributions
very similar to the proportional odds model and, consequently, it is
appropriate to interpret $\shiftparm$ as log-odds ratios. Larger values of
$\gamma_1$ indicate a more variable distribution and thus the choice $\pZ =
\text{expit}$ might be questionable.

%%% this takes too long because we have 1531 independent clusters...
%We obtain
%<<mtram-CHFLS-Polr, cache = TRUE>>=
%CHFLS_Polr <- Polr(orgasm ~ AincomeSD + AheightSD + RAdurationSD + 
%                   RageSD + edudiffSD + wealthdiffSD + Redu + 
%                   Rhealth + Rhappy + Region, data = orgAcc)
%logLik(CHFLS_Polr)
%orgAcc$ID <- factor(1:nrow(orgAcc))
%CHFLS_mtram <- mtram(CHFLS_Polr, ~ (1 | ID), 
%                     data = orgAcc)
%logLik(CHFLS_mtram)
%coef(CHFLS_mtram)
%@
%and from $\hat{\gamma}_1 = \Sexpr{round(coef(CHFLS_mtram)["gamma1"], 3)}$ 
%can conclude that the proportional odds model is appropriate here.

%%% start of static output
We obtain
\begin{Schunk}
\begin{Sinput}
R> CHFLS_Polr <- Polr(orgasm ~ AincomeSD + AheightSD + RAdurationSD + 
+                     RageSD + edudiffSD + wealthdiffSD + Redu + 
+                     Rhealth + Rhappy + Region, data = orgAcc)
R> logLik(CHFLS_Polr)
\end{Sinput}
\begin{Soutput}
'log Lik.' -1852.615 (df=27)
\end{Soutput}
\begin{Sinput}
R> orgAcc$ID <- factor(1:nrow(orgAcc))
R> CHFLS_mtram <- mtram(CHFLS_Polr, ~ (1 | ID), 
+                       data = orgAcc)
R> logLik(CHFLS_mtram)
\end{Sinput}
\begin{Soutput}
'log Lik.' -1852.829 (df=28)
\end{Soutput}
\begin{Sinput}
R> coef(CHFLS_mtram)
\end{Sinput}
\begin{Soutput}
         orgasm1          orgasm2          orgasm3          orgasm4 
     -3.00585271      -1.38078485       1.20377067       3.14108781 
       AincomeSD        AheightSD     RAdurationSD           RageSD 
      0.02603981      -0.02217060       0.07529263      -0.35310380 
       edudiffSD     wealthdiffSD         Redujcol        Reduupmid 
     -0.17413324      -0.03428024       0.13682888       0.17668815 
      Redulowmid      Reduprimary     Redunoschool  Rhealthnot good 
     -0.42694677      -0.95912533      -1.80236615       1.36152834 
     Rhealthfair      Rhealthgood Rhealthexcellent    Rhappynot too 
      1.69707889       1.83806835       1.86386980       0.29162606 
Rhappyrelatively       Rhappyvery  RegionNortheast      RegionNorth 
      0.75856526       1.02374623       0.40651849       0.20369845 
   RegionInlandS   RegionCoastalE   RegionCoastalS           gamma1 
      0.49050565       0.20059228       0.58338272       0.01678467 
\end{Soutput}
\end{Schunk}
and from $\hat{\gamma}_1 = 0.017$ 
can conclude that the proportional odds model is appropriate here.
%%% end of static output


\bibliography{mlt,packages}


%<<mtram-mlt-coef, echo = FALSE, results = "hide">>= 
%### print coefs for regression tests
%objs <- ls()
%mltobj <- objs[grep("^mlt_", objs)]
%sapply(mltobj, function(m) eval(parse(text = paste("coef(", m, ")"))))
%@

<<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")
#pkg("lme4")
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")
@

\end{document}

