{'fraction_non_alphanumeric': 0.03568938469583406, 'fraction_numerical': 0.016384870141188775, 'mean_word_length': 4.47708830548926, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Analysis of Augmented Mixed Beta Models with Skew-Normal Random Effects

Many studies in different areas include data in the form of rates or proportions that should be analyzed. The data may also accept values zero and one. Augmented beta regression models are an appropriate choice for continuous response variables in the closed unit interval [0, 1]. The data in this model are based on a combination of three distributions, degenerate distribution at 0 and 1, and a beta density in (0, 1). The random effects are usually added to the model for accommodating the data structures as well as correlation impacts. In most of these models, the random effects are generally assumed to be normally distributed, while this assumption is frequently violated in applied studies. In this paper, the augmented mixed beta regression model with skew-normal distributed random effects is presented. A Bayesian approach is adopted for parameter estimation using Markov Chain Monte Carlo method. The proposed model is applied to analyze a real data set from Labor Force Survey.


Introduction
Beta regression models are commonly used to model variables that assume values in the open interval (0, 1). They are based on the assumption that the dependent variable follows a beta distribution, also its mean is related to a set of explanatory variables through a linear predictor with unknown coefficients and a suitable link function, (Figueroa-Z et al., 2013). The beta distribution is flexible for modeling proportions and rates since its density has different shapes depending on the values of the two parameters. Many studies have been done for observed outcomes containing zero and one which are not in the support of the beta distribution. Inflated beta distributions were explained as a mixed continuous-discrete distribution to model proportion values that contain zeros and ones by Ospina and Ferrari (2010). Inflation happens when the probability mass of some points exceeds that of a proposed distribution. For example, zero-inflated counts or zero inflated binomial regression models that were introduced by Lachenbruch (2002) and Hall and Berenhaut (2002), respectively. In beta distribution which is a continues one, the model with extra zeros and ones is preferred to be called as augmented points instead of inflated points at zero and one (Galvis et al., 2014). They were explained augmented beta regression by mixed continuous-discrete distribution to model proportion values containing zeros and ones. Verkuilen and Smithson (2012) extended generalized linear mixed models for beta distributed response variables. Mixed beta regression model under the Bayesian approach while including fixed and random effects are discussed by (Figueroa-Z et al., 2013). Random effects models are commonly used to model correlated data. In longitudinal data, these effects are usually assumed to be normally distributed. Some studies showed that the normality assumption is inappropriate when data exhibit non-normal behavior (Zhang and Davidian, 2001). Thus, some extension of statistical models with remarkable flexibility in the distributional assumptions are required that can readily adapt to the non normality behavior of the random effects. Azzalini (1985) presented a flexible family of skew-normal distributions that contains the normal family as a special case. Studies on the multivariate skew-normal distributions are also presented by (Azzalini and Dalla Valle, 1996). The skew-normal linear mixed models are presented by Arellano-Valle et al. (2005), Lachos et al. (2010) and Lin (2008).
The augmented mixed beta models was defined under normality assumption for random effects which is an extension of linear mixed effect model.
A commonly used approach for analyzing proportion cluster data such as longitudinal proportion is augmented mixed beta models. In this paper, we have propose using a skew-normal random effects instead of the normal one to allow more flexibility in the augmented mixed beta regression models. So we consider a zero and one augmented Beta regression model with a skewnormal distributional assumption for random effects. This model is amenable to the Bayesian implementation and the MCMC scheme described for sampling from the posterior distribution. We have used the available software WinBUGS and R package for fitting the models, where the criteria of Log Pseudo Marginal Likelihood (LPML) have been used for model comparison.
The motivating data example in this paper comes from a Labor Force Survey (LFS) data to determine the proportion of employed persons in the household. For each household, this variable is measured four times and since these data are considered as longitudinal data. We need to investigate random effects in the model for considering the correlation among the response variable during the time. In the application Section, we intend to compare the performance of the augmented mixed beta regression model under normal and skew-normal random effects and we conclude that the augmented mixed beta regression under skew-normal random effects is a better fit compared to other one.
The rest of the paper proceeds as follows. After a brief introduction to the beta regression model and the augmented mixed beta regression model under normal distribution for random effects, Section 2 introduces the augmented mixed Beta models with skew-normal random effects. Section 3 develops the Bayesian estimation scheme. In Section 4 some Bayesian criteria have been presented for model selection. In Section 5, some simulation studies are carried out to examine the performance of the proposed model. Section 6 contains an application using an Iranian LFS data set and concluding remarks are given in Section 7.

Statistical Model and Bayesian Inference
In this section, the beta regression model is introduced for responses lying in (0, 1). Then an augmented mixed beta regression model with normal random effects is presented for modeling the observations in the closed interval [0, 1]. Next, we utilize skew-normal distributions for random effects in the augmented mixed beta regression model as a tool for robust modeling under a Bayesian approach.

Beta Regression Model
The beta distribution is often applied for modeling continuous data that are confined in the interval (0, 1). A random variable Y follows a beta distribution parameterized in terms of its mean and precision parameter ϕ is given by where 0 < y < 1, 0 < µ < 1 and . A beta regression model can be considered under a generalized linear model (GLM) framework by linking the subject-specific mean µ i and covariate vector X i (Cribari-Neto and Zeilis, 2010). This is given as g(µ i ) = X T i β, where β is the vector of regression parameters.

Augmented Mixed Beta Regression Models
In some studies the data may lie in [0, 1]. For modeling such data the augmented beta regression is defined by (Galvis et al., 2014) that comprises of a mixture of three distributions, with two degenerated distributions at 0 and 1, and a beta distribution. Thus the density of Y i , i = 1, . . . , n, is given by: where p 0 = P (Y i = 0), p 1 = P (Y i = 1) and f (y i ; µ i , ϕ) is given in (1). We denote the density of Y i by AB(p 0 , p 1 , µϕ, (1 − µ)ϕ) which using indicator functions would be written as follows: (3) where I(·) denotes the indicator function. The mean of Y i and its variance can be derived as: The augmented mixed beta regression model is defined as follows. Let Y 1 , . . . , Y n be independent continuous random vectors, where Y i = (y i1 , . . . , y in i ) represents an observed response vector for a sample unit i, with components, y ij in the interval [0, 1]. The regression model is obtained over a suitably transformed µ i , with the following structure by (Galvis et al., 2014).
where X i is the design matrix of dimension p×n i corresponding to the vector of fixed effects β = (β 1 , . . . , β p ) T , and Z i is the design matrix of dimension q × n i corresponding to the vector of random effects

Augmented Mixed Beta Models with Skew-Normal Random Effects
The normality of the random effects is a routine assumption for the GLM framework, but it may be unrealistic. We have utilized skew-normal distribution as a tool for robust modeling under a Bayesian approach. A random vector b follows a multivariate skew-normal distribution with location vector µ, covariance matrix Σ b and skewness parameter vector λ, if its pdf is given by: where ϕ p (·; µ, Σ b ) is the density function of the q-variate normal distribution and Φ(·) is the cumulative distribution function of the standard univariate normal distribution (Azzalini and Dalla Valle, 1996). Now consider the augmented mixed beta model with a single random intercept, where the random effects has pdf as in (5) . A particularly useful and easy interpretable link function for proportional data is the logit link function (Ferrari and Cribari-Neto, 2004), therefore the logit link function is used.
For the logit link function we can consider the following structure:

Model fitting Using MCMC Sampling
We use the Bayesian methodology using MCMC sampling for the Augmented mixed beta regression model with skew-normal distributional assumptions for random effects since a full Bayesian approach is difficult to implement for such models. One MCMC type approach, that needs only the characteristics of the conditional posterior distribution of each parameter, is the Gibbs sampler.
For convenience, we suppose that elements of θ in the model are independent and are defined by θ = (β, p 0 , p 1 , ϕ, σ 2 b , λ). The primary goal here is to estimate θ. Considering the random effects design matrix Z i to be an identity matrix, (y 1 , x 1 ), . . . , (y n , x n ) are the observed sample for n subjects, with y i as the response vector for subject i. Considering the fixed effects as β, a Bayesian specification of the model requires consideration of the prior distribution for all parameters. We attempt to assign weakly informative distribution for the parameters. Also, the hyper parameters are selected such that they lead to the low informative prior distributions. We adopt the following prior specifications: where HN denotes the half-Normal distribution, and IΓ denotes the inverse gamma distribution. Let b = (b 1 , b 2 , . . . , b n ), the complete conditional likelihood of the outcome y is given by: where θ = (β, p 0 , p 1 , ϕ, σ 2 b , λ) denotes the unknown model parameters vector. For convenience, we use the stochastic representation of the skew-normal and def = meaning 'distributed as', with independent Z 0 and Z 1 by Arellano-Valle et al. (2005).
The joint posterior distribution using the Bayes Theorem is given by . The posterior distribution (9) is analytically intractable, but the MCMC methods such as the Gibbs sampler and Metropolis-Hastings algorithm can be used to draw samples. The Gibbs sampler works by drawing samples iteratively through full conditional distributions derived from (9). Since some of the full conditional distributions are not available, the Metropolis-Hastings algorithm is applied. We implement the Bayesian methodology using MCMC techniques for the beta Regression model with the skew-normal distributional assumption for random effects. A key feature of this model, which allows writing WinBugs codes, is that it can be formulated in a flexible hierarchical representation as follows: The relevant MCMC steps were implemented using the R2WinBUGS package which connects the R with the WinBUGS software.

Bayesian Model Selection
Many criteria have been defined for model selection in the Bayesian inference. One of the popular criterion, which is usually used is the conditional predictive ordinate (CPO) statistic. Let the vector Y be the data set, Y (−i) , i = 1, · · · , n denote all data except the ith observation and π(Θ|Y (−i) ) denotes the posterior distribution of the model parameter vector Θ given where m is the number of iterations after a burn-in period (Gelfand and Dey, 1994) The larger values of LPML indicates the better fit for the model. Two other criteria EAIC and EBIC (Carlin and Louis, 2008) can also be used for model selection. Let θ (k) be the MCMC posterior sample generated at the iteration k of the algorithm, η is the number of parameters, n is the number of observations and D = 1 Then EAIC = −2D + 2η, and EBIC = −2D + 2η log n. Unlike LPML, the smaller values of EAIC and EBIC indicate a better fit for the model.
In order to assess the convergence of MCMC algorithm, we have used the Gelman and Rubin's diagnostic test, (Gelman and Rubin, 1992). This method requires two or more parallel chains (the number of which is denoted by m ′ ) to be generated, each with different starting values. For assessing convergence of individual model parameters, θ, this convergence diagnostic test is proposed as a univariate statistic, referred to as the potential scale reduction factor (PSRF). Calculation of this ststistic is based on the last n ′ samples in each of m ′ parallel chains. This factor is calculated by: where (B/n ′ ) is the between-chain variance and W is the within-chain variance and are calculated by: As chains converge to a target distribution, the estimated varianceV = (1 − 1 n ′ )W + 1 n ′ B should be almost equivalent to the within-chain variance, so PSRF should approximately equal one.

Simulation Study
Some simulation studies are illustrated to assess the performance of the proposed models. In all cases, we generate 100 data sets with sample size n = 100 with n i = 5, i = 1, . . . , 100. To generate the dependent variable with p 0 = p 1 = 0.1, first we draw a Bernoulli sequence with success probability (1 − p 0 − p 1 ) that is allocated to y ij ∈ (0, 1). Then, conditional on this draw, we carry out another Bernoulli with equal success probability for one-occurences. In this way, the sequence of y ij ∈ (0, 1), y ij = 0 and y ij = 1 with p 0 = 0.1, p 1 = 0.1 and (1 − p 0 − p 1 ) = 0.8 are generated. We took the augmented mixed beta regression model with location parameter µ ij . For this purpose, µ ij is generated as: logit(µ ij ) = β 1 +β 2 x 2ij +β 3 x 3ij +b i where i = 1, . . . , n, j = 1, . . . , 5. The values of the covariates are generated from a standard normal distribution and regression parameters are fixed at: β 1 = 0.5, β 2 = 0.4, and β 3 = 0.6. The parameter ϕ is given 40. The random effects is generated by using Formula 8, when Z 0 and Z 1 are independent random variables that have standard normal distributions. Also we consider σ 2 b = 1. Three data sets are generated by assuming λ ∈ {0, 2, 5}. The Bayesian analysis is specified under priors (7). For each model, two parallel chains in different initial points were run, with 800,000 iterations for each chain. The results are presented considering the last 500,000 iterations. In addition, to avoid correlation problems in the generated chains, the thin value is considered 100.
In Table 1, random effects are assumed to follow normal and skew-normal distributions. Three data sets were generated by assuming λ ∈ {0, 2, 5}, representing behaviors from weak to strong skewness. To compare normal and skew-normal distributions for random effects, we calculate two criteria, the average of bias and the standard deviation (SD) that were computed as follow: whereθ i is the estimation of θ from the ith simulated data set.
In Table 1, the random effects are generated from skew-normal distribution for three cases, λ =0,2,5. Which λ =0 is the normal one. Then, the normal and the skew-normal distributions are allocated to random effects. For λ =2,5 by comparing biase's and SD's, the skew-normal for random effects is better than of the normal distribution especially when λ=5. When λ =0 by comparing the biase's and the SD's of parameter estimates, the results show that although the normal one has a little better performance but the skew-normal distribution is also admitted for random effects. But when λ > 0, allocation of a normal distribution to random effects is not desirable and it declines by increasing λ. Generally, our evaluation suggests that the skewnormal distribution for random effects is better than normal distribution. In other words, allocating skew-normal distribution to the random effects provide better estimates for the model parameters.
In addition, since the Bayesian analysis is based on prior distributions, λ ∼ HN (0, 10 k ), we performed the Bayesian inference with different values of k. The results corresponding to the posterior analysis for λ=5 and k = 2, . . . , 5 are reported in Table 2. Changing values of k for k ⩾ 2 shows that results are not sensitive to the values of k. Also, we allowed σ 2 b ∼  IΓ(a, a) and ϕ ∼ Γ(a, a), where a ∈ {0.01, 0.1}. The results show that the parameter estimates are not sensitive for values of a. In order to see how stable estimates are, two parallel chains with different initial values should be used. The convergence of the MCMC samples of parameters were verified by Gelman and Rubin test. Furthermore, for each parameter, we checked that the convergence is achieved for each chain.

Application
In this section, we apply our proposed model for analysis of LFS data that is a quarterly survey of households for measuring the economically active population conducted by the Statistical Centre of Iran. The LFS consists of probabilistic samples of household units that produces quarterly labor force and related estimates for all members of private settled households whose place of usual residence was located in Iran at the time of the enumeration. The sample rotation follows a 2-2-2 pattern. Therefore, each unit is considered as one household. A household unit is interviewed in two consecutive quarters and not in the sample for the next two quarters. It returns to the sample again for the next two quarters.
The LFS interview is divided into two parts: 1-household and demographic information, 2-labor force information. The second part for each 10 years or older person in the household is obtained. One of the main goals of labor force information is to classify persons as employed, unemployed or inactive. So the proportion of employed persons in the household can be calculated. We consider a panel data on 159 households who had been interviewed in the first and second quarters of 2013 and also the same quarters of 2014 in Tehran, the capital of Iran. Therefore, each household is interviewed four times, consisting 636 observations of 159 households. The definitions and concepts used in LFS were consistent with international recommendations of the International Labor Organisation in 1998. Figure 1 displays the histogram of the proportion of employed persons in the household based on 159 households for each quarter separately that drop in the closed interval [0, 1], where 0 and 1 respectively, represent households with no employed member and households with all members employed. Therefore, we have used the augmented mixed Beta regression model for which random effects have two different normal and skew-normal distributions. Based on the augmented mixed beta regression model, the prior distributions in (7) and logit link function (6), we considered the following structure to this panel data: The density of Y ij , i = 1, . . . , 159, j = 1, 2, 3, 4 is given by: (1−y ij ) (1−µ ij )ϕ−1 . In this model, y ij is the proportion of employed individuals in the household, p 0 is the probability of households with no employed member, p 1 is the probability of households with all members employed, x 1ij is the number of persons that have 10 years or older in the household and x 2ij is the household size. Also µ ij is the mean of the response variable, ϕ is the precision parameter and b i is a random effects which considers the potential correlation between responses over time. Hence significance of its variance shows the existence of correlated responses over time.
It should be noted that multicolinearity for these variables is investigated and the problem was not seen.
Considering prior densities in (7), we generated two parallel independent runs of Gibbs sampler chain with size 800,000 iterations for each parameter, discarding the first 500,000 iterations and avoiding correlation problems, we considered the thin value 100. To monitor the convergence of the Gibbs samples we used the methods presented by (Gelman and Rubin, 1992). Table 3 shows the comparison among the two different models by using the selection criteria considered in Section 4. Notice that, the model with skew-normal distribution for random effects outperforms the one with normal distribution for all criteria. Table 4 gives the parameter estimates (Est.) and standard deviations (SD) obtained from fitting models with random effects under skew-normal and normal distributions to the LF S data. The results show that the number of 10 year old or older persons in the household has a positive effect to explain the proportion of employed persons in the household. Reversely, the household size variable with a negative coefficient is included in the model. The estimate of precision parameter,φ = 43.092. The skewness parameter λ = 0.547, p 0 = 0.309 and p 1 = 0.038.
To monitor Markov chain convergence, we checked the PSRF of all pa- A sensitivity analysis was performed on the prior assumptions and initial values. We allowed σ 2 b ∼ IΓ(a, a) and ϕ ∼ Γ(a, a), where a ∈ {0.01, 0.1}. The results in Table 5 show that the parameter estimations are not very sensitive on the selected values of a. We also checked the sensitivity of the posterior estimates of β on different priors. The results in Figure 2 show that they are not sensitive on t-student priros with degrees of freedom (df) greater than 5.

Conclusion
The augmented Beta regression is utilizable to model data that are seen on [0, 1), (0, 1] or [0, 1]. In this paper, we have considered the augmented Beta regression model with random effects under the skew-normal distributional assumption. Since the model was complicated, the Bayesian method was used through the Markov chain Monte Carlo algorithm. All the Bayesian computation were performed using accessible software WinBUGS. To show the performance of our proposed model, some intensive simulation studies have been done. In the application section, by Labor Force Survey data set, the proportion of employed persons in the household has been analyzed where the augmented mixed Beta with skew-normal random effects have been used. Since random effects on real data is generally not symmetric, the proposed model can be more appropriate model. In this paper we used a random effects model under skew-normal distribution to analyzed the longitudinal proportion data, which is the most famous approach for considering association among repeated measures but one can consider marginal or transition model for analyzing such data as another future work.