Modelling multiple outcomes in repeated measures studies: Comparing aesthetic eyelid surgery techniques

Abstract

We propose a multivariate regression model to deal with multiple outcomes along with repeated measures in the context of longitudinal data analysis. Our model allows for flexible and interpretable modelling of the covariance structure within outcomes by using a linear combination of known matrices, while the generalized Kronecker product is employed to take into account the correlation between outcomes. We present maximum likelihood estimation along with extensions of the classical multivariate analysis of variance and multiple comparison hypothesis tests to deal with multivariate longitudinal data. The model and the associated multivariate hypothesis test are motivated by a prospective study conducted to compare three aesthetic eyelid surgery techniques, namely blepharoplasty, endoscopic forehead lift and endoscopic forehead lift associated with blepharoplasty. The effect of the techniques was assessed using measurements of a horizontal line through pupil centre and then three vertical lines, which go in direction to lateral canthus, middle pupil and medial canthus to the top of the brow. In this study, 30 female patients were randomly divided into three groups. Preoperative measurements were compared with postoperative measurements taken 30 days, 90 days and 10 years after the surgery. The presented multivariate model provided a better fit than its univariate counterpart. The results showed that the three surgery techniques tend to increase all considered outcomes in a long-term perspective, that is, from preoperative to 10 years postoperative evaluations. The only exception was for the outcome lateral eyebrow, for which the blepharoplasty had no significant effect.

Keywords

blepharoplasty brow lift frontplasty endoscopic multivariate hypothesis test multiple outcomes

1 Introduction

The periorbital region has received great attention in the field of aesthetic surgery. It can be seen through the evolution of rejuvenation techniques developed for this region (Ramirez (1994); Paul (2001)). On the evaluation of the patient's periorbital region, we should analyse the eyelid excess tissues, such as skin and adipocyte fat, as well as the eyebrow ptosis (McKinney et al., 1991). Based on this analysis, three diagnostics are possible: (a) the patient presents excess of eyelid tissues and eyebrow in the correct position; (b) the patient presents eyebrow ptosis without an excess of eyelid tissues; and (c) the patient presents both eyebrow ptosis and an excess of eyelid tissues. Thus, superior third face surgeries are recommended such as blepharoplasty (Naik et al., 2009), endoscopic forehead lift (Isse, 1994) or even a combination of both.

Often in biomedical studies, multiple outcomes are taken over time for a group of patients. Consequently, the joint analysis of multiple outcomes in biomedical studies has been of increased interest in the medical and statistical literature. Recent contributions include the analysis of ordinal and continuous outcomes (Grigorova and Gueorguieva, 2016), multivariate repeated measures and time to event data in crossover trials (Liu and Li, 2016), joint analysis of longitudinal and survival data (Andrinopoulou and Rizopoulos, 2016; Wen et al., 2016; Hagar et al., 2016) and multivariate disease mapping using linear models of coregionalization (MacNab, 2016; Botella-Rocamora et al., 2015; Martinez-Beneito, 2013).

Similarly to the aforementioned articles, we are interested in the analysis of multiple outcomes in the context of longitudinal data analysis. The study we shall describe in Section 2 was conducted to compare the effect of three aesthetic eyelid surgery techniques, namely blepharoplasty, endoscopic forehead lift and endoscopic forehead lift with blepharoplasty. In this study, 30 female patients were randomly divided into three groups. Thus, each group was assigned 10 patients. To evaluate the effect of each surgery technique, three measurements (medial, central and lateral) were made: one vertical measurement from the top of eyebrow to a horizontal line (through pupil centre) in medial portion, one vertical height from the top of eyebrow to the middle pupil and one from the top of eyebrow to the lateral portion of horizontal pupil centreline. Preoperative measurements were compared to postoperative measurements taken $30$ days, $90$ days and $10$ years after the surgery. These measurements were taken over time characterize the longitudinal structure of the experiment. It is important to highlight that the measurements are not evenly spaced over time. Thus, popular approaches to deal with longitudinal data as the autoregressive models are not suitable for this study (Zeger et al., 1988). Furthermore, three outcomes, lateral measure, middle pupil measure and medial measure are taken at each patient for both right and left eyes, so it is expected that these measurements are also correlated, which in turn characterize the multiple outcomes structure of the experiment.

The main data analysis goal is to assess the effect of the surgery techniques for all outcomes in a long-term perspective, that is, $10$ years after the surgery, which implies that multivariate hypothesis tests are demanded. Furthermore, our goal is to compare all possible contrasts between the measures over time and possibly interactive effects of them with the surgery techniques. Thus, procedures for multiple comparison hypothesis testing and multivariate extensions of them are required to reach such a goal. In the cases in which blepharoplasty is applied only, we expect a little decrease or no difference in the eyebrow height in relation to the pupil. It is due to the resection of the eyelid skin. On the other hand, in cases of endoscopic forehead lift, we expect an increase in the eyebrow height in relation to the pupil, which is due to the reposition and the lift of the frontal region (Graf et al., 2008). Finally, when combining endoscopic forehead lift and blepharoplasty, we expect to keep the increase in the eyebrow height, since the superior vertical vector is stronger than inferior vertical vector (Codner et al., 2010).

The main goal of this article is to propose a multivariate regression model suitable to deal with the combination of repeated measures and multiple outcomes and applied such a model to analyse the dataset described in Section 2. Furthermore, we present maximum likelihood estimation and propose extensions of the classical multivariate analysis of variance (MANOVA) along with multiple comparison techniques to deal with multivariate longitudinal data.

When dealing with multiple outcomes in longitudinal studies, a general issue is how to specify the joint covariance matrix to take into account the within outcome correlation structure induced by the measures taken over time, as well as the correlation introduced by the multiple outcomes. In this article, we adopt the multivariate covariance generalized linear models (McGLM) framework (Bonat and Jørgensen, 2016), which provides flexible and interpretable modelling of the covariance structure. The within outcomes covariance matrix is specified for each marginal outcome using a linear combination of known matrices, while the joint covariance matrix is specified using the generalized Kronecker product (Martinez-Beneito, 2013; Bonat and Jørgensen, 2016). An advantage of this specification is that the univariate counterparts are easily identified as special cases of the multivariate model. Therefore, orthodox measures of goodness-of-fit as the log-likelihood value, Akaike and Bayesian information criteria can be used to compare the fit of the multivariate model with its univariate counterparts.

Given the recent developments in the McGLMs framework, the main contributions of this article are as follows (a) introducing a suitable specification of the McGLMs to deal with multiple outcomes in the context of longitudinal data; (b) presenting maximum likelihood estimation for the model parameters in contrast to the estimating function approach often used in the context of McGLMs; (c) proposing an extension of the standard MANOVA where we can explicitly model the covariance structure to take into account the correlation induced by the longitudinal structure; (d) extending standard multiple comparison test procedures to deal with multiple outcomes and longitudinal data; (e) providing R code for the multivariate tests; and (f) analysing the eyelid dataset.

We present the dataset in Section 2. Section 3 presents the model and discusses its properties and components. Section 4 discusses the maximum likelihood estimation for model parameters. In Section 5, we present extensions of the MANOVA test and multivariate multiple comparison hypothesis tests. In Section 6, we apply the model to the data and present the main results. Finally, the main results are discussed in Section 7, including some directions for future investigations. The dataset and R code are available as supplementary material.

Figure 1:

Three measurements of the vertical height from the lateral canthus, middle pupil and medial canthus to the top of the brow (from ear to nose)

2 Data collection

The case study analysed in this article uses data collected in the Clinic's Hospital of the Paraná Federal University, Curitiba, Brazil. The study was conducted to assess the effect of three eyelid aesthetic surgery techniques, namely blepharoplasty, endoscopic forehead lift and a combination of both techniques. Thirty female patients were randomly assigned into three groups, that is, 10 patients for each surgery technique. A trained surgeon did the surgical procedures during the period from March to December in $2005$ . In order to assess the effect of the techniques, we collected three measurements of vertical lines, which go in direction to lateral canthus, middle pupil and medial canthus to the top of the brow until a horizontal line through pupil centre, as shown in the Figure 1.

The measurements were collected using a digital photo in the preoperative, $30$ days, $90$ days and $10$ years after the surgery. To obtain the measurements, the patients were positioned seated in front of the camera with the chin supported on a tripod and the camera (Nikon® $4.500$ , $3.0$ mega-pixels) was positioned at the same height $42$ cm away from the patient. The photometric parameters of the camera were fitted equally for all patients. The software MIRROR® $6.0$ was used to convert the measurements from pixels to millimeters. Hence, we have four observations per patient for each eye and 30 patients, resulting in the total of $240$ records with $27$ missing data, which were considered as missing at random.

Box plots in Figure 2 suggested for the blepharoplasty cases, there is a slight decrease in the six response variables from the preoperative to the $30$ days postoperative evaluation. However, after that, the values of the three outcomes return to values similar to the ones obtained in the preoperative evaluation. On the other hand, for endoscopic forehead lift and endoscopic forehead lift with blepharoplasty, the outcomes increased from the preoperative to the $30$ and $90$ days postoperative evaluations. Furthermore, the slight decrease from the $90$ days postoperative to the $10$ years postoperative evaluation suggests that the effects of these surgery techniques are lasting. In general, the six response variables present a similar pattern for all surgery techniques considered. Finally, the symmetry of the box plots suggests that the Gaussian distribution is a suitable choice for all outcomes.

Figure 2:

Boxplots by response variables, surgery techniques, evaluation times and eyes

3 Multivariate regression models for longitudinal data

In this section, we shall propose a multivariate regression model for dealing with multiple continuous outcomes in the context of longitudinal studies. Consider the case of $R$ response variables and $N$ observations. Let $Y = (Y_{1}^{⊤}, \dots, Y_{R}^{⊤})^{⊤}$ be the $RN \times 1$ stacked vector of response variables. Similarly, let $μ = (μ_{1}^{⊤}, \dots, μ_{R}^{⊤})^{⊤}$ denote the $RN \times 1$ corresponding stacked vector of expected values and $μ_{r} = X_{r} β_{r}$ denotes the $N \times 1$ vector of expected values for the response variable $r = 1, \dots, R$ . Furthermore, $X_{r}$ denote an $N \times k_{r}$ design matrix composed by the values of $k$ covariates and the index $r$ indicates that we can have different covariates for each response variable. Consequently, $β_{r}$ is a $k_{r} \times 1$ regression parameter vector.

Let $Ω_{r}$ denote the $N \times N$ covariance matrix within the outcome $r$ . Similarly, let $Ω_{b}$ be the $R \times R$ correlation matrix between outcomes. Following the general specification of McGLMs, we propose to model the expected vector and covariance matrix of a multivariate Gaussian distribution by the components

\begin{array}{l} E (Y) = μ \\ var (Y) = Σ = Ω R \begin{matrix} G \\ \otimes \end{matrix} Ω b, \end{array}

(3.1)

where

Ω_{R}_{\otimes}^{G} Ω_{b} = Bdiag ({\tilde{Ω}}_{1}, \dots, {\tilde{Ω}}_{R}) (Ω_{b} \otimes I) Bdiag ({\tilde{Ω}}_{1}^{T}, \dots, {\tilde{Ω}}_{R}^{T})

is the generalized Kronecker product (Martinez-Beneito, 2013). The matrix ${\tilde{Ω}}_{r}$ denotes the lower triangular matrix of the Cholesky decomposition of $Ω_{r}$ . The operator $Bdiag$ denotes a block diagonal matrix and $I$ denotes an $N \times N$ identity matrix.

The key point to deal with multiple response variables in longitudinal studies is the specification of the covariance matrix within outcomes $Ω_{r}$ . Following the ideas of Anderson (1973) and Pourahmadi (2000) and consolidated in Bonat and Jørgensen (2016) and Bonat et al. (2017), we model the within covariance matrix as a linear combination of known matrices, that is,

Ω (τ_{r}) = τ_{r 0} Z_{r 0} + \dots + τ_{rD} Z_{rD},

(3.2)

where $Z_{rd}$ with $d = 0, \dots, D$ are known matrices reflecting the covariance structure within the response variable $r$ , and $τ_{r} = (τ_{r 0}, \dots, τ_{rD})$ is a $(D + 1) \times 1$ parameter vector. Bonat et al. (2017) and Petterle et al. (2019) following the arguments of Demidenko (2013) showed that popular approaches to deal with longitudinal autocorrelated data, as the exchangeable, moving average and first order autoregressive models are covariance linear models, that is, have the form of (3.2).

An important issue that appears when dealing with general covariance models as the one presented in Equation 3.2 is the definition of the parameter space for the parameter vector $τ$ . Unfortunately, there is no easy way to define the parameter space for the $τ$ parameters by restricting each of its components. Thus, we took an alternative approach and define generically such a set as $Θ_{τ} = {τ | Σ > 0}$ where the notation $Σ > 0$ denotes the set of positive definite matrices. In Section 4, we shall present a modification of the Newton scoring algorithm, where we introduce a tuning constant parameter to guarantee that we iterate only over suitable values for the $τ$ components, that is, in $τ$ candidates that the covariance matrix is positive definite. In what follows, we discuss some possibilities for the specification of the matrix linear predictor for the eyelid dataset presented in Section 2.

It is important to highlight that the matrix linear predictor is specified for each response variable. Thus, suppose without loss of generality that $r = 1$ . Denote $y_{ij}$ an observation $j = 1, \dots, J_{i}$ within the patient $i = 1, \dots, I$ . In particular, for the eyelid dataset, we have $i = 1, \dots, 30$ and $j = 1, \dots, 4$ , since we have four records taken over time for each eye of $30$ patients. In this article, since the observations over time are not evenly spaced, we focus on the exchangeable, and inverse of Euclidean distance covariance models to describe the longitudinal structure within response variables. Furthermore, we considered the possibility of variance heterogeneity over time by allowing the variance vary over time in a regression model fashion.

The components of the matrix linear predictor for the outcome $r$ and for a specific patient $i$ are given by

Ω_{r}^{(i)} (τ) = τ_{r 0} [\begin{matrix} 1 & . & . & . \\ . & 1 & . & . \\ . & . & 1 & . \\ . & . & . & 1 \end{matrix}] + τ_{r 1} [\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \end{matrix}],

in the exchangeable case. Similarly, in the inverse of Euclidean distance case, we have

Ω_{r}^{(i)} (τ) = τ_{r 0} [\begin{matrix} 1 & . & . & . \\ . & 1 & . & . \\ . & . & 1 & . \\ . & . & . & 1 \end{matrix}] + τ_{r 1} [\begin{matrix} 0 & 1 / d_{12} & 1 / d_{13} & 1 / d_{14} \\ 1 / d_{12} & 0 & 1 / d_{23} & 1 / d_{24} \\ 1 / d_{13} & 1 / d_{23} & 0 & 1 / d_{34} \\ 1 / d_{14} & 1 / d_{24} & 1 / d_{34} & 0 \end{matrix}],

where $d_{ij}$ denote the Euclidean distance between observations at time $i$ and $j$ .

Finally, to deal with variance heterogeneity in a regression model fashion, we specify the diagonal of the matrix linear predictor components as

Ω_{r}^{(i)} (τ) = τ_{r 0} [\begin{matrix} 1 & . & . & . \\ . & 1 & . & . \\ . & . & 1 & . \\ . & . & . & 1 \end{matrix}] + τ_{r 1} [\begin{matrix} . & . & . & . \\ . & 1 & . & . \\ . & . & . & . \\ . & . & . & . \end{matrix}] τ_{r 2} [\begin{matrix} . & . & . & . \\ . & . & . & . \\ . & . & 1 & . \\ . & . & . & . \end{matrix}] τ_{r 3} [\begin{matrix} . & . & . & . \\ . & . & . & . \\ . & . & . & . \\ . & . & . & 1 \end{matrix}] .

To obtain the components of the matrix linear predictor for the entire response variable vector $Y$ , we assume independent patients. Therefore, the covariance structure within the response variable $r$ is given by

Ω_{r} (τ) = Bdiag (Ω_{r}^{(1)} (τ), \dots, Ω_{r}^{(I)} (τ)),

(3.3)

where the operator $Bdiag$ denotes a block diagonal matrix. It is important to note that the aforementioned structure for the matrix linear predictor can be combined to specify more flexible models.

The correlation matrix $Ω_{b}$ in the eyelid case is given by

Ω_{b} = [\begin{matrix} 1 & ρ_{12} & ρ_{13} & ρ_{14} & ρ_{15} & ρ_{16} \\ ρ_{12} & 1 & ρ_{23} & ρ_{24} & ρ_{25} & ρ_{26} \\ ρ_{13} & ρ_{23} & 1 & ρ_{34} & ρ_{35} & ρ_{36} \\ ρ_{14} & ρ_{24} & ρ_{34} & 1 & ρ_{45} & ρ_{46} \\ ρ_{15} & ρ_{25} & ρ_{35} & ρ_{45} & 1 & ρ_{56} \\ ρ_{16} & ρ_{26} & ρ_{36} & ρ_{46} & ρ_{56} & 1 \end{matrix}],

where denotes the correlation between the outcomes $r$ and $r^{'}$ . An important property of the model is that if we fix the correlation parameters at zero, it corresponds to fit separate models for each response variable. This property allows us to compare the fit of the multivariate model with its univariate counterparts using orthodox measures of goodness-of-fit as the log-likelihood value, Akaike and Bayesian information criteria. The model proposed in this section is a special case of the McGLMs obtained by using identity link function, constant variance function and identity covariance link function and can be easily fitted by the maximum likelihood method, as explained in Section 4. 4 Estimation and inference In this section, we describe the likelihood approach used to estimate the model parameters. We divide the set of parameters into two subsets $θ = (β^{⊤}, λ^{⊤})^{⊤}$ . In this notation, $β = (β_{1}^{⊤}, \dots, β_{R}^{⊤})^{⊤}$ denote a $P \times 1$ vector containing all regression parameters. Similarly, we let $λ = (τ_{1}^{⊤}, \dots, τ_{R}^{⊤}, ρ_{1}, \dots, ρ_{R (R - 1) / 2})^{⊤}$ be a $Q \times 1$ vector of all parameters involved in the covariance structure. We use the convention to stack the correlation parameters $ρ_{{rr}^{'}}$ by columns. For the stacked vector of observed values $y$ , the Gaussian log-likelihood function for $θ$ is given by

L (θ; y) = - \frac{NR}{2} log (2 π) - \frac{1}{2} log | Σ (λ) | - \frac{1}{2} r (β)^{⊤} Σ (λ)^{- 1} r (β),

(4.1)

where $r (β) = (y - μ (β))$ . We use the notation $μ (β)$ and $Σ (λ)$ to emphasize that the vector of expected values and the covariance matrix depending on the regression and covariance parameter vectors, respectively. The maximum likelihood estimator is obtained by the maximization of the log-likelihood function in equation (4.1) with respect to the parameter vector $θ$ . It is easy to show using matrix calculus (Wand, 2002) that the score function for $β$ is given by

U_{β} (θ; y) = D^{⊤} Σ (λ)^{- 1} r (β),

(4.2)

where $D = \frac{\partial μ (β)}{\partial β} = Bdiag (X_{1}, \dots, X_{R})$ . By solving equation (4.2), we obtain a closed-form expression for the maximum likelihood estimator of the regression coefficients, whose general form is

\hat{β} = (D^{⊤} Σ (λ)^{- 1} D)^{- 1} (D^{⊤} Σ (λ)^{- 1} y) .

Similarly, the score function for $λ$ is given by

U_{λ} (θ; y) = - \frac{1}{2} {Σ (λ)^{- 1} - Σ (λ)^{- 1} r (β) r (β)^{⊤} Σ (λ)^{- 1}} \frac{\partial Σ (λ)}{\partial λ} .

(4.3)

By using similar calculations, we can show that the Fisher information matrix for $β$ is given by

F_{β} = D^{⊤} Σ (λ)^{- 1} D .

In a similar way, the $(i, j)$ entry of the Fisher information matrix for the components $λ_{i}$ and $λ_{j}$ is

F_{λ_{i} λ_{j}} = \frac{1}{2} tr (W_{λ_{i}} Σ (λ) W_{λ_{j}} Σ (λ)),

(4.4)

where $W_{λ_{i}} = - \partial Σ (λ)^{- 1} / \partial λ_{i}$ .

The cross-entry $(i, j)$ of the Fisher information matrix between $β$ and $λ$ is given by

F_{β_{i}, λ_{j}} = E [\frac{\partial U_{β} (θ; y)}{\partial λ_{j}}] = 0 .

Let $F_{λ}$ denote the $Q \times Q$ Fisher information matrix of the vector $λ$ . Thus, the joint Fisher information matrix for the parameter vector $θ$ has the form

F_{θ} = (\begin{matrix} F_{β} & 0 \\ 0 & F_{λ} \end{matrix}) .

(4.5)

Let $\hat{θ} = ({\hat{β}}^{⊤}, {\hat{λ}}^{⊤})^{⊤}$ be the maximum likelihood estimator of $θ$ . Thus, we have that the asymptotic distribution of $\hat{θ}$ is

\hat{θ} \sim N (θ, F_{θ}^{- 1}) .

Given the block diagonal structure of equation (4.5), the asymptotic variance of the maximum likelihood estimator of $\hat{β}$ and $\hat{λ}$ are given by $F_{β}^{- 1}$ and $F_{λ}^{- 1}$ respectively. Note that equations (4.3) and (4.4) depend on the derivative of the covariance matrix with respect to the $λ$ parameters. Such derivatives depend on the Cholesky decomposition of $Ω_{r}$ , which although tedious can be done, for details see (Bonat and Jørgensen, 2016). Finally, closed-form expression for the maximum likelihood estimator of $λ$ is in general non-available, thus we employed the following Newton scoring iterative algorithm

λ^{i + 1} = λ^{i} - α F_{λ}^{- 1} U_{λ} (\tilde{θ}; y),

where $\tilde{θ} = ({\hat{β}}^{⊤}, λ^{⊤})^{⊤}$ and $α$ is a tuning parameter introduced to control the step length and guarantee that we iterate only over $τ$ candidates where the covariance matrix is positive definite. The algorithm presented in this section can easily be implemented in R through the mcglm package (Bonat, 2018), whose code is available in the supplementary material.

5 Multivariate hypothesis test for longitudinal data

The classical MANOVA is a huge topic in itself (Smith et al., 1962; Berndt and Savin, 1977). In this section, we explore the Wald test to propose an extension of the classical MANOVA test to deal with multiple responses in the context of longitudinal data. The main difference from our multivariate test to the classical MANOVA is that under our approach, we can model the covariance structure to explicitly take into account the correlation induced by the longitudinal design, while in the classical MANOVA the covariance matrix is completely non-specified, that is, we have to estimate all of its parameters.

In general, under the model in equation (3.1), the general linear hypothesis may be stated as

H_{0} : L β = 0,

(5.1)

where $L = G \otimes F$ . The $G$ matrix has dimension $R \times R$ , while the $F$ matrix has dimension $s \times p$ , where $s$ is the number of linear constraints and $p$ is the number of regression coefficients for each response variable $r$ . In the context of multivariate hypothesis tests, we assume that the linear predictor is composed of the same set of covariates for all response variables. Thus, the $L$ matrix has dimension $(sR \times P)$ . The $G$ matrix is used to state between responses hypotheses. On the other hand, the $F$ matrix is used to state between treatments hypotheses. The alternative hypothesis may be stated in the form,

H_{1} : L β = n,

where $n$ is not the null vector but nothing else is specified about its elements (Smith et al., 1962).

To assess the linear hypothesis specified in equation (5.1), we use the Wald statistics given by

W_{s} = (L β)^{⊤} (L F_{β}^{- 1} L^{⊤})^{- 1} (L β),

which under the null hypothesis is asymptotically chi-square distributed with $sR$ degrees of freedom. It is interesting to note that this simple construction allows us to perform hypothesis test for all response variables as well as between combination of them and all possible contrasts between treatment levels. In this sense, the Wald statistics may be used to perform multivariate multiple hypothesis tests (Bretz et al., 2016).

To give a hint of the approach flexibility, we consider a simple bivariate case and a treatment with three levels ( $A$ , $B$ and $C$ ). In this case, the linear predictor for each response variable is given by

μ_{r} = β_{r 0} + β_{r 1} [X = B] + β_{r 2} [X = C],

where $β_{r 0}$ is the intercept and $β_{r 1}$ and $β_{r 2}$ represent the difference from the treatment $A$ to $B$ and $A$ to $C$ , respectively. In this parametrization, the joint design matrix $D$ has the form,

D = (\begin{matrix} 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 1 & \cdot & \cdot \\ \cdot & \cdot & \cdot & 1 & 1 & \cdot \\ \cdot & \cdot & \cdot & 1 & \cdot & 1 \end{matrix}) .

The null hypothesis of the classical MANOVA test is stated by

H_{0} : β_{r 1} = β_{r 2} = 0, for r = 1, 2 .

Equivalently, in a matrix notation

H_{0} : (\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}) (\begin{matrix} β_{11} \\ β_{12} \\ β_{21} \\ β_{22} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \end{matrix}),

where it is easy to see that

G = (\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}) and F = (\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) .

Given the $L$ matrix, we compute the Wald statistics and perform the test interpretation. It is important to highlight that although the covariance matrix $Σ$ does not appear explicitly in the Wald statistics, it has a crucial impact on the Wald statistics because the Fisher information matrix $F_{β}$ depends on it. Consequently, the presented test has the classical MANOVA as a special case and a lot of new ones obtained by different specifications of the $Ω_{r}$ matrices.

The $F$ matrix can be used to specify contrast between treatment levels. Suppose, we are interested to evaluate the significance of all combinations between the treatment levels. Thus, we have the contrasts $A - B$ , $A - C$ and $B - C$ . The $F$ matrix associated to these contrasts are

F_{AB} = (\begin{matrix} 0 & - 1 & 0 \end{matrix}), F_{AC} = (\begin{matrix} 0 & 0 & - 1 \end{matrix}) and F_{BC} = (\begin{matrix} 0 & 1 & - 1 \end{matrix}) .

Finally, we can use similar ideas for joint testing the covariance parameters. For example, to assess the significance of the longitudinal effect over all response variables.

6 Results

In this section, we apply the model proposed in Section 3 to analyse the eyelid dataset. In this application, for composing the linear predictor we have two covariates surgery techniques (technique) which is a factor with three levels (blepharoplasty, endoscopic forehead lift and endoscopic forehead lift with blepharoplasty) and evaluation times (time) which we consider as a factor with four levels (preoperative, $30$ days, $90$ days and $10$ years postoperative). Thus, the linear predictor for each response variable is composed of the main and interactive effects of technique and time covariates.

To specify the matrix linear predictor for each response variable, we consider the structures: variance heterogeneity, exchangeable and the inverse of Euclidean distance as described in Section 3. For comparison purposes, we also fitted a model considering independent observations. It is obtained by specifying the matrix linear predictor as an identity matrix. Finally, we fitted all the aforementioned models fixing the correlation parameters at zero, which correspond to fit separated models for each outcome. Table 1 presents the maximized log-likelihood values, degrees of freedom (df) along with the Akaike (AIC) and Bayesian (BIC) information criteria for the univariate and multivariate models fitted using each of the within covariance models.

Table 1:
Maximized log-likelihood values (logLik), degrees of freedom (df), Akaike (AIC) and Bayesian (BIC) information criterion by fitted models

Models Univariate Multivariate

logLik df AIC BIC logLik df AIC BIC

Independent $- 1 527.33$ $78$ $3 210.66$ $3 557.42$ $- 1 030.77$ $93$ $2 247.54$ $2 660.99$

Variance het. $- 1 515.13$ $96$ $3 222.26$ $3 649.04$ $- 1 022.25$ $111$ $2 266.50$ $2 759.98$

Euclidean dist. $- 1 411.65$ $84$ $2 991.30$ $3 364.74$ $- 969.20$ $99$ $2 136.40$ $2 576.52$

Exchangeable $- 1 289.83$ $84$ $2 747.66$ $3 121.10$ $- 942.39$ $99$ $2 082.78$ $2 522.90$

Models	Univariate	Multivariate
Independent	$- 1 527.33$	$78$	$3 210.66$	$3 557.42$	$- 1 030.77$	$93$	$2 247.54$	$2 660.99$
Variance het.	$- 1 515.13$	$96$	$3 222.26$	$3 649.04$	$- 1 022.25$	$111$	$2 266.50$	$2 759.98$
Euclidean dist.	$- 1 411.65$	$84$	$2 991.30$	$3 364.74$	$- 969.20$	$99$	$2 136.40$	$2 576.52$
Exchangeable	$- 1 289.83$	$84$	$2 747.66$	$3 121.10$	$- 942.39$	$99$	$2 082.78$	$2 522.90$

The results in Table 1 show that for all measures of goodness-of-fit the models improved from the independent to the inverse of Euclidean distance and exchangeable covariance models in both univariate and multivariate specifications. It shows the relevancy to model the within covariance structure. Variance heterogeneity does not seem concerning. Thus, based on the log-likelihood value we conclude that the exchangeable structure provides the best fit. Such a conclusion is corroborated by the AIC and BIC, which show that the model specified using the exchangeable covariance structure provides the best balance between model complexity and goodness-of-fit.

For all within covariance structures, the multivariate model presented measures of goodness-of-fit better than its univariate counterparts. It highlights the importance of considering the multivariate model as proposed in this article. Hence, we opted to continue the data analysis by analysing the results of the multivariate model fitted using the exchangeable within covariance structure.

The results of the MANOVA test presented in Table 2 show that all effects considered in the model are highly significant.

Table 2:

Wald statistics $(W_{s})$ , degrees of freedom ( $df)$ and $p$ -values; Overall effects

Effects	$df$	$W_{s}$	$p$ -value
Technique	$12$	$36.49$	$< 0.001$
Time	$18$	$78.96$	$< 0.001$
Technique:Time	$36$	$89.57$	$< 0.001$

The MANOVA test showed the overall significance of the effects. Then, we extend the analysis by performing the test for each response variable. Results in Table 3 show that many of the effects composing the linear predictor are highly significant ( $p$ -value $< 0.001$ ), with exception of the main time effect for the response variable lateral canthus in both eyes and the response variable middle pupil for the right eye. However, since the interaction effect between time and technique was significant, we opted to keep this effect in the model.

Table 3:

Wald statistics ( $W_{s}$ ), degrees of freedom (df) and $p$ -values: Response variable specific effects

Effects	df	Lateral canthus		Middle pupil		Medial canthus
Effects	df	$W_{s}$	$p$ -value	$W_{s}$	$p$ -value	$W_{s}$	$p$ -value
Right eye
Technique	$2$	$20.00$	$< 0.01$	$21.92$	$< 0.01$	$10.59$	$< 0.01$
Time	$3$	$0.76$	$0.85$	$5.40$	$0.14$	$33.87$	$< 0.01$
Technique:Time	$6$	$49.75$	$< 0.01$	$51.71$	$< 0.01$	$41.57$	$< 0.01$
Left eye
Technique	$2$	$21.27$	$< 0.01$	$18.43$	$< 0.01$	$9.50$	$0.01$
Time	$3$	$2.93$	$0.40$	$9.11$	$0.03$	$26.66$	$< 0.01$
Technique:Time	$6$	$56.93$	$< 0.01$	$57.87$	$< 0.01$	$40.56$	$< 0.01$

Table 4 shows parameter estimates, standard errors and Z-values for the components of the within covariance structure. The coefficients $τ_{0}$ and $τ_{1}$ are associated with the variance and covariance between observations over time. The exchangeable covariance structure implies that the intraclass correlation does not change over time. Thus, based on the estimates presented in Table 4 the intraclass correlation for the response variables lateral canthus, middle pupil and medial canthus for each eye are presented in Table 5.

Table 4:

Parameter estimates (Est), standard errors (SE) and Z-statistics (Z) for the components of the matrix linear predictor for each response variable

Parameters	Lateral canthus			Middle pupil			Medial canthus
Parameters	Est	SE	Z	Est	SE	Z	Est	SE	Z
Eye right
$τ_{\cdot 0}$	$1.86$	$0.26$	$7.14$	$2.17$	$0.30$	$7.12$	$2.48$	$0.34$	$7.22$
$τ_{\cdot 1}$	$2.54$	$0.45$	$5.62$	$2.63$	$0.48$	$5.52$	$3.44$	$0.61$	$5.67$
Eye left
$τ_{\cdot 0}$	$1.82$	$0.27$	$6.75$	$2.05$	$0.28$	$7.23$	$2.84$	$0.39$	$7.19$
$τ_{\cdot 1}$	$1.53$	$0.31$	$4.88$	$2.25$	$0.40$	$5.64$	$3.62$	$0.66$	$5.46$

Table 5:

Parameter estimates (Est), standard errors (SE) and Z-statistics (Z) for the intraclass correlation

Parameters	Eye	Lateral canthus			Middle pupil			Medial canthus
Parameters	Eye	Est	SE	Z	Est	SE	Z	Est	SE	Z
$ρ_{(times)}$	Right	$0.58$	$0.05$	$11.01$	$0.55$	$0.05$	$11.03$	$0.58$	$0.05$	$11.73$
$ρ_{(times)}$	Left	$0.46$	$0.06$	$7.94$	$0.52$	$0.05$	$10.38$	$0.56$	$0.05$	$10.90$

It is interesting to note that the intraclass correlation is larger in the right eye than in the left eye for all response variables. The estimates and standard errors for the entries of the $Ω_{b}$ matrix were as follows:

{\hat{Ω}}_{b} = [\begin{matrix} 1 \\ 0.85 (0.17) & 1 \\ 0.69 (0.17) & 0.84 (0.19) & 1 \\ 0.68 (0.14) & 0.67 (0.12) & 0.58 (0.13) & 1 \\ 0.61 (0.13) & 0.70 (0.13) & 0.69 (0.15) & 0.86 (0.14) & 1 \\ 0.54 (0.12) & 0.70 (0.13) & 0.83 (0.18) & 0.65 (0.12) & 0.83 (0.13) & 1 \end{matrix}] .

All correlations are significantly different from zero and substantial in size. Furthermore, these correlations take into account the effect of all covariates and the longitudinal structure.

Finally, to address the main goal of this data analysis which is to assess the overall effect of the three eyelid aesthetic surgery techniques Table 6 presents all contrasts between evaluation times and surgery techniques, along with the associated Wald statistics and $p$ -values. As usual in multiple comparison test, the $p$ -values were computed using the Bonferroni correction.

Table 6:

Wald statistics $(W_{s})$ , degrees of freedom ( $df$ ) and $p$ -values; Overall effects

Effects	df	Blepharoplasty		Frontplasty		Frontplasty + Blepharoplasty
Effects	df	$W_{s}$	$p$ -value	$W_{s}$	$p$ -value	$W_{s}$	$p$ -value
Pre–30 days	$6$	$6.36$	$1.00$	$40.68$	$< 0.01$	$37.06$	$< 0.01$
Pre–90 days	$6$	$24.12$	$< 0.01$	$86.00$	$< 0.01$	$112.26$	$< 0.01$
Pre–10 years	$6$	$35.55$	$< 0.01$	$25.13$	$< 0.01$	$58.40$	$< 0.01$
30–90 days	$6$	$38.61$	$< 0.01$	$23.59$	$< 0.01$	$38.03$	$< 0.01$
30–10 years	$6$	$49.07$	$< 0.01$	$5.36$	$1.00$	$18.79$	$0.03$
90–10 years	$6$	$3.07$	$1.00$	$7.72$	$1.00$	$9.07$	$1.00$

Concerning the effects presented in Table 6 is important to note that the differences from preoperative to $90$ days postoperative measures evaluate the surgery's effect, while the differences from $90$ days to $10$ years postoperative measures evaluate whether the effects are lasting. Results in Table 6 show that the three surgery techniques present overall effects, that is, $p$ -values $< 0.01$ for the differences from preoperative to $90$ days postoperative measures while being lasting, that is, no difference from $90$ days to $10$ years postoperative measures.

In order to further explore the results, Tables 7 and 8 presents the estimated differences and the associated standard errors, Z-statistics and $p$ -values for each response variable, for the right and left eyes, respectively. Note that, in this case, we have one degree of freedom for all tests. Thus, the square-root of the Wald statistics is standard Gaussian distributed.

Table 7:

Estimates (Est), standard errors (SE), Z-statistics (Z) and $p$ -values for the difference between evaluation times according to surgery techniques and response variables; Right eye

Difference	Lateral canthus				Middle pupil				Medial canthus
	Est	SE	Z	$p$ -value	Est	SE	Z	$p$ -value	Est	SE	Z	$p$ -value
Blepharoplasty
Pre–30 days	$0.52$	$0.61$	$0.85$	$1.00$	$0.95$	$0.66$	$1.44$	$0.90$	$1.20$	$0.70$	$1.70$	$0.53$
Pre–90 days	$0.30$	$0.63$	$0.48$	$1.00$	$- 0.44$	$0.68$	$- 0.64$	$1.00$	$- 2.08$	$0.73$	$- 2.86$	$0.03$
Pre–10 years	$0.37$	$0.66$	$0.56$	$1.00$	$- 0.40$	$0.71$	$- 0.56$	$1.00$	$- 2.55$	$0.76$	$- 3.36$	$< 0.01$
30–90 days	$- 0.22$	$0.63$	$- 0.34$	$1.00$	$- 1.39$	$0.68$	$- 2.03$	$0.25$	$- 3.28$	$0.73$	$- 4.50$	$< 0.01$
30–10 Years	$- 0.15$	$0.66$	$- 0.23$	$1.00$	$- 1.35$	$0.71$	$- 1.90$	$0.34$	$- 3.75$	$0.76$	$- 4.95$	$< 0.01$
90–10 Years	$0.07$	$0.68$	$0.10$	$1.00$	$0.04$	$0.73$	$0.05$	$1.00$	$- 0.47$	$0.78$	$- 0.60$	$1.00$
Frontplasty endoscopic
Pre–30 days	$- 3.50$	$0.64$	$- 5.44$	$< 0.01$	$- 3.74$	$0.69$	$- 5.39$	$< 0.01$	$- 3.34$	$0.74$	$- 4.50$	$< 0.01$
Pre–90 days	$- 4.74$	$0.64$	$- 7.38$	$< 0.01$	$- 5.86$	$0.69$	$- 8.43$	$< 0.01$	$- 6.51$	$0.74$	$- 8.77$	$< 0.01$
Pre–10 years	$- 2.95$	$0.86$	$- 3.43$	$< 0.01$	$- 3.54$	$0.93$	$- 3.81$	$< 0.01$	$- 4.32$	$0.99$	$- 4.34$	$< 0.01$
30–90 days	$- 1.24$	$0.64$	$- 1.94$	$0.32$	$- 2.11$	$0.69$	$- 3.04$	$0.01$	$- 3.17$	$0.74$	$- 4.27$	$< 0.01$
30–10 Years	$0.55$	$0.86$	$0.64$	$1.00$	$0.20$	$0.93$	$0.22$	$1.00$	$- 0.97$	$0.99$	$- 0.98$	$1.00$
90–10 Years	$1.79$	$0.86$	$2.08$	$0.22$	$2.32$	$0.93$	$2.49$	$0.08$	$2.20$	$0.99$	$2.21$	$0.16$
Frontplasty endoscopic associated with blepharoplasty
Pre–30 days	$- 3.15$	$0.61$	$- 5.16$	$< 0.01$	$- 3.58$	$0.66$	$- 5.43$	$< 0.01$	$- 3.59$	$0.70$	$- 5.10$	$< 0.01$
Pre–90 days	$- 4.78$	$0.63$	$- 7.55$	$< 0.01$	$- 6.01$	$0.68$	$- 8.79$	$< 0.01$	$- 7.45$	$0.73$	$- 10.20$	$< 0.01$
Pre–10 years	$- 3.60$	$0.66$	$- 5.47$	$< 0.01$	$- 4.21$	$0.71$	$- 5.92$	$< 0.01$	$- 5.44$	$0.76$	$- 7.16$	$< 0.01$
30–90 days	$- 1.63$	$0.63$	$- 2.58$	$0.06$	$- 2.43$	$0.68$	$- 3.55$	$< 0.01$	$- 3.86$	$0.73$	$- 5.28$	$< 0.01$
30–10 Years	$- 0.45$	$0.66$	$- 0.68$	$1.00$	$- 0.63$	$0.71$	$- 0.88$	$1.00$	$- 1.85$	$0.76$	$- 2.43$	$0.09$
90–10 Years	$1.18$	$0.67$	$1.77$	$0.46$	$1.80$	$0.72$	$2.50$	$0.08$	$2.01$	$0.77$	$2.61$	$0.05$

Table 8:

Estimates (Est), standard errors (SE), Z-statistics (Z) and $p$ -values for the difference between evaluation times according to surgery techniques and response variables ’ Left eye

Difference	Lateral canthus				Middle pupil				Medial canthus
	Est	SE	Z	$p$ -value	Est	SE	Z	$p$ -value	Est	SE	Z	$p$ -value
Blepharoplasty
Pre–30 days	$0.98$	$0.60$	$1.62$	$0.63$	$1.48$	$0.64$	$2.31$	$0.12$	$1.57$	$0.75$	$2.08$	$0.22$
Pre–90 days	$0.62$	$0.62$	$0.99$	$1.00$	$- 0.23$	$0.66$	$- 0.34$	$1.00$	$- 1.52$	$0.78$	$- 1.95$	$0.31$
Pre–10 Years	$0.80$	$0.65$	$1.23$	$1.00$	$- 0.12$	$0.69$	$- 0.18$	$1.00$	$- 2.18$	$0.81$	$- 2.69$	$0.04$
30–90 days	$- 0.36$	$0.62$	$- 0.58$	$1.00$	$- 1.71$	$0.66$	$- 2.58$	$0.06$	$- 3.09$	$0.78$	$- 3.96$	$< 0.01$
30–10 Years	$- 0.18$	$0.65$	$- 0.28$	$1.00$	$- 1.60$	$0.69$	$- 2.33$	$0.12$	$- 3.75$	$0.81$	$- 4.62$	$< 0.01$
90–10 Years	$0.18$	$0.67$	$0.27$	$1.00$	$0.10$	$0.71$	$0.15$	$1.00$	$- 0.66$	$0.84$	$- 0.79$	$1.00$
Frontplasty endoscopic
Pre–30 days	$- 3.63$	$0.64$	$- 5.71$	$< 0.01$	$- 3.59$	$0.67$	$- 5.32$	$< 0.01$	$- 2.96$	$0.79$	$- 3.72$	$< 0.01$
Pre–90 days	$- 4.61$	$0.64$	$- 7.25$	$< 0.01$	$- 5.37$	$0.67$	$- 7.96$	$< 0.01$	$- 6.38$	$0.79$	$- 8.03$	$< 0.01$
Pre–10 Years	$- 3.57$	$0.85$	$- 4.21$	$< 0.01$	$- 4.05$	$0.90$	$- 4.50$	$< 0.01$	$- 4.14$	$1.06$	$- 3.89$	$< 0.01$
30–90 days	$- 0.98$	$0.64$	$- 1.54$	$0.75$	$- 1.78$	$0.67$	$- 2.64$	$0.05$	$- 3.42$	$0.79$	$- 4.31$	$< 0.01$
30–10 Years	$0.07$	$0.85$	$0.08$	$1.00$	$- 0.46$	$0.90$	$- 0.51$	$1.00$	$- 1.18$	$1.06$	$- 1.11$	$1.00$
90–10 Years	$1.04$	$0.85$	$1.23$	$1.00$	$1.32$	$0.90$	$1.46$	$0.86$	$2.24$	$1.06$	$2.11$	$0.21$
Frontplasty endoscopic associated with blepharoplasty
Pre–30 days	$- 3.33$	$0.60$	$- 5.52$	$< 0.01$	$- 3.31$	$0.64$	$- 5.18$	$< 0.01$	$- 3.44$	$0.75$	$- 4.56$	$< 0.01$
Pre–90 days	$- 4.47$	$0.63$	$- 7.15$	$< 0.01$	$- 5.92$	$0.68$	$- 8.93$	$< 0.01$	$- 7.33$	$0.78$	$- 9.38$	$< 0.01$
Pre–10 Years	$- 3.18$	$0.65$	$- 4.90$	$< 0.01$	$- 4.34$	$0.69$	$- 6.30$	$< 0.01$	$- 5.71$	$0.81$	$- 7.02$	$< 0.01$
30–90 days	$- 1.14$	$0.63$	$- 1.83$	$0.41$	$- 2.61$	$0.66$	$- 3.94$	$< 0.01$	$- 3.89$	$0.78$	$- 4.98$	$< 0.01$
30–10 Years	$0.15$	$0.65$	$0.23$	$1.00$	$- 1.03$	$0.69$	$- 1.49$	$0.81$	$- 2.27$	$0.81$	$- 2.79$	$0.03$
90–10 Years	$1.29$	$0.66$	$1.95$	$0.31$	$1.58$	$0.70$	$2.26$	$0.14$	$1.62$	$0.83$	$1.97$	$0.29$

Regarding the blepharoplasty technique, results in Table 7 show that this technique has no effect on the outcomes Lateral canthus and Middle pupil for both eyes. For the outcome Medial canthus the difference from preoperative to $90$ days postoperative measures was significant ( $p$ -value $= 0.03$ ) for the right eye, however, such a difference was not significant for the left eye ( $p$ -value $= 0.31$ ). Comparing the preoperative measures with $10$ years postoperative measures, we detected a significant effect for both eyes at $5 %$ confidence level.

The dynamic for the techniques endoscopic forehead lift and endoscopic forehead lift with blepharoplasty is similar. In general, we have a significant effect of the surgery techniques, that is, the differences from preoperative to $90$ days postoperative measures are highly significant for all response variables and both eyes, while the differences from $90$ days to $10$ years postoperative measures are not significantly.

7 Discussion

We presented a statistical model to deal with multiple outcomes along with repeated measures in the context of longitudinal data analysis. The model was motivated for an experiment concerning the comparison of three aesthetic surgery techniques in relation to three outcomes, namely, lateral canthus, middle pupil and medial canthus. For modelling the correlation induced by the measures taken over time, we compared the fit of the variance heterogeneity, exchangeable and inverse of the Euclidean distance covariance structures. By using standard measures of goodness-of-fit, we concluded that the multivariate model combined with the exchangeable within covariance structure provided the best balance between goodness-of-fit and complexity. Furthermore, based on the presented multivariate model we extend the classical MANOVA and multiple comparison hypothesis test to deal with multivariate longitudinal data. Although, the model specification is simple and the fitting method amounts to finding the root for a set of non-linear equations, it also involves computational demanding tasks such as the Cholesky decomposition and matrix multiplication. Thus, in the case of large datasets, consisting of many outcomes and long follow-up, the computational cost may prevent routine use of the proposed model.

In our study the measurements were taken based on the horizontal line passing through the pupil centre, not in medial and lateral canthus (eye canthus) because these measurements change when blepharoplasty and endoscopic forehead lift are performed. So, to avoid a bias, a horizontal fixed line through pupil centre was determined. Lateral measurements are our focus in surgery evaluation, since they tend to increase. That is what we look for when surgery is made; especially in cases that endoscopic forehead technique has been used. Our results, showed that both the endoscopic forehead as well as the endoscopic forehead associated with blepharoplasty have a strong and significant effect on the Lateral canthus measures, even in a long-term perspective ( $10$ years after the surgery). Lift forehead techniques have a stronger superior vertical vector, especially at lateral part of eyebrow, larger than inferior vertical vector, which is found in blepharoplasty techniques, in which a larger part of skin is resected at the lateral aspect of the upper eyelid.

In cases of both techniques have been performed (endoscopic forehead and blepharoplasty), a lateral elevation of eyebrow was noted but it tends to be less than the elevation of lateral eyebrow endoscopic forehead surgery. We can explain this difference by a little inferior vertical vector after skin resection in blepharoplasty. Medial and central measurements tend to increase, especially in cases which endoscopic forehead technique has been used (associated or not to blepharoplasty). When forehead lift techniques are applied, a major lateral undermining (ligaments detach) is made, more than in the central part of the forehead. These findings have been noted in all patients, with a high tendency to increase lateral eyebrow elevation in cases in which forehead lift technique have been used. A long-term evaluation showed long lasting results and a lateral eyebrow elevation tendency in cases of endoscopic forehead lift (associated or not to blepharoplasty). When blepharoplasty was performed, as expected, vertical vectors have no significant elevation.

An additional feature of the multivariate model is the possibility to consider the six outcomes in one joint model and thus compute the correlation between outcomes. Although, the correlation between outcomes is not the main interest in our data analysis. It is still important to consider such a structure because it can potentially increase the power of the hypotheses tests associated with the regression coefficients and consequently increase the power of our analysis. Thus, we let as a future work to determine the gain in terms of power obtained by using the multivariate model and further investigate how it can be used to compute the sample size required to reach some determined power in the hypotheses tests, as discussed for example in (Li and McKeague, 2013). Another interesting extension consists of to allow more flexible models to deal with the longitudinal structure that are not linear covariance models as for instance models based on correlation functions as the exponential or spherical and the impact of such structures in the power of the multivariate hypothesis tests.

Supplemental material

Supplementary materials for this article are available from http://www.statmod.org/smij/archive.html and http://www.leg.ufpr.br/doku.php/publications:papercompanions:plastica.

Footnotes

Declaration of conflicting interests

Funding

The authors received no financial support for the research, authorship and/or publication of this article.

References

Anderson

(1973) Asymptotically efficient estimation of covariance matrices with linear structure. The Annals of Statistics , 1, 135–141.

Andrinopoulou

E-R

Rizopoulos

(2016) Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures. Statistics in Medicine , 35, 4813–4823.

Berndt

Savin

(1977) Conflict among criteria for testing hypotheses in the multivariate linear regression model. Econometrica , 45, 1263–1277.

Bonat

(2018) Multiple response regression models in R: The mcglm package. Journal of Statistical Software , 85, 1–30.

Bonat

Jørgensen

(2016) Multivariate covariance generalized linear models. Journal of the Royal Statistical Society: Series C (Applied Statistics) , 65, 649–675.

Bonat

Olivero

Grande-Verga

Farfan

(2017) Modelling the covariance structure in marginal multivariate count models. Journal of Agricultural, Biological and Environmental Statistics , 22, 446–464.

Botella-Rocamora

Martinez-Beneito

Banerjee

(2015) A unifying modeling framework for highly multivariate disease mapping. Statistics in Medicine , 34, 1548–1559.

Bretz

Hothorn

Westfall

(2016). Multiple Comparisons Using R . Boca Raton, FL: CRC Press.

Codner

Kikkawa

Korn

Pacella

(2010) Blepharoplasty and brow lift. Plastic and Reconstructive Surgery , 126, 1–17.

10.

Demidenko

(2013) Mixed Models: Theory and Applications with R . Hoboken, NJ: Wiley.

11.

Graf

Tolazzi

Mansur

Teixeira

(2008) Endoscopic periosteal brow lift: Evaluation and follow-up of eyebrow height. Plastic and Reconstructive Surgery , 121, 609–616.

12.

Grigorova

Gueorguieva

(2016) Correlated probit analysis of repeatedly measured ordinal and continuous outcomes with application to the health and retirement study. Statistics in Medicine , 35, 4202–4225.

13.

Hagar

Harvey

Beckett

(2016) A multivariate cure model for left-censored and right-censored data with application to colorectal cancer screening patterns. Statistics in Medicine , 35, 3347–3367.

14.

Isse

(1994) Endoscopic facial rejuvenation: Endoforehead, the functional lift. Case reports. Aesthetic Plastic Surgery , 18, 21–29.

15.

McKeague

(2013) Power and sample size calculation for generalized estimating equations via local asymptotics. Statistica Sinica , 23, 231–250.

16.

Liu

(2016). A Bayesian model for joint analysis of multivariate repeated measures and time to event data in crossover trials. Statistical Methods in Medical Research , 25, 2180–2192.

17.

MacNab

(2016) Linear models of coregionalization for multivariate lattice data: Order-dependent and order-free cMCARs. Statistical Methods in Medical Research , 25, 1118–1144.

18.

Martinez-Beneito

(2013). A general modelling framework for multivariate disease mapping. Biometrika , 100, 539–553.

19.

McKinney

Mossie

Zukowski

(1991) Criteria for the forehead lift. Aesthetic Plastic Surgery , 15, 141–147.

20.

Naik

Honavar

S G

Das

Desai

Dhepe

(2009) Blepharoplasty: An overview. Journal of Cutaneous and Aesthetic Surgery , 2, 6–11.

21.

Paul

(2001) The evolution of the brow lift in aesthetic plastic surgery. Plastic and Reconstructive Surgery , 108, 1409–1422.

22.

Petterle

Bonat

Scarpin

(2019) Quasi-beta longitudinal regression model applied to water quality index data. Journal of Agricultural, Biological and Environmental Statistics , 24, 346–368.

23.

Pourahmadi

(2000) Maximum likelihood estimation of generalised linear models for multivariate normal covariance matrix. Biometrika , 87, 425–435.

24.

Ramirez

(1994). Endoscopic full facelift. Aesthetic Plastic Surgery , 18, 363–371.

25.

Smith

Gnanadesikan

Hughes

(1962) Multivariate analysis of variance (MANOVA). Biometrics , 18, 22–41.

26.

Wand

(2002) Vector differential calculus in statistics. The American Statistician , 56, 55–62.

27.

Wen

Huang

Frankowski

Cormier

Pisters

(2016) A Bayesian multivariate joint frailty model for disease recurrences and survival. Statistics in Medicine , 35, 4794–4812.

28.

Zeger

Liang

K-Y

Albert

(1988) Models for longitudinal data: A generalized estimating equation approach. Biometrics , 44, 1049–1060.