A robust mixed-effects parametric quantile regression model for continuous proportions: Quantifying the constraints to vitality in cushion plants

There is no literature on outlier-robust parametric mixed-effects quantile regression models for continuous proportion data as an alternative to systematically identifying and eliminating outliers. To fill this gap, we formulate a robust method by extending the recently proposed fixed-effects quantile regression model based on the heavy-tailed Johnson-t distribution for continuous proportion data to the mixed-effects modeling context, using a Bayesian approach. Our proposed method is motivated by and used to model the extreme quantiles of the vitality of cushion plants to provide insights into the ecology of the system in which the plants are dominant. We conducted a simulation study to assess the new method’s performance and robustness to outliers. We show that the new model has good accuracy and confidence interval coverage properties and is remarkably robust to outliers.

rent approach for modeling the quantiles of correlated continuous proportions when outliers are present in the data.

K E Y W O R D S
Bayesian, continuous proportions, cushion plants, mixed-effects, outliers, quantile regression, sub-Antarctic

Need for robust quantile regression models
Traditional regression approaches focus on modeling the relationship between the average response and a set of covariates.As an alternative to modeling the mean response, regression curves can be fitted to selected parts of the response variable's distribution through quantile regression (Cade & Noon, 2003;Koenker & Bassett Jr, 1978), thereby modeling the relationship between the outcome variable and covariates for any portion (e.g., extreme quantiles) of the probability distribution.For example, Wei, Kehm, Goldberg, and Terry (2019)'s study showed that certain risk factors have a limited impact on adult BMI's lower quantiles but are significantly associated with the mean BMI.Thus, quantile regression allows quantifying the strength of the association between the response and the set of covariates that may go unnoticed when only the mean response is considered, as conventional techniques do.Quantile regression can be performed using either parametric or nonparametric methods (Min & Kim, 2004;Yirga, Melesse, Mwambi, & Ayele, 2021).This paper focuses on modeling the quantiles of correlated data bound to the unit interval using random effects.The proposed model is applied to an exemplar dataset from an ecological study.Ideally, one would consider using nonparametric quantile regression methods that do not make any distributional assumptions about the response variable or random effects.However, no nonparametric methods for implementing random-effects models for bounded outcomes are available in the literature.We choose a fully parametric Bayesian approach for the current paper by implementing the models in JAGS (Plummer, 2003).Unlike nonparametric quantile regression models, a significant disadvantage of parametric quantile regression models is that the inference about their parameters can be sensitive to outliers.Our application dataset contains significant outliers; therefore, we particularly focus on models that are robust to outliers.As software, we choose JAGS because it is convenient, user-friendly, relatively fast, and can handle complex models with many parameters (such as random-effects models).Finally, the choice of the Bayesian approach is motivated by the fact that Bayesian inference does not rely on asymptotic theory as frequentist methods generally do, and therefore Bayesian confidence intervals often show better coverage than their frequentist counterparts.
Ecological and clinical research commonly encounters extreme or unusual observations in the data, that is, "outliers" (Benhadi-Marín, 2018).An outlier is a data point that differs considerably from other data points.Such an outlier is called influential if the important features of typical analyses would be altered if they were to be retained or deleted from the dataset (Begashaw & Yohannes, 2020).Several different mechanisms can result in outliers in the data, such as sampling errors (Begashaw & Yohannes, 2020).In contrast, natural outliers that are not due to a sampling error may also emerge.There has been considerable discussion about handling outliers in a variety of research fields (Benhadi-Marín, 2018;Kwak & Kim, 2017;Leys, Klein, Dominicy, & Ley, 2018).One commonly implemented technique is the systematic identification and elimination of outliers.However, deleting outliers may fail to compensate for the uncertainty in the exclusion process and, consequently, may result in underestimated SEs of estimates (Lange, Little, & Taylor, 1989).An alternative to excluding outliers from the data is to employ robust regression techniques that downweigh the influence of outliers on statistical inference.
In the context of parametric regression modeling, one way to achieve robustness to outliers is by assuming heavy-tailed distributions for response variables.For example, replacing the conventional normal distribution for the response with the t-distribution yields inference about the mean outcome robust to outliers (Lange et al., 1989).Alternatively, the modeling of the mean response can be substituted by that of the median, which is considered a more robust measure of central tendency when the data exhibit skewness and contain outliers (Burger & Lesaffre, 2021).However, in parametric quantile regression (Burger & Lesaffre, 2021;Cancho, Bazán, & Dey, 2020), the estimates of the regression coefficients, even for the median, can be prone to outliers if the underlying distribution does not accommodate skewness or heavy tails.Therefore, in the presence of outliers, heavy-tailed distributions can be considered for robust parametric quantile regression modeling, similar to robustly regressing the mean.
For the robust modeling of the average of continuous proportion data, the rectangular beta and flexible beta regression models were proposed by Bayes, Bazán, and García (2012) and Migliorati, di Brisco, & Ongaro, 2018, respectively, to replace the conventional beta model (when outliers are present).di Brisco and Migliorati (2020), for example, modeled EQ-VAS scores, a patient-reported outcome ranging from 0% to 100%, in a Parkinson's disease (PD) longitudinal study.Since the PD dataset of di Brisco and Migliorati (2020) contains outliers (e.g., due to a small group of outliers in EQ-VAS scores near 0), they used the mixed-effects augmented flexible beta model to model the mean EQ-VAS scores over time robustly.
For robust modeling of the quantiles of bounded data, the recently proposed heavy-tailed unit-interval distributions, namely the power normal-logistic distribution (Cancho et al., 2020) and the Johnson-t distribution (Lemonte & Moreno-Arenas, 2020), can serve as alternatives to the Kumaraswamy distribution.The robust logistic quantile regression model of Galarza, Zhang, and Lachos (2020) can also model the quantiles of bounded data in the presence of outliers.

Objectives and outline
To the best of our knowledge, the Kumaraswamy mixed model of Bayes, Bazán, and de Castro (2017) is the only parametric mixed-effects quantile regression approach available in the literature for bounded data.The Kumaraswamy distribution is available in R-INLA (Lindgren & Rue, 2015) and can be implemented as a random-effects model; Flores, Prates, Bazán, and Bolfarine (2021) applied it in a spatial modeling context.However, the Kumaraswamy distribution consists of two parameters and lacks flexibility regarding its tails (mainly, it cannot accommodate heavy tails).Hence, estimates of the regression coefficients from the Kumaraswamy model may be prone to outliers in the data.We note that there is no literature on parametric mixed-effects quantile regression models for bounded data that can accommodate outliers.This paper proposes a parametric mixed-effects quantile regression model for bounded outcomes that is robust to outliers in the data.We extend the recently proposed fixed-effects quantile regression model based on the Johnson-t distribution of Lemonte and Moreno-Arenas (2020) for bounded (continuous proportion) outcomes to the mixed-effects modeling context.Hence, we provide a robust quantile regression method for hierarchically structured data using a mixed-effects approach.We compare the mixed-effects Johnson-t and Kumaraswamy models to assess the suitability of our model.That is, we contrast the current model, namely, the Kumaraswamy model, for hierarchically structured bounded data, with our proposed robust method, namely, the Johnson-t model.We also consider robust and nonrobust models for the mean outcome, respectively, based on the conventional and rectangular beta distributions.We limit our study to the rectangular beta model for robustly modeling the mean (alternatives include the model of Migliorati et al., 2018 (mentioned earlier)).
The paper is organized as follows: Section 2 describes the ecological dataset that motivates our proposed methodology.Section 3 introduces the nonrobust and robust Bayesian mixed-effects models for bounded outcome data, namely the beta, rectangular beta, Kumaraswamy, and Johnson-t models.Section 4 applies the mixed-effects models to the ecological dataset.Section 5 presents simulation studies to investigate the robustness of the Kumaraswamy and Johnson-t models to outliers (data contamination), as well as assess the performance of the Johnson-t model.Finally, Section 6 presents a discussion of the results and findings of the paper.

Cushion plant vitality
Raath-Krüger et al. ( 2022) compiled a long-term dataset of repeated measures to examine the impact of the grass species Agrostis magellanica on the cushion-forming plant, Azorella selago, using sub-Antarctic Marion Island as a model system.These two species are the dominant vascular plants in the sub-Antarctic, and A. magellanica is the most common vascular plant species growing on A. selago (Huntley, 1972).Because of its cushion growth form, A. selago can modify the local microenvironment (e.g., ameliorate temperature conditions (McGeoch, le Roux, Hugo, & Nyakatya, 2008;Nyakatya & McGeoch, 2008)) and therefore has the potential to positively impact species associated with it, particularly in cold, wind-exposed areas where the cushion plant is commonly found.

Objectives
Our proposed statistical method is motivated by modeling the extreme quantiles of the vitality (i.e., dead stem cover) of A. selago, which may provide several key insights into the ecology of the system in which A. selago cushion plants are dominant.The ecological research questions we aim to address are as follows: 1. Given that A. selago vitality could be negatively affected by the cover of A. magellanica and other vascular and nonvascular plant species (see le Roux et al., 2005;Owen, 1995;Raath-Krüger et al., 2022), we ask: Is the limit to A. selago vitality constrained by the cover of vascular and nonvascular plant species growing on A. selago?We model the 0.95th quantile of dead stem cover to address this research question.Here, the 0.95th quantile quantifies the upper extreme quantile of cushion plant vitality, representing the unhealthiest cushion plants.2. Second, given that (i) the western and eastern aspects of Marion island are abiotically different (see, e.g., Goddard, Craig, Schoombie, & le Roux, 2022), resulting in plant populations across the two aspects experiencing different abiotic stressors and (ii) plant species on Marion Island are exposed to increasingly stressful abiotic conditions with increasing altitude (see le Roux, 2008), we ask: Is the limit to A. selago vitality imposed by altitude and aspect?Similar to the above, we model the 0.95th quantile of dead stem cover to address this research question.3. Given that environmental conditions may have different effects on the upper or lower limits of an organism's health (broadly in line with, for example, Wei et al., 2019), we ask: Does the effect of vascular and nonvascular plant species, altitude, and aspect on A. selago vary across plants of different vitality?As a result, we model a broad range of quantiles of dead stem cover, particularly the 0.1th, 0.25th, 0.5th, 0.75th, and 0.95th quantiles, to address this research question.Here, the 0.1th quantile quantifies the vitality of the cushion plants' lower extreme quantile, representing the healthiest cushion plants.
We also model the mean dead stem cover to compare the median and mean fits for completeness sake.Since bounded data tends to be skewed, the median may be a better measure of central tendency than the mean.

Model covariates
We test the long-term constraints on cushion plant vitality by modeling the mean and the quantiles of the final A. selago dead stem cover in response to (i) initial A. selago size, (ii) initial A. magellanica cover on A. selago, (iii) initial combined cover of other vascular plants and mosses on A. selago, (iv) altitude, and (v) aspect.Note that the dead stem cover on A. selago denotes the proportion of black and grey parts on the cushion plant and is therefore considered a continuous proportion outcome.We also include the initial A. selago dead stem cover in the model as a predictor variable because we expect A. selago individuals with greater dead stem cover in the initial year of the survey to have increasingly more dead stem cover in the final year.In order to account for the data's spatial structure, we include a random effect for "plot" in the model.The plots represent specific sampling sites situated on the eastern and western aspects of Marion Island at mid and high altitudes, from which A. selago individuals were surveyed.The data from the low-altitude sites of Raath-Krüger et al. ( 2022) have been excluded since the low and mid altitudes sites exhibited similar patterns.We specify the covariates for modeling the continuous proportion of dead stem cover on Azorella individual j of plot i measured in 2016 as follows: • DS ij , Azo ij , Agr ij , and CO ij denote the covariates corresponding to Azorella individual j of plot i taken in 2003 (i.e., initial measurements): DS ij is the Azorella dead stem cover (%), Azo ij the Azorella size (cm 2 ; log-transformed), Agr ij the Agrostis cover on Azorella (%), and CO ij is the combined cover of other vascular plants and mosses on Azorella (%).
• The indicator variables Mid i and West i are assigned as follows: Mid i = 1 if plot i is situated at a mid-altitude, and Mid i = 0 otherwise (i.e., high altitudes); West i = 1 if plot i is situated on the western side of the island, and West i = 0 otherwise (i.e., eastern side).
The linear predictor corresponding to the final dead stem cover on Azorella individual j of plot i is written as: where  0 ,  1 , … ,  6 are the corresponding regression coefficients (fixed effects), and u i is the random intercept of plot i.

The dataset and outliers
In their original analysis, Raath-Krüger et al. ( 2022) addressed one of their particular research questions by modeling the mean final dead stem cover in response to the covariates above.The data analysis considered the fit of a beta regression model to the final dead stem cover based on the logit link function (i.e., modeling the mean outcome) using the R package glmmTMB (Brooks et al., 2017).Raath-Krüger et al. (2022) excluded A. selago individuals that died (i.e., data points considered as significant outliers due to catastrophic loss of biomass) over the 13-year observation period from their analysis.In contrast, in the present paper, we opt for drawing inferences about the quantiles of dead stem cover for this dataset while accommodating and safeguarding against outliers by employing robust techniques without excluding specific data points (outliers) from the analysis.
Figure 1 shows the final versus initial A. selago dead stem cover by plot, altitude, and aspect, distinguishing between the data included and excluded from the original analysis.For certain plots, it is evident that outliers (relative to the relationship between the final vs. initial dead stem cover) are primarily associated with the previously excluded data.Aside from the previously excluded data due to catastrophic loss, we note that dead stem cover on a few plants (included in the previous analysis of Raath-Krüger et al., 2022) changed considerably during the monitoring period, thus also resulting in data outliers.Figure S1a,b in Section A of Data S1 gives a photograph example of the dead stem cover that remained similar on the vast majority of individual plants over the 13-year monitoring period, not resulting in an outlier.In contrast, Figures S1c,d is an example of an individual plant that yielded an outlier; in this case, the dead stem cover on the individual plant increased from about 1% to 50% during the monitoring period.
The dataset we consider consists of eight plots and 308 cushion plants for the current publication, including the previously excluded data points of the original analysis.Hence, we model the outliers in the current paper instead of excluding them.

Descriptive analysis
Summary statistics of the final dead stem cover are presented in Table 1.The quantiles of the final dead stem cover of the island's western side are higher than those of the eastern side, and the quantiles at mid-altitudes are higher than those at high altitudes.The preliminary investigation of the data, particularly the 0.75th and 0.95th quantiles, suggests that the limit to the cushion plants' vitality may be imposed by altitude and aspect.However, a regression analysis considering the covariates, as mentioned earlier, needs to be performed to make formal conclusions about the constraints to the vitality of cushion plants in this system.

MIXED-EFFECTS MODELS FOR CONTINUOUS PROPORTIONS
This section formulates the robust mixed-effects regression model for the quantiles of bounded responses, namely the Johnson-t model and its nonrobust competitor, the Kumaraswamy model of Bayes et al. (2017).In addition, we consider the rectangular beta model of Bayes et al. (2012) and its nonrobust competitor, namely the conventional beta model of Ferrari and Cribari-Neto (2004).The latter two models fit the mean as a function of the covariates, and we are interested in comparing the results of these models with the results of the two quantile regression models.
For the four models under consideration in this paper (namely beta, rectangular beta, Kumaraswamy, and Johnson-t; see Sections 3.1-3.4),we use the logit link function to model the mean and quantiles as a function of covariates.
Suppose that y ij is the bounded outcome for cluster i = 1, … , I and observation j = 1, … , J i .Let  and u i denote fixed and cluster-specific random effects vectors, and x ij and z ij the covariate vectors, respectively.Assume the u i follow a multivariate normal distribution with mean 0 and d-dimensional unstructured covariance matrix , such that u i | ∼ N d (0, ).
Details on the Bayesian specification of the candidate models are presented in Section 3.5.

Beta regression model
The probability density function of the reparameterized beta regression model of Ferrari and Cribari-Neto (2004) for a given bounded outcome 0 < y ij < 1 is: where ) is the conditional mean of y ij given the ith random effect under the beta model, and  > 0 is the precision parameter of the beta distribution.
It is noted that the beta distribution does not accommodate heavy-tailed data (Bayes et al., 2012), and therefore, the estimation of the mean parameter is prone to data outliers.Moreover, the beta distribution's quantile function is not available in closed form, making the beta distribution not appropriate for parametric quantile regression modeling.

Rectangular beta regression model
The probability density function of the rectangular beta regression model of Bayes et al. (2012) for a given bounded outcome 0 < y ij < 1 is: ) where ) is the conditional mean of y ij given the ith random effect under the x a 1 a 2 −1 (1 − x) a 2 (1−a1)−1 , and  > 0 and 0 ≤  ≤ 1 are, respectively, the precision and shape parameters of the rectangular beta distribution.
The parameter  governs the tail of the rectangular beta distribution: larger values of  yield heavier tails, making the distribution more robust to outliers than the conventional beta distribution (Bayes et al., 2012).The probability density function of the rectangular beta model reduces to that of the conventional beta model when  = 0.
Even though the rectangular beta model accommodates heavy-tailed data and gross outliers, its quantile function cannot be expressed analytically.Therefore, this model is not suitable for quantile regression modeling.

Kumaraswamy regression model
The probability density function of the Kumaraswamy regression model of Bayes et al. (2017) for a given bounded outcome 0 < y ij < 1 is: , ) is the qth conditional quantile of y ij given the ith random effect under the Kumaraswamy model, and  > 0 is the precision parameter of the Kumaraswamy distribution.The Kumaraswamy model is suitable for quantile regression because the quantile function of the Kumaraswamy distribution can be expressed analytically.However, similar to the beta distribution, the Kumaraswamy distribution has only a location and precision parameter but no shape parameter, as the rectangular beta distribution does.Hence, the Kumaraswamy model may be prone to outliers in the data.

Johnson-t regression model
The probability density function of the Johnson-t regression model of Lemonte and Moreno-Arenas (2020) for a given bounded outcome 0 < y ij < 1 is: ) Γ ( where ) is the qth conditional quantile of y ij given the ith random effect under the Johnson-t model,  > 0 and  > 0 are respectively the dispersion parameter and degrees of freedom of the Johnson-t distribution, and t  (q) is the qth quantile of the conventional t-distribution with degrees of freedom . Figure S2 in Section A of Data S1 shows examples of the Johnson-t distribution's probability density function for various values of q and .The parameter  governs the tail of the Johnson-t distribution: smaller values of  yield heavier tails (Lemonte & Moreno-Arenas, 2020).The Johnson-t distribution reduces to the Johnson-normal distribution (Johnson, 1949) for infinite degrees of freedom ( → ∞).Hence, the Johnson-t distribution accommodates heavy-tailed data (i.e., outliers).
The Johnson-t distribution accommodates heavy-tailed data (i.e., outliers) and is suitable for quantile regression given its closed-form quantile function.

Bayesian specification
The prior distributions are specified in such a way as to assure vagueness about prior belief on the model parameters (i.e., weakly informative priors).
A normal prior distribution, namely Normal (0, 10,000), is specified for each component of the vector of fixed effects (i.e., ) for each regression model.
The precision parameter of the beta, rectangular beta, and Kumaraswamy models (i.e., ) and the dispersion parameter of the Johnson-t model (i.e., ) are assigned a gamma prior distribution, namely Gamma (0.0001, 0.0001).
The exponential distribution is often used as a prior distribution for the degrees of freedom of the t-distribution.However, the exponential distribution has a relatively light tail, and therefore, the degrees of freedom's posterior distribution may be unduly influenced by the exponential prior distribution, in some cases, yielding confidence interval coverage very far below the nominal value (Simpson, Rue, Riebler, Martins, & Sørbye, 2017).Therefore, as an alternative to the widely used exponential prior distribution, the Johnson-t distribution's degrees of freedom (i.e., ) are assigned the hierarchical prior distribution of Juárez and Steel (2010), which is more heavy-tailed relative to the exponential and gamma distributions.In particular, the hierarchical prior distribution of  is expressed as a mixture of an exponential distribution with rate parameter 1 for , namely Exp (1), and a gamma distribution with shape parameter 2 and rate parameter , namely Gamma (2, ).From the law of total probability (integrating out ), the resulting probability density function of  is written as P () = 2 (1+) 3 .The shape parameter of the rectangular beta distribution (i.e., ) is assigned a uniform prior distribution, namely Uniform (0, 1).
We specify the matrix-generalized half-t (MGH-t) prior distribution of Huang and Wand (2013) for the variance-covariance matrix (i.e., ) for each model as a more appropriate alternative to the conventional inverse Wishart distribution.The MGH-t prior distribution of  is expressed as a mixture representation of Gamma (0.5, 0.25) for the diagonal entries of diagonal matrix  = diag ( 1 , … ,  z , … ,  d ), and a Wishart distribution with inverse scale matrix 4 and degrees of freedom d + 1, namely Wishart (4, d + 1) (Burger, Schall, Ferreira, & Chen, 2020).This mixture representation results in the specification of the half-t prior distribution with location parameter 0, scale parameter 4, and two degrees of freedom, namely t (0, 4, 2) T (0, ∞), for the SD terms in , and the uniform prior distribution, namely Uniform (−1, 1), for the correlation terms in .From the law of total probability, the set of nuisance parameters  integrated out results in the MGH-t prior distribution, namely: where  > 0, and (  −1 ) zz is the zth diagonal entry of  −1 .The MCMC Gibbs sampling algorithm can draw samples from the joint posterior distribution of the model parameters obtained by forming the product of all likelihoods and prior distributions (Gelfand & Smith, 1990).Software such as JAGS (Plummer, 2003) can be employed to carry out the Gibbs sampling procedure.
For each model, 150,000 samples were simulated from the joint posterior distribution for 15 parallel chains.Among those 150,000 samples (per chain), the initial 10,000 samples were discarded (burn-in).The convergence of posterior samples was checked using trace plots and Brooks-Gelman-Rubin statistics (Brooks & Gelman, 1998).We used a thinning factor of 25 to reduce autocorrelation among the samples.Ultimately, we obtained 84,000 posterior samples in total for each model parameter (hence, K = 84,000).
We reported the posterior distributions' mean and highest posterior density (HPD) intervals as point and interval estimates of the model parameters.The 95% HPD intervals are constructed by finding the shortest interval covering 95% of posterior samples and thus are more informative than the classic symmetric interval.

DATA ANALYSIS
This section examines the fit of the four regression models to the cushion plant data.Comparing the fit of the Kumaraswamy model with that of the Johnson-t model illustrates the impact of outliers on the two model fits.Furthermore, a comparison of the parameter estimates of the beta and rectangular beta distribution, on the one hand, and the Kumaraswamy and Johnson-t distribution, on the other hand, illustrates the difference in modeling the mean compared to the quantiles of the bounded data.We then discriminate between models formally, assess model adequacy, and address the ecological research questions based on the results from the preferred mean and quantile models.

Model implementation
We report regression fits of the mean and the q ∈ {0.1, 0.25, 0.5, 0.75, 0.95} quantiles of dead stem cover to answer the ecological research questions outlined in Section 2. The extreme quantiles of the cushion plant vitality quantified by the 0.1th and 0.95th quantiles of dead stem cover are of primary interest.We fitted the mixed-effects beta, rectangular beta, Kumaraswamy, and Johnson-t regression models in Section 3 to the proportion of dead stem cover in 2016 (see Section 2).We are primarily interested in evaluating the robustness of the Johnson-t and Kumaraswamy models to the outliers present in our dataset.
As per Equation ( 1), the terms in the linear predictor x ′ ij  + z ′ ij u i (Section 3) are defined as follows: ) ′ and z ij = 1 are the covariate vectors, and  = ( 0 ,  1 , … ,  6 ) ′ and u i = u i are respectively the vectors of fixed and random effects.Furthermore,  is the unstructured covariance matrix of u i .Since the model contains a single random effect, its variance component is expressed in scalar rather than vector form (i.e.,  =  2 ).The models were fitted using JAGS (Plummer, 2003) via the R package jagsUI (Kellner, 2021).We ran our models on a desktop computer with a 3.00 GHz Intel® CoreTM i9-10980XE processor and 64 GB installed memory (RAM).

Regression fits and model comparison
This section presents the mean and quantile regression fits derived from the candidate regression models.In addition, we calculated the deviance information criterion (DIC) statistic marginalized over the models' random effects to compare the candidate models, that is, the marginal DIC (mDIC) of Quintero and Lesaffre (2018).The mDIC statistic is a more appropriate model comparison tool than the conventional DIC statistic of Spiegelhalter, Best, Carlin, and van der Linde (2002) (i.e., conditional on the random effects) when population-average inferences are of interest, which is the case for our application.Details on the calculation of the mDIC are presented in Appendix A. We first compare the fits of the nonrobust and robust quantile models (i.e., the Kumaraswamy and Johnson-t models).We do this by studying the posterior estimates (PEs) and 95% highest posterior density (HPD) intervals for the regression coefficients.The quantile regression coefficients' PEs and 95% HPD intervals calculated from the Kumaraswamy and Johnson-t models are presented in Figures 2 and S3 in Section A of Data S1.In addition, the complete set of quantile model parameters, including the model comparison statistic, are presented in Table S1 in Section B of Data S1.From Figures 2 and S3, and Table S1, we observe the following: • The mDIC strongly favors the robust quantile model over the nonrobust model.
• The degrees of freedom estimate under the Johnson-t model confirms our finding that the data are heavy-tailed ( ν ≈ 3), implying that the robust model better fits the data than the nonrobust model.
• The regression coefficients' 95% HPD intervals under the nonrobust quantile model are generally wider than those under the robust model.Therefore, the robust model produces shorter 95% HPD intervals for the regression coefficients.• The effect of the covariates on the quantiles generally differs considerably between the robust and nonrobust models.Under the robust model, the effect of the covariates on the range of quantiles is similar.In contrast, the covariate effects under the nonrobust model differ considerably among the range of quantiles.
Secondly, we compare the fits of the mean and quantile models; we restrict the comparison to the robust models (i.e., the rectangular beta and Johnson-t models).Figure S4 presents the mean and quantile regression coefficients' PEs and 95% HPD intervals calculated from the rectangular beta and Johnson-t models.From Figure S4, we observe that the effect of some covariates on the mean differs somewhat from that of the quantiles.
Thirdly, we compare the fits of the nonrobust and robust mean and median models (i.e., the beta and rectangular beta models for the mean and the Kumaraswamy and Johnson-t models for the median).The PEs and 95% HPD intervals for the regression coefficients of the mean and median computed from the beta, rectangular beta, Kumaraswamy, and Johnson-t models are presented in Figure S5.In addition, the complete set of mean and median model parameters, including the model comparison statistic, are presented in Table S2.From Figure S5 and Table S2, we observe that the robust mean and median models are favored over the nonrobust models and that the PEs and 95% HPD intervals under the nonrobust and robust models differ somewhat.The estimates and confidence intervals of the two central tendency measures' regression coefficients (i.e., the mean and the median) are somewhat different, probably due to skewness in the data.

Model adequacy
We determine the influence of observations on model fits using leave-one-out cross-validation to assess model adequacy (Wang & Luo, 2016): in particular, we evaluate the effect of a data point when absent and present in the dataset using the Kullback-Leibler (K-L) divergence (see Section B.1 of Appendix B).We also assess the model residuals and empirical predictive coverage as the goodness of fit measures based on the posterior predictive distribution (see Sections B.2 and B.3 of Appendix B). Figure S6 in Section A of Data S1 presents the K-L divergence estimates calculated from the beta and rectangular beta models, whereas Figure 3 presents those from the Kumaraswamy and Johnson-t models.Many observations considerably affect the estimates of the regression coefficients of mean and quantiles under the nonrobust models (influential outliers), whereas fewer observations significantly affect these estimates under the robust models.
Figures S7 and S8 show the residual diagnostics calculated from the Kumaraswamy and Johnson-t models.Under the Kumaraswamy model, the scaled residuals do not vary uniformly between 0 and 1 for all quantiles, implying that the model does not fit the data well.Furthermore, the quantile regression fits of the scaled residuals against the corresponding fitted values deviate considerably from the expected quantiles.In contrast, the residual diagnostics suggest that the Johnson-t model fits the data well.
Figure S9 presents the candidate models' empirical predictive coverage and the average length of predictive intervals.The coverage of the predictive intervals from the robust mean and quantile models is close to the nominal value, whereas the predictive coverage under the nonrobust models is too high.The robust mean and quantile models yield narrower predictive intervals than the nonrobust models.
Per the current statistical software implementing the Kumaraswamy model, R-INLA (Lindgren & Rue, 2015), we checked the model fits using predictive integral transforms (PITs).The PITs presented in Figure S10 suggest that the Kumaraswamy model fits the data poorly for all investigated quantiles (potentially due to outliers).Therefore, the Kumaraswamy model is not adequate for our dataset.

Ecological findings
Since the Johnson-t model clearly performs better than the Kumaraswamy model, as shown in the previous section, we quantify the constraints to vitality in cushion plants according to the Johnson-t model.We exponentiate the regression coefficients and interpret them as odds ratios, analogous to logistic regression modeling of binomial data.a Predictor: The mean' odds ratios are derived from the rectangular beta distribution.b The odds ratios are calculated according to a 5% increase in the initial cover of vascular and nonvascular plant species growing on the cushion plants (i.e.,  3 and  4 ).c The quantiles' odds ratios are derived from the Johnson-t distribution.
and quantile models' odds ratio estimates, that is, exp , calculated from the rectangular beta and Johnson-t models.The estimates of the odds ratios suggest the following: 1.A 5% increase in the initial Agrostis magellanica cover results in a 2% increase in the odds of all quantiles of final dead stem cover (see  3 ).Similarly, a 5% increase in the initial cover of other vascular plants and mosses results in an approximately 5% decrease in the odds of all quantiles of final dead stem cover (see  4 ).Therefore, despite evidence for A. selago altering the population structure, biomass, reproductive output, cover, and abundance of A. magellanica and other vascular and nonvascular plants compared to surrounding areas where A. selago is absent (le Roux, Shaw, & Chown, 2013;Raath-Krüger et al., 2021), the extreme quantiles of vitality in A. selago are not constrained by the cover of vascular plants and mosses.2. The odds of all quantiles of final dead stem cover are about 64% higher for mid-altitude than high-altitude sites (see  5 ).Similarly, the odds of all quantiles of final dead stem cover are about 57% higher for the western aspect of Marion Island than the eastern aspect (see  6 ).Therefore, the limit to A. selago vitality is constrained by altitude and aspect: the odds ratios reveal that, across all quantiles, dead stem cover is higher at the mid-altitude than at high-altitude sites and on the western aspect of the island than the eastern aspect.However, the difference in the extreme quantiles of dead stem cover between mid versus high-altitude sites may depend on A. magellanica cover, and therefore, additional covariates may be necessary for the model to account for such interactions.3. The effect of vascular and nonvascular plant species, altitude, and aspect on A. selago does not vary across plants of different vitality.Therefore, it suggests that neither upper nor lower A. selago vitality is constrained by environmental conditions (altitude, aspect, and plant cover) differently from the median vitality, a robust central tendency measure.

Model performance
We assessed the performance of the Johnson-t regression model in Section 3 in a simulation study under two distributional assumptions, namely data generated from the Johnson-t distribution with heavy and light tails, each under three quantile levels, namely the first decile, the median, and the 0.95th quantile (i.e., q ∈ {0.1, 0.5, 0.95}; hence a total of six parameter scenarios).The parameter values were chosen to partly reflect the ecology study's typical data in Section 2. For all scenarios, the sample size and the values of the fixed effects, covariance matrix, and dispersion parameter were chosen as I = 10, J = 5, , and  = 1.5.Note that the model is specified as a random intercept-slope mixed-effects model instead of a random intercepts model as in our application.The degrees of freedom were chosen as (i)  = 5 (i.e., the heavy-tailed scenario) and (ii)  = 30 (i.e., the light-tailed scenario).Hence, the following six scenarios were considered in total: (i) q = 0.1 and  = 5, (ii) q = 0.1 and  = 30, (iii) q = 0.5 and  = 5, (iv) q = 0.5 and  = 30, (v) q = 0.95 and  = 5, (vi) q = 0.95 and  = 30.We simulated the t ij from a unit interval uniform distribution, namely Uniform (0, 1).We fitted the Johnson-t model to 500 simulated datasets for each of the six scenarios.
The bias of the PEs was calculated (see details in Appendix C), as was the coverage probability of the associated 95% HPD intervals.We used the autorun.jagsfunction of the runjags package (Denwood, 2016) to guarantee convergence of the posterior samples for each fitted dataset.The corresponding results are presented in Table 3. From Table 3, we observe the following: • Under each parameter scenario, the bias of the estimates of the fixed effects (i.e.,  0 ,  1 ) and of the dispersion parameter (i.e., ) is small.In contrast, the bias of the estimates of the variance components (i.e.,  2 1 ,  2 2 ,  12 ) is large.• The estimates of the degrees of freedom (i.e., ) are relatively unbiased for low degrees of freedom and considerably biased for high degrees of freedom.
• Under each parameter scenario, the coverage probability of the HPD interval of the fixed effects (i.e.,  0 ,  1 ) is close to or slightly higher than the nominal value.In contrast, the HPD intervals of the dispersion parameter (i.e., ) and variance components (i.e.,  2 1 ,  2 2 ,  12 ) are too conservative.
• When the data are heavy-tailed, the coverage probability of the HPD interval of the degrees of freedom (i.e., ) is slightly higher than the nominal value.In contrast, when the data are light-tailed, the degrees of freedom's HPD interval coverage probability is somewhat lower than the nominal value.
In summary, the estimates and HPD interval coverage of the fixed effects, which are the main parameters of interest, are acceptable under all parameter scenarios.However, the HPD intervals of the variance components are generally too conservative, similar to the findings of Burger et al. (2020) and Burger and Lesaffre (2021), in particular, under small variance components.The HPD interval coverage of the degrees of freedom is slightly lenient when the data are light-tailed; in contrast, it is acceptable when the data are heavy-tailed.

Data contamination
We performed a simulation study to investigate the robustness of the Kumaraswamy and Johnson-t regression models to outliers.Datasets were simulated from the two models where we chose the model parameter values for I, J, x ij , z ij ,  0 ,  1 ,  2 1 ,  2 2 , and  12 as per Section 5.1.We chose  = 30, and  = 3 and  = 1.5 for data simulated from the Kumaraswamy and Johnson-t models, respectively.The selection of the degrees of freedom, precision, and dispersion parameter values yields comparable datasets, ensuring a sensible comparison between the two candidate models.The datasets were randomly contaminated by replacing y ij with data simulated from the Uniform (0.975, 1) distribution at a rate of 5%, thereby contaminating the data with outliers close to the upper bound of the parameter space (i.e., 1).Hence, our contamination scheme is based on a mixture of the Kumaraswamy/Johnson-t distribution and the uniform distribution.The candidate models were fitted to both the uncontaminated and contaminated versions of the simulated datasets.Under each scenario, we considered two quantile levels: the median and the 0.95th quantile (i.e., q ∈ {0.5, 0.95}; hence a total of four parameter scenarios for each model).We fitted the models to 500 simulated datasets for each scenario.
The bias and root mean square error (RMSE) of the PEs were calculated (see details in Appendix C), as were the average length and empirical coverage probability of the associated 95% HPD intervals.The corresponding results are presented in Table 4 (fixed effects only).From Table 4, we observe the following: • Both models perform well under no contamination, as judged by the accuracy characteristic (i.e., bias) and HPD interval coverage.• Under "contamination" relative to "no contamination": -The Kumaraswamy model's median coefficient estimates are more biased and less precise, most notably the fixed intercept term; in contrast, the bias under the Johnson-t model is small.-The 0.95th quantile's coefficient estimates are considerably more biased and less precise under both models.However, the lack of precision of the estimates of the intercept term is more prominent for the Kumaraswamy model.-The HPD interval of the median slope term is extremely conservative under the Kumaraswamy model, whereas the coverage of the HPD interval of the median slope term from the Johnson-t model is acceptable.-The 0.95th quantile's HPD interval of the intercept term is extremely lenient under the Kumaraswamy model; in contrast, the HPD interval's lack of coverage is considerably less problematic under the Johnson-t model.-The increase in the HPD intervals' average length under the Johnson-t model is generally considerably less than that of the Kumaraswamy model.
In summary, the data contamination simulation study suggests that the Johnson-t model is considerably more robust to outliers than the Kumaraswamy model.

DISCUSSION
The currently available approach for modeling the quantiles of hierarchically structured continuous proportion data is the Kumaraswamy model of Bayes et al. (2017).However, our application dataset contains significant outliers, and the Kumaraswamy model is not adequate for modeling heavy-tailed data (i.e., containing outliers).We, therefore, considered a robust model to accommodate outliers by extending the fixed effects Johnson-t model of Lemonte and Moreno-Arenas (2020).According to the mDIC statistic, the robust models (i.e., the rectangular beta and Johnson-t models) fit the cushion plant dataset better than the nonrobust models (i.e., the beta and Kumaraswamy models).Furthermore, the quantile model fits differ considerably between the nonrobust and robust models.Under the robust models, the covariate effects on the mean and median differ somewhat; therefore, one should carefully consider which measure of central tendency to report (i.e., mean vs. median) for bounded data.Based on the model adequacy checks, it is clear that the outliers in the data have a substantial effect on the quantile fits, and thus, the ecological research questions were addressed using the Johnson-t model instead of the Kumaraswamy model.We used a Bayesian implementation of the models considered for this manuscript.We chose the specification of the MGH-t prior for the variance components (Huang & Wand, 2013) over the conventional Wishart prior since the latter may yield too lenient confidence intervals than the former.Overly lenient confidence intervals for the variance components can be much more problematic for inferences about fixed effects than extremely conservative confidence intervals for the variance components (Burger et al., 2020).
The data contamination simulation study suggests that the Johnson-t model for modeling the quantiles of continuous proportions is remarkably robust to outliers, whereas the Kumaraswamy model is susceptible to outliers, especially when modeling the extreme quantiles.
The HPD intervals of the variance components are generally very conservative; however, the coverage for the fixed effects is satisfactory.The degrees of freedom estimates are biased when there are no outliers in the data.However, accurate estimation of the t-distribution's degrees of freedom is known to be challenging: see Fonseca, Ferreira, and Migon (2008) for a detailed explanation; furthermore, in the current context, precise estimation of the degrees of freedom is not required.Instead, we need to determine whether the degrees of freedom are (i) small, which leads to a model robust to outliers, or (ii) large, suggesting that outliers are not an issue.Overall, the simulation study suggests that the proposed robust model has good accuracy and confidence interval coverage properties.
The Jeffreys prior of Juárez and Steel (2010) for the degrees of freedom can be used as an alternative to the hierarchical prior as a sensitivity analysis.However, the trigamma function is not available in JAGS; the Stan software (Carpenter et al., 2017) can alternatively be used to specify the Jeffreys prior as it contains the trigamma function.
It should be noted that, in some cases, excluding outliers because they deviate considerably from the other observations may misrepresent vital ecological processes, ultimately leading to misleading conclusions.From a statistical perspective, it seems preferable to carry out an analysis that is robust to outliers rather than an analysis that is preceded by outlier removal.
The robust regression models introduced can be extended to model the precision and dispersion parameters as a function of covariates.
Should inferences about the quantiles of the continuous proportions on the population level be of interest (i.e., as opposed to conditional on the random effects), additional computationally expensive steps to integrate (marginalize) over the distribution of the random effects are needed for our proposed model.
In conclusion, our study demonstrated that the proposed Johnson-t model is an appropriate robust alternative to the current approach, the Kumaraswamy model, for modeling the quantiles of correlated continuous proportions when outliers are present in the data.
Cushion plant dataset: final versus initial Azorella selago dead stem cover by plot, altitude, and aspect.a Plots are situated either at a mid or high altitude on the western or eastern side of the island.b The red triangles represent the data originally excluded from Raath-Krüger et al. (2022)'s analysis; in contrast, the blue dots represent the data originally included in the study ofRaath-Krüger et al. (2022) Cushion plant dataset: quantile regression coefficients' PEs and 95% HPD intervals computed from the Kumaraswamy and Johnson-t models by parameter, model, and quantile.HPD, highest posterior density.PE, posterior estimate.a Predictor:  1 = Initial dead stem cover (%),  2 = Initial Azorella selago size (cm 2 ),  3 = Initial Agrostis cover (%),  4 = Initial cover of other (%),  5 = Mid versus high,  6 = West versus east.b The red dots and solid lines/shaded bands represent the quantile regression coefficients' PEs and 95% HPD intervals under the Kumaraswamy model, whereas the blue triangles and dash lines/shaded bands represent those under the Johnson-t model; the black dotted line denotes the reference line at zero.c The logit link function was used to model the final A. selago dead stem cover's quantiles as a function of covariates.
Cushion plant dataset: descriptive statistics of final Azorella selago dead stem cover by variable.A. selago individuals were surveyed at sampling sites situated on the eastern and western aspects of Marion Island at mid and high altitudes.
TA B L E 1 a Table 2 summarizes the mean Cushion plant dataset: mean and quantile regression models' odds ratio estimates computed from the rectangular beta and Johnson-t regression model.
F I G U R E 3 Cushion plant dataset: Kullback-Leibler divergence estimates computed from the Kumaraswamy and Johnson-t regression models by quantile, observation, and model.aThe blue lines represent data points that are not considered outliers, whereas the red lines represent influential observations that significantly affect the quantile regression coefficients' estimates (i.e., outliers).TA B L E 2 Simulation study: performance of the Johnson-t regression model.a The logit link function was used to model the response variable's quantiles as a function of covariates.
a b Coverage of 95% highest posterior density intervals (%).

a Johnson-t a Quantile Rate b Parameter Value Bias RMSE Coverage c Length d Value Bias RMSE Coverage c
Simulation study: robustness of the Kumaraswamy and Johnson-t regression models.
TA B L E 4

4 Continued Kumaraswamy a Johnson-t a Quantile Rate b Parameter Value Bias RMSE Coverage c Length d Value Bias RMSE Coverage c
Abbreviations: HPD, Highest posterior density; NA, not applicable.a The logit link function was used to model the response variable's quantiles as a function of covariates.b Contamination rate.c 95% HPD interval coverage (%).d 95% HPD interval average length.