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

\documentclass[article,nojss,shortnames]{jss}\usepackage[]{graphicx}\usepackage[]{color}

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

\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}
}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}


\footnote{Please cite this document as: Torsten Hothorn (2022)
Some Applications of Marginally Interpretable Linear 
Transformation Models for Clustered Observations.
\textsf{R} package vignette version 0.7-0, 
URL \url{https://CRAN.R-project.org/package=tram}.
Running this vignette takes too long for CRAN; you can access the code by
\code{demo("mtram", 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
\begin{Schunk}
\begin{Sinput}
R> install.packages("tram")
R> vignette("mtram", package = "tram")
\end{Sinput}
\end{Schunk}

\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]
\begin{Schunk}


{\centering \includegraphics{mtram-sleep-plot-1} 

}

\end{Schunk}
\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}:
%
\begin{Schunk}
\begin{Sinput}
R> library("lme4")
R> sleep_lmer <- lmer(Reaction ~ Days + (Days | Subject), 
+                     data = sleepstudy, REML = FALSE)
\end{Sinput}
\end{Schunk}
%
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
\begin{Schunk}
\begin{Sinput}
R> sleep_LM <- Lm(Reaction ~ Days, data = sleepstudy)
R> sleep_LMmer <- mtram(sleep_LM, ~ (Days | Subject), data = sleepstudy)
\end{Sinput}
\end{Schunk}
%
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
\begin{Schunk}
\begin{Sinput}
R> logLik(sleep_lmer)
\end{Sinput}
\begin{Soutput}
'log Lik.' -875.9697 (df=6)
\end{Soutput}
\begin{Sinput}
R> logLik(sleep_LMmer)
\end{Sinput}
\begin{Soutput}
'log Lik.' -875.9697 (df=6)
\end{Soutput}
\end{Schunk}
Looking at the model coefficients, the two procedures lead to almost
identical inverse residual standard deviations
\begin{Schunk}
\begin{Sinput}
R> (sdinv <- 1 / summary(sleep_lmer)$sigma)
\end{Sinput}
\begin{Soutput}
[1] 0.03907485
\end{Soutput}
\begin{Sinput}
R> coef(sleep_LMmer)["Reaction"]
\end{Sinput}
\begin{Soutput}
  Reaction 
0.03907741 
\end{Soutput}
\end{Schunk}
and fixed effects (the slope can be interpreted as inverse coefficient of
variation)
\begin{Schunk}
\begin{Sinput}
R> fixef(sleep_lmer) * c(-1, 1) * sdinv
\end{Sinput}
\begin{Soutput}
(Intercept)        Days 
 -9.8236175   0.4090077 
\end{Soutput}
\begin{Sinput}
R> coef(sleep_LMmer)[c("(Intercept)", "Days")]
\end{Sinput}
\begin{Soutput}
(Intercept)        Days 
 -9.8243917   0.4089265 
\end{Soutput}
\end{Schunk}
The random-effect parameters $\varparm$ are also reasonably close
\begin{Schunk}
\begin{Sinput}
R> sleep_lmer@theta
\end{Sinput}
\begin{Soutput}
[1] 0.92919061 0.01816575 0.22264321
\end{Soutput}
\begin{Sinput}
R> coef(sleep_LMmer)[-(1:3)]
\end{Sinput}
\begin{Soutput}
    gamma1     gamma2     gamma3 
0.92901066 0.01843056 0.22280431 
\end{Soutput}
\end{Schunk}
Consequently, the variance-covariance and correlation matrices
\begin{Schunk}
\begin{Sinput}
R> sleep_LMmer$G * (1 / sdinv)^2
\end{Sinput}
\begin{Soutput}
2 x 2 sparse Matrix of class "dsCMatrix"
                      
[1,] 565.2580 11.21410
[2,]  11.2141 32.73513
\end{Soutput}
\begin{Sinput}
R> cov2cor(sleep_LMmer$G * (1 / sdinv)^2)
\end{Sinput}
\begin{Soutput}
2 x 2 sparse Matrix of class "dsCMatrix"
                          
[1,] 1.00000000 0.08243925
[2,] 0.08243925 1.00000000
\end{Soutput}
\begin{Sinput}
R> unclass(VarCorr(sleep_lmer))$Subject
\end{Sinput}
\begin{Soutput}
            (Intercept)     Days
(Intercept)   565.47697 11.05512
Days           11.05512 32.68179
attr(,"stddev")
(Intercept)        Days 
  23.779760    5.716799 
attr(,"correlation")
            (Intercept)       Days
(Intercept)  1.00000000 0.08132109
Days         0.08132109 1.00000000
\end{Soutput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> library("survival")
R> sleepstudy$Reaction_I <- with(sleepstudy, Surv(Reaction - 20, Reaction + 20, 
+                                                 type = "interval2"))
R> sleepstudy$Reaction_I[1:5]
\end{Sinput}
\begin{Soutput}
[1] [229.5600, 269.5600] [238.7047, 278.7047] [230.8006, 270.8006]
[4] [301.4398, 341.4398] [336.8519, 376.8519]
\end{Soutput}
\end{Schunk}
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:
\begin{Schunk}
\begin{Sinput}
R> sleep_LM_I <- Lm(Reaction_I ~ Days, data = sleepstudy)
R> sleep_LMmer_I <- mtram(sleep_LM_I, ~ (Days | Subject), data = sleepstudy)
\end{Sinput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> logLik(sleep_LMmer_I)
\end{Sinput}
\begin{Soutput}
'log Lik.' -214.9675 (df=6)
\end{Soutput}
\begin{Sinput}
R> coef(sleep_LMmer_I)
\end{Sinput}
\begin{Soutput}
(Intercept)  Reaction_I        Days      gamma1      gamma2      gamma3 
-9.78770607  0.03900116  0.41633415  0.83398952  0.07584130  0.19038611 
\end{Soutput}
\begin{Sinput}
R> coef(sleep_LMmer)
\end{Sinput}
\begin{Soutput}
(Intercept)    Reaction        Days      gamma1      gamma2      gamma3 
-9.82439168  0.03907741  0.40892652  0.92901066  0.01843056  0.22280431 
\end{Soutput}
\end{Schunk}

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
\begin{Schunk}
\begin{Sinput}
R> sleep_BC <- BoxCox(Reaction ~ Days, data = sleepstudy)
R> sleep_BCmer <- mtram(sleep_BC, ~ (Days | Subject), data = sleepstudy, 
+                       Hessian = TRUE)
R> logLik(sleep_BCmer)
\end{Sinput}
\begin{Soutput}
'log Lik.' -859.5455 (df=11)
\end{Soutput}
\end{Schunk}
%
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]
\begin{Schunk}


{\centering \includegraphics{mtram-sleep_BoxCoxPlot-1} 

}

\end{Schunk}
\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]
\begin{Schunk}


{\centering \includegraphics{mtram-sleep_marginal-1} 

}

\end{Schunk}
\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 
\begin{Schunk}
\begin{Sinput}
R> cov2cor(sleep_BCmer$G)
\end{Sinput}
\begin{Soutput}
2 x 2 sparse Matrix of class "dsCMatrix"
                          
[1,]  1.0000000 -0.1946629
[2,] -0.1946629  1.0000000
\end{Soutput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> VC <- solve(sleep_BCmer$Hessian)
R> idx <- (nrow(VC) - 2):nrow(VC)
R> Rcoef <- rmvnorm(1000, mean = coef(sleep_BCmer), sigma = VC)[,idx]
R> 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])
+  })
\end{Sinput}
\end{Schunk}
The $95\%$ confidence intervals
\begin{Schunk}
\begin{Sinput}
R> ### variance random intercept
R> quantile(ret[1,], c(.025, .5, .975))
\end{Sinput}
\begin{Soutput}
     2.5%       50%     97.5% 
0.9127821 2.5713595 5.2493469 
\end{Soutput}
\begin{Sinput}
R> ### variance random slope
R> quantile(ret[2,], c(.025, .5, .975))
\end{Sinput}
\begin{Soutput}
      2.5%        50%      97.5% 
0.01890987 0.05348231 0.10594879 
\end{Soutput}
\begin{Sinput}
R> ### correlation random intercept / random slope
R> quantile(ret[3,], c(.025, .5, .975))
\end{Sinput}
\begin{Soutput}
      2.5%        50%      97.5% 
-0.6193527 -0.1883314  0.4689778 
\end{Soutput}
\end{Schunk}
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]
\begin{Schunk}


{\centering \includegraphics{mtram-toenail-plot-1} 

}

\end{Schunk}
\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:
\begin{Schunk}
\begin{Sinput}
R> ### Laplace
R> toenail_glmer_RI_1 <- 
+      glmer(outcome ~ treatment * time + (1 | patientID),
+            data = toenail, family = binomial(link = "probit"), 
+            nAGQ = 1)
R> summary(toenail_glmer_RI_1)
\end{Sinput}
\begin{Soutput}
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( probit )
Formula: outcome ~ treatment * time + (1 | patientID)
   Data: toenail

     AIC      BIC   logLik deviance df.resid 
  1279.0   1306.8   -634.5   1269.0     1898 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.507 -0.017 -0.004  0.000 54.046 

Random effects:
 Groups    Name        Variance Std.Dev.
 patientID (Intercept) 20.68    4.548   
Number of obs: 1903, groups:  patientID, 289

Fixed effects:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -3.39650    0.22091 -15.375   <2e-16 ***
treatmentterbinafine      -0.01532    0.25359  -0.060   0.9518    
time                      -0.21749    0.02256  -9.639   <2e-16 ***
treatmentterbinafine:time -0.07155    0.03425  -2.089   0.0367 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) trtmnt time  
trtmnttrbnf -0.593              
time        -0.009  0.102       
trtmnttrbn:  0.093 -0.143 -0.629
\end{Soutput}
\begin{Sinput}
R> toenail_glmer_RI_1@theta
\end{Sinput}
\begin{Soutput}
[1] 4.547891
\end{Soutput}
\begin{Sinput}
R> ### Adaptive Gaussian Quadrature
R> toenail_glmer_RI_2 <- 
+      glmer(outcome ~ treatment * time + (1 | patientID),
+            data = toenail, family = binomial(link = "probit"), 
+            nAGQ = 20)
R> summary(toenail_glmer_RI_2)
\end{Sinput}
\begin{Soutput}
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( probit )
Formula: outcome ~ treatment * time + (1 | patientID)
   Data: toenail

     AIC      BIC   logLik deviance df.resid 
  1277.8   1305.6   -633.9   1267.8     1898 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.847 -0.189 -0.078 -0.001 33.997 

Random effects:
 Groups    Name        Variance Std.Dev.
 patientID (Intercept) 4.485    2.118   
Number of obs: 1903, groups:  patientID, 289

Fixed effects:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -0.93061    0.23176  -4.015 5.93e-05 ***
treatmentterbinafine      -0.07609    0.30921  -0.246   0.8056    
time                      -0.19074    0.02059  -9.263  < 2e-16 ***
treatmentterbinafine:time -0.06419    0.03099  -2.071   0.0383 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) trtmnt time  
trtmnttrbnf -0.655              
time        -0.186  0.212       
trtmnttrbn:  0.193 -0.287 -0.611
\end{Soutput}
\begin{Sinput}
R> toenail_glmer_RI_2@theta
\end{Sinput}
\begin{Soutput}
[1] 2.117846
\end{Soutput}
\end{Schunk}

%%%% 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)
\begin{Schunk}
\begin{Sinput}
R> m <- ctm(as.basis(~ outcome, data = toenail), 
+           shifting = ~ treatment * time, 
+           data = toenail, todistr = "Normal")
R> toenail_probit <- mlt(m, data = toenail, 
+                        fixed = c("outcomemoderate or severe" = 0))
R> toenail_mtram_RI <- 
+      mtram(toenail_probit, ~ (1 | patientID), 
+            data = toenail, Hessian = TRUE)
R> logLik(toenail_mtram_RI)
\end{Sinput}
\begin{Soutput}
'log Lik.' -633.9638 (df=5)
\end{Soutput}
\begin{Sinput}
R> coef(toenail_mtram_RI)
\end{Sinput}
\begin{Soutput}
              (Intercept)      treatmentterbinafine 
               0.92947317                0.07699415 
                     time treatmentterbinafine:time 
               0.19056726                0.06355500 
                   gamma1 
               2.11448400 
\end{Soutput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> vcov(toenail_glmer_RI_2)
\end{Sinput}
\begin{Soutput}
4 x 4 Matrix of class "dpoMatrix"
                            (Intercept) treatmentterbinafine          time
(Intercept)                0.0537124797         -0.046953444 -0.0008877052
treatmentterbinafine      -0.0469534440          0.095609463  0.0013522344
time                      -0.0008877052          0.001352234  0.0004239968
treatmentterbinafine:time  0.0013871416         -0.002754672 -0.0003896726
                          treatmentterbinafine:time
(Intercept)                            0.0013871416
treatmentterbinafine                  -0.0027546716
time                                  -0.0003896726
treatmentterbinafine:time              0.0009603451
\end{Soutput}
\begin{Sinput}
R> solve(toenail_mtram_RI$Hessian)[1:4, 1:4]
\end{Sinput}
\begin{Soutput}
              [,1]         [,2]          [,3]          [,4]
[1,]  0.0535097431 -0.046828788 -0.0008863008  0.0013728404
[2,] -0.0468287884  0.095401757  0.0013454564 -0.0027201727
[3,] -0.0008863008  0.001345456  0.0004223674 -0.0003889295
[4,]  0.0013728404 -0.002720173 -0.0003889295  0.0009479493
\end{Soutput}
\end{Schunk}
%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}
\begin{Schunk}
\begin{Sinput}
R> toenail_glmer_RS <- 
+      glmer(outcome ~ treatment * time + (1 + time | patientID),
+            data = toenail, family = binomial(link = "probit"))
R> summary(toenail_glmer_RS)
\end{Sinput}
\begin{Soutput}
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( probit )
Formula: outcome ~ treatment * time + (1 + time | patientID)
   Data: toenail

     AIC      BIC   logLik deviance df.resid 
   985.8   1024.7   -485.9    971.8     1896 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.85421 -0.00210 -0.00037  0.00000  2.35828 

Random effects:
 Groups    Name        Variance Std.Dev. Corr 
 patientID (Intercept) 118.433  10.883        
           time          3.305   1.818   -0.90
Number of obs: 1903, groups:  patientID, 289

Fixed effects:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -4.30119    0.26361 -16.316   <2e-16 ***
treatmentterbinafine       0.05419    0.34652   0.156   0.8757    
time                      -0.06792    0.07847  -0.866   0.3867    
treatmentterbinafine:time -0.23478    0.13885  -1.691   0.0909 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) trtmnt time  
trtmnttrbnf -0.662              
time        -0.453  0.342       
trtmnttrbn:  0.270 -0.438 -0.335
\end{Soutput}
\begin{Sinput}
R> toenail_glmer_RS@theta
\end{Sinput}
\begin{Soutput}
[1] 10.8826820 -1.6359584  0.7930834
\end{Soutput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> toenail_mtram_RS <- 
+      mtram(toenail_probit, ~ (1 + time | patientID), 
+            data = toenail)
R> logLik(toenail_mtram_RS)
\end{Sinput}
\begin{Soutput}
'log Lik.' -545.1164 (df=7)
\end{Soutput}
\begin{Sinput}
R> coef(toenail_mtram_RS)
\end{Sinput}
\begin{Soutput}
              (Intercept)      treatmentterbinafine 
                1.5765323                -0.2666843 
                     time treatmentterbinafine:time 
                0.5323919                 0.1842506 
                   gamma1                    gamma2 
                5.2172371                -0.3723897 
                   gamma3 
                0.5285640 
\end{Soutput}
\end{Schunk}
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}:
\begin{Schunk}
\begin{Sinput}
R> m <- ctm(as.basis(~ outcome, data = toenail), 
+           shifting = ~ treatment * time, 
+           data = toenail, todistr = "Logistic")
R> toenail_logit <- mlt(m, data = toenail, 
+                       fixed = c("outcomemoderate or severe" = 0))
R> toenail_mtram_logit <- mtram(toenail_logit, ~ (1 | patientID), 
+                                data = toenail)
R> toenail_mtram_logit_s <- mtram(toenail_logit, ~ (1 | patientID), 
+                                 data = toenail, standardise = TRUE, 
+                                 Hessian = TRUE)
\end{Sinput}
\end{Schunk}
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
\begin{Schunk}
\begin{Sinput}
R> tmp <- toenail_logit
R> cf <- coef(tmp)
R> cf <- cf[names(cf) != "outcomemoderate or severe"]
R> sdrf <- rev(coef(toenail_mtram_logit_s))[1]
R> cf <- coef(toenail_mtram_logit_s)[names(cf)] / sqrt(sdrf^2 + 1)
R> cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
R> coef(tmp) <- cf
R> time <- 0:180/10
R> treatment <- sort(unique(toenail$treatment))
R> nd <- expand.grid(time = time, treatment = treatment)
R> nd$prob_logit_s <- predict(tmp, newdata = nd, type = "distribution")[1,]
R> nd$odds <- exp(predict(tmp, newdata = nd, type = "trafo")[1,])
\end{Sinput}
\end{Schunk}

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]
\begin{Schunk}


{\centering \includegraphics{mtram-toenail_OR_2-1} 

}

\end{Schunk}
\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)
\begin{Schunk}
\begin{Sinput}
R> tmp <- toenail_logit
R> cf <- coef(tmp)
R> cf <- cf[names(cf) != "outcomemoderate or severe"]
R> sdrf <- rev(coef(toenail_mtram_logit))[1]
R> cf <- coef(toenail_mtram_logit)[names(cf)] 
R> cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
R> coef(tmp) <- cf
R> pr <- predict(tmp, newdata = nd, type = "distribution")[1,]
R> nd$prob_logit <- pnorm(qnorm(pr) / sdrf)
\end{Sinput}
\end{Schunk}
and (logit)
\begin{Schunk}
\begin{Sinput}
R> tmp <- toenail_probit
R> cf <- coef(tmp)
R> cf <- cf[names(cf) != "outcomemoderate or severe"]
R> sdrf <- rev(coef(toenail_mtram_RI))[1]
R> cf <- coef(toenail_mtram_RI)[names(cf)] / sqrt(sdrf^2 + 1)
R> cf <- c(cf[1], "outcomemoderate or severe" = 0, cf[-1])
R> coef(tmp) <- cf
R> nd$prob_probit <- predict(tmp, newdata = nd, type = "distribution")[1,]
\end{Sinput}
\end{Schunk}
The marginal time-dependent probabilities obtained from all three models are
very similar as shown in Figure~\ref{fig:toenailprob}.

\begin{figure}[t]
\begin{Schunk}


{\centering \includegraphics{mtram-toenail_probplot-1} 

}

\end{Schunk}
\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]
\begin{Schunk}


{\centering \includegraphics{mtram-neck_plot-1} 

}

\end{Schunk}
\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*}

\begin{Schunk}
\begin{Sinput}
R> library("ordinalCont")
R> neck_ocm <- ocm(vas ~ laser * time + (1 | id), data = pain_df, 
+                  scale = c(0, 1))
\end{Sinput}
\end{Schunk}
The results
\begin{Schunk}
\begin{Sinput}
R> summary(neck_ocm)
\end{Sinput}
\begin{Soutput}
Call:
ocm(formula = vas ~ laser * time + (1 | id), data = pain_df, 
    scale = c(0, 1))

Random effects:
         Name Variance Std.Dev.
 Intercept|id    5.755    2.399

Coefficients:
                         Estimate   StdErr t.value   p.value    
laserActive              -2.07922  0.65055 -3.1961  0.001918 ** 
time7 weeks              -0.60366  0.35744 -1.6889  0.094689 .  
time12 weeks             -0.23804  0.36365 -0.6546  0.514395    
laserActive:time7 weeks   4.40817  0.56073  7.8615 7.604e-12 ***
laserActive:time12 weeks  3.38593  0.53925  6.2790 1.159e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{Soutput}
\end{Schunk}
suggest that there is a difference at baseline; the pain distribution of 
subjects in the placebo group on the odds scale is only 
$13\%$ of the odds in the active 
group for any cut-off $\ry$:
\begin{Schunk}
\begin{Sinput}
R> exp(cbind(coef(neck_ocm)[2:6], confint(neck_ocm)[2:6,]))
\end{Sinput}
\begin{Soutput}
                                          2.5 %      97.5 %
laserActive               0.1250278  0.03493482   0.4474608
time7 weeks               0.5468040  0.27137954   1.1017581
time12 weeks              0.7881704  0.38643700   1.6075391
laserActive:time7 weeks  82.1194073 27.36208405 246.4577275
laserActive:time12 weeks 29.5454666 10.26785879  85.0162253
\end{Soutput}
\end{Schunk}
In contrast, there seems to be a very large treatment effect (at week 7, the
odds in the placebo group is $1$ 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
\begin{Schunk}
\begin{Sinput}
R> neck_Colr <- Colr(vas ~ laser * time, data = pain_df, 
+                    bounds = c(0, 1), support = c(0, 1),
+                    extrapolate = TRUE)
R> neck_Colrmer <- mtram(neck_Colr, ~ (1 | id), data = pain_df, 
+                        Hessian = TRUE)
R> logLik(neck_Colrmer)
\end{Sinput}
\begin{Soutput}
'log Lik.' 74.35908 (df=13)
\end{Soutput}
\end{Schunk}
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]
\begin{Schunk}


{\centering \includegraphics{mtram-neck_Colr_distr-1} 

}

\end{Schunk}
\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}.



\begin{figure}
\begin{Schunk}


{\centering \includegraphics{mtram-CAO-plot-1} 

}

\end{Schunk}
\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}).

\begin{Schunk}
\begin{Sinput}
R> ### convert "exact" event dates to interval-censoring (+/- one day)
R> tmp <- CAOsurv$iDFS
R> exact <- tmp[,3] == 1 
R> tmp[exact,2] <- tmp[exact,1] + 2
R> tmp[exact,1] <- pmax(tmp[exact,1] - 2, 0)
R> tmp[exact,3] <- 3
R> CAOsurv$iDFS2 <- tmp
\end{Sinput}
\end{Schunk}

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:
\begin{Schunk}
\begin{Sinput}
R> CAO_SR <- Survreg(iDFS2 ~ randarm, data = CAOsurv)
R> CAO_SR_mtram_s <- mtram(CAO_SR, ~ (1 | Block), data = CAOsurv,
+                          standardise = TRUE, Hessian = TRUE)
R> logLik(CAO_SR_mtram_s)
\end{Sinput}
\begin{Soutput}
'log Lik.' -2081.542 (df=4)
\end{Soutput}
\begin{Sinput}
R> (cf <- coef(CAO_SR_mtram_s))
\end{Sinput}
\begin{Soutput}
              (Intercept)                log(iDFS2) 
               -6.2990054                 0.7412855 
randarm5-FU + Oxaliplatin                    gamma1 
                0.2328600                 0.1683613 
\end{Soutput}
\begin{Sinput}
R> (OR <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
\end{Sinput}
\begin{Soutput}
randarm5-FU + Oxaliplatin 
                 0.794829 
\end{Soutput}
\end{Schunk}
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 $0.79$ as
\begin{Schunk}
\begin{Sinput}
R> S <- solve(CAO_SR_mtram_s$Hessian)
R> sqrt(diag(S))
\end{Sinput}
\begin{Soutput}
[1] 0.29019989 0.03872268 0.10722445 0.12433065
\end{Soutput}
\begin{Sinput}
R> rbeta <- rmvnorm(10000, mean = coef(CAO_SR_mtram_s), 
+                   sigma = S)
R> s <- rbeta[,ncol(rbeta)]
R> rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1)
R> quantile(exp(-rbeta[, ncol(rbeta)]), prob = c(.025, .5, .975))
\end{Sinput}
\begin{Soutput}
     2.5%       50%     97.5% 
0.6478684 0.7953985 0.9747696 
\end{Soutput}
\end{Schunk}

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:
\begin{Schunk}
\begin{Sinput}
R> CAO_SR_2 <- Survreg(iDFS2 | 0 + strat_n:strat_t ~ randarm, data = CAOsurv)
R> CAO_SR_2_mtram_s <- mtram(CAO_SR_2, ~ (1 | Block), data = CAOsurv,
+                            standardise = TRUE, Hessian  = TRUE)
R> logLik(CAO_SR_2_mtram_s)
\end{Sinput}
\begin{Soutput}
'log Lik.' -2067.797 (df=10)
\end{Soutput}
\begin{Sinput}
R> (cf <- coef(CAO_SR_2_mtram_s))
\end{Sinput}
\begin{Soutput}
(Intercept):strat_ncN0:strat_tcT1-3  log(iDFS2):strat_ncN0:strat_tcT1-3 
                         -7.8833653                           0.9584499 
(Intercept):strat_ncN+:strat_tcT1-3  log(iDFS2):strat_ncN+:strat_tcT1-3 
                         -6.2225174                           0.7198965 
  (Intercept):strat_ncN0:strat_tcT4    log(iDFS2):strat_ncN0:strat_tcT4 
                         -3.0467542                           0.3711277 
  (Intercept):strat_ncN+:strat_tcT4    log(iDFS2):strat_ncN+:strat_tcT4 
                         -4.8207089                           0.6214653 
          randarm5-FU + Oxaliplatin                              gamma1 
                          0.2240023                           0.1474685 
\end{Soutput}
\begin{Sinput}
R> (OR_2 <- exp(-cf["randarm5-FU + Oxaliplatin"] / sqrt(cf["gamma1"]^2 + 1)))
\end{Sinput}
\begin{Soutput}
randarm5-FU + Oxaliplatin 
                0.8012313 
\end{Soutput}
\end{Schunk}
The corresponding confidence interval for the marginal treatment effect is
then
\begin{Schunk}
\begin{Soutput}
 [1] 0.68882415 0.09402924 0.34618756 0.04633063 1.01843933 0.13874107
 [7] 0.68657931 0.09534206 0.10731700 0.13455556
\end{Soutput}
\begin{Soutput}
     2.5%       50%     97.5% 
0.6528915 0.8044045 0.9856882 
\end{Soutput}
\end{Schunk}
%%% 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
\begin{Schunk}
\begin{Sinput}
R> library("coxme")
R> m <- coxme(DFS ~ randarm + (1 | Block), data = CAOsurv)
R> summary(m)
\end{Sinput}
\begin{Soutput}
Cox mixed-effects model fit by maximum likelihood
  Data: CAOsurv
  events, n = 357, 1236
  Iterations= 22 91 
                    NULL Integrated    Fitted
Log-likelihood -2432.971  -2430.475 -2414.482

                  Chisq    df         p  AIC    BIC
Integrated loglik  4.99  2.00 0.0823850 0.99  -6.76
 Penalized loglik 36.98 16.72 0.0029833 3.53 -61.33

Model:  DFS ~ randarm + (1 | Block) 
Fixed coefficients
                                coef exp(coef)  se(coef)     z    p
randarm5-FU + Oxaliplatin -0.2310483 0.7937011 0.1067215 -2.16 0.03

Random effects
 Group Variable  Std Dev    Variance  
 Block Intercept 0.21741797 0.04727057
\end{Soutput}
\begin{Sinput}
R> sd <- sqrt(diag(vcov(m)))
R> exp(coef(m) + c(-1, 0, 1) * qnorm(.975) * sd)
\end{Sinput}
\begin{Soutput}
[1] 0.6438957 0.7937011 0.9783596
\end{Soutput}
\end{Schunk}


\section{Assessment of Unexplained Variability}



\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 = 1531$ 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 = 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-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}

