{'fraction_non_alphanumeric': 0.03772402556711368, 'fraction_numerical': 0.01558257091531938, 'mean_word_length': 4.506786289394985, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Clustering Airbnb Reviews

A clustering approach is developed for Boston Airbnb reviews, in the English language, collected since 2008. This approach is based on a mixture of latent variables model, which provides an appealing framework for handling clustered binary data. In the broader context of social science applications (e.g., voting data, web reviews, and survey data), extremely large numbers of variables rule out the use of a mixture of latent trait models. A penalized mixture of latent traits approach is developed to reduce the number of parameters and identify variables that are not informative for clustering. The introduction of component-specific rate parameters avoids the over-penalization that can occur when inferring a shared rate parameter on clustered data. A variational expectation-maximization algorithm is developed and provides closed-form estimates for model parameters; this is in contrast to an intensive search over the rate parameters via a model selection criterion. This approach is important for a whole class of applications, but the focus herein is the Boston Airbnb reviews data.


Introduction
Since 2008, 65,275 guests have provided detailed English comments on the Boston Airbnb website. The goal is to separate these guests into meaningful clusters to better understand the market and, from the point of view of the proprietor, the main driving forces behind positive and negative reviews. Ignoring rarely used terms, there are a total of 278 words used in the reviews (details in Section 4). The frequency with which these 278 words are used can be visualized via a word cloud (Figure 1), and the objective of the analysis herein is to cluster similar reviews into meaningful groups. Clustering high-dimensional (sparse) binary data is a challenging problem and, because of the lack of a suitable method, an approach is introduced herein based on a mixture of latent variables model.
Broadly speaking, cluster analysis is the organization of a data set into meaningful groups and mixture model-based clustering is recently receiving broad interest. According to this approach, data are clustered using some assumed mixture modelling structure, and a mixture model is a convex linear combination of a finite number of component distributions.
Extensive details on mixture model-based clustering are given by Fraley and Raftery (2002) and McNicholas (2016a,b). Popular clustering methods for quantitative data are based on the Gaussian mixture model, which assume that each mixture component is represented by a multivariate Gaussian probability (e.g., Wolfe, 1965;Banfield and Raftery, 1993;Celeux and Govaert, 1995). Recent work on the analysis of clustered binary data based on mixtures of latent trait models include Muthen et al. (2006), Vermunt (2007), Browne and McNicholas (2012), and Gollini and Murphy (2014). Browne and McNicholas (2012) introduce a mixture of latent variables models for the model-based clustering of data with mixed type, and a data set comprising only binary variables fits within their modelling framework as a special case. Within their parameter estimation procedure, Browne and McNicholas (2012) draw on the deterministic annealing approach of Zhou and Lange (2010). This approach focuses on increasing the chance of finding the global maximum; however, Gauss-Hermite quadrature is required to approximate the likelihood. A mixture of item response models (Muthen et al., 2006;Vermunt, 2007) uses a probit structure, and numerical integration is required to compute the likelihood. Thus, it can be difficult to apply to large heterogeneous data sets in practice. A similar approach has also been discussed by Cagnone and Viroli (2012), who use Gauss-Hermite quadrature to approximate the likelihood. In addition, they assume a semiparametric distributional form for the latent variables by adding extra parameters to the model. Gollini and Murphy (2014) propose a mixture of latent trait analyzers (MLTA) for model-based clustering of binary data, wherein a categorical latent variable identifies clusters of observations and a latent trait is used to accommodate within-cluster dependency. They consider a lower bound approximation to the log-likelihood. This approach is easy to implement and converges quickly in comparison with other numerical approximations to the likelihood.
A problem that arises with high-dimensional binary data is the large number of parameters; consequently, there has recently been an increasing interest in penalized latent variable models for binary data (see Houseman et al., 2007;DeSantis et al., 2008, for examples). Houseman et al. (2007) propose a penalized item response theory model with univariate traits and penalized the item-response slopes with ridge penalties. However, it does not take into account the potential group structure of the data and Gauss-Hermite quadrature is required to approximate the likelihood. DeSantis et al. (2008) develop a penalized latent class model to facilitate analysis of high-dimensional ordinal data. A ridge penalty is introduced to the feature-based parameterization of class-specific response probabilities to stabilize maximum likelihood estimation. Both methods adopt a shared rate parameter among variables and require a model selection criterion, such as the Bayesian information criterion (BIC; Schwarz, 1978) to choose it.
For these reasons, a penalized mixture of latent trait models (PMLTM) is proposed for clustered binary data. The data are assumed to have been generated by a mixture of latent trait models (Gollini and Murphy, 2014) and we shrink the slope parameters using a gamma-Laplace penalty function. The PMLTM model enables us to encourage sparsity in estimating the slope parameters, thus reduces the number of free parameters considerably and achieves automatic variable selection. Moreover, the component-specific independent rate parameter avoids the over-penalization that can occur when inferring a shared rate parameter on clustered data. The newly developed variational expectation-maximization (VEM) algorithm (Tipping, 1999;Gollini and Murphy, 2014) provides closed form estimates for model parameters and avoids intensive searches of the rate parameters through model selection criterion, e.g., the Bayesian information criterion (BIC; Schwarz, 1978).
2 Model-Based Clustering via Penalized Mixture of Latent Trait Models

Overview
In this paper, we assume that each observation x i , i = 1, . . . , n, comes from one of the G components and we use z i = (z i1 , . . . , z iG ) to identify the component membership, where z ig = 1 if observation i is in component g and z ig = 0 otherwise. The conditional distribution of X i in component g is a latent trait model and takes the form, where and the response function for each categorical variable in each component is where α mg and w mg are the model parameters and the multivariate latent variable Y i ∼ MVN(0, I D ). Under these conditions, the approach represents a generalization of the latent class model where observations are not necessarily conditionally independent given the group memberships. In fact, the observations within groups are modelled using a latent trait analysis model and thus dependence is accommodated. This model is known as mixture of latent trait analyzers (MLTA; see Gollini and Murphy, 2014).

Penalized MLTA Models via Non-Convex Penalties
A potential drawback of the MLTA for high-dimensional data is its large number of parameters. In particular, the model in (1)  To reduce the number of free parameters, we propose a penalized log-likelihood of the form where l(Θ) is the log-likelihood of (1) and C(Θ) is a penalty term. Similar to the LASSO penalty for regression (Tibshirani, 1996), we propose use of heavy-tailed and sparsity-inducing independent Laplace prior for each coefficient w mg . To account for uncertainty about the appropriate level of variable-specific regularization, each Laplace rate parameter λ jg is left unknown with a gamma hyperprior. Thus, for s, r, λ mg > 0. However, available cross-validation (e.g., via solution paths) and fully Bayesian (i.e., through Monte-Carlo marginalization) methods for estimating w mg under unknown λ mg are prohibitively expensive. A novel algorithm is proposed for finding posterior mode estimates of the slope parameters, i.e., maximum a posteriori (MAP) estimation, while treating λ mg as missing data via an EM algorithm. The MAP inference with fixed λ mg is equivalent to likelihood maximization under an L 1 -penalty in the LASSO estimation and λ mg ∼ Gamma(s, r) leads us to a a non-convex penalty ( Figure 2)

Motivation for Gamma-Laplace Penalties
One unique aspect of our approach is the use of independent gamma-Laplace priors for each slope parameter w mg . The Laplace prior for w mg encourages sparsity in w mg through a sharp density spike at w mg = 0 and MAP inference with fixed λ mg is equivalent to likelihood maximization under an L 1 penalty in the LASSO estimation and selection procedure of Tibshirani (1996). In the Bayesian inference for LASSO regression, conjugate gamma hyerpriors are a common choice for the rate parameter λ (e.g., Park and Casella, 2008;Yuan and Wei, 2014). However, that independent rate parameter λ mg is thought to provide a better representation of prior utility, and it avoids the over-penalization that can occur when inferring a shared rate parameter on clustered data. As detailed in Section 2.2, our approach yields an estimation procedure that corresponds to likelihood maximization under a specific non-convex penalty that can be seen as a reparametrization of the "log-penalty" described in Mazumder et al. (2012). Similar to the standard LASSO, singularity at zero in C(w mg ) causes some coefficients to be set to zero.
However, unlike the LASSO, the gamma-Laplace has gradient , which disappears as D d=1 |w dmg | → ∞, leading to the property of unbiasedness for large coefficients (Fan and Li, 2001). Commonly, the rate parameter λ is selected using crossvalidation or an information criterion such as the BIC. However, our independent λ mg would require searches of impossibly massive dimension. Moreover, cross-validation is an estimation technique that is sensitive to the data sample on which it is applied. That said, one may wish to use cross-validation to choose s or r in the hyperprior and, because results are less sensitive to these parameters than they are to a fixed penalty, a small grid of search locations should suffice.

Interpretation of the Model Parameters
The model parameters can be interpreted exactly as for MLTA and item response models.
In the finite mixture model, η g is the proportion of observations in the gth component. The characteristics of component g are determined by the parameters α mg and w mg . In particular, the intercept α mg has a direct effect on the probability of a positive response to the variable m given by an individual in group g, through the relationship The value π mg (0) is the probability that the median individual in group g has a positive response for the variable m. However, when the data set has very low percentage of positive responses (e.g., text data), the value of π mg (0) can be very low for all items across all components. Thus we use the slope parameters to characterize each component in Section 4. The slope parameters w mg are known as discrimination parameters in item response theory. The larger the value of w dmg , the greater the effect of factor y d on the probability of a positive response to item m in group g. The quantity w dmg can be used to calculate the correlation coefficient between the observed item x i and the multivariate latent variable Y i . In the latent trait case, the slope parameters cannot be interpreted as correlation coefficients, because they are not bounded by 0 and 1. However, it is possible to transform the loadings so that they can be interpreted as correlation coefficients in exactly the same way as in factor analysis. The standardized w dmg is given by The purpose of the Laplace prior for w mg is to encourage sparsity in w mg , therefore identifying non-informative variables for each component. When the mth row of the slope parameter matrix in gth component is zero everywhere (w 1mg = w 2mg = · · · = w Dmg = 0), then the corresponding variable is not informative. In addition, w * dmg = 0 (or, equivalently, w dmg = 0) indicates that item m is independent from latent trait y d in component g.

Model Identifiability
The identifiability of our model depends on the identifiability of the latent trait part as well as the identifiability of the mixture model. The identifiability of finite mixture models has been discussed by several authors (e.g., McLachlan and Peel, 2000). Knott and Bartholomew (1999) discuss the model identifiability issue in the latent trait analysis. One necessary condition for model identifiability is that the number of the free parameters to be estimated not exceed the number of possible data patterns. However, this condition is not sufficient because the actual information in a data set can be less depending on the size of the data set. As with the mixture factor analysis model, the solution is not unique when d > 1, i.e., an orthogonal rotation of the latent variable coupled with corresponding rotation of the estimated slope parameters w mg leaves the log posterior unchanged. However, the model rotates the slope parameters automatically by shrinking the slope parameter through a penalized likelihood method. Jaakkola and Jordan (2000) introduced a variational approximation for the predictive likelihood in a Bayesian logistic regression model and also briefly considered the "dual" problem, which is closely related to the latent trait model. It obtains a closed form approximation to the posterior distribution of the parameters within a Bayesian framework. Their method is based on a lower bound variational approximation of the logistic function

Variational Approximation
It can be approximated by the exponential of a quadratic form involving variational parameters ξ ig = (ξ i1g , ..., ξ iM g ), where ξ img = 0 for all m = 1, ..., M . Now, the lower bound of each term in the log-likelihood is given by wherẽ This approximation has the property thatp(x im |y i , z ig = 1, ξ img ) ≤ p(x im |y i , z ig = 1, ξ img ) with equality when |ξ img | = A img . At first glance, the derivation of slope parameters w mg is a challenging task, due to the fact that the penalization term is not differentiable at w mg . To this end, we write and exploit the concavity of this square-root. In particular, with equality if and only if w mg = w mg . Using these inequalities, we can obtain a surrogate function which can be used to obtain parameter estimates (Section 3.2).

Parameter Estimation
The VEM algorithm is a natural approach for MAP estimation when data are incomplete, and it is used for the PMLTM. log p(x i |θ g , y i , z ig ) + log p(θ g |z ig , λ g ).
(6) VE-step. We compute the expected value of the complete-data log posterior in the VE-step using the expected values of the missing data in log p(Θ|x, y, z, λ). We require the expectations .

Compute the location vector and the covariance matrix forp(y
The expected value of the independent multidimensional rate parameter can be written where s, r > 0 are predetermined shape and rate parameters of the gamma hyperprior. M-step 1. The first M-step on the (t + 1) iteration is to optimize the variational parameter ξ img in order to make the approximationp(x i |z ) as close as possible to p(x i |z ig = 1): M-step 2. Optimize the parameters w mg and α g in order to increase the log (complete-date) posterior log p(w mg , α g |x, y (t+1) , λ We adopt a numerically more convenient form of the update for w mg : Dmg | 1/2 . This avoids estimating |w −1(t) dmg |, some of which are expected to go to zero.
The M-step 3. Estimate η g : Compute the log posterior. Obtain the lower bound of the log (complete-data) posterior at the expansion point ξ ig : Convergence criterion. The convergence of our VEM algorithm is determined using a criterion based on the Aitken acceleration (Aitken, 1926), i.e., |l and l (t) is the log posterior at iteration t. See Böhning et al. (1994) for details.

Model Selection
We use the Bayesian information criterion (BIC; Schwarz, 1978) as a criterion for model selection, i.e., BIC = −2l + k log n, where l is the maximized log-likelihood, k is the number of free parameters to be estimated in the model, and n is the number of observations. The presence of a penalty term reduces the number of free parameters of the slope parameter matrix as the effective degrees of freedom of W equals the number of nonzero terms in the loading parameter matrix. Other methods for computing the effective degrees of freedom of W are also available (e.g., Hunter and Li, 2005). Within the framework of MLTA models, the number of components G and the dimension of the latent variable Y (i.e., d) need to be determined. When defined as in (7), models with lower values of BIC are preferable. The BIC value could be overestimated using the variational approximation of log-likelihood, which is always less than or equal to the true value. For model selection purposes, we calculate maximum the log-posterior using Gauss-Hermite quadrature after convergence is attained.
For high-dimensional binary data, particularly when the number of observations n is not very large relative to their dimension m, it is common to have a large number of patterns with small observed frequency. Accordingly, we cannot use a χ 2 test to check the goodness of the model fit. In the simulated examples in Section 3.5, where the true classes are known, the adjusted Rand index (ARI; Hubert and Arabie, 1985) can be used to assess model performance. The ARI is the corrected-for-chance version of the Rand index (Rand, 1971). The general form of the ARI is index − expected index maximum index − expected index , which is bounded above by 1, and has expected value 0 under random classification.
In real examples, such as the analysis of the Boston Airbnb reviews in Section 4, the analysis of the clusters in the selected model can be used to interpret the model.

Selection of Programming Languages
When fitting the PMLTM model using R the task becomes increasing burdensome as the number of items becomes large. Therefore, we implement our algorithm in two scripting languages, R and Python, and compare their performance. Python is an elegant open-source language that has become popular in the scientific community. We use the Numpy library for matrix operations and Scipy.stats library for the use of probability distributions and statistical functions.

Simulation Study
A simulation study is performed to illustrate the proposed PMLTM model. A set of 100 samples of n = 500 observations is generated from a PMLTG model with a two-component mixture (G = 2, π 1 = π 2 = 0.5). The latent variable is generated from a Gaussian distribution, i.e., Y ∼ N (0, 1). Table 1 reports the slope parameters w used to generate a set of M = 10 observed variables. For each sample, the value of the gamma hyperparameters, i.e., (s, r), are selected from {(0.1, 0.5), (0.5, 0.5), (1, 0.5), (2, 0.5)}. Table 2 shows the BIC and ARI values averaged on the 100 samples for each pair (s, r). As shown in Table 2, on average, the BIC has a minimum and ARI has a maximum when (s, r) = (1, 0.5). Therefore, we use (s, r) = (1, 0.5) for the analysis of the Boston Airbnb Reviews. A comparison of computing times between R and Python, based on these simulated data, is given in Appendix A. Table 1: Component-specific slope parameters.

Boston Airbnb Reviews
This data set includes detailed English comments on the Airbnb website in the Boston area from 65,275 guests (i.e., n = 65275) since 2008. We perform some pre-processing of the text data (i.e., converting the text to lower case, removing numbers and punctuation, removing stop words, and stemming). These basic transforms are available within the R package tm (Feinerer and Hornik, 2015). We then create a matrix with each comment as a row and each word as a column. If a word is mentioned in a comment, the response for the corresponding cell is coded as 1, and otherwise is 0. The term matrix contains 43,584 words but most of them are infrequently used, i.e., so-called "sparse terms". Sparse terms that appear in less than 2% of all reviews are not of interest and so are removed. At the end of this preprocessing step, the term matrix consists of 278 words (i.e., M = 278) and the word cloud in Figure 1 provides a quick visual overview of the frequency of the words in the final term matrix.
The PMLTM model is fitted to these data for D = 1, . . . , 5 and G = 1, . . . , 5. We run all models in both R and Python. Using R, it takes more than 24 hours to obtain results while Python takes only 106 minutes. The minimum BIC occurs at the three-component, two-dimensional (i.e., two latent traits) PMLTM model. The BIC value is 869, 496. The clusters (components) for the selected model (G = 3, d = 2) are summarized in Table 3. Average sentiment scores for each cluster are calculated using the built-in Python library nlkt (natural language toolkit). The sentiment of each comment -positive, negative, or neutral -is presented using a score ranging from 0 to 1. Because each comment could contain positives and negatives at the same time, a compound score is presented as well. Each compound score in on [−1, 1], where −1 corresponds to an overall unpleasant tone and 1 is an overall pleasant tone. Cluster 1 consists of comments with a positive tone overall, as indicated by the compound score of 0.95. The average positivity score in Cluster 1 is 0.30. Cluster 3 consists of mainly positive comments as well, but the average compound score is lower in Cluster 3 when compared to Cluster 1. Cluster 2 is a small group that consists of comments that have a slightly negative tone overall. It is worth noting that the average negativity score is higher than the average positivity score in Cluster 2. Table 4 shows the high-loading words for each latent trait in each cluster. We note that the first latent trait Y 1 is concerned with the property (e.g., location, condition, etc.) whereas the second latent trait Y 2 is concerned with the host. The comments in Cluster 1 and Cluster 3 both have an overall positive tone. However, comments in Cluster 1 are more specific about the listings and the hosts. Moreover, the

Summary
An MLTA approach that encourages sparsity in estimating the slope parameters -thus reducing the number of free parameters considerably -is introduced for clustering the Boston Airbnb data. The component-specific rate parameters avoid the over-penalization that can occur when inferring a shared rate parameter on clustered data. The PMLTM model retains the ability to investigate the dependence between variables while clustering with the added advantage of being able to model very high-dimensional binary data (e.g., text data). Applying the PMLTM model to the Boston Airbnb reviews data showed that the method scales to far larger datasets than any existing model-based clustering methods for binary data. The results for these data reveal two latent traits and three clusters of reviews. The latent traits can be interpreted as concerning the property and the host, respectively. One cluster contains highly positive reviews, another contains positive reviews, and the other contains reviews that are not positive, i.e., moderate and negative reviews. Table 6 shows a comparison, based on the simulated data from Section 3.5, of the average run time over 100 loops of the VE-step and M-steps using R and Python (G = 2, D = 2, n = 100). Python runs approximately 103 times faster than R for the VE-step and 190 times faster for the M-steps.