NEJSDS The New England Journal of Statistics in Data Science 2693-71662693-7166 New England Statistical Society NEJSDS68 10.51387/24-NEJSDS68 Statistical Methodology Methodology Article Linear Mixed-effects Models for Censored Data with Serial Correlation Errors Using the Multivariate Student’s t-distribution ZhongKelinkelin.zhong@uconn.edu OlivariRommy C.rco1@de.ufpe.br GarayAldo M.agaray@de.ufpe.br∗ LachosVíctor H.hlachos@uconn.edu Department of Statistics, University of Connecticut, Storrs, USA. E-mail address: kelin.zhong@uconn.edu Department of Statistics, Federal University of Pernambuco, Brazil. E-mail address: rco1@de.ufpe.br Department of Statistics, Federal University of Pernambuco, Prof. Moraes Rego, 1235 Cidade Universitaria, Recife, PE, CEP: 50670-901, Brazil. E-mail address: agaray@de.ufpe.br Department of Statistics, University of Connecticut, Storrs, USA. E-mail address: hlachos@uconn.edu Corresponding author. 2025307202431119134672024 © 2025 New England Statistical Society2025 Open access article under the CC BY license.

The purpose of this paper is to develop a practical framework for the analysis of the linear mixed-effects models for censored or missing data with serial correlation errors, using the multivariate Student’s t-distribution, being a flexible alternative to the use of the corresponding normal distribution. We propose an efficient ECM algorithm for computing the maximum likelihood estimates for these models with standard errors of the fixed effects and likelihood function as a by-product. This algorithm uses closed-form expressions at the E-step, which relies on formulas for the mean and variance of a truncated multivariate Student’s t-distribution. In order to illustrate the usefulness of the proposed new methodology, artificial and a real dataset are analyzed. The proposed algorithm and methods are implemented in the R package ARpLMEC.

 AR(p) models Censored data Damped exponential correlation EM algorithm Linear mixed effects models Student’s t-distribution The research of Aldo M. Garay was supported by Grant APQ-0950-1.02/22 from Fundação de Amparo à Ciência e Tecnologia do Estado de Pernambuco FACEPE - Brazil and by Grant 441476/2023-6 from National Council for Scientific and Technological Development – CNPq – Brazil. Victor H. Lachos acknowledges the partial financial support from the Office of the Vice President for Research and UConn - CLAS’s Summer Research Funding Initiative 2023. Introduction

Linear mixed effects (LME) models are frequently used to analyze repeated measures data because they model flexibly the within-subject correlation often present in this type of data. Usually, for mathematical convenience, it is assumed that both random effect and error term follow normal distributions (N-LME). These restrictive assumptions, however, may result in a lack of robustness against departures from the normal distribution and invalid statistical inferences, especially when the data show heavy tails. For instance, substantial bias in the maximum likelihood (ML) estimates of regression parameters can result when the random-effects distribution is misspecified [11].

Several researchers have investigated alternative distributions for errors in LME models. For example, [28] propose a robust hierarchical linear mixed model in which the random effects and the within-subject errors have a multivariate Student’s t-distribution (t-LME). [18] and [16] developed some additional tools for t-LME models from likelihood-based and Bayesian perspectives. However, in longitudinal studies, such as those on environmental pollution and infectious diseases, measurements of some variables may be subjected to certain threshold values below or above which the measurements are not quantifiable. For instance, viral load measures the amount of actively replicating virus and, depending upon the diagnostic assays used, its measurement may be subjected to detection limits, below or above, which are not quantifiable. For this kind of data, [30] proposed the LME model for censored responses (LMEC), and an exact EM-type algorithm is proposed. Recently, [21] and [23] have considered the Student’s t-distribution in the context of LMEC models (t-LMEC), including influence diagnostics analyses, under different perturbation schemes. A common feature of these classes of LMECs is that the error terms are unconditionally independent. However, in longitudinal studies, repeated measures are collected over time, and hence, the error term tends to be serially correlated.

There are several proposals in the literature that account for the time dependence in longitudinal data. For instance, [9] studied longitudinal data with random effects and AR(1) errors under the Bayesian paradigm, [4] extended the classical random effects model to deal with time dependence, leaving the random effects distribution unspecified, [13] proposed a robust structure for a censored linear model based on the multivariate Student’s t-distribution, considering a damped exponential correlation (DEC) structure. More recently, [26] have proposed a full likelihood-based approach for the Gaussian N-LMEC modeling with autoregressive correlation of order p (AR(p)) errors (AR(p)-LMEC), including the implementation of the EM algorithm for maximum likelihood (ML) estimation.

Even though some proposals have been made to deal with the problem of serial correlation among the observations in LMEC models, to the best of our knowledge, there are no studies of t-LMEC with serially correlated error structures, such as DEC or AR(p). Thus, in this paper we propose the t-LMEC modeling considering some useful correlation structures (CSs). We develop a full likelihood-based treatment, including the implementation of a computationally efficient estimation method via the EM algorithm with the likelihood function, predictions of unobservable values of the response, and the asymptotic standard errors as byproducts. The model developed here is an extension of those previously presented by [21], [22] and [26] for the analysis of mixed-effects models with censored responses and HIV data.

The proposed algorithm and methods are implemented in the R package ARpLMEC [27], available on CRAN repository.

The rest of the paper is organized as follows: Section 2 discusses some preliminary results related to the multivariate Student’s t-distribution and its truncated version. Some of its key properties are also presented. In Section 3, the related likelihood-based inference is presented, including estimation of the standard errors for the regression parameters. The application of the proposed method is presented in Sections 4 and 5 through some simulation studies and the analysis of data from a real HIV case study, respectively. Finally, Section 6 concludes with a short discussion of issues raised by this study and some possible directions for future research.

 Preliminaries

We begin our exposition by defining the notation and presenting some basic concepts which are used throughout the development of our methodology. We denote a random variable by an upper-case letter and its realization by the corresponding lower case and use boldface letters for vectors and matrices. I p represents a p × p identity matrix and A ⊤ denotes the transpose of A matrix. Thus, for a = ( a 1 , … , a p ) ⊤ and b = ( b 1 , … , b p ) ⊤ we have that, if the Borel set in R p has the form: A = { ( y 1 , … , y p ) ∈ R p : a 1 ≤ y 1 ≤ b 1 , … , a p ≤ y p ≤ b p } = { y ∈ R p : a ≤ y ≤ b } then, we use the shorthand notations { Y ∈ A } = { a ≤ Y ≤ b } and ∫ a b f ( y ) d y = ∫ a 1 b 1 … ∫ a p b p f ( y 1 , … , y p ) d y p … d y 1 .

The Multivariate Student’s <italic>t</italic>-distribution

Y ∼ t p ( μ , Σ , ν ) denotes a random vector Y following a p-variate Student’s-t distribution, with location vector μ, positive-definite scale-covariance matrix Σ and degrees of freedom ν and t p ( y ∣ μ , Σ , ν ) denotes its probability density function (pdf).

The cumulative distribution function (cdf) of Y on [ a , b ] is denoted by: L p ( a , b ; μ , Σ , ν ) = ∫ a b t p ( y | μ , Σ , ν ) d y .

An important property of the random vector Y is that it can be written as a scale mixture of the MVN (multivariate normal) random vector coupled with a positive random variable, i.e., Y = μ + U − 1 / 2 Z , where Z ∼ N p ( 0 , Σ ) independent of U ∼ Gamma ( ν / 2 , ν / 2 ), where N p ( μ , Σ ) and Gamma ( a , b ) denote the p-variate normal with mean vector μ and variance-covariance matrix Σ and the gamma distribution with mean a / b, respectively.

 The Multivariate Truncated Student’s <italic>t</italic>-distribution

A p-dimensional random vector Y is said to follow a doubly truncated Student’s t-distribution with location vector μ, scale-covariance matrix Σ and degrees of freedom ν over the truncation region A, defined in Eq. (2.1), denoted by Y ∼ T t p ( μ , Σ , ν ; A ), if it has the pdf: T t p ( y | μ , Σ , ν ; A ) = t p ( y | μ , Σ , ν ) L p ( a , b ; μ , Σ , ν ) , a ≤ y ≤ b , where L p ( a , b ; μ , Σ , ν ) as defined in Eq. (2.2).

The cdf of Y evaluated at the region { a ≤ y ≤ b } is: T T p ( y | μ , Σ , ν ; A ) = 1 L p ( a , b ; μ , Σ , ν ) ∫ a y t p ( x | μ , Σ , ν ) d x = L p ( a , y ; μ , Σ , ν ) L p ( a , b ; μ , Σ , ν ) .

Some propositions and properties related to the p-variate Student’s t-distribution and the marginal and conditional moments of the first two moments of TMVT distributions under a double truncation, which are useful for our theoretical developments, can be found in Appendix A and B, respectively.

 Linear Mixed-effects with Censored Response Model Specification

We proceed as in [21, 28], by considering a generalization of the classical N-LME model in which the random terms are assumed to follow a Student’s-t distribution as follows: Y i = X i β + Z i b i + ϵ i , with ( b i , ϵ i ) ⊤ ∼ t n i + q ( 0 , D i a g ( D , Ω i ) , ν ) . The subscript “i” is the subject index; Diag(A,B) is a block diagonal matrix whose elements are the matrices A and B. Y i = ( Y i 1 , … , Y i n i ) ⊤ is a n i × 1 vector of observed continuous responses for sample unit “i”, X i is the n i × s design matrix corresponding to the fixed effects β, which is the s × 1 vector of population-averaged regression coefficients. Z i is the n i × q design matrix corresponding to the q × 1 vector of random effects b i . ϵ i is the n i × 1 vector of random errors, the dispersion matrix D = D ( α ) depends on unknown and reduced parameters α. The correlation structure of the error vector is assumed to be Ω i = σ 2 R i , where R i = R i ( ϕ ), for ϕ = ( ϕ 1 , … , ϕ p ) ⊤ is a n i × n i matrix, that incorporates a time-dependence structure.

Note that b i and ϵ i are uncorrelated, once Cov ( b i , ϵ i ) = E [ b i ϵ i ⊤ ] = E [ E ( b i ϵ i ⊤ | U i ) ] = 0, where U i is a scalar generated from Gamma ( ν / 2 , ν / 2 ). Classical inference on the parameter vector θ = ( β ⊤ , σ 2 , α ⊤ , ϕ ⊤ , ν ) ⊤ is based on the marginal distribution of Y i , thus Y i ∼ ind . t n i ( μ i , Σ i , ν ), for i = 1 , … , n, where μ i = X i β and Σ i = σ 2 R i + Z i D Z i ⊤ . The estimates from the multivariate t-LME are more robust against outliers than those based on the standard LME; in this sense, [28] showed by simulation study that the t-LME substantially outperforms the normal or standard LME when outliers are present in the data. This issue has also been discussed by [21] in the context of censored mixed-effects models.

In this study, we follow [13, 15, 21, 23, 26] by considering non-informative censoring observations (or censoring at random), i.e. we assume that the response Y i j is not fully observed for all i , j. Thus, let ( V i , C i ) be the observed data for the i-th subject, where V i represents the vector of uncensored readings ( V i j = V 0 i ) or censoring interval ( V 1 i j , V 2 i j ) and C i is the vector of censoring indicators, such that: C i j = 1 if V 1 i j ≤ Y i j ≤ V 2 i j , 0 if Y i j = V 0 i , for all i ∈ { 1 , … , n } and j ∈ { 1 , … , n i }, i.e., C i j = 1 if Y i j is located within a specific interval. Note that for a right-censored observation V 2 i j = + ∞ and for a left-censored observation V 1 i j = − ∞. Moreover, missing at random (MAR) observations can be handled by considering V 1 i j = − ∞ and V 2 i j = + ∞. The model defined in Eqs. (3.1)–(3.3) is henceforth called the t-LMEC model.

 Within-subject Dependence Structures

In order to enable some flexibility when modeling the error covariance, we consider essentially three dependence structures: unconditionally independent (UNC), the continuous-type autoregressive process of order p and the damped exponential correlation, which will be discussed next.

Unconditional independence

The most common and simplest approach is to assume that the error terms are Unconditionally independent (UNC), i.e. we have R i = I n i , for each i = 1 , … , n. In practice, the linear mixed models considering UNC errors is very frequent, for example it has been considered by Matos et al. [21]. It will be denoted as UNC-t-LMEC.

However, in longitudinal studies, repeated measures are collected over time and hence the error term might be serially correlated.

In order to account for the within-subject serial correlation, we consider other two general structures.

Autoregressive dependence of order p

In this case, we propose R i as a structured AR(p) dependence matrix [8]. Specifically, R i = R i ( ϕ ) = 1 1 − ϕ 1 ρ 1 − ⋯ − ϕ p ρ p [ ρ | r − s | ] , where r , s = 1 , … , n i and ρ 1 , … , ρ p are the theoretical autocorrelations of the process, and thereby they are functions of autoregressive parameters ϕ = ( ϕ 1 , … , ϕ p ) ⊤ and satisfy the Yule–Walker equations [7], i.e., ρ k = ϕ 1 ρ k − 1 + ⋯ + ϕ p ρ k − p , ρ 0 = 1 , k = 1 , … , p . In addition, the roots of 1 − ϕ 1 B − ϕ 2 B 2 − ⋯ − ϕ p B p = 0 must lie outside the unit circle to ensure stationarity of the AR(p) model.

Following [3], the autoregressive process can be reparameterized using a one-to-one, continuous and differentiable transformation in order to simplify the conditions for stationarity. For details on the estimation of the autoregressive coefficients, we refer to [29], see also [17, 16] for details on estimation of LME model with AR(p) dependence.

The model formulated in Eqs. (3.1)–(3.3) with Ω i = σ 2 R i , where R i is given by Eq. (3.4), i = 1 , … , n, will be denoted AR(p)-t-LMEC. In order to accommodate situations in which measurements are taken irregularly over discrete time, we modify R i by computing it for a regular range of time and then suppressing the line and column regarding the position from the missing measurements.

Damped exponential correlation

In this case, following [25], we propose to structure R i as a damped exponential correlation (DEC) matrix, as follows: R i = R i ( ϕ , t i ) = [ ϕ 1 | t i j − t i k | ϕ 2 ] , 0 ≤ ϕ 1 < 1 , ϕ 2 ≥ 0 , where j , k = 1 , … , n i , for i = 1 , … , n and ϕ = ( ϕ 1 , ϕ 2 ) ⊤ . The correlation parameter ϕ 1 describes the autocorrelation between observations separated by the absolute length of two time points, and the damping parameter ϕ 2 allows the acceleration of the exponential decay of the autocorrelation function defining a continuous-time autoregressive (AR) model. It is important to stress that considering the DEC structure it is possible to obtain different correlation structures. For example:

if ϕ 2 = 0, then R i reduces to the compound symmetry correlation structure (DEC-SYM);

if ϕ 2 = 1, then R i reduces to the DEC-AR(1) correlation structure;

if 0 < ϕ 2 < 1, then R i generates a decay rate slower than the DEC-AR(1) structure;

if ϕ 2 > 1, then R i generates a decay rate faster than the DEC-AR(1) structure; and

if ϕ 2 → ∞, then R i converges to the correlation matrix of a moving-average of order 1 (DEC-MA(1)).

 The model presented in Eqs. (3.1)–(3.3) considering Ω i = σ 2 R i , where R i is given by Eq. (3.5), i = 1 , … , n, will be denoted DEC-t-LMEC.

 The Likelihood Function

Let Y i = ( Y i o , Y i c ) ⊤ where Y i o and Y i c represent the n i o and n i c vectors of observed outcomes and censored observations for subject i, respectively; with n i = n i o + n i c , such that C i j = 0 for all elements in Y i o and 1 for all elements in Y i c . After reordering, Y i , V i , μ i , and Σ i can be partitioned as follows: Y i = vec ( Y i o , Y i c ) , V i = vec ( V i o , V i c ) , μ i ⊤ = ( μ i o , μ i c ) and Σ i = Σ i o o Σ i o c Σ i c o Σ i c c , where vec ( . ) denotes the function which stacks vectors or matrices of the same number of columns.

Using Proposition 1 (Appendix A), we have that Y i o ∼ t n i o ( μ i o , Σ i o o , ν ) , and Y i c | Y i o = y i o ∼ t n i c ( μ i c o , S i c o , ν + n i o ) , where μ i o = X i o β and μ i c = X i c β , μ i c o = μ i c + Σ i c o Σ i o o − 1 ( y i o − μ i o ) , S i c o = ( ν + δ 2 ( y i o ) ν + n i o ) S i , S i = Σ i c c − Σ i c o Σ i o o − 1 Σ i o c and δ 2 ( y i o ) = ( y i o − μ i o ) ⊤ Σ i o o − 1 ( y i o − μ i o ) . Let θ = ( β ⊤ , σ 2 , α ⊤ , ϕ ⊤ , ν ) ⊤ be the parameters vector, as presented by [21], the likelihood function for subject i is given by: L i ( θ ) = f ( y i | θ ) = f ( V i | C i , θ ) = f ( y i o | θ ) P ( V 1 i c ≤ Y i c ≤ V 2 i c | V i o , θ ) = t n i o ( V i o ; μ i o , Σ i o o , ν ) L n i c ( V 1 i c , V 2 i c ; μ i c o , S i c o , ν + n i o ) , where L p ( a , b ; μ , Σ , ν ) is as defined in Eq. (2.2), which can be easily evaluated by using the R package MomTrunc.

The log-likelihood function for the observed data is given by ℓ ( θ | y ) = ∑ i = 1 n log L i ( θ ), and the estimates obtained by maximizing the log-likelihood function ℓ ( θ | y ) are the ML estimates.

 The EM Algorithm

In order to obtain the ML estimation of the parameters in the t-LMEC model, presented in Section 3.1, we implemented the ECM algorithm, proposed by [24]. This algorithm preserves the properties of simplicity and stability of EM algorithm, however replace the M-step with a sequence of conditional maximization (CM) steps. Based in property of multivariate Student’s t-distribution, the t-LMEC model can be expressed in the following hierarchical model: Y i | b i , u i ∼ ind . N n i ( μ i , u i − 1 Ω i ) , b i | u i ∼ ind . N q ( 0 , u i − 1 D ) , u i ∼ iid . G a m m a ( ν 2 , ν 2 ) . Thus, it is possible to apply the ECM algorithm by assuming that Y = ( Y 1 ⊤ , … , Y n ⊤ ), b = ( b 1 ⊤ , … , b n ⊤ ) and u = ( u 1 , … , u n ) ⊤ are hypothetical missing variables, and augmenting with the observed variables ( V , C ) where V = vec ( V 1 , … , V n ) and C = vec ( C 1 , … , C n ).

Considering Y c = ( C ⊤ , V ⊤ , y ⊤ , b ⊤ , u ⊤ ) ⊤ and θ = ( β ⊤ , σ 2 , α ⊤ , ϕ ⊤ , ν ) ⊤ , the complete-data log-likelihood, is given by ℓ c ( θ | y c ) = ∑ i = 1 n − n i 2 log σ 2 − 1 2 log ( | R i | ) − u i 2 σ 2 ( y i − μ i − Z i b i ) ⊤ R i − 1 ( y i − μ i − Z i b i ) − 1 2 log | D | − u i 2 b i ⊤ D − 1 b i + log h ( u i | ν ) + C , where C is a constant that does not depend on the θ. h ( u i | ν ) represents the pdf of a G a m m a ( ν / 2 , ν / 2 ) distribution.

Given the estimate θ = θ ˆ ( k ) , at iteration “ ( k )” of the algorithm, in the E-step computes the conditional expectation of the complete-data log-likelihood function, given by: Q ( θ | θ ˆ ( k ) ) = E [ ℓ c ( θ | y c ) | V , C , θ ˆ ( k ) ] = ∑ i = 1 n Q 1 i ( β , σ 2 , ϕ | θ ˆ ( k ) ) + ∑ i = 1 n Q 2 i ( α | θ ˆ ( k ) ) , where Q 1 i ( β , σ 2 , ϕ | θ ˆ ( k ) ) = − n i 2 log σ 2 − 1 2 log ( | R i | ) − 1 2 σ 2 [ a ˆ i ( k ) − 2 μ i ⊤ R i − 1 ( u i y i ˆ ( k ) − Z i u i b i ˆ ( k ) ) + u i ˆ ( k ) μ i ⊤ R i − 1 μ i ] and Q 2 i ( α | θ ˆ ( k ) ) = − 1 2 log | D | − 1 2 tr ( u i b i b i ⊤ ˆ ( k ) D − 1 ) , with a ˆ i ( k ) = tr ( u i y i y i ⊤ ˆ ( k ) R i − 1 − 2 u i y i b i ⊤ ˆ ( k ) Z i ⊤ R i − 1 + u i b i b i ⊤ ˆ ( k ) Z i ⊤ R i − 1 Z i ) , u i b i ˆ ( k ) = E [ u i b i | V i , C i , θ ˆ ( k ) ] = φ i ( u i y i ˆ ( k ) − u i ˆ ( k ) μ i ) , u i b i b i ⊤ ˆ ( k ) = E [ u i b i b i ⊤ | V i , C i , θ ˆ ( k ) ] = Λ i + φ i ( u i y i y i ⊤ ˆ ( k ) − 2 u i y i ˆ ( k ) μ i + u i ˆ ( k ) μ i μ i ⊤ ) φ i ⊤ , u i y i b i ⊤ ˆ ( k ) = E [ u i y i b i ⊤ | V i , C i , θ ˆ ( k ) ] = ( u i y i y i ⊤ ˆ ( k ) − u i y i ˆ ( k ) μ i ⊤ ) φ i ⊤ , where Λ i = ( D − 1 + Z i ⊤ R i − 1 Z i / σ 2 ) − 1 and φ i = Λ i Z i ⊤ R i − 1 / σ 2 .

In conditional maximization (CM) steps, we update θ ˆ ( k + 1 ) by maximizing Q ( θ | θ ˆ ( k ) ) over θ, as follows: β ˆ ( k + 1 ) = ( ∑ i = 1 n u i ˆ ( k ) X i ⊤ R ˆ i − 1 ( k ) X i ) − 1 ∑ i = 1 n X i ⊤ R ˆ i − 1 ( k ) × ( u i y i ˆ ( k ) − Z i u i b i ˆ ( k ) ) , σ 2 ˆ ( k + 1 ) = 1 N ∑ i = 1 n a ˆ i ( k ) − 2 μ ˆ i ( k + 1 ) ⊤ R i − 1 × ( u i y i ˆ ( k ) − Z i u i b i ˆ ( k ) ) + u i ˆ ( k ) μ ˆ i ( k + 1 ) ⊤ R i − 1 μ ˆ i ( k + 1 ) , D ˆ ( k + 1 ) = 1 n ∑ i = 1 n u i b i b i ⊤ ˆ ( k ) , ϕ ˆ ( k + 1 ) = arg max ϕ ( − 1 2 log ( | R i | ) − 1 2 σ 2 ˆ ( k + 1 ) [ a ˆ i ( k ) − 2 μ ˆ i ( k + 1 ) ⊤ R i − 1 × ( u i y i ˆ ( k ) − Z i u i b i ˆ ( k ) ) + u i ˆ ( k ) μ ˆ i ( k + 1 ) ⊤ R i − 1 μ ˆ i ( k + 1 ) ] ) , ν ˆ ( k + 1 ) = arg max ν { ∑ i = 1 n log L n i c ( V 1 i c , V 2 i c ; μ i c o ( k + 1 ) , S i c o ( k + 1 ) , ν + n i o ) + ∑ i = 1 n log t n i o ( V i o ; μ i o ( k + 1 ) , Σ i o o ( k + 1 ) , ν ) } , where N = ∑ i = 1 n n i .

The algorithm is iterated until a suitable convergence rule is satisfied. In this case, the process is iterated until some distance between two successive evaluations of the currently penalized log-likelihood ℓ p ( θ , λ ), such as | ℓ ( θ ˆ ( k + 1 ) ) / ℓ ( θ ˆ ( k ) ) − 1 | becomes small enough. For example, we adopted ϵ = 10 − 6 .

As proposed by [26], a set of starting values may be achieved by computing β ˆ ( 0 ) , σ 2 ˆ ( 0 ) , D ˆ ( 0 ) and ϕ ˆ ( 0 ) , as the solution of the normal linear mixed-effects model, using the R package nlme.

On the other hand, it is important to stress that, from Eqs. (3.9)–(3.12), the E-step reduces to the computation of u i y i y i ⊤ ˆ = E [ u i y i y i ⊤ | V i , C i , θ ] , u i y i ˆ = E [ u i y i | V i , C i , θ ] and u i ˆ = E [ u i | V i , C i , θ ] .

These expected values can be determined, using the marginal and conditional moments of the first two moments of the TMVT distribution under a double truncation, presented in Propositions 3 and 4 (Appendix B). Thus:

If the i-th subject has only non-censored components, then: u i ˆ = ( ν + n i ν + δ 2 ( y i ) ) , u i y i ˆ = u i ˆ y i , u i y i y i ⊤ ˆ = u i ˆ y i y i ⊤ , where δ 2 ( y i ) = ( y i − μ i ) ⊤ Σ i − 1 ( y i − μ i ).

If the i-th subject has only censored components then, we have: u i ˆ = L n i ( V 1 i , V 2 i ; μ i , Σ i ∗ , ν + 2 ) L n i ( V 1 i , V 2 i ; μ i , Σ i , ν ) , u i y i ˆ = u i ˆ E ( W i ) , u i y i y i ⊤ ˆ = u i ˆ E ( W i W i ⊤ ) , where W ∼ T t n i ( μ i , Σ i ∗ , ν + 2 ; ( V 1 i , V 2 i ) ), μ i = X i β, Σ i ∗ = ν ν + 2 Σ i , Σ i = Z i D Z i ⊤ + Ω i .

If the i-th subject has censored and uncensored components and given that Y i | V i , C i , Y i | V i , C i , Y i o and Y i c | V i , C i , Y i o are equivalent process, then, we have: u i ˆ = ( n i o + ν ν + δ 2 ( y i o ) ) × L n i c ( V 1 i c , V 2 i c ; μ i c o , S ˜ i c o , ν + n i o + 2 ) L n i c ( V 1 i c , V 2 i c ; μ i c o , S i c o , ν + n i o ) , u i y i ˆ = vec ( u i ˆ y i o , u i ˆ E [ W i ] ) , u i y i y i ⊤ ˆ = u i ˆ y i o y i o ⊤ u i ˆ y i o E ⊤ [ W i ] u i ˆ E [ W i ] y i o ⊤ u i ˆ E [ W i W i ⊤ ] , where W i ∼ T t n i c ( μ i c o , S ˜ i c o , ν + n i o + 2 ; ( V 1 i c , V 2 i c ) ), S ˜ i c o = ( ν + δ 2 ( y i o ) ν + n i o + 2 ) S i and S i , S i c o and μ i c o are defined in Section 3.3.

Formulas for E [ W ] and E [ W W ⊤ ], where W ∼ T t p ( μ , Σ , ν ; A ), have been recently developed in the R package relliptical [31], using a slice sampling algorithm with Gibbs sampler steps.

 The Expected Information Matrix

As developed by [21] and [26], we used the [20]’s technique to obtain an asymptotic approximation for the variances of the fixed effects in our t-MLE model. This approximation is given by: I s ( β | y ) = I c ( β | y ) − I m ( β | y ) , where I s ( β | y ) represents the information about β in the observed data y, I c ( β | y ) is the conditional expectation of the complete-data information and I m ( β | y ) is the missing information. Therefore, the approximated covariance-matrix of β ˆ is given by: Cov ˆ ( β ˆ ) ≈ I s − 1 ( β ) | θ ˆ , where I s ( β ) = ∑ i = 1 n ( ν + n i ν + n i + 2 ) X i ⊤ Σ i − 1 X i − X i ⊤ Σ i − 1 ( ( ν + n i + 2 ν + n i ) E 2 − E 1 ) Σ i − 1 X i , where E 1 = ( u i y i ˆ − u i ˆ μ i ) ( u i y i ˆ − u i ˆ μ i ) ⊤ and E 2 = ( u i 2 y i y i ⊤ ˆ − u i 2 y i ˆ μ i ⊤ − μ i u i 2 y i ⊤ ˆ + u i 2 ˆ μ i μ i ⊤ ) . Note that E 1 depend on the computation of u i ˆ, u i y i ˆ that can be obtained in Section 3.4 and E 2 depend on the following quantities u i 2 ˆ = E [ ( ν + n i ν + δ 2 ( y i ) ) 2 | V i , C i , θ ] , u i 2 y i ˆ = E [ ( ν + n i ν + δ 2 ( y i ) ) 2 y i | V i , C i , θ ] and u i 2 y i y i ⊤ ˆ = E [ ( ν + n i ν + δ 2 ( y i ) ) 2 y i y i ⊤ | V i , C i , θ ] . These expected values can be determined in closed form using Propositions 3 and 4, as follows:

If the i-th subject has only non-censored components, then, u i 2 ˆ = ( ν + n i ν + δ 2 ( y i ) ) 2 , u i 2 y i ˆ = u i 2 ˆ y i , u i 2 y i y i ⊤ ˆ = u i 2 ˆ y i y i ⊤ , where δ 2 ( y i ) = ( y i − μ i ) ⊤ Σ i − 1 ( y i − μ i ).

If the i-th subject has only censored components then u i 2 ˆ = c p ( ν , 2 ) L n i ( V 1 i , V 2 i ; μ i , Σ i ∗ , ν + 4 ) L n i ( V 1 i , V 2 i ; μ i , Σ i , ν ) , u i 2 y i ˆ = u i 2 ˆ E [ W i ] , u i 2 y i y i ⊤ ˆ = u i 2 ˆ E [ W i W i ⊤ ] , where c p ( ν , 2 ) = ( n i + ν ) ( ν + 2 ) ν ( n i + ν + 2 ) , W i ∼ T t n i ( μ i , Σ i ∗ , ν + 4; ( V 1 i , V 2 i ) ), Σ i ∗ = ν ν + 4 Σ i , μ i = X i β, Σ i = Z i D Z i ⊤ + Ω i .

If the i-th subject has censored and uncensored components and given that { Y i | V i , C i }, { Y i | V i , C i , Y i o } and { Y i c | V i , C i , Y i o } are equivalent process, we have: u i 2 ˆ = d p ( n i o , ν , 2 ) ( ν + δ 2 ( y i o ) ) 2 × L n i c ( V 1 i c , V 2 i c ; μ i c o , S ˜ i c o , ν + n i o + 4 ) L n i c ( V 1 i c , V 2 i c ; μ i c o , S i c o , ν + n i o ) , u i 2 y i ˆ = vec ( u i 2 ˆ y i o , u i 2 ˆ E [ W i ] ) , u i 2 y i y i ⊤ ˆ = u i 2 ˆ y i o y i o ⊤ u i 2 ˆ y i o E ⊤ [ W i ] u i 2 ˆ E [ W i ] y i o ⊤ u i 2 ˆ E [ W i W i ⊤ ] , where d p ( n i o , ν , 2 ) = ( ν + n i ) ( n i o + ν + 2 ) ( n i o + ν ) n i + ν + 2 , W i ∼ T t n i c ( μ i c o , S ˜ i c o , ν + n i o + 4 ; ( V 1 i c , V 2 i c ) ), S ˜ i c o = ( ν + δ 2 ( y i o ) ν + n i o + 4 ) S i and S i , S i c o and μ i c o are presented in Section 3.3.

Simulation Studies

In this Section, we conduct two simulation studies in order to analyze the performance of our proposed methods. The first simulation study shows the asymptotic behavior of the ML estimators developed in Section 3.4 and the second one examines the consistency of approximate the standard error (SE) presented in Section 3.5. We consider the AR(p)-LMEC model, defined in Eqs. (3.1)–(3.3) with left censored cases and n i = 7 measurements for each subject. As recommended by [26], for both simulation schemes, the initial set of values to generate the sample were β = ( 1 , 2 , 1 ) ⊤ , σ 2 = 0.5. For the dispersion matrices of the random error and random effects we assume an autoregressive dependence of order p = 2 with ϕ = ( 0.6 , − 0.2 ) ⊤ and D = 0.490 0.001 0.001 0.002 , respectively. The design matrices for each individual are given by X i = ( X i 1 , X i 2 , X i 3 ) and Z i = ( 1 7 , Z i 2 ), with X i k = ( X i k 1 , … , X i k 7 ) ⊤ and Z i 2 = ( Z i 21 , … , Z i 27 ) ⊤ where 1 m represents a 1 × m vector of ones. Each vector, X i k j and Z i 2 j are generated independently of a uniform distribution, thus X i k j ∼ U ( 0 , 1 ) and Z i 2 j ∼ U ( 1 , 3 ), i = 1 , … , n, k = 2 , 3 and j = 1 , … , 7. These sets of values are fixed for all the replications. It is important to stress that, for the two simulation, we consider the parameter ν as fixed at 4. Thus, the parameter vector is θ = ( β ⊤ , σ 2 , α ⊤ , ϕ ⊤ ) ⊤ .

Simulation Study 1

We generated R = 100 datasets from the AR(2)-t-LMEC model, considering each one of the combinations, of different sample sizes n ∈ { 50 , 100 , 200 , 400 , 600 } and censoring levels l % ∈ { 0 % , 5 % , 20 % , 40 % }. In each combination, were fitted the AR(2)-t-LMEC and AR(2)-N-LMEC models. Thus, for both models, we computed the bias ( B i a s) and mean squared error ( M S E) of the estimates obtained. These measures are defined by: B i a s ( θ i ˆ ) = 1 R ∑ j = 1 R ( θ i ˆ ( j ) − θ i ) and M S E ( θ i ˆ ) = 1 R ∑ j = 1 R ( θ i ˆ ( j ) − θ i ) 2 , where θ i ˆ ( j ) denotes the ML estimate of θ i , for the j-th replication, j = 1 , … , R.

Simulation study 1. Average Bias of parameter estimates under the AR(2)-t-LMEC model.

Simulation study 1. Average MSE of parameter estimates under the AR(2)-t-LMEC model.

Simulation study 1. Average Bias of parameter estimates under the AR(2)-N-LMEC model.

Simulation study 1. Average MSE of parameter estimates under the AR(2)-N-LMEC model.

From Figures 1–2, we observe as a general rule that, when the distribution’s assumption is correct, the B i a s and M S E tend to zero when the sample size increases, indicating that the ML estimates based on the proposed EM-type algorithm do provide good asymptotic properties, under the t-LMEC model. However, from Figures 3–4, we noted that when the distribution’s assumption is incorrect, the ML estimates do not provide good asymptotic properties; for example, the Bias of ϕ ˆ 1 under the AR(2)-N-LMEC model, departs from zero when the sample size increases.

 Simulation Study 2

In this simulation study, we consider sample sizes n ∈ { 40 , 100 } and generated R = 100 Monte Carlo samples from the AR(2)-N-LMEC and AR(2)-t-LMEC models, with different censoring levels l % ∈ { 0 % , 5 % , 15 % , 25 % , 40 % }. We analyzed the S E for the β = ( β 1 , β 2 , β 3 ) parameter vector (MC-SE), the averages values of the standard errors computed, using the empirical information matrix (MC-IM-SE) described in Section 3.5, and the percentage of times in the R samples, that the 95 % confidence interval contained the true parameter values (MC-COV), assuming asymptotic normality. These measures are defined by: MC-IM-SE ( β i ˆ ) = 1 R ∑ j = 1 R S E ( β i ) ˆ ( j ) and MC-SE ( β i ˆ ) = 1 R − 1 ( ∑ j = 1 R ( β i ˆ ( j ) ) 2 − ( 1 R ∑ j = 1 R β i ˆ ( j ) ) 2 ) , where for j = 1 , … , R, we have that β i ˆ ( j ) and S E ( β i ) ( j ) ˆ represent the ML estimates of parameter β i and the SE estimate of β i ˆ ( j ) , respectively.

Standard errors of parameter estimates (MC-SE), average values of the standard errors (MC-IM-SE) and MC-COV, under t-MLE model with n = 40.

Censoring AR(2)-N-LMEC
Levels Criteria β 1 β 2 β 3
0% MC-IM-SE 0.2106 0.2298 0.2160
MC-SE 0.1648 0.1791 0.1604
MC-COV 0.98 0.97 0.98
5% MC-IM-SE 0.2220 0.2151 0.2066
MC-SE 0.1840 0.1601 0.1770
MC-COV 0.97 0.96 0.96
15% MC-IM-SE 0.2136 0.2155 0.2100
MC-SE 0.1538 0.1661 0.1611
MC-COV 1 0.96 0.97
25% MC-IM-SE 0.2177 0.2161 0.2076
MC-SE 0.2009 0.1708 0.1840
MC-COV 0.96 0.98 0.93
40% MC-IM-SE 0.2569 0.2559 0.2551
MC-SE 0.2050 0.2440 0.2124
MC-COV 0.96 0.95 0.96
Censoring AR(2)-t-LMEC
Levels Criteria β 1 β 2 β 3
0% MC-IM-SE 0.1469 0.1600 0.1467
MC-SE 0.1366 0.1450 0.1338
MC-COV 0.95 0.98 0.96
5% MC-IM-SE 0.1654 0.1575 0.1567
MC-SE 0.1543 0.1418 0.1492
MC-COV 0.94 0.95 0.97
15% MC-IM-SE 0.1834 0.1811 0.1743
MC-SE 0.1443 0.1471 0.1315
MC-COV 0.98 0.97 0.98
25% MC-IM-SE 0.2053 0.1973 0.2004
MC-SE 0.1636 0.1451 0.1439
MC-COV 0.98 0.97 0.96
40% MC-IM-SE 0.2525 0.2545 0.2503
MC-SE 0.1690 0.1774 0.1781
MC-COV 0.98 0.94 0.93

Standard errors of parameter estimates (MC-SE), average values of the standard errors (MC-IM-SE) and MC-COV, under t-MLE model with n = 100.

Censoring AR(2)-N-LMEC
Levels Criteria β 1 β 2 β 3
0% MC-IM-SE 0.1267 0.1296 0.1231
MC-SE 0.1159 0.1240 0.1242
MC-COV 0.94 0.96 0.95
5% MC-IM-SE 0.1144 0.1151 0.1241
MC-SE 0.1040 0.1063 0.1201
MC-COV 0.95 0.97 0.95
15% MC-IM-SE 0.1149 0.1240 0.1208
MC-SE 0.1003 0.1052 0.1062
MC-COV 0.97 0.97 0.99
25% MC-IM-SE 0.1259 0.1257 0.1245
MC-SE 0.1103 0.1315 0.1151
MC-COV 0.99 0.95 0.96
40% MC-IM-SE 0.1446 0.1412 0.1371
MC-SE 0.1144 0.1260 0.1271
MC-COV 0.97 0.98 0.96
Censoring AR(2)-t-LMEC
Levels Criteria β 1 β 2 β 3
0% MC-IM-SE 0.0915 0.0960 0.0913
MC-SE 0.0938 0.0923 0.0882
MC-COV 0.97 0.95 0.96
5% MC-IM-SE 0.1012 0.1006 0.1041
MC-SE 0.0938 0.0962 0.0939
MC-COV 0.96 0.97 0.97
15% MC-IM-SE 0.1096 0.1102 0.1145
MC-SE 0.1003 0.0978 0.0941
MC-COV 0.97 0.97 0.97
25% MC-IM-SE 0.1257 0.1272 0.1249
MC-SE 0.0932 0.1133 0.1021
MC-COV 0.99 0.96 0.96
40% MC-IM-SE 0.1495 0.1478 0.1554
MC-SE 0.0907 0.0968 0.1085
MC-COV 0.98 0.99 0.98

From Tables 1 and 2, we notice that the MC-IM-SE and MC-SE for the β parameter vector, from the t-LMEC model are smaller than the corresponding MC-IM-SE and MC-SE, from the N-LMEC model. Besides, the MC-COV under both models are similar.

 Application

In this Section, we study the viral load data of patients infected with the human immunodeficiency virus type 1 (HIV-1), from the clinical trial “A5055”, previously analyzed by [32]. This clinical trial consists of 44 HIV-1 infected patients that were treated with one of the two powerful antiretroviral therapies: (i) IDV 800 mg and RTV 200 mg, twice a day every 12 hs. and (ii) IDV 400 mg and RTV 400 mg, twice a day every 12 hs. (see [33] for more details). Were measured the number of RNA copies in the blood plasma (viral load) and the immunological marker of T cells of differentiation groups 4 (CD4), at different times of the treat. It is important to note that the viral loads were not necessarily measured at equally spaced intervals of days.

Following [33], we investigate the longitudinal evolution of RNA viral load (in log-base-10 scale) of the 44 patients, during antiretroviral treatment. We focus on analyzing only 42 HIV-1 infected patients; two patients were excluded from the analysis because both of them dropped out very early ( # 4 and # 8 with 3 and 1 measurements, respectively).

The trajectories of transformed RNA viral load along with the detection limit by plotting dotted red lines are shown in Figure 5. The lower detection limit for RNA viral load is 50 copies / ml, thus 106 out of 312 viral load measurements, below the detection limit, are considered to be censoring (approximately 34%). We start by presenting the results obtained by fitting the LMEC model defined in (3.1) and (3.3), as follows: Y i = X i β + Z i b i + ϵ i , with

Y i = ( Y i 1 , … , Y i n i ) ⊤ ,

X i = ( 1 n i ⊤ , t i ⊤ , Treat i × 1 n i ⊤ , CD 4 1 / 2 × 1 n i ⊤ , Treat i × t i ⊤ ),

Z i = ( 1 n i ⊤ , t i ⊤ ), β = ( β 0 , … , β 4 ) ⊤ , b i = ( b i 0 , b i 1 ) ⊤ ,

ϵ i = ( ϵ 1 , … , ϵ n i ) ⊤ , for j = 1 , … , n i and i = 1 , … , 42,

 where

1 m represents a 1 × m vector of ones and t i = ( t i 1 , … , t i n i ).

Y i j is the is the log10-transformation of the viral load of subject “i” at time t i j (days).

Treat i is a treatment indicator (i.e. 0 for the treatments 1 and 1 for the treatments 2).

CD 4 i j 1 / 2 represents the square root of CD4 counts of patient i, measured at time t i j .

The random intercept and random slope are represent by b i 0 and b i 1 , respectively.

 For comparison purposes, we fitted the N-LMEC and the t-LMEC counterpart with the structure dependency errors discussed in Section 3.2. Thus,

Assuming the N-LMEC model.

We made a preliminary analysis using unconditional independence (UNC) in the random effects and random errors, i.e., the model proposed by [30]. The scatter plot of the estimated random intersections ( b 0 effect) and random slopes ( b 1 effect), presented in Figure 6(a), suggests that there is a possible linear relationship between both of them. We can suspect that there is a dependency structure for the random effect. From Figure 6(b), we observe that, apparently, the behavior of the marginal residuals is not totally random, highlighting a trend. Finally, in Figure 6(c)–(e), we show the QQ-plots for the b 0 and b 1 effect and the standard marginal residuals. These figures clearly indicate that the normality assumption for random effects and within-subject errors may be unrealistic. Moreover, the variogram plot of the standard marginal residuals, presented in Figure 7, indicates a long-term autocorrelation, showing a possible intra-individual correlation over time.

Tables 3 presents the results for the “A5055” dataset considering the N-LMEC model, under the different correlation structures. We observed that only the estimation of σ 2 parameters is affected by the choice of the correlation structure.

Assuming the t-LMEC model.

In order to address the possible serial correlation among within-subject errors, we fit several models: i) the unconditionally independent (UNC), ii) the damped exponential correlation (DEC), and iii) the continuous-type autoregressive process (AR) of order p ∈ { 1 , 2 , 3 } cases. It is important to stress that considering the DEC structure it is possible to obtain different correlation structures, for instance, the symmetry (SYM) and DEC-AR(1).

From Table 4, we observed that there are no significant differences in the estimation of the fixed and random effects parameters under the different correlation structures. On the other hand, the estimation of the ν and σ 2 parameters are affected by the choice of the correlation structure.

Note from Figure 8(a)–(b) that the inclusion of a correlation structure improves the model fit. The QQ plots for the b 0 , b 1 effect, and standard marginal residuals, under the t-LMEC model with DEC-AR(1) correlation structure, presented in Figure 8(c)–(e) show a significant improvement, with respect to fit obtained under the N-LMEC model. On the other hand, in Figure 9(a)–(b), we present the simulated envelope plots of Marginal residuals, from fitting the DEC(AR)-N-LMEC and DEC(AR)-t-LMEC models. We observed that some marginal residuals under the DEC(AR)-N-LMEC model are outside the boundary of the envelope, indicating the inadequacy of this model. However, the marginal residuals corresponding to the DEC(AR)-t-LMEC model are within the boundary of the envelope. This shows that although we consider the DEC(AR) correlation structure, a correct specification of the data distribution assumption is essential. Finally, as suggested by [19], we performed the likelihood ratio test to examine the appropriateness of using DEC(AR) dependency structure errors versus UNC errors under the t-LMEC model, obtaining a p-value of 0.0058. As expected, this result indicates the existence of strong dependence among within-patient errors across occasions under the t-LMEC model.

 On the other hand, Table 5 shows the summaries of the log-likelihood, AIC and BIC criteria for t-LMEC and N-LMEC models, respectively, under the different autocorrelation structures. We observe that t-LMEC consistently outperforms the normal counterpart in all cases. In particular, these criteria indicate a preference of the unspecified correlation structure (DEC(AR)), i.e., the estimated correlation structure of the error is given by: Ω ˆ i = σ 2 ˆ R ˆ i = σ 2 ˆ [ ϕ ˆ 1 | t i j − t i k | ] = 0.51 [ 0.899 | t i j − t i k | ] . Thus, by the previous analysis, we conclude that the t-LMEC model with DEC(AR) correlation structure for the error, is the most appropriate model for these dataset.

A5055 Data. Trajectories of log 10 ( R N A ) for 42 HIV-1 infected patients who were randomized to one of two treatment regimens.

A5055 Data. Scatter diagrams of the random effect and marginal errors and Q-Q plots of estimates of random intercepts and slopes and marginal residuals for log10(RNA) using the N-LMEC model with conditional independence in the random effects and the random errors.

A5055 Data. Variogram from the marginal residuals for log10(RNA) using the N-LMEC model with conditional independence in the random effects and the random errors.

A5055 Data. Summary of parameter estimates and the standard errors of fixed effects (in parentheses) under the N-LMEC model, with different correlation structures for the errors.

Correlation structures
Parameter estimates Damped exponential AR(p)
UNC DEC DEC-AR(1) SYM AR(1) AR(2) AR(3)
β 0 4.383 4.379 4.378 4.388 4.384 4.382 4.287
(0.533) (0.002) (0.002) (0.010) (0.540) (0.553) (0.587)
β 1 −0.004 −0.004 −0.004 −0.004 −0.004 −0.004 −0.003
(0.009) (0.008) (0.008) (0.008) (0.009) (0.009) (0.009)
β 2 0.305 0.304 0.303 0.304 0.307 0.308 0.361
(0.309) (0.058) (0.058) (0.013) (0.310) (0.309) (0.316)
β 3 −0.111 −0.111 −0.111 −0.112 −0.111 −0.111 −0.113
(0.031) (0.010) (0.010) (0.008) (0.032) (0.032) (0.034)
β 4 −0.003 −0.003 −0.003 −0.003 −0.003 −0.003 −0.003
(0.005) (0.004) (0.004) (0.004) (0.005) (0.005) (0.006)
σ 2 0.733 0.733 0.732 0.765 0.427 0.417 0.193
(–) (–) (–) (–) (–) (–) (–)
α 11 0.383 0.382 0.381 0.352 0.382 0.384 0.294
(–) (–) (–) (–) (–) (–) (–)
α 12 −0.004 −0.004 −0.004 −0.004 −0.004 −0.004 −0.003
(–) (–) (–) (–) (–) (–) (–)
α 22 >0.001 >0.001 >0.001 >0.001 >0.001 >0.001 >0.001
(–) (–) (–) (–) (–) (–) (–)
ϕ 1 – 0.998 1.000 0.000 0.652 0.762 0.273
(–) (–) (–) (–) (–) (–) (–)
ϕ 2 – 0.106 1 0 – −0.484 0.102
(–) (–) (–) (–) (–) (–) (–)
ϕ 1 – – – – – – −0.820
(–) (–) (–) (–) (–) (–) (–)

A5055 Data. Summary of parameter estimates and the standard errors of fixed effects (in parentheses)under the t-LMEC model, with different correlation structures for the errors.

Correlation structures
Parameter estimates Damped exponential AR(p)
UNC DEC DEC-AR(1) SYM AR(1) AR(2) AR(3)
β 0 3.995 3.998 4.051 4.001 3.988 3.987 4.080
(0.604) (0.509) (0.531) (0.513) (0.607) (0.606) (0.643)
β 1 −0.005 −0.005 −0.004 −0.005 −0.005 −0.005 −0.004
(0.008) (0.007) (0.007) (0.007) (0.008) (0.008) (0.009)
β 2 0.350 0.350 0.385 0.345 0.354 0.355 0.353
(0.309) (0.250) (0.264) (0.253) (0.308) (0.309) (0.323)
β 3 −0.086 −0.086 −0.092 −0.086 −0.086 −0.086 −0.094
(0.020) (0.020) (0.020) (0.020) (0.020) (0.020) (0.021)
β 4 −0.003 −0.003 −0.004 −0.004 −0.004 −0.004 −0.004
(0.005) (0.005) (0.005) (0.005) (0.005) (0.006) (0.006)
σ 2 0.473 0.476 0.510 0.507 0.286 0.282 0.175
(–) (–) (–) (–) (–) (–) (–)
α 11 0.220 0.222 0.210 0.194 0.218 0.218 0.222
(–) (–) (–) (–) (–) (–) (–)
α 12 −0.002 −0.002 −0.002 −0.002 −0.002 −0.002 −0.002
(–) (–) (–) (–) (–) (–) (–)
α 22 >0.001 >0.001 >0.001 >0.001 >0.001 >0.001 >0.001
(–) (–) (–) (–) (–) (–) (–)
ϕ 1 – 0.105 0.899 0.066 0.634 0.731 0.276
(–) (–) (–) (–) (–) (–) (–)
ϕ 2 – 0.998 1 0 – −0.421 0.109
(–) (–) (–) (–) (–) (–) (–)
ϕ 3 – – – – – – −0.789
(–) (–) (–) (–) (–) (–) (–)
ν 3.382 3.407 3.953 3.304 3.247 3.296 4.915
(–) (–) (–) (–) (–) (–) (–)

A5055 Data. Summary of selection criteria for the t-LMEC and N-LMEC model under different correlation structures.

Correlation structures
Model Criteria Damped exponential AR(p)
UNC DEC DEC-AR(1) SYM AR(1) AR(2) AR(3)
t-LMEC loglik −353.840 −353.877 −349.323 −353.696 −353.439 −353.456 −352.021
AIC 729.681 731.755 720.646 729.393 728.878 730.913 730.042
BIC 770.854 776.671 761.819 770.566 770.0517 775.829 778.701
N-LMEC loglik −362.571 −362.643 −362.671 −362.569 −362.343 −362.385 −355.407
AIC 743.142 747.285 745.342 745.139 744.686 746.771 734.815
BIC 776.829 788.458 782.773 782.569 782.116 787.944 779.731

All the computational procedures were realized using a 64-bit Windows environment on a notebook machine with a 1.99 GHz Intel Core i7 processor with 12GB of RAM, the EM algorithm took approximately 9.8 min to converge by each model.

A5055 Data. Scatter diagrams of the random effects and marginal errors and Q-Q plots of estimates of random effects and marginal residuals for log10(RNA) using the DEC(AR)-t-LMEC model.

A5055 Data. Simulated envelope of the Marginal residuals from the fitted models: (a) DEC(AR)-N-LMEC model and (b) DEC(AR)-t-LMEC model.

 Conclusions

In this paper, we proposed a robust linear mixed effect model to analyze longitudinal censored or missing data under the multivariate Student’s t-distribution with within-subject correlation, considering some useful dependence structures. In practical implementation, the EM algorithm is used to obtain the ML estimates of the model parameters. The model proposed in this manuscript is an extension of [21], [22] and [26] by developing some additional tools for robust inferences in practical data analysis.

Two simulation studies were performed in order to evaluate the proposed model. The simulation studies validate our method’s performance and indicate an efficiency gain of the t-LMEC model over the N-LMEC model when data present heavy tails. A real data set, previously analyzed with a conditional independent structure, is analyzed. The results show the t-LMEC model with DEC-AR(1) correlation structure fits the data better. The proposed methods were implemented as part of the new R package ARpLMEC, which is available for download at the CRAN repository. Moreover, we have prepared a Rmarkdown file which contains the codes to reproduce the results of the application and it is available with this manuscript via GitHub https://github.com/hlachos/tlmec/.

Although the t-LMEC model considered here has shown great flexibility for modeling symmetric data with evidence of heavier tails than the normal distributions, its robustness against outliers can be seriously affected by the presence of skewness. Thus, it is of interest to generalize the t-LMEC model by considering a more flexible family of distributions, such as the scale mixtures of skew-normal (SMSN) distributions [14], to accommodate the censoring, skewness and heavy tails simultaneously. It is also of interest to develop effective algorithms for non-linear censored mixed effect (NLMEC) models or semiparametric structures in order to accommodate substantial non-linearity trends observed in the data.

Another possible to explore is to study models that estimate the conditional heteroscedasticity of dataset, as ARCH [12], GARCH [5, 6] models, among others. Finally, another interesting avenue for research is to propose methods within a unified framework in LMEC/NLMEC models for checking random-effects distribution following [1] for generalized linear mixed models and [10] for non-linear mixed models. An in-depth investigation of such extensions is beyond the scope of the present paper, but these are interesting topics for further research.

 <italic>p</italic>-variate Student’s-<italic>t</italic>

The following properties of the p-variate Student’s t-distribution are useful for our theoretical developments. We start with the marginal-conditional decomposition of a p-variate Student’s t random vector. The proof of the following propositions can be found in [2].

Let Y ∼ t p ( μ , Σ , ν ) partitioned as Y = ( Y 1 ⊤ , Y 2 ⊤ ) ⊤ with d i m ( Y 1 ) = p 1 , d i m ( Y 2 ) = p 2 , where p 1 + p 2 = p. Let μ = ( μ 1 ⊤ , μ 2 ⊤ ) ⊤ and Σ = Σ 11 Σ 12 Σ 21 Σ 22 be the corresponding partitions of μ and Σ, respectively. Then, we have

Y 1 ∼ t p 1 ( μ 1 , Σ 11 , ν ); and

The conditional distribution of Y 2 ∣ ( Y 1 = y 1 ) is given by Y 2 ∣ ( Y 1 = Y 1 ) ∼ t p 2 ( μ 2.1 , Σ ˜ 22.1 , ν + p 1 ) , where μ 2.1 = μ 2 + Σ 21 Σ 11 − 1 ( y 1 − μ 1 ) and Σ ˜ 22.1 = ( ν + δ 1 ν + p 1 ) Σ 22.1 with δ 1 = ( y 1 − μ 1 ) ⊤ Σ 11 − 1 ( y 1 − μ 1 ) and Σ 22.1 = Σ 22 − Σ 21 Σ 11 − 1 Σ 12 .

Let Y ∼ t p ( μ , Σ , ν ). Then for any fixed vector b ∈ R m and matrix A ∈ R m × p of full rank we get V = b + A Y ∼ t m ( b + A μ , A Σ A ⊤ , ν ) .

TMVT Distributions

The following propositions are related to the marginal and conditional moments of the first two moments of the TMVT distributions under a double truncation. The proof are similar to those given in [21]. In what follows, we shall use the notation Y ( 0 ) = 1, Y ( 1 ) = Y, Y ( 2 ) = Y Y ⊤ and W ∼ T t p ( μ , Σ , ν ; ( a , b ) ) stands for a p-variate doubly truncated Student’s-t distribution on ( a , b ) ∈ R p .

If Y ∼ T t p ( μ , Σ , ν ; ( a , b ) ) then it follows that E [ ( ν + p ν + δ ( Y ) ) r Y ( k ) ] = c p ( ν , r ) L p ( a , b ; μ , Σ ∗ , ν + 2 r ) L p ( a , b ; μ , Σ , ν ) E [ W ( k ) ] , where c p ( ν , r ) = ( ν + p ν ) r Γ ( p + ν 2 ) Γ ( ν + 2 r 2 ) Γ ( ν 2 ) Γ ( p + ν + 2 r 2 ) , Σ ∗ = ν ν + 2 r Σ and ν + 2 r > 0 with W ∼ T t p ( μ , Σ ∗ , ν + 2 r ; ( a , b ) ).

Having established the formula on the k-order moment of Y, we provide an explicit formula for the conditional moments with respect to a two-component partition of Y.

Let Y ∼ T t p ( μ , Σ , ν ; ( a , b ) ). Consider the partition Y = ( Y 1 ⊤ , Y 2 ⊤ ) ⊤ with dim ( Y 1 ) = p 1 , dim ( Y 2 ) = p 2 , p 1 + p 2 = p, and the corresponding partitions: a = ( a 1 ⊤ , a 2 ⊤ ) ⊤ , b = ( b 1 ⊤ , b 2 ⊤ ) ⊤ , μ = ( μ 1 ⊤ , μ 2 ⊤ ) ⊤ and Σ = Σ 11 Σ 12 Σ 21 Σ 22 . Then, E [ ( ν + p ν + δ ( Y ) ) r Y 2 ( k ) ∣ Y 1 ] = d p ( p 1 , ν , r ) ( ν + δ ( Y 1 ) ) r × L p 2 ( a 2 , b 2 ; μ 2.1 , Σ ˜ 22.1 ∗ , ν + p 1 + 2 r ) L p 2 ( a 2 , b 2 ; μ 2.1 , Σ ˜ 22.1 , ν + p 1 ) × E [ W 2 ( k ) ] , for ν + p 1 + 2 r > 0, with δ ( Y 1 ) = δ ( Y 1 ; μ 1 , Σ 11 ), Σ ˜ 22.1 ∗ = ( ν + δ 1 ν + 2 r + p 1 ) Σ 22.1 and d p ( p 1 , ν , r ) = ( ν + p ) r Γ ( p + ν 2 ) Γ ( p 1 + ν + 2 r 2 ) Γ ( p 1 + ν 2 ) Γ ( p + ν + 2 r 2 ) , where Σ 22.1 is defined as in Proposition 1. Moreover, W 2 ∼ T t p 2 ( μ 2.1 , Σ ˜ 22.1 ∗ , ν + p 1 + 2 r ; ( a 2 , b 2 ) ).

 Acknowledgements

The authors thank the editor and anonymous referees for their valuable comments and suggestions, which significantly contributed to the improvement of this paper.

 References Abad, A. A., Litière, S. and Molenberghs, G. (2010). Testing for misspecification in generalized linear mixed models. Biostatistics 11(4) 771–786. Arellano-Valle, R. B. and Bolfarine, H. (1995). On some characterizations of the t-distribution. Statistics & Probability Letters 25(1) 79–85. https://doi.org/10.1016/0167-7152(94)00208-P. MR1364821 Barndorff-Nielsen, O. and Schou, G. (1973). On the parametrization of autoregressive models by partial autocorrelations. Journal of Multivariate Analysis 3(4) 408–419. https://doi.org/10.1016/0047-259X(73)90030-4. MR0343510 Bartolucci, F. and Farcomeni, A. (2009). A multivariate extension of the dynamic logit model for longitudinal data based on a latent Markov heterogeneity structure. Journal of the American Statistical Association 104(486) 816–831. https://doi.org/10.1198/jasa.2009.0107. MR2751454 Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31(3) 307–327. https://doi.org/10.1016/0304-4076(86)90063-1. MR0853051 Bollerslev, T. (1987). A conditionally heteroskedastic time series model for speculative prices and rates of return. The Review of Economics and Statistics 69(3) 542–547. Box, G. E. P. and Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco. MR0436499 Box, G. E. P., Jenkins, G. M., Reinsel, G. C. and Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control. John Wiley & Sons. MR3379415 Chi, E. M. and Reinsel, G. C. (1989). Models for longitudinal data with random effects and AR(1) errors. Journal of the American Statistical Association 84(406) 452–459. MR1010333 Drikvandi, R. (2020). Nonlinear mixed-effects models with misspecified random-effects distribution. Pharmaceutical Statistics 19(3) 187–201. Drikvandi, R., Verbeke, G. and Molenberghs, G. (2017). Diagnosing misspecification of the random-effects distribution in mixed models. Biometrics 73(1) 63–71. https://doi.org/10.1111/biom.12551. MR3632352 Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica 50(4) 987–1007. https://doi.org/10.2307/1912773. MR0666121 Garay, A. M., Castro, L. M., Leskow, J. and Lachos, V. H. (2017). Censored linear regression models for irregularly observed longitudinal data using the multivariate-t distribution. Statistical Methods in Medical Research 26(2) 542–566. https://doi.org/10.1177/0962280214551191. MR3635923 Lachos, V. H., Ghosh, P. and Arellano-Valle, R. B. (2010). Likelihood based inference for skew-normal independent linear mixed models. Statistica Sinica 20(1) 303–322. MR2640696 Lachos, V. H., A. Matos, L., Castro, L. M. and Chen, M.-H. (2019). Flexible longitudinal linear mixed models for multiple censored responses data. Statistics in Medicine 38(6) 1074–1102. https://doi.org/10.1002/sim.8017. MR3916716 Lin, T. I. and Lee, J. C. (2007). Bayesian analysis of hierarchical linear mixed modeling using the multivariate t distribution. Journal of Statistical Planning and Inference 137(2) 484–495. https://doi.org/10.1016/j.jspi.2005.12.010. MR2298952 Lin, T. I. and Lee, J. C. (2003). On modelling data from degradation sample paths over time. Australian & New Zealand Journal of Statistics 45(3) 257–270. https://doi.org/10.1111/1467-842X.00282. MR1999510 Lin, T. I. and Lee, J. C. (2006). A robust approach to t linear mixed models applied to multiple sclerosis data. Statistics in Medicine 25(8) 1397–1412. https://doi.org/10.1002/sim.2384. MR2226794 Lin, T.-I. and Wang, W.-L. (2020). Multivariate-t linear mixed models with censored responses, intermittent missing values and heavy tails. Statistical Methods in Medical Research 29(5) 1288–1304. https://doi.org/10.1177/0962280219857103. MR4097145 Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological) 44(2) 226–233. MR0676213 Matos, L. A., Prates, M. O., Chen, M. H. and Lachos, V. H. (2013). Likelihood-based inference for mixed-effects models with censored response using the multivariate-t distribution. Statistica Sinica 23 1323–1342. MR3114716 Matos, L. A., Castro, L. M. and Lachos, V. H. (2016). Censored mixed-effects models for irregularly observed repeated measures with applications to HIV viral loads. Test 25(4) 627–653. https://doi.org/10.1007/s11749-016-0486-2. MR3554408 Matos, L. A., Bandyopadhyay, D., Castro, L. M. and Lachos, V. H. (2015). Influence assessment in censored mixed-effects models using the multivariate Student’s t distribution. Journal of Multivariate Analysis 141 104–117. https://doi.org/10.1016/j.jmva.2015.06.014. MR3390061 Meng, X. L. and Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika 81 633–648. https://doi.org/10.1093/biomet/80.2.267. MR1243503 Muñoz, A., Carey, V., Schouten, J. P., Segal, M. and Rosner, B. (1992). A parametric family of correlation structures for the analysis of longitudinal data. Biometrics 48 733–742. Olivari, R. C., Garay, A. M., Lachos, V. H. and Matos, L. A. (2021). Mixed-effects models for censored data with autoregressive errors. Journal of Biopharmaceutical Statistics 31(3) 273–294. Olivari, R. C., Zhong, K., Garay, A. M. and Lachos, V. H. (2022). ARpLMEC: Fitting autoregressive censored mixed-effects models. R package version 2.3. https://CRAN.R-project.org/package=ARpLMEC. Pinheiro, J. C., Liu, C. H. and Wu, Y. N. (2001). Efficient algorithms for robust estimation in linear mixed-effects models using a multivariate t-distribution. Journal of Computational and Graphical Statistics 10 249–276. https://doi.org/10.1198/10618600152628059. MR1939700 Schumacher, F. L., Lachos, V. H. and Dey, D. K. (2017). Censored regression models with autoregressive errors: A likelihood-based perspective. Canadian Journal of Statistics 45(4) 375–392. https://doi.org/10.1002/cjs.11338. MR3729976 Vaida, F. and Liu, L. (2009). Fast implementation for normal mixed effects models with censored response. Journal of Computational and Graphical Statistics 18(4) 797–817. https://doi.org/10.1198/jcgs.2009.07130. MR2750442 Valeriano, K., Matos, L. A. and Morales, C. G. (2022). relliptical: The truncated elliptical family of distributions. R package version 1.1.0. https://CRAN.R-project.org/package=relliptical. Wang, W.-L. (2013). Multivariate t linear mixed models for irregularly observed multiple repeated measures with missing outcomes. Biometrical Journal 55(4) 554–571. https://doi.org/10.1002/bimj.201200001. MR3079990 Wang, W.-L., Lin, T.-I. and Lachos, V. H. (2018). Extending multivariate-t linear mixed models for multiple longitudinal data with censored responses and heavy tails. Statistical Methods in Medical Research 27(1) 48–64. https://doi.org/10.1177/0962280215620229. MR3745654