A new, development version, of package <a target="_blank" href="https://cran.r-project.org/package=JMbayes">JMbayes</a> has been rolled-out on the dedicated <a target="_blank" href="https://github.com/drizopoulos/JMbayes">GitHub repo</a>. The major addition in this version is a set of new functions that can fit multivariate joint models for multiple longitudinal outcomes and a time-to-event. As with other GitHub packages, it can be easily installed using the <span><span>install_github()</span></span> function from package <a target="_blank" href="https://cran.r-project.org/package=devtools">devtools</a>. The new function that fits multivariate joint models is called <span><span>mvJointModelBayes() </span></span>and has a very similar syntax as the <span><span>jointModelBayes()</span></span> function described in a <a target="_blank" href="http://iprogn.blogspot.nl/2014_09_01_archive.html">previous post</a>. However, contrary to <span><span>jointModelBayes() </span></span>that is entirely written in R, the main bulk of computations of <span><span>mvJointModelBayes() </span></span>are based on C++ code building upon the excellent <a target="_blank" href="https://cran.r-project.org/package=Rcpp">Rcpp</a> and <a target="_blank" href="https://cran.r-project.org/package=RcppArmadillo">RcppArmadillo</a> packages.<br><br>To make the connection between the two functions from a practical viewpoint more clear, let's fit first a univariate joint model using <span><span>mvJointModelBayes()</span></span>. As for <span><span>jointModelBayes()</span></span>, we first need to separately fit a mixed-effects model for the longitudinal outcome, and a Cox model for the survival outcome. However, contrary to <span><span>jointModelBayes()</span></span>, which required the mixed model to be fitted with function <span><span>lme()</span></span> from the <a target="_blank" href="https://cran.r-project.org/package=nlme">nlme</a> package, for <span><span>mvJointModelBayes()</span></span> the mixed model needs to be fitted with the new function <span><span>mvlgmer()</span></span>. This follows the syntax of <span><span>lmer()</span></span> from package <a target="_blank" href="https://cran.r-project.org/package=lme4">lme4</a>. For example, the code below fits a linear mixed model for the longitudinal biomarker serum bilirubin from the Mayo Clinic PBC data set:<br><a target="_blank" href="https://github.com/drizopoulos/JMbayes"> </a><br><a target="_blank" href="https://github.com/drizopoulos/JMbayes"> </a><br><br>The main arguments of this function are: <span><span>formulas</span></span> a list of lme4-like formulas (a formula per outcome), <span><span>data</span></span> a data.frame that contains all the variables specified in <span><span>formulas</span></span> (NAs allowed), and <span><span>families</span></span> a list of family objects specifying the type of each outcome (currently the gaussian, binomial and poisson families are allowed). Following, we fit a Cox model using function <span><span>coxph()</span></span> from the <a target="_blank" href="https://cran.r-project.org/package=survival">survival</a> package; argument <span><span>model</span></span> of <span><span>coxph()</span></span> needs to be set to <span><span>TRUE</span></span>:<br><br><br>Then the univariate joint model is fitted by the following simple call to <span><span>mvJointModelBayes()</span></span>:<br><br><br>To fit a multivariate joint model for multiple longitudinal outcomes, the only thing that needs to be adapted from the procedure described above, is the call to <span><span>mvglmer()</span></span>. For example, to fit a bivariate joint model for serum bilirubin and presence of spiders, we use the syntax:<br><br><br>In addition, similarly to <span><span>jointModelBayes()</span></span>, <span><span>mvJointModelBayes() </span></span>allows to specify different functional forms for the longitudinal outcomes that are included in the Cox model.  As an example, we extend model <span><span>JMFit2</span></span> by including the value and slope term for bilirubin, and the area term for spiders (in the log-odds scale). To include these terms into the model we specify the <span><span>Formulas</span></span> argument. This is specified in a similar manner as the <span><span>derivForms</span></span> argument of <span><span>jointModelBayes()</span></span>. In particular, it should be a list of lists. Each component of the outer list should have as name the name of the corresponding outcome variable. Then in the inner list we need to specify four components, namely, the <span><span>fixed</span></span> &amp; <span><span>random</span></span> R formulas specifying the fixed and random effects part of the term to be included, and <span><span>indFixed</span></span> &amp; <span><span>indRandom</span></span> integer indices specifying which of the original fixed and random effects are involved in the calculation of the new term. In the inner list you can also optionally specify a name for the term you want to include. A couple of notes:<br><ol><li>For terms not specified in the <span><span>Formulas</span></span> list, the default value functional form is used. </li><li>If for a particular outcome you want to include both the value functional form and an extra term, then you need to specify that in the <span><span>Formulas</span></span> using two entries. To include the value functional form you only need to set the corresponding to <span><span>&quot;value&quot;</span></span>, and for the second term to specify the inner list. See example below on how to include the value and slope for serum bilirubin (for example, if the list below the entry <span><span>&quot;log(serBilir)&quot; = &quot;value&quot;</span></span> was not given, then only the slope term would have been included in the survival submodel).</li></ol><br>We further extend the previous model by including the interactions terms between the terms specified in <span><span>Formulas</span></span> and the randomized treatment indicator <span><span>drug</span></span>. The names specified in the list that defined the interaction factors should match with the names of the output from <span><span>JMFit3</span>. </span>Because of the many association parameters we have, we place a shrinkage prior on the association coefficients 'alpha'. In particular, if we have K association parameters, we assume that alpha_k ~ N(0, tau * phi_k), k = 1, ..., K. The precision parameters tau and phi_k are given Gamma priors. Parameter tau is global shrinkage parameter, and phi_k a specific shrinkage parameter per alpha coefficient:<br><br><br><br>Future updates of the package will include additional features, including among others<br><ul><li>a method for the generic functions <span><span>survfitJM()</span></span> and <span><span>predict()</span></span> to compute dynamic predictions for multivariate longitudinal data setting for both the longitudinal and survival outcomes.</li><li>allow for competing risks.</li><li>make the C++ code (even) faster.</li></ul>