# Introduction of mlm4omics

Irene SL Zeng

#### 15 May 2019

#### Package

mlm4omics 1.3.0

# Contents

    library(knitr)
    library(rmarkdown)

    ##
    ## Attaching package: 'rmarkdown'

    ## The following objects are masked from 'package:BiocStyle':
    ##
    ##     html_document, md_document, pdf_document

    knitr::opts_chunk$set(echo = TRUE)

# 1 Introduction

## 1.1 Summary

High throughput mass-spectrometry clinical proteomics studies curate abundance data with hierarchical structure. The non-random missingness of the measurements from a vast amount of information also adds complexity in data analysis. We propose a multivariate 2 level model to analyze protein abundances and to handle abundance-dependent Missingness within a Bayesian framework. The proposed model enables the variance decomposition at different levels of the data hierarchy and provides central shrinkage of protein-level estimates for proteins with small numbers of observations. A logistic regression model with informative prior of the regression coefficient is used to model the missing probability. The non-random missing and censored values are treated as unknown parameters in a Bayesian model. Hamiltonian MC/No-U-Turn Sampling algorithm is created to derive the posterior distribution of the unknown parameters. Hamiltonian MC is demonstrated to gain more efficiency for these high-dimensional correlated data.

## 1.2 The proposed model

We use an inference model within a multivariate-multilevel framework to analyze the protein abundance data that are curated in clinical proteomics studies. These models decompose variations from different sources, such as experimental factors, proteins’ biological features and biological samples’ physiological properties. It is also known that in the ion intensities data which is used to construct protein’s relative abundance, its missing values depend on the unobserved abundance values and the probability of missing associated with some observable experiment factors. In addition to this, due to the detection limit of a device, we also observe left-censoring values in the output from many studies. We model non-random missingness as either left-censored or completely missing using one of the proposed likelihood-based method of Little and Rubin.
Little (1993), Little and Rubin (1987) proposed three likelihood-based methods (selection model, pattern-mixture model and shared-parameter model) for Missing Not At Random (MNAR). Selection model, in which the missing data mechanism and parameters of the inferential model are conditionally estimated given the hypothetical completed data, is an appropriate approach given the abundance-dependent missingness. The hypothetical completed data comprises the missing and observed values. The known physical properties of peptides and mass-spectrometers provide the auxiliary information for the missing probability and missing values. We use the selection model factorization and include the left-censored missingness under a Bayesian framework. In two case studies, we estimate the likelihood of missingness using auxiliary variable mass-to-charge ratio (m/z) and the intensity values.

## 1.3 The computing algorithm

Our proposed method utilizes Hamiltonian MC/Non-U-Turn Sampling for the posterior distribution. Compared to the Gibb sampling, Hamiltonian MC (HMC) avoids the random walk by introducing the leapfrog function. It provides an alternative to approximating the solution on the continuous time scale from the solutions on the discrete time scale, with a specified step-size. The logarithmic posterior probability function was simulated by one of the paired partial differentiated equations, namely the Hamiltonian. Larger moving steps are generated from the leap frog scheme, and this helps to improve the convergence compared to the random walk. It has been shown to have higher efficiency in sampling high-dimensional correlated multivariate distributions. RStan is a new tool recently developed to implement HMC modeling for Bayesian data analysis, of which the posterior distribution can be sampled using the No-U-Turn Sampling method -an extension of HMC. No-U-Turn sampling implements a recursive algorithm that will enable auto-tuning of the couple parameters in HMC: numbers of leap-frog steps and the discrete step. We based on Rstan library function stan() to set up the program for our proposed model. The R library “mlm4omics” is written to implement posterior samplings using HMC.

# 2 The mlm4omics R package

The functions included in the pacakge “mlm4omics” are:

1. setinitialvalue(): this function generates initiating values for the parameters in the computing.

2. mlmc(): this function provides Rstan program to estimate the posterior distribution of parameters in a linear regression model with multivariate normal distributed respondents , for example the log (peptide ion intensity) from mass-spec, and explanatory variables such as a protein identification. Other explanatory variables can be experimental variables and clinical design variables. The protein identification is an ID to link the peptides to its parent protein. The probability of missing, missing and censored values of the respondent variable are estimated in the algorithm through a logistic regression and a truncated normal distribution for the censored values. The probability of missing can be fitted via explanatory variables such as the abundance and mass-to-charge ratio of the observed ion intensity. The censored values are modelled through defining a truncated normal as its prior distribution. mlmc() function also provides option (respond_dep_missing=TRUE/FALSE) to choose if the missing probability is depending on the magnitude of missing value.

3. mlmm(): this function is similar to mlmc(), but it is written for data without any censored respondent values in the regression model. The input variables in this function are similar to those in mlmm() but without terms for censored missing.

# 3 Installation

There are two ways to install the package. Before install the package, please use the instructions from https://github.com/stan-dev/rstan/wiki to install the rstan package firstly.

## 3.1 Use bioconductor as installation source

    if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("mlm4omics")
    library(mlm4omics)

## 3.2 Use github as installation source

    devtools::install_github("superOmics/mlm4omics")
    library(mlm4omics)

# 4 Instructions and Examples of using main functions for regressions

The following section provides examples using a simulated dataset: pdata. The examples include:

1. Examples with no informative prior to call the compiled model; and the default priors are used.

2. Example with informative priors.

When set up the formulas and parameters for using the mlmm() and mlmc() functions, some instructions are as follows:

1.The formula_subject will need to have at least one variable. For example, to include subject-id for subject and this will provide the mean value for each subject.

2.The missing and censor values in the response variable require Different indicators, the missing and censor indictors are mutually exclusive.

3.The number of iterations are better to set to be >1000. The default burn-in is half of the set iterations.

## 4.1 Read data and check data

    library(mlm4omics)

    ## Loading required package: Rcpp

    ## Registered S3 methods overwritten by 'ggplot2':
    ##   method         from
    ##   [.quosures     rlang
    ##   c.quosures     rlang
    ##   print.quosures rlang

    data(pdata,package="mlm4omics")
    pdata[1:2,]

    ##    var1       var2 treatment miss censor geneid sid
    ## 1    NA 1.63230970         0    1      0      1   1
    ## 2 0.002 0.06299626         0    0      1      2   2

### 4.1.1 Check for missing and censor indictors

Check if there is overlapped definition for missing and censor.
The definition for missing and censored are differently defined. The following code is to check if the indicators for missing and censored are mutually exclusive.

    table(pdata$miss, pdata$censor)

    ##
    ##      0  1
    ##   0 48 17
    ##   1 27  8

If the indictors coded for missing and censored are overlapped, the following codes convert missing indictor from 1 to 0 where censor =1.

    n=dim(pdata)[1]
    for (i in seq_len(n))
    if (pdata$miss[i]==1 && pdata$censor[i]==1) pdata$miss[i]=0

After that, we can set the formula for response, probability of missing and subject as the followings:

### 4.1.2 Set the formula for response, missingness and subject

    formula_completed=var1~var2+treatment;
    formula_missing=miss~var2;
    formula_censor=censor~1;
    formula_subject=~treatment;
    response_censorlim=0.002

Where var1 is the response variable, miss and censor are the indicator variables for missing and censoring respectively.

## 4.2 Set the formula and fit the regression model

The formula for subject requires at least one variable. For example, to include subject-id for subject and this will provide the mean value for each subject.

### 4.2.1 R examples for mlmc()

#### 4.2.1.1 Example 1a. Use default prior and call complied codes

When no priors information available, we can have posterior samples based on the compiled codes, given the default flat prior of precision matrix of the protein level covariates (a second level unit in the multilevel model), and prior for regression coefficients of the logistic regression predicting the missing probability.

    model1 <- mlmc(formula_completed = var1~var2+treatment,
    formula_missing = miss~var2,
    formula_censor = censor~1, formula_subject = ~sid,
    pdata = pdata, response_censorlim = 0.002, respond_dep_missing = TRUE,
    pidname = "geneid", sidname = "sid",
    iterno = 100, chains = 2, savefile = TRUE)

    ##
    ## SAMPLING FOR MODEL 'mlmc_code' NOW (CHAIN 1).
    ## Chain 1:
    ## Chain 1: Gradient evaluation took 0 seconds
    ## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
    ## Chain 1: Adjust your expectations accordingly!
    ## Chain 1:
    ## Chain 1:
    ## Chain 1: WARNING: There aren't enough warmup iterations to fit the
    ## Chain 1:          three stages of adaptation as currently configured.
    ## Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
    ## Chain 1:          the given number of warmup iterations:
    ## Chain 1:            init_buffer = 7
    ## Chain 1:            adapt_window = 38
    ## Chain 1:            term_buffer = 5
    ## Chain 1:
    ## Chain 1: Iteration:  1 / 100 [  1%]  (Warmup)
    ## Chain 1: Iteration: 10 / 100 [ 10%]  (Warmup)
    ## Chain 1: Iteration: 20 / 100 [ 20%]  (Warmup)
    ## Chain 1: Iteration: 30 / 100 [ 30%]  (Warmup)
    ## Chain 1: Iteration: 40 / 100 [ 40%]  (Warmup)
    ## Chain 1: Iteration: 50 / 100 [ 50%]  (Warmup)
    ## Chain 1: Iteration: 51 / 100 [ 51%]  (Sampling)
    ## Chain 1: Iteration: 60 / 100 [ 60%]  (Sampling)
    ## Chain 1: Iteration: 70 / 100 [ 70%]  (Sampling)
    ## Chain 1: Iteration: 80 / 100 [ 80%]  (Sampling)
    ## Chain 1: Iteration: 90 / 100 [ 90%]  (Sampling)
    ## Chain 1: Iteration: 100 / 100 [100%]  (Sampling)
    ## Chain 1:
    ## Chain 1:  Elapsed Time: 4.344 seconds (Warm-up)
    ## Chain 1:                4.438 seconds (Sampling)
    ## Chain 1:                8.782 seconds (Total)
    ## Chain 1:
    ##
    ## SAMPLING FOR MODEL 'mlmc_code' NOW (CHAIN 2).
    ## Chain 2:
    ## Chain 2: Gradient evaluation took 0 seconds
    ## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
    ## Chain 2: Adjust your expectations accordingly!
    ## Chain 2:
    ## Chain 2:
    ## Chain 2: WARNING: There aren't enough warmup iterations to fit the
    ## Chain 2:          three stages of adaptation as currently configured.
    ## Chain 2:          Reducing each adaptation stage to 15%/75%/10% of
    ## Chain 2:          the given number of warmup iterations:
    ## Chain 2:            init_buffer = 7
    ## Chain 2:            adapt_window = 38
    ## Chain 2:            term_buffer = 5
    ## Chain 2:
    ## Chain 2: Iteration:  1 / 100 [  1%]  (Warmup)
    ## Chain 2: Iteration: 10 / 100 [ 10%]  (Warmup)
    ## Chain 2: Iteration: 20 / 100 [ 20%]  (Warmup)
    ## Chain 2: Iteration: 30 / 100 [ 30%]  (Warmup)
    ## Chain 2: Iteration: 40 / 100 [ 40%]  (Warmup)
    ## Chain 2: Iteration: 50 / 100 [ 50%]  (Warmup)
    ## Chain 2: Iteration: 51 / 100 [ 51%]  (Sampling)
    ## Chain 2: Iteration: 60 / 100 [ 60%]  (Sampling)
    ## Chain 2: Iteration: 70 / 100 [ 70%]  (Sampling)
    ## Chain 2: Iteration: 80 / 100 [ 80%]  (Sampling)
    ## Chain 2: Iteration: 90 / 100 [ 90%]  (Sampling)
    ## Chain 2: Iteration: 100 / 100 [100%]  (Sampling)
    ## Chain 2:
    ## Chain 2:  Elapsed Time: 0.5 seconds (Warm-up)
    ## Chain 2:                0.453 seconds (Sampling)
    ## Chain 2:                0.953 seconds (Total)
    ## Chain 2:

    ## Warning: There were 89 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
    ## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

    ## Warning: There were 11 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
    ## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

    ## Warning: Examine the pairs() plot to diagnose sampling problems

    ## Inference for Stan model: mlmc_code.
    ## 2 chains, each with iter=100; warmup=50; thin=1;
    ## post-warmup draws per chain=50, total post-warmup draws=100.
    ##
    ##                   mean se_mean    sd    2.5%     25%     50%    75%  97.5%
    ## U[1,1]            1.92    1.58  1.88    0.20    0.35    0.73   3.20   5.49
    ## U[1,2]            0.92    0.30  0.31    0.42    0.60    0.95   1.22   1.30
    ## U[1,3]           -0.11    0.77  0.98   -2.15   -0.96    0.46   0.69   0.76
    ## U[2,1]           -1.77    3.33  3.95  -10.47   -3.93    0.03   1.56   1.67
    ## U[2,2]            0.11    1.03  1.09   -1.25   -0.93    0.03   1.00   1.79
    ## U[2,3]            0.73    0.16  0.61   -0.46    0.44    0.61   0.78   2.03
    ## U[3,1]            1.08    0.19  0.68   -0.42    0.93    1.29   1.47   2.26
    ## U[3,2]           -0.01    1.19  1.21   -1.40   -1.29    0.17   1.18   1.25
    ## U[3,3]            0.18    0.38  0.88   -1.65   -0.03    0.19   0.42   2.03
    ## U[4,1]           -0.72    1.27  2.04   -5.38   -2.30    0.45   0.71   0.88
    ## U[4,2]           -0.06    1.05  1.10   -1.26   -1.11   -0.56   1.03   1.60
    ## U[4,3]            0.59    0.44  0.92   -0.74   -0.23    0.82   1.25   2.88
    ## U[5,1]           -0.55    1.77  1.86   -3.89   -2.15    0.16   1.25   1.49
    ## U[5,2]           -0.44    1.66  1.67   -2.45   -2.03   -0.22   1.22   1.25
    ## U[5,3]            0.18    0.98  1.31   -2.47   -0.86    0.90   1.17   1.63
    ## U[6,1]           -0.89    2.47  2.53   -4.09   -3.59   -0.42   1.73   1.95
    ## U[6,2]           -0.78    2.17  2.20   -3.28   -3.00   -0.87   1.32   1.97
    ## U[6,3]            0.76    0.55  0.77   -1.05    0.49    0.92   1.28   1.87
    ## U[7,1]           -0.63    2.75  3.17   -7.62   -2.47    0.68   2.28   2.47
    ## U[7,2]           -0.11    1.30  1.35   -1.93   -1.73    0.38   1.18   1.25
    ## U[7,3]           -0.03    0.53  0.76   -1.01   -0.74   -0.12   0.52   1.47
    ## U[8,1]            1.86    1.10  1.24    0.45    0.89    1.12   3.11   4.33
    ## U[8,2]            0.14    1.00  1.01   -1.20   -0.87    0.40   1.13   1.20
    ## U[8,3]            0.93    0.68  1.14   -1.52    0.50    1.34   1.72   2.12
    ## U[9,1]            1.43    0.91  1.89   -1.17    0.50    0.95   1.68   5.61
    ## U[9,2]           -0.53    1.69  1.70   -2.41   -2.19   -0.54   1.15   1.24
    ## U[9,3]            0.35    0.21  0.68   -0.80    0.19    0.29   0.48   2.43
    ## U[10,1]          -0.19    1.40  1.53   -3.10   -1.30    0.52   1.27   1.54
    ## U[10,2]           0.03    1.18  1.19   -1.43   -1.16    0.15   1.21   1.25
    ## U[10,3]          -0.27    0.49  1.00   -2.36   -0.84    0.00   0.32   1.69
    ## beta2[1,1]       -1.00    0.71  0.97   -1.94   -1.71   -1.65  -0.12   0.86
    ## beta2[1,2]        0.33    1.02  1.24   -1.61   -1.08    1.24   1.31   1.54
    ## beta2[2,1]       -0.49    0.25  0.75   -2.30   -0.67   -0.40  -0.07   0.51
    ## beta2[2,2]       -0.94    1.10  1.16   -2.17   -1.99   -1.53   0.32   0.99
    ## beta2[3,1]        0.62    0.21  0.57   -0.28    0.30    0.44   0.92   1.69
    ## beta2[3,2]        0.12    0.08  0.24   -0.36    0.02    0.17   0.29   0.50
    ## beta2[4,1]        1.30    0.96  1.01   -0.54    0.55    1.28   2.27   2.54
    ## beta2[4,2]       -0.23    0.45  0.69   -1.03   -0.71   -0.53   0.27   1.12
    ## beta2[5,1]        1.12    0.35  0.57   -0.62    0.98    1.22   1.55   1.70
    ## beta2[5,2]        0.37    0.18  0.25    0.12    0.19    0.30   0.46   0.99
    ## beta2[6,1]        0.13    0.40  0.63   -0.37   -0.31   -0.14   0.32   1.76
    ## beta2[6,2]        0.76    0.53  0.59    0.04    0.24    0.76   1.18   1.87
    ## beta2[7,1]       -0.27    0.58  0.70   -0.95   -0.81   -0.65   0.25   1.29
    ## beta2[7,2]        0.42    0.12  0.28   -0.10    0.29    0.47   0.56   0.95
    ## beta2[8,1]       -0.01    0.24  0.40   -0.45   -0.29   -0.22   0.25   0.83
    ## beta2[8,2]        0.47    0.17  0.22    0.02    0.28    0.52   0.65   0.74
    ## beta2[9,1]        0.17    1.45  1.72   -1.72   -1.23   -1.11   1.88   3.45
    ## beta2[9,2]        0.33    0.19  0.33   -0.28    0.14    0.39   0.51   0.98
    ## beta2[10,1]       0.57    0.70  0.84   -0.50   -0.09    0.25   1.44   2.21
    ## beta2[10,2]       0.45    0.05  0.12    0.17    0.40    0.48   0.54   0.60
    ## alpha[1]         -0.08    0.20  0.41   -0.72   -0.35   -0.13   0.23   0.79
    ## alpha[2]         -0.92    0.07  0.37   -1.51   -1.14   -1.03  -0.78  -0.09
    ## alpha_response    0.04    0.05  0.13   -0.30    0.02    0.07   0.11   0.20
    ## lp__           -120.65   88.51 89.20 -225.09 -207.10 -120.90 -33.12 -22.00
    ##                n_eff  Rhat
    ## U[1,1]             1  3.35
    ## U[1,2]             1  3.42
    ## U[1,3]             2  2.13
    ## U[2,1]             1  2.78
    ## U[2,2]             1  5.92
    ## U[2,3]            15  1.07
    ## U[3,1]            12  1.11
    ## U[3,2]             1 10.44
    ## U[3,3]             5  2.33
    ## U[4,1]             3  3.00
    ## U[4,2]             1  3.35
    ## U[4,3]             4  1.28
    ## U[5,1]             1  4.33
    ## U[5,2]             1 15.43
    ## U[5,3]             2  2.59
    ## U[6,1]             1  5.40
    ## U[6,2]             1  7.85
    ## U[6,3]             2  1.58
    ## U[7,1]             1  2.79
    ## U[7,2]             1  5.28
    ## U[7,3]             2  1.74
    ## U[8,1]             1  2.21
    ## U[8,2]             1 10.42
    ## U[8,3]             3  2.87
    ## U[9,1]             4  2.22
    ## U[9,2]             1 23.32
    ## U[9,3]            11  1.20
    ## U[10,1]            1  3.53
    ## U[10,2]            1 16.40
    ## U[10,3]            4  1.97
    ## beta2[1,1]         2  2.07
    ## beta2[1,2]         1  2.08
    ## beta2[2,1]         9  1.15
    ## beta2[2,2]         1  5.56
    ## beta2[3,1]         7  1.27
    ## beta2[3,2]         8  1.22
    ## beta2[4,1]         1  4.54
    ## beta2[4,2]         2  2.88
    ## beta2[5,1]         3  1.73
    ## beta2[5,2]         2  2.28
    ## beta2[6,1]         2  1.58
    ## beta2[6,2]         1  4.00
    ## beta2[7,1]         1  2.03
    ## beta2[7,2]         5  1.86
    ## beta2[8,1]         3  1.35
    ## beta2[8,2]         2  2.26
    ## beta2[9,1]         1  2.42
    ## beta2[9,2]         3  2.24
    ## beta2[10,1]        1  1.85
    ## beta2[10,2]        5  1.44
    ## alpha[1]           4  1.36
    ## alpha[2]          24  1.05
    ## alpha_response     6  1.21
    ## lp__               1 17.27
    ##
    ## Samples were drawn using NUTS(diag_e) at Wed May 15 00:43:56 2019.
    ## For each parameter, n_eff is a crude measure of effective sample size,
    ## and Rhat is the potential scale reduction factor on split chains (at
    ## convergence, Rhat=1).

#### 4.2.1.2 Example 1b-1c. Use a provided informative priors

Example 1b is to set priors for regression coefficients in the logistic regression model for missingness. alpha_prior is a prior for coefficients (including intercept) in the logistic regression to predict missing probability; it is a vector with a length of (1+number of predictors).

Example 1c is to Use an informative prior for precision matrix in the regression for completed data. prec_prior is the precision matrix prior of regression coefficients -intercept and predictor (var2) in the completed data formula.

    #Example 1b
    model1 <- mlmc(formula_completed=var1~var2+treatment,
    formula_missing=miss~var2,
    formula_censor=censor~1, formula_subject=~sid,
    pdata=pdata, response_censorlim=0.002,
    respond_dep_missing=TRUE, pidname="geneid", sidname="sid",
    alpha_prior <- c(0,0.001), iterno=100, chains=2,
    savefile=TRUE)

    #Example 1c
    prec_example <- matrix(c(0.01,0.001,0.001,0.01),nrow=2,ncol=2)

    model3 <- mlmc(formula_completed=var1~var2,
    formula_missing=miss~var2,
    formula_censor=censor~1, formula_subject=~sid+treatment,
    pdata=pdata, response_censorlim=0.002, respond_dep_missing=TRUE,
    pidname="geneid", sidname="sid", prec_prior=prec_example,
    iterno=1000, chains=2, savefile=TRUE)

### 4.2.2 R examples for mlmm()

Formulas for completed data, subject and missingness all require at least one variable.

#### 4.2.2.1 Example 2a. Use default prior

    model5 <- mlmm(formula_completed=var1~var2+treatment,
    formula_missing=miss~var2,
    formula_subject=~sid, pdata=pdata,
    respond_dep_missing=FALSE, pidname="geneid", sidname="sid",
    iterno=1000, chains=2, savefile=FALSE)

#### 4.2.2.2 Example 2b-2c. Provide informative priors

The user define priors include:

1. prec_prior : the precision matrix of the first-level parameters (var2 and treatment);

2. alpha_prior: coefficients for intercept and var2 to predict probability of having missing value in logistic regression.

    prec_example <- matrix(c(0.01,0.001,0.001,0.001,0.01,0.001,0.001,0.001,0.01),
    nrow=3, ncol=3)

    #Example 2b use both priors
    model4 <- mlmm(formula_completed=var1~var2+treatment,
    formula_missing=miss~var2,
    formula_subject=~sid, pdata=pdata,
    respond_dep_missing=FALSE, pidname="geneid", sidname="sid",
    alpha_prior <- c(0,0.001), prec_prior=prec_example,
    iterno=100, chains=2, savefile=TRUE)

    #Example 2c. Use alpha prior
    model5 <- mlmm(formula_completed=var1~var2+treatment,
    formula_missing=miss~var2,
    formula_subject=~sid+treatment, pdata=pdata,
    respond_dep_missing=FALSE, pidname="geneid",alpha_prior <- c(0,0.001),
    sidname="sid", iterno=100, chains=2, savefile=TRUE)

## 4.3 How to make plots of the parameters and to diagnose convergence

If the estimate converges, it will have high number of efficient samples and the rhat value will close to 1.

The following codes are to:

1. plot posterior parameter using the outsummary results from multiple chains.

2. plot the trajectory of posterior sampling including those values generated from burn-in iterations.

### 4.3.1 plot posterior parameters

    summaryreader <- read.csv(file = file.path(getwd(),"outsummary.csv"),
    header=TRUE, sep=",", skip=0)

    iterno <- dim(summaryreader)[1]; burnin <- iterno/2

    U.1.1 <- rowMeans(matrix(c(summaryreader$chain.1.U.1.1,
    summaryreader$chain.2.U.1.1), nrow=iterno, ncol=2))[burnin:iterno]

    meanU <- mean(U.1.1)
    qU <- quantile(U.1.1,p <- seq(0, 1, by=0.025))
    scale <- seq(0, 1, by=0.025)

    plot(scale, qU, pch=19, ylab="quantiles of estimate", xlab="quantiles")

    segments(0,qU[names(qU)=="50%"],1,qU[names(qU)=="50%"],lwd=2,col="red")
    segments(0,qU[names(qU)=="2.5%"],1,qU[names(qU)=="2.5%"],lty=2,lwd=2,col="red")
    segments(0,qU[names(qU)=="97.5%"],1,qU[names(qU)=="97.5%"],lty=2,lwd=2,
    col="red")

    legend(0.5,qU[names(qU)=="50%"],"median",cex=0.8,bty="n")
    legend(0.03,qU[names(qU)=="2.5%"],"2.5%",cex=0.8,bty="n")
    legend(0.90,qU[names(qU)=="97.5%"],"97.5%",cex=0.8,bty="n")

    qU

    ##       0%     2.5%       5%     7.5%      10%    12.5%      15%    17.5%
    ## 1.654635 1.715877 1.762414 1.786891 1.794677 1.815199 1.903362 1.970699
    ##      20%    22.5%      25%    27.5%      30%    32.5%      35%    37.5%
    ## 2.024153 2.042247 2.182003 2.504250 2.569184 2.572293 2.585412 2.602249
    ##      40%    42.5%      45%    47.5%      50%    52.5%      55%    57.5%
    ## 2.621565 2.624096 2.625682 2.625851 2.626606 2.628631 2.635147 2.637798
    ##      60%    62.5%      65%    67.5%      70%    72.5%      75%    77.5%
    ## 2.637873 2.697287 2.747229 2.782961 2.814567 2.839814 2.843753 2.850625
    ##      80%    82.5%      85%    87.5%      90%    92.5%      95%    97.5%
    ## 2.859453 2.862589 2.866191 2.870493 2.895748 2.924094 2.943853 2.984407
    ##     100%
    ## 3.038825

### 4.3.2 plot trajectory of the posterior samples

    sample1reader <- read.csv(file <- file.path(getwd(),"samples_1.csv"),
    header=TRUE, sep=",", skip=25)

    sample2reader <- read.csv(file <- file.path(getwd(),"samples_2.csv"),
    header=TRUE, sep=",", skip=25)

    #plot variable U.1.1 - the intercept of first unit
    trajectory_length <- dim(sample1reader)[1]

    plot(seq(1, trajectory_length, by=1), sample1reader$U.1.1, xlab="trajectory
    number", ylab="U.1.1", type="n",
    ylim = c(min(sample1reader$U.1.1, sample2reader$U.1.1, na.rm=TRUE),
    max(sample1reader$U.1.1, sample2reader$U.1.1, na.rm=TRUE)))

    trajectory <- seq(1, trajectory_length, by=1)

    lines(trajectory, sample1reader$U.1.1)
    lines(trajectory, sample2reader$U.1.1, col="red")

# 5 sessionInfo

All of the outputs in this vignette are produced under the following conditions:

    sessionInfo()

    ## R version 3.6.0 Patched (2019-05-13 r76493)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] mlm4omics_1.3.0  Rcpp_1.0.1       rmarkdown_1.12   knitr_1.22
    ## [5] BiocStyle_2.13.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] compiler_3.6.0     pillar_1.4.0       BiocManager_1.30.4
    ##  [4] plyr_1.8.4         prettyunits_1.0.2  tools_3.6.0
    ##  [7] pkgbuild_1.0.3     digest_0.6.18      evaluate_0.13
    ## [10] tibble_2.1.1       gtable_0.3.0       lattice_0.20-38
    ## [13] pkgconfig_2.0.2    rlang_0.3.4        Matrix_1.2-17
    ## [16] cli_1.1.0          parallel_3.6.0     yaml_2.2.0
    ## [19] xfun_0.7           loo_2.1.0          gridExtra_2.3
    ## [22] stringr_1.4.0      dplyr_0.8.1        stats4_3.6.0
    ## [25] grid_3.6.0         tidyselect_0.2.5   inline_0.3.15
    ## [28] glue_1.3.1         R6_2.4.0           processx_3.3.1
    ## [31] bookdown_0.10      rstan_2.18.2       callr_3.2.0
    ## [34] ggplot2_3.1.1      purrr_0.3.2        magrittr_1.5
    ## [37] codetools_0.2-16   rstantools_1.5.1   matrixStats_0.54.0
    ## [40] StanHeaders_2.18.1 ps_1.3.0           scales_1.0.0
    ## [43] htmltools_0.3.6    MASS_7.3-51.4      assertthat_0.2.1
    ## [46] colorspace_1.4-1   stringi_1.4.3      lazyeval_0.2.2
    ## [49] munsell_0.5.0      crayon_1.3.4

# 6 References

1. Hoffman, M.D. and Gelman, A. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of machine learning research 2011;12.
2. Little, R.J.A. Pattern-mixture models for multivariate incomplete data. J. Am. Statist. Assoc. 1993;88:125-134.
3. Little, R.J.A. and Rubin, D.B. Statistical Analysis with Missing Data. 1987. Neal, R.M. MCMC using Hamiltonian dynamics. In: S. Brooks, e.a., editor, Handbook of Markov Chain Monte Carlo. 2011.
4. Zeng,I.S. PhD thesis. Statistical methods in clinical proteomic studies 2007-2014 -A protein concerto-. Research Space: The University of Auckland.
5. Zeng, I.S. Topics in Study Design and Analysis for Multistage Clinical Proteomics In: Jung, K., editor, Statistical Analysis in Proteomics New York: Humana Press: Springer; 2015.
6. Hrydziuszko, O., and Viant, M. R.. Missing values in mass spectrometry based metabolomics: an undervalued step in the data processing pipeline. Metabolomics 2012;8,S161S174.
7. Zeng, I.S., Lumley, T., Ruggiero, K. and Middleditch, M. A Bayesian approach to multivariate and multilevel modelling with non-random missingness for hierarchical clinical proteomics data. In: https://www.biorxiv.org/content/early/2017/06/21/153049; 2017.