% % Copyright (C) 2004-2005 Friedrich Leisch % $Id: flexmix-intro.Rnw 5187 2020-06-25 17:59:39Z gruen $ % \documentclass[nojss]{jss} \title{FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in \proglang{R}} \Plaintitle{FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R} \Shorttitle{FlexMix: Finite Mixture Models in \proglang{R}} \author{Friedrich Leisch\\Universit\"at f\"ur Bodenkultur Wien} \Plainauthor{Friedrich Leisch} \Address{ Friedrich Leisch\\ Institut f\"ur Angewandte Statistik und EDV\\ Universit\"at f\"ur Bodenkultur Wien\\ Peter Jordan Stra\ss{}e 82\\ 1190 Wien, Austria\\ E-mail: \email{Friedrich.Leisch@boku.ac.at} } \usepackage[utf8]{inputenc} \usepackage{listings} \newcommand{\R}{\proglang{R}} <>= suppressWarnings(RNGversion("3.5.0")) set.seed(1504) options(width=70, prompt = "R> ", continue = "+ ", useFancyQuotes = FALSE) grDevices::ps.options(family="Times") library("graphics") library("flexmix") data("NPreg") @ \Abstract{ This article was originally published as \cite{flexmix:Leisch:2004a} in the \emph{Journal of Statistical Software}. FlexMix implements a general framework for fitting discrete mixtures of regression models in the \R{} statistical computing environment: three variants of the EM algorithm can be used for parameter estimation, regressors and responses may be multivariate with arbitrary dimension, data may be grouped, e.g., to account for multiple observations per individual, the usual formula interface of the \proglang{S} language is used for convenient model specification, and a modular concept of driver functions allows to interface many different types of regression models. Existing drivers implement mixtures of standard linear models, generalized linear models and model-based clustering. FlexMix provides the E-step and all data handling, while the M-step can be supplied by the user to easily define new models. } \Keywords{\proglang{R}, finite mixture models, model based clustering, latent class regression} \Plainkeywords{R, finite mixture models, model based clustering, latent class regression} \Volume{11} \Issue{8} \Month{October} \Year{2004} \Submitdate{2004-04-19} \Acceptdate{2004-10-18} %%\usepackage{Sweave} %% already provided by jss.cls %%\VignetteIndexEntry{FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R} %%\VignetteDepends{flexmix} %%\VignetteKeywords{R, finite mixture models, model based clustering, latent class regression} %%\VignettePackage{flexmix} \begin{document} \section{Introduction} \label{sec:introduction} Finite mixture models have been used for more than 100 years, but have seen a real boost in popularity over the last decade due to the tremendous increase in available computing power. The areas of application of mixture models range from biology and medicine to physics, economics and marketing. On the one hand these models can be applied to data where observations originate from various groups and the group affiliations are not known, and on the other hand to provide approximations for multi-modal distributions \citep{flexmix:Everitt+Hand:1981,flexmix:Titterington+Smith+Makov:1985,flexmix:McLachlan+Peel:2000}. In the 1990s finite mixture models have been extended by mixing standard linear regression models as well as generalized linear models \citep{flexmix:Wedel+DeSarbo:1995}. An important area of application of mixture models is market segmentation \citep{flexmix:Wedel+Kamakura:2001}, where finite mixture models replace more traditional cluster analysis and cluster-wise regression techniques as state of the art. Finite mixture models with a fixed number of components are usually estimated with the expectation-maximization (EM) algorithm within a maximum likelihood framework \citep{flexmix:Dempster+Laird+Rubin:1977} and with MCMC sampling \citep{flexmix:Diebolt+Robert:1994} within a Bayesian framework. \newpage The \R{} environment for statistical computing \citep{flexmix:R-Core:2004} features several packages for finite mixture models, including \pkg{mclust} for mixtures of multivariate Gaussian distributions \citep{flexmix:Fraley+Raftery:2002,flexmix:Fraley+Raftery:2002a}, \pkg{fpc} for mixtures of linear regression models \citep{flexmix:Hennig:2000} and \pkg{mmlcr} for mixed-mode latent class regression \citep{flexmix:Buyske:2003}. There are three main reasons why we have chosen to write yet another software package for EM estimation of mixture models: \begin{itemize} \item The existing implementations did not cover all cases we needed for our own research (mainly marketing applications). \item While all \R{} packages mentioned above are open source and hence can be extended by the user by modifying the source code, we wanted an implementation where extensibility is a main design principle to enable rapid prototyping of new mixture models. \item We include a sampling-based variant of the EM-algorithm for models where weighted maximum likelihood estimation is not available. FlexMix has a clean interface between E- and M-step such that variations of both are easy to combine. \end{itemize} This paper is organized as follows: First we introduce the mathematical models for latent class regression in Section~\ref{sec:latent-class-regr} and shortly discuss parameter estimation and identifiability. Section~\ref{sec:using-flexmix} demonstrates how to use FlexMix to fit models with the standard driver for generalized linear models. Finally, Section~\ref{sec:extending-flexmix} shows how to extend FlexMix by writing new drivers using the well-known model-based clustering procedure as an example. \section{Latent class regression} \label{sec:latent-class-regr} Consider finite mixture models with $K$ components of form \begin{equation}\label{eq:1} h(y|x,\psi) = \sum_{k = 1}^K \pi_k f(y|x,\theta_k) \end{equation} \begin{displaymath} \pi_k \geq 0, \quad \sum_{k = 1}^K \pi_k = 1 \end{displaymath} where $y$ is a (possibly multivariate) dependent variable with conditional density $h$, $x$ is a vector of independent variables, $\pi_k$ is the prior probability of component $k$, $\theta_k$ is the component specific parameter vector for the density function $f$, and $\psi=(\pi_1,,\ldots,\pi_K,\theta_1',\ldots,\theta_K')'$ is the vector of all parameters. If $f$ is a univariate normal density with component-specific mean $\beta_k'x$ and variance $\sigma^2_k$, we have $\theta_k = (\beta_k', \sigma_k^2)'$ and Equation~(\ref{eq:1}) describes a mixture of standard linear regression models, also called \emph{latent class regression} or \emph{cluster-wise regression} \citep{flexmix:DeSarbo+Cron:1988}. If $f$ is a member of the exponential family, we get a mixture of generalized linear models \citep{flexmix:Wedel+DeSarbo:1995}, known as \emph{GLIMMIX} models in the marketing literature \citep{flexmix:Wedel+Kamakura:2001}. For multivariate normal $f$ and $x\equiv1$ we get a mixture of Gaussians without a regression part, also known as \emph{model-based clustering}. The posterior probability that observation $(x,y)$ belongs to class $j$ is given by \begin{equation}\label{eq:3} \Prob(j|x, y, \psi) = \frac{\pi_j f(y | x, \theta_j)}{\sum_k \pi_k f(y | x, \theta_k)} \end{equation} The posterior probabilities can be used to segment data by assigning each observation to the class with maximum posterior probability. In the following we will refer to $f(\cdot|\cdot, \theta_k)$ as \emph{mixture components} or \emph{classes}, and the groups in the data induced by these components as \emph{clusters}. \subsection{Parameter estimation} \label{sec:parameter-estimation} The log-likelihood of a sample of $N$ observations $\{(x_1,y_1),\ldots,(x_N,y_N)\}$ is given by \begin{equation}\label{eq:4} \log L = \sum_{n=1}^N \log h(y_n|x_n,\psi) = \sum_{n=1}^N \log\left(\sum_{k = 1}^K \pi_kf(y_n|x_n,\theta_k) \right) \end{equation} and can usually not be maximized directly. The most popular method for maximum likelihood estimation of the parameter vector $\psi$ is the iterative EM algorithm \citep{flexmix:Dempster+Laird+Rubin:1977}: \begin{description} \item[Estimate] the posterior class probabilities for each observation \begin{displaymath} \hat p_{nk} = \Prob(k|x_n, y_n, \hat \psi) \end{displaymath} using Equation~(\ref{eq:3}) and derive the prior class probabilities as \begin{displaymath} \hat\pi_k = \frac1N \sum_{n=1}^N \hat p_{nk} \end{displaymath} \item[Maximize] the log-likelihood for each component separately using the posterior probabilities as weights \begin{equation}\label{eq:2} \max_{\theta_k} \sum_{n=1}^N \hat p_{nk} \log f(y_n | x_n, \theta_k) \end{equation} \end{description} The E- and M-steps are repeated until the likelihood improvement falls under a pre-specified threshold or a maximum number of iterations is reached. The EM algorithm cannot be used for mixture models only, but rather provides a general framework for fitting models on incomplete data. Suppose we augment each observation $(x_n,y_n)$ with an unobserved multinomial variable $z_n = (z_{n1},\ldots,z_{nK})$, where $z_{nk}=1$ if $(x_n,y_n)$ belongs to class $k$ and $z_{nk}=0$ otherwise. The EM algorithm can be shown to maximize the likelihood on the ``complete data'' $(x_n,y_n,z_n)$; the $z_n$ encode the missing class information. If the $z_n$ were known, maximum likelihood estimation of all parameters would be easy, as we could separate the data set into the $K$ classes and estimate the parameters $\theta_k$ for each class independently from the other classes. If the weighted likelihood estimation in Equation~(\ref{eq:2}) is infeasible for analytical, computational, or other reasons, then we have to resort to approximations of the true EM procedure by assigning the observations to disjoint classes and do unweighted estimation within the groups: \begin{displaymath} \max_{\theta_k} \sum_{n: z_{nk=1}} \log f(y_n | x_n, \theta_k) \end{displaymath} This corresponds to allow only 0 and 1 as weights. Possible ways of assigning the data into the $K$ classes are \begin{itemize} \item \textbf{hard} \label{hard} assignment to the class with maximum posterior probability $p_{nk}$, the resulting procedure is called maximizing the \emph{classification likelihood} by \cite{flexmix:Fraley+Raftery:2002}. Another idea is to do \item \textbf{random} assignment to classes with probabilities $p_{nk}$, which is similar to the sampling techniques used in Bayesian estimation (although for the $z_n$ only). \end{itemize} Well known limitations of the EM algorithm include that convergence can be slow and is to a local maximum of the likelihood surface only. There can also be numerical instabilities at the margin of parameter space, and if a component gets to contain only a few observations during the iterations, parameter estimation in the respective component may be problematic. E.g., the likelihood of Gaussians increases without bounds for $\sigma^2\to 0$. As a result, numerous variations of the basic EM algorithm described above exist, most of them exploiting features of special cases for $f$. \subsection{Identifiability} \label{sec:identifiability} An open question is still identifiability of many mixture models. A comprehensive overview of this topic is beyond the scope of this paper, however, users of mixture models should be aware of the problem: \begin{description} \item[Relabelling of components:] Mixture models are only identifiable up to a permutation of the component labels. For EM-based approaches this only affects interpretation of results, but is no problem for parameter estimation itself. \item[Overfitting:] If a component is empty or two or more components have the same parameters, the data generating process can be represented by a smaller model with fewer components. This kind of unidentifiability can be avoided by requiring that the prior weights $\pi_k$ are not equal to zero and that the component specific parameters are different. \item[Generic unidentifiability:] It has been shown that mixtures of univariate normal, gamma, exponential, Cauchy and Poisson distributions are identifiable, while mixtures of discrete or continuous uniform distributions are not identifiable. A special case is the class of mixtures of binomial and multinomial distributions which are only identifiable if the number of components is limited with respect to, e.g., the number of observations per person. See \cite{flexmix:Everitt+Hand:1981}, \cite{flexmix:Titterington+Smith+Makov:1985}, \cite{flexmix:Grun:2002} and references therein for details. \end{description} FlexMix tries to avoid overfitting because of vanishing prior probabilities by automatically removing components where the prior $\pi_k$ falls below a user-specified threshold. Automated diagnostics for generic identifiability are currently under investigation. Relabelling of components is in some cases more of a nuisance than a real problem (``component 2 of the first run may be component 3 in the second run''), more serious are interactions of component relabelling and categorical predictor variables, see \cite{flexmix:Grun+Leisch:2004} for a discussion and how bootstrapping can be used to assess identifiability of mixture models. \pagebreak[4] \section{Using FlexMix} \label{sec:using-flexmix} \SweaveOpts{width=12,height=8,eps=FALSE,keep.source=TRUE} The standard M-step \texttt{FLXMRglm()} of FlexMix is an interface to R's generalized linear modelling facilities (the \texttt{glm()} function). As a simple example we use artificial data with two latent classes of size \Sexpr{nrow(NPreg)/2} each: \begin{center} \begin{tabular}{ll} Class~1: & $ y = 5x+\epsilon$\\ Class~2: & $ y = 15+10x-x^2+\epsilon$\\ \end{tabular} \end{center} with $\epsilon\sim N(0,9)$ and prior class probabilities $\pi_1=\pi_2=0.5$, see the left panel of Figure~\ref{fig:npreg}. We can fit this model in \R{} using the commands <<>>= library("flexmix") data("NPreg") m1 <- flexmix(yn ~ x + I(x^2), data = NPreg, k = 2) m1 @ and get a first look at the estimated parameters of mixture component~1 by <<>>= parameters(m1, component = 1) @ and <<>>= parameters(m1, component = 2) @ for component~2. The paramter estimates of both components are close to the true values. A cross-tabulation of true classes and cluster memberships can be obtained by <<>>= table(NPreg$class, clusters(m1)) @ The summary method <<>>= summary(m1) @ gives the estimated prior probabilities $\hat\pi_k$, the number of observations assigned to the corresponding clusters, the number of observations where $p_{nk}>\delta$ (with a default of $\delta=10^{-4}$), and the ratio of the latter two numbers. For well-seperated components, a large proportion of observations with non-vanishing posteriors $p_{nk}$ should also be assigned to the corresponding cluster, giving a ratio close to 1. For our example data the ratios of both components are approximately 0.7, indicating the overlap of the classes at the cross-section of line and parabola. \begin{figure}[htbp] \centering <>= par(mfrow=c(1,2)) plot(yn~x, col=class, pch=class, data=NPreg) plot(yp~x, col=class, pch=class, data=NPreg) @ \caption{Standard regression example (left) and Poisson regression (right).} \label{fig:npreg} \end{figure} Histograms or rootograms of the posterior class probabilities can be used to visually assess the cluster structure \citep{flexmix:Tantrum+Murua+Stuetzle:2003}, this is now the default plot method for \texttt{"flexmix"} objects \citep{flexmix:Leisch:2004}. Rootograms are very similar to histograms, the only difference is that the height of the bars correspond to square roots of counts rather than the counts themselves, hence low counts are more visible and peaks less emphasized. \begin{figure}[htbp] \centering <>= print(plot(m1)) @ \caption{The plot method for \texttt{"flexmix"} objects, here obtained by \texttt{plot(m1)}, shows rootograms of the posterior class probabilities.} \label{fig:root1} \end{figure} Usually in each component a lot of observations have posteriors close to zero, resulting in a high count for the corresponing bin in the rootogram which obscures the information in the other bins. To avoid this problem, all probabilities with a posterior below a threshold are ignored (we again use $10^{-4}$). A peak at probability 1 indicates that a mixture component is well seperated from the other components, while no peak at 1 and/or significant mass in the middle of the unit interval indicates overlap with other components. In our simple example the components are medium well separated, see Figure~\ref{fig:root1}. Tests for significance of regression coefficients can be obtained by <<>>= rm1 <- refit(m1) summary(rm1) @ Function \texttt{refit()} fits weighted generalized linear models to each component using the standard \R{} function \texttt{glm()} and the posterior probabilities as weights, see \texttt{help("refit")} for details. The data set \texttt{NPreg} also includes a response from a generalized linear model with a Poisson distribution and exponential link function. The two classes of size \Sexpr{nrow(NPreg)/2} each have parameters \begin{center} \begin{tabular}{ll} Class~1: & $ \mu_1 = 2 - 0.2x$\\ Class~2: & $ \mu_2 = 1 + 0.1x$\\ \end{tabular} \end{center} and given $x$ the response $y$ in group $k$ has a Poisson distribution with mean $e^{\mu_k}$, see the right panel of Figure~\ref{fig:npreg}. The model can be estimated using <>= options(width=55) @ <<>>= m2 <- flexmix(yp ~ x, data = NPreg, k = 2, model = FLXMRglm(family = "poisson")) summary(m2) @ <>= options(width=65) @ \begin{figure}[htbp] \centering <>= print(plot(m2)) @ \caption{\texttt{plot(m2)}} \label{fig:root2} \end{figure} Both the summary table and the rootograms in Figure~\ref{fig:root2} clearly show that the clusters of the Poisson response have much more overlap. For our simple low-dimensional example data the overlap of the classes is obvious by looking at scatterplots of the data. For data in higher dimensions this is not an option. The rootograms and summary tables for \texttt{"flexmix"} objects work off the densities or posterior probabilities of the observations and thus do not depend on the dimensionality of the input space. While we use simple 2-dimensional examples to demonstrate the techniques, they can easily be used on high-dimensional data sets or models with complicated covariate structures. \subsection{Multiple independent responses} \label{sec:mult-indep-resp} If the response $y=(y_1,\ldots,y_D)'$ is $D$-dimensional and the $y_d$ are mutually independent the mixture density in Equation~(\ref{eq:1}) can be written as \begin{eqnarray*} h(y|x,\psi) &=& \sum_{k = 1}^K \pi_k f(y|x,\theta_k)\\ &=& \sum_{k = 1}^K \pi_k \prod_{d=1}^D f_d(y|x,\theta_{kd}) \end{eqnarray*} To specify such models in FlexMix we pass it a list of models, where each list element corresponds to one $f_d$, and each can have a different set of dependent and independent variables. To use the Gaussian and Poisson responses of data \texttt{NPreg} simultaneously, we use the model specification \begin{Sinput} > m3 = flexmix(~x, data=NPreg, k=2, + model=list(FLXMRglm(yn~.+I(x^2)), + FLXMRglm(yp~., family="poisson"))) \end{Sinput} <>= m3 <- flexmix(~ x, data = NPreg, k = 2, model=list(FLXMRglm(yn ~ . + I(x^2)), FLXMRglm(yp ~ ., family = "poisson"))) @ Note that now three model formulas are involved: An overall formula as first argument to function \texttt{flexmix()} and one formula per response. The latter ones are interpreted relative to the overall formula such that common predictors have to be specified only once, see \texttt{help("update.formula")} for details on the syntax. The basic principle is that the dots get replaced by the respective terms from the overall formula. The rootograms show that the posteriors of the two-response model are shifted towards 0 and 1 (compared with either of the two univariate models), the clusters are now well-separated. \begin{figure}[htbp] \centering <>= print(plot(m3)) @ \caption{\texttt{plot(m3)}} \label{fig:root3} \end{figure} \subsection{Repeated measurements} \label{sec:repe-meas} If the data are repeated measurements on $M$ individuals, and we have $N_m$ observations from individual $m$, then the log-likelihood in Equation~(\ref{eq:4}) can be written as \begin{displaymath} \log L = \sum_{m=1}^M \sum_{n=1}^{N_m} \log h(y_{mn}|x_{mn},\psi), \qquad \sum_{m=1}^M N_m = N \end{displaymath} and the posterior probability that individual $m$ belongs to class $j$ is given by \begin{displaymath} \Prob(j|m) = \frac{\pi_j \prod_{n=1}^{N_m} f(y_{mn} | x_{mn}, \theta_j)}{\sum_k \pi_k \prod_{n=1}^{N_m} f(y_{mn} | x_{mn}, \theta_k)} \end{displaymath} where $(x_{mn}, y_{mn})$ is the $n$-th observation from individual $m$. As an example, assume that the data in \texttt{NPreg} are not 200 independent observations, but 4 measurements each from 50 persons such that $\forall m: N_m=4$. Column \texttt{id2} of the data frame encodes such a grouping and can easily be used in FlexMix: <<>>= m4 <- flexmix(yn ~ x + I(x^2) | id2, data = NPreg, k = 2) summary(m4) @ Note that convergence of the EM algorithm is much faster with grouping and the two clusters are now perfectly separated. \subsection{Control of the EM algorithm} \label{sec:control-em-algorithm} Details of the EM algorithm can be tuned using the \texttt{control} argument of function \texttt{flexmix()}. E.g., to use a maximum number of 15 iterations, report the log-likelihood at every 3rd step and use hard assignment of observations to clusters (cf. page~\pageref{hard}) the call is <<>>= m5 <- flexmix(yn ~ x + I(x^2), data = NPreg, k = 2, control = list(iter.max = 15, verbose = 3, classify = "hard")) @ Another control parameter (\texttt{minprior}, see below for an example) is the minimum prior probability components are enforced to have, components falling below this threshold (the current default is 0.05) are removed during EM iteration to avoid numerical instabilities for components containing only a few observations. Using a minimum prior of 0 disables component removal. \subsection{Automated model search} In real applications the number of components is unknown and has to be estimated. Tuning the minimum prior parameter allows for simplistic model selection, which works surprisingly well in some situations: <<>>= m6 <- flexmix(yp ~ x + I(x^2), data = NPreg, k = 4, control = list(minprior = 0.2)) m6 @ Although we started with four components, the algorithm converged at the correct two component solution. A better approach is to fit models with an increasing number of components and compare them using AIC or BIC. As the EM algorithm converges only to the next local maximum of the likelihood, it should be run repeatedly using different starting values. The function \texttt{stepFlexmix()} can be used to repeatedly fit models, e.g., <<>>= m7 <- stepFlexmix(yp ~ x + I(x^2), data = NPreg, control = list(verbose = 0), k = 1:5, nrep = 5) @ runs \texttt{flexmix()} 5 times for $k=1,2,\ldots,5$ components, totalling in 25 runs. It returns a list with the best solution found for each number of components, each list element is simply an object of class \texttt{"flexmix"}. To find the best model we can use <<>>= getModel(m7, "BIC") @ and choose the number of components minimizing the BIC. \section{Extending FlexMix} \label{sec:extending-flexmix} One of the main design principles of FlexMix was extensibility, users can provide their own M-step for rapid prototyping of new mixture models. FlexMix was written using S4 classes and methods \citep{flexmix:Chambers:1998} as implemented in \R{} package \pkg{methods}. The central classes for writing M-steps are \texttt{"FLXM"} and \texttt{"FLXcomponent"}. Class \texttt{"FLXM"} specifies how the model is fitted using the following slots: \begin{description} \item[fit:] A \texttt{function(x,y,w)} returning an object of class \texttt{"FLXcomponent"}. \item[defineComponent:] Expression or function constructing the object of class \texttt{"FLXcomponent"}. \item[weighted:] Logical, specifies if the model may be fitted using weighted likelihoods. If \texttt{FALSE}, only hard and random classification are allowed (and hard classification becomes the default). \item[formula:] Formula relative to the overall model formula, default is \verb|.~.| \item[name:] A character string describing the model, this is only used for print output. \end{description} The remaining slots of class \texttt{"FLXM"} are used internally by FlexMix to hold data, etc. and omitted here, because they are not needed to write an M-step driver. The most important slot doing all the work is \texttt{fit} holding a function performing the maximum likelihood estimation described in Equation~(\ref{eq:2}). The \texttt{fit()} function returns an object of class \texttt{"FLXcomponent"} which holds a fitted component using the slots: \begin{description} \item[logLik:] A \texttt{function(x,y)} returning the log-likelihood for observations in matrices \texttt{x} and \texttt{y}. \item[predict:] A \texttt{function(x)} predicting \texttt{y} given \texttt{x}. \item[df:] The degrees of freedom used by the component, i.e., the number of estimated parameters. \item[parameters:] An optional list containing model parameters. \end{description} In a nutshell class \texttt{"FLXM"} describes an \emph{unfitted} model, whereas class \texttt{"FLXcomponent"} holds a \emph{fitted} model. \lstset{frame=trbl,basicstyle=\small\tt,stepnumber=5,numbers=left} \begin{figure}[tb] \centering \begin{minipage}{0.94\textwidth} \lstinputlisting{mymclust.R} \end{minipage} \caption{M-step for model-based clustering: \texttt{mymclust} is a simplified version of the standard FlexMix driver \texttt{FLXmclust}.} \label{fig:mymclust.R} \end{figure} \subsection{Writing an M-step driver} \label{sec:writing-an-m} Figure~\ref{fig:mymclust.R} shows an example driver for model-based clustering. We use function \texttt{dmvnorm()} from package \pkg{mvtnorm} for calculation of multivariate Gaussian densities. In line~5 we create a new \texttt{"FLXMC"} object named \texttt{retval}, which is also the return value of the driver. Class \texttt{"FLXMC"} extends \texttt{"FLXM"} and is used for model-based clustering. It contains an additional slot with the name of the distribution used. All drivers should take a formula as their first argument, this formula is directly passed on to \texttt{retval}. In most cases authors of new FlexMix drivers need not worry about formula parsing etc., this is done by \texttt{flexmix} itself. In addition we have to declare whether the driver can do weighted ML estimation (\texttt{weighted=TRUE}) and give a name to our model. The remainder of the driver creates a \texttt{fit()} function, which takes regressors \texttt{x}, response \texttt{y} and weights \texttt{w}. For multivariate Gaussians the maximum likelihood estimates correspond to mean and covariance matrix, the standard R function \texttt{cov.wt()} returns a list containing estimates of the weighted covariance matrix and the mean for given data. Our simple example performs clustering without a regression part, hence $x$ is ignored. If \texttt{y} has $D$ columns, we estimate $D$ parameters for the mean and $D(D-1)/2$ parameters for the covariance matrix, giving a total of $(3D+D^2)/2$ parameters (line~11). As an additional feature we allow the user to specify whether the covariance matrix is assumed to be diagonal or a full matrix. For \texttt{diagonal=TRUE} we use only the main diagonal of the covariance matrix (line~14) and the number of parameters is reduced to $2D$. In addition to parameter estimates, \texttt{flexmix()} needs a function calculating the log-likelihood of given data $x$ and $y$, which in our example is the log-density of a multivariate Gaussian. In addition we have to provide a function predicting $y$ given $x$, in our example simply the mean of the Gaussian. Finally we create a new \texttt{"FLXcomponent"} as return value of function \texttt{fit()}. Note that our internal functions \texttt{fit()}, \texttt{logLik()} and \texttt{predict()} take only \texttt{x}, \texttt{y} and \texttt{w} as arguments, but none of the model-specific parameters like means and covariances, although they use them of course. \R{} uses \emph{lexical scoping} rules for finding free variables \citep{flexmix:Gentleman+Ihaka:2000}, hence it searches for them first in the environment where a function is defined. E.g., the \texttt{fit()} function uses the variable \texttt{diagonal} in line~24, and finds it in the environment where the function itself was defined, which is the body of function \texttt{mymclust()}. Function \texttt{logLik()} uses the list \texttt{para} in lines~8 and 9, and uses the one found in the body of \texttt{defineComponent()}. Function \texttt{flexmix()} on the other hand never sees the model parameters, all it uses are function calls of form \texttt{fit(x,y,w)} or \texttt{logLik(x,y)}, which are exactly the same for all kinds of mixture models. In fact, it would not be necessary to even store the component parameters in the \texttt{"FLXcomponent"} object, they are there only for convenience such that users can easily extract and use them after \texttt{flexmix()} has finished. Lexical scope allows to write clean interfaces in a very elegant way, the driver abstracts all model details from the FlexMix main engine. \subsection{Example: Using the driver} \label{sec:example:-model-based} \SweaveOpts{width=12,height=6,eps=FALSE} <>= library("flexmix") set.seed(1504) options(width=60) grDevices::ps.options(family="Times") suppressMessages(require("ellipse")) suppressMessages(require("mvtnorm")) source("mymclust.R") @ As a simple example we use the four 2-dimensional Gaussian clusters from data set \texttt{Nclus}. Fitting a wrong model with diagonal covariance matrix is done by <<>>= data("Nclus") m1 <- flexmix(Nclus ~ 1, k = 4, model = mymclust()) summary(m1) @ The result can be seen in the left panel of Figure~\ref{fig:ell}, the result is ``wrong'' because we forced the ellipses to be parallel to the axes. The overlap between three of the four clusters can also be inferred from the low ratio statistics in the summary table (around 0.5 for components 1, 3 and 4), while the much better separated upper left cluster has a much higher ratio of 0.85. Using the correct model with a full covariance matrix can be done by setting \texttt{diagonal=FALSE} in the call to our driver \texttt{mymclust()}: <<>>= m2 <- flexmix(Nclus ~ 1, k = 4, model = mymclust(diagonal = FALSE)) summary(m2) @ \begin{figure}[htbp] \centering <>= par(mfrow=1:2) plotEll(m1, Nclus) plotEll(m2, Nclus) @ \caption{Fitting a mixture model with diagonal covariance matrix (left) and full covariance matrix (right).} \label{fig:ell} \end{figure} \pagebreak[4] \section{Summary and outlook} \label{sec:summary} The primary goal of FlexMix is extensibility, this makes the package ideal for rapid development of new mixture models. There is no intent to replace packages implementing more specialized mixture models like \pkg{mclust} for mixtures of Gaussians, FlexMix should rather be seen as a complement to those. By interfacing R's facilities for generalized linear models, FlexMix allows the user to estimate complex latent class regression models. Using lexical scope to resolve model-specific parameters hides all model details from the programming interface, FlexMix can in principle fit almost arbitrary finite mixture models for which the EM algorithm is applicable. The downside of this is that FlexMix can in principle fit almost arbitrary finite mixture models, even models where no proper theoretical results for model identification etc.\ are available. We are currently working on a toolset for diagnostic checks on mixture models to test necessary identifiability conditions for those cases where results are available. We also want to implement newer variations of the classic EM algorithm, especially for faster convergence. Another plan is to have an interactive version of the rootograms using \texttt{iPlots} \citep{flexmix:Urbanek+Theus:2003} such that the user can explore the relations between mixture components, possibly linked to background variables. Other planned extensions include covariates for the prior probabilities and to allow to mix different distributions for components, e.g., to include a Poisson point process for background noise. \section*{Computational details} <>= SI <- sessionInfo() pkgs <- paste(sapply(c(SI$otherPkgs, SI$loadedOnly), function(x) paste("\\\\pkg{", x$Package, "} ", x$Version, sep = "")), collapse = ", ") @ All computations and graphics in this paper have been done using \proglang{R} version \Sexpr{getRversion()} with the packages \Sexpr{pkgs}. \section*{Acknowledgments} This research was supported by the Austrian Science Foundation (FWF) under grant SFB\#010 (`Adaptive Information Systems and Modeling in Economics and Management Science'). Bettina Gr\"un has modified the original version to include and reflect the changes of the package. \bibliography{flexmix} \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: