Fitting the Cusp Catastrophe in R : A cusp -Package Primer.

This vignette for the cusp package for R is a altered version of (Grasman, van der Maas, and Wagenmakers 2009) published in the Journal of Statistical Software . Of the seven elementary catastrophes in catastrophe theory, the \cusp" model is the most widely applied. Most applications are however qualitative. Quantitative techniques for catastrophe modeling have been developed, but so far the limited availability of (cid:13)exible software has hindered quantitative assessment. We present a package that implements and extends the method of (Cobb and Watson 1980; Cobb, Koppstein, and Chen 1983), and makes it easy to quantitatively (cid:12)t and compare di(cid:11)erent cusp catastrophe models in a statistically principled way. After a short introduction to the cusp catastrophe, we demonstrate the package with two instructive examples.

Stochastic formulations of catastrophe theory have been found, and statistical methods have been developed that allow quantitative comparison of catastrophe models with data (Cobb and Ragade 1978;Cobb and Watson 1980;Cobb 1981;Cobb et al. 1983;Guastello 1982;Oliva, Desarbo, Day, and Jedidi 1987;Lange, Oliva, and McDade 2000;Wagenmakers, Molenaar, Grasman, Hartelman, and van der Maas 2005a;Guastello 1992). Of these methods, the maximum likelihood approach of Cobb and Watson (1980); Cobb et al. (1983) is arguably most appealing for reasons that we discuss in the Appendix A.
In this paper we present an add-on package for the statistical computing environment R (R Development Core Team 2009) which implements the method of (Cobb et al. 1983), and extends it in a number of ways. In particular, the approach of Oliva et al. (1987) is adopted to allow for a behavioral variable that is embedded in multivariate response space. Furthermore, the implementation incorporates many of the suggestions of Hartelman (1997), and indeed intends to succeed and update the "cuspfit" FORTRAN program of that author. The package is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/ package=cusp

Cusp catastrophe
Consider a (deterministic) dynamical system that obeys equations of motion of the form In this equation y(t) represents the system's state variable(s), and c represents one or multiple (control) parameter(s) who's value(s) determine the specific structure of the system. If the system state y is at a point where ∂V (y; c)/dx = 0 the system is in equilibrium. If the system is at a non-equilibrium point, the system will move towards an equilibrium point where the function V (y; c) acquires a minimum with respect to y. These equilibrium points are stable equilibrium points because the system will return to such a point after a small perturbation to the system's state. Equilibrium points that correspond to maxima of V (y; c) are unstable equilibrium points because a perturbation of the system's state will cause the system to move away from the equilibrium point towards a stable equilibrium point. In the physical world no system is fully isolated and there are always forces impinging on a system. It is therefore highly unlikely to encounter a system in an unstable equilibrium state. Equilibrium points that correspond neither to maxima nor to minima of V (y; c), at which the Hessian matrix (∂ 2 V (y)/∂y i ∂y j ) has eigenvalues equal to zero, are called degenerate equilibrium points, and these are the points at which a system can give rise to unexpected bifurcations in its equilibrium states when the control variables of the system are changed.
or "universal unfoldings", in only one or two canonical state variables. These universal unfoldings are called the "elementary catastrophes". The canonical state variables are obtained by smooth transformations of the original state variables. The elementary catastrophes constitute the different families of catastrophe models.
In the biological and behavioral sciences, the so-called cusp catastrophe model has been applied most frequently, as it is the simplest of catastrophe models that exhibits discontinuous transitions in equilibrium states as control parameters are varied. The canonical form of the potential function for the cusp catastrophe is − V (y; α, β) = αy + 1 2 β y 2 − 1 4 y 4 . (2) Its equilibrium points, as a function of the control parameters α and β, are solutions to the equation α + β y − y 3 = 0.
This equation has one solution if δ = 27α − 4β 3 , which is known as Cardan's discriminant, is greater than zero, and has three solutions if δ < 0. These solutions are depicted in Figure 1 as a two dimensional surface living in three dimensional space, the floor of which is the two dimensional (α, β) coordinate system called the "control plane". The set of values of α and β for which δ = 0 demarcates the bifurcation set, the cross hatched cusp shaped region on the floor in Figure 1. From a regression perspective, the cusp equilibrium surface may be conceived of as response surface, the height of which predicts the value of the dependent variable y given the values of the control variables. This response surface has the peculiar property however, that for some values of the control variables α and β the surface predicts two possible values instead of one. In addition, this response surface has the unusual characteristic that it "antipredicts" an intermediate value for these values of the control variables-that is, the surface predicts that certain state values, viz. unstable equilibrium states, should not occur (Cobb 1980). As indicated, the dependent variable y is not necessarily an observed quantity that characterizes the system under study, but is in fact a canonical variable that in general depends on a number of actually measurable dependent variables. Similarly, the control coordinates α and β represent canonical coordinates that are dependent on actual measured or controlled independent variables. The α coordinate is called the "normal" or "asymmetry" coordinate, while the β coordinate is called the "bifurcation" or "splitting" coordinate (Stewart and Peregoy 1983).

Cusp catastrophe as a model
In evaluating the cusp as a model for data there are two complementary approaches. The first approach evaluates whether certain qualitative phenomena occur in the system under consideration. In the second approach, a parameterized cusp is fitted to the data. Gilmore (1993) derived a number of qualitative behavioral characteristics of the cusp model; the so-called catastrophe flags. Among the more prominent are sudden jumps in the value of the (canonical) state variables; hysteresis-i.e., memory for the path through the phase space of the system; and multi-modality-i.e., the simultaneously presence of multiple preferred states. Verification of the presence of these flags constitute one important stage in gathering evidence for the presence of a cusp catastrophe in the system under scrutiny. For extensive discussions of these qualitative flags we refer to Gilmore (1993) (see also Stewart and Peregoy 1983;van der Maas and Molenaar 1992;van der Maas, Kolstein, and van der Pligt 2003). Most applications of catastrophe theory in general, and the cusp catastrophe in particular, have focused entirely on this qualitative verification. An alternative approach is constituted by a quantitative evaluation of an actual match between a cusp catastrophe model and the data using statistical fitting procedures. This is indispensable for sound verification that a cusp catastrophe does a better job at describing the data than other conceivable models for the data.
A difficulty that arises when constructing empirically testable catastrophe models is the fact that catastrophe theory applies to deterministic systems as described by equation (1). Being inherently deterministic, catastrophe theory cannot be applied directly to systems that are subject to random influences which is commonly the case for real physical systems, especially in the biological and behavioral sciences.
To bridge the gap between determinism of catastrophe theory and applications in stochastic environments, Loren Cobb and his colleagues (Cobb and Ragade 1978;Cobb 1980;Cobb and Watson 1980;Cobb and Zacks 1985) proposed to turn catastrophe theory into a stochastic catastrophe theory by adding to equation (1) a (white noise) Wiener process, dW (t), with variance σ 2 , and to treat the resulting equation as a stochastic differential equation (SDE): This SDE is then associated with a probability density that describes the distribution of the system's states on any moment in time, which may be expressed as Here ψ is a normalizing constant, and λ merely determines the origin of scale of the state variable. In this stochastic context β is called the bifurcation factor, as it determines the number of modes of the density function, while α is called "asymmetry" factor as it determines the direction of the skew of the density (the density is symmetric if α = 0 and becomes left or right skewed depending on the sign of α; Cobb 1980). The function is implemented in the R package described below as dcusp(). Figure 2 displays the density for different regions of the control plane.

Estimation methods
As indicated earlier, the state variable y of the cusp is a "canonical variable". This means that it is a (generally unknown) smooth transformation of the actual state variables of the system. If we have a set of measured dependent variables Y 1 , Y 2 , . . . , Y p , to a first order approximation we may say where w 0 , w 1 , . . . , w p are the first order coefficients of a polynomial approximation to the "true" smooth transformation. Similarly, the parameters α and β are "canonical variates" in the sense that they are (generally unknown) smooth transformations of actual control variates. Again to first approximation, for experimental parameters or measured independent variables X 1 , . . . , X q , we may write α = a 0 + a 1 X 1 + a 2 X 2 + a q X q ,  Fitting the cusp model to empirical data then reduces to estimating the parameters w 0 , w 1 , . . ., w p , a 0 , . . . , a q , b 0 , . . . , b q .
On the basis of equation (4) and the associated stochastic differential equation discussed earlier, Cobb (1981) and Cobb et al. (1983) respectively developed method of moments estimators and maximum likelihood estimators. The maximum likelihood method of Cobb (Cobb and Watson 1980;Cobb and Zacks 1985) has not been commonly used however. Two important reasons for this are instability of Cobb's software for fitting the cusp density, and difficulties in its use (see van der Maas et al. 2003 for a discussion; see the appendix for a discussion of other estimation methods that have been proposed in the literature).
Cobb's methods assume that the state variable is directly accessible to measurement. As argued by Oliva et al. (1987), in the behavioral sciences both dependent (state) as well as independent (control) variables are more often than not constructs which cannot be easily measured directly. It is therefore important to incorporate equation (5) so that the state variable that adheres to the cusp catastrophe may be embedded in the linear space spanned by a set of dependent variables.
In the cusp package that we describe below, we use the maximum likelihood approach of Cobb and Watson (1980), augmented with the subspace fitting method (Equation 5) proposed by (Oliva et al. 1987)

Package description
The cusp package is a package for the statistical computing language and environment R (R Development Core Team 2009). The package contains a number of functions for use with catastrophe modeling, including utility functions to generate observations from the cusp density, to evaluate the cusp density and cumulative distribution function, to fit the cusp catastrophe, to evaluate the model fit, and to display the results. The core user interface functions of the package are listed in Table 1.

Model specification
As is the with many model fitting routines R, the fitting routines in the cusp package allow the user to specify models in terms of dependent and independent variables in a compact symbolic form. Basic to the formation of such models is the~(tilde) operator. An expression of the form y~model signifies that the response y is modeled by a linear predictor that is specified in model. In the cusp package, the dependent variables are however always y, α, and β, in accordance with Equations (8-10) discussed below. Model formulas are used to generate an appropriate design matrix. It should be mentioned that R makes a strict distinction between "factors" and variables. The former are study design factors that have an enumerable number of levels (e.g., education level, treatment level, etc.), whereas the former are continuous variates (e.g., age, heart rate, response time). Factors and variates are (appropriately) treated differently in the constructing the design matrix.

Cost function (likelihood function)
At the heart of the package is the fitting routine that performs maximum likelihood estimation of all the parameters in Equations (5-7). That is, for observed dependent variables Generates a graphical display of the fit, including the estimated locations on the cusp control surface for each observation, nonparametric density estimates for different regions of the control surface, and a residuals plot. An example is provided in Figure 3. The graph is fully customizable. plot(fit) cusp3d Generates a three dimensional display of the cusp equilibrium surface on which the estimated states are displayed. The graphic is fully customizable. Figure 5 displays an example. cusp3d(fit) cusp3d.surface Generates a three dimensional display of the cusp equilibrium surface. Figure 1 was generated with it.
Y i1 , Y i2 , . . . , Y ip , and independent variables X i1 , X i2 , . . . , X iq , for subjects i = 1, . . . , n, the distribution of If we collect the observed dependent variables in the n × where 1 n is an n-vector whose entries all equal 1, and furthermore collect the coefficients in the vectors w = (w 0 , w 1 , . . . , w p ) , a = (a 0 , a 1 , . . . , a q ) , and b = (b 0 , b 1 , . . . , b q ) , where denotes transpose, we can write these equations succinctly as where y = (y 1 , . . . , y n ) , α = (α 1 , . . . , α n ) , and β = (β 1 , . . . , β n ) . Therefore, the negative log-likelihood for a sample of observed values ( Note that compared to equation (4), here we have absorbed the location and scale parameters λ and σ into the coefficients w 0 , w 1 , . . . , w p . It should be noted furthermore that there is an ambiguity in the equations above: The signs of the a j and w j may all be switched without affecting the value negative log-likelihood function. This is the case, because the quadratic and quartic terms in y i do not determine the sign of y i (i.e., y 2 i = (−y i ) 2 and y 4 i = (−y i ) 4 ), while for the linear term α i y i = (−α i )(−y i ) holds. Thus the sign of y i (and hence, the sign for the w j 's) can be switched without affecting the value of (11) if the sign of α i , and hence the signs of the a j 's, is switched as well. The estimates of a 1 , . . . , a q , w 1 , . . . , w p are therefore identifiable up to a change in sign.
The cusp routine of the package minimizes L with respect to the parameters w 0 , . . ., w p , a 0 , . . . , a p , b 0 , . . . , b q . The computational load of evaluating L is mainly burdened by the calculation of the normalization constants ψ i , which has to be carried out numerically. To this end, we use an adaptive quadrature routine to minimize the computational cost. To speed up the calculations substantively, this is carried out in linked C code. Speed is especially important from the user experience perspective, if different models are to be tried and compared.
Internally, the data are standardized using a QR decomposition. This is both for stability of the estimation algorithm, as well as handling collinear predictors in the design matrix. Only coefficients for each dimension of the column space of the design matrix are estimated. The estimated coefficients are related back to the actual dependent variables. If the design matrix does not have full column rank, coefficients for some independent variables-those that are fully explained by others-are not determined. For which of the independent variables coefficients are estimated in such collinear cases depends partly on the order of their occurrence in the formula.
Note that, because the models for α i 's and β i 's are specified separately, different sets of independent variables can figure in each, hence allowing for confirmatory data analysis. This has been considered a deficiency the method of Cobb and Watson (Stewart and Peregoy 1983;Alexander, Herbert, DeShon, and Hanges 1992), although this option has been present in Hartelman's (1997) program.

Optimization algorithm
The negative log-likelihood function is minimized using one of the built-in optimization routines of R. By default the limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm with bounds (L-BFGS-B; Zhu, Byrd, Lu, and Nocedal 1997) is used by the cusp function. Other methods provided by R can be used, but in our experience L-BFGS-B works best. The main problem is that the cusp density quickly decreases beyond numerical precision for |α|, |β| > 5, and hence some boundary restrictions are helpful. Boundaries (both lower and upper) default to −10 and 10, which is sufficiently large in combination with the internal standardization of the data in al the cases we have encountered. Different boundaries may be specified, but our experience with earlier applications of the method is that the values of |α| and |β| rarely exceed 3.

Starting point
The starting values used by default turn out to often render convergence of the optimization algorithm without problems. When convergence problems arise a warning is issued, and this may be interpreted as an invitation to the user to provide alternative starting values. In our experience convergence to a proper minimum of the negative log-likelihood can be achieved most of the time by providing alternative starting values. A strong indication of improper convergence is a warning about NaN's that is issued when summary statistics are computed. If this happens, the model should be refit using a different set of starting values.

Statistical evaluation of model fit
To assess the correspondence of data with the predictions made by the cusp catastrophe, a number of diagnostic tools have been suggested.

Maxwell vs delay convention and R 2
First of all, Cobb (1998; see also Stewart and Peregoy 1983;Hartelman 1997 have suggested a pseudo-R 2 as a measure of explained variance defined by This is clearly inspired by the familiar relation between the squared (multiple) correlation coefficient and the "explained variance" from ordinary regression. However, the term "error variance" needs to be defined, which is not trivial for the cusp catastrophe model.
As indicated previously, the cusp catastrophe is not an ordinary regression model. Rather, it is an implicit regression model of an irregular type. 1 As such, unlike ordinary regression models, for a given set of values of the independent variables the model may predict multiple values for the dependent variable. In ordinary regression the predicted value is usually the expected value of the dependent variable given the values of the independent variables. In the case of the cusp density, for certain values of α and β the cusp density is bimodal, and the expected value of this density lies in a region of low probability between the two modi. That is, the mathematical expectation of the cusp density is a value that in itself is relatively unlikely to be observed. Two alternatives for the expected value as the predicted value can be used, which are closely related to a similar problem regarding interpretation conventions in deterministic catastrophe theory: One can choose the mode of the density closest to the state values as the predicted value, or one can use the mode at which the density is highest. The former is known as the delay convention, the latter as the Maxwell convention. Although in the physical sciences both have their uses, Watson (1980, Cobb 1998) suggest to use the delay convention (Stewart and Peregoy 1983). Both conventions are provided by the package, the default being the delay convention.
Hence, in the pseudo-R 2 statistic defined above, the error variance is defined as the variance of the differences between the observed (or estimated) states and the mode of the distribution that is closest to this value. It should be noted however that this pseudo-R 2 can become negative if many of the α i 's deviate from 0 in the same direction-in that case the distribution is strongly skewed, and deviation from the mode is on average larger than deviation from the mean. Negative pseudo-R 2 's are thus perfectly legitimate for the cusp density, which limits its value for model fit assessment. Alternatives to the pseudo-R 2 are discussed below.

AIC, BIC and logistic curve
In addition to the pseudo-R 2 statistic, to establish convincingly the presence of a cusp catastrophe, Cobb gives three guidelines for evaluating the model fit (Cobb 1998;Hartelman 1997): First, the fit of the cusp should be substantially better than multiple linear regression-i.e., its likelihood should be significantly higher than that of the ordinary regression model. Second, any of the coefficients w 1 , . . . , w p should deviate significantly from zero (w 0 does not have to), as well as at least one of the a j 's or b j 's. Thirdly, at least 10% of the (α i , β i ) pairs should lie within the bifurcation region. The former two guidelines can be assessed with the summary function of the packages; the latter guideline can be assessed with the plot function; both of these are detailed in the Examples section below.
A problem arises when the general case of equation (8) with more than one dependent variable is used: The linear regression model with which to compare the cusp model is not uniquely defined. The most natural approach seems to be linear subspace regression: The first canonical variate between the X ji 's and the Y ki 's is used as the predicted value, and the first canonical correlation is used for calculating the explained and residual variance. For models in which (8) where y i , α i , β i are defined in Equations (8-10), and the i 's are zero mean random disturbances. (One may assume i to be normally distributed, but this is not necessary; see e.g., Seber and Wild 1989.) The rationale for the logistic function is that this function does not posses degenerate critical points, while it does have the possibility to model arbitrarily steep changes in the (canonical) state variable as a function of minute changes in an independent variable, thus mimicking "sudden" transitions of the cusp (Hartelman 1997). The summary function in the package provides an option to fit this logistic curve to the data. Because the cusp density and the logistic functions are not nested models, the fit cannot be assessed on the basis of differences in the likelihood, and one has to resort to other indicators such as AIC and BIC. These fit indices are both computed by the summary function, as well as an AIC corrected for small sample sizes (AICc; Burnham and Anderson 2004).
Instead of the 10% guideline, Hartelman (1997) proposes to require that the AIC and BIC indicate a better fit for the cusp density than for the logistic curve. Wagenmakers, van der Maas, and Molenaar (2005b) suggest to use the BIC of both models to compute approximation of the posterior odds for the cusp relative to the logistic curve, assuming equal prior probabilities for both models.

Examples
For illustration purposes, we provide three example analyses with the package. The first two examples entail data that have been analyzed with cusp catastrophe methods before and have been published elsewhere.

Example I
The first example is taken from van der Maas et al. (2003), and concerns attitudinal response transitions with respect to the statement "The government must force companies to let their workers benefit from the profit as much as the shareholders do". Some 3000 Dutch respondents indicated their level of agreement with this statement on a 5 point scale (1 = totally agree, 5 = totally disagree). As a normal factor political orientation (measures on a 10 point scale from 1 = left wing to 10 = right wing) was used. As a bifurcation factor the total score on a 12 item "political involvement" scale was used. The theoretical social psychological details are discussed in (van der Maas et al. 2003). The data thus consist of a table with three columns: Orient, Involv, and Attitude. We use the same subset of 1387 cases that was also analyzed in van der Maas et al. (2003). In that paper an extensive comparison of models was done. These data are made available in the package and can be loaded using the data("attitudes") command. Here we only present the analysis for the best fitting model (  solution quickly for this example. Both were loaded with a data call. Note the use of formula's discussed earlier: The first formula specifies that the state variable y i is determined by the variable Attitude for which a location and scaling parameter (w 0 and w 1 in Equation 8) have to be estimated. The second formula states that the normal factor α i is determined by the variables Orient and Involv for which regression coefficients and an intercept must be estimated. Similarly the third formula which specifies that the bifurcation factor is determined by the variable Involv for which an intercept and a regression coefficient have to be estimated. When the cusp function returns, the estimates can be printed to the console window of R by typing print(fit).More informative however is a summary of the parameters and the associated statistics, which is obtained with the statement

R> summary(fit, logist = TRUE)
where logist is set to TRUE to compare the cusp to the logistic curve fit, in addition to a comparison with a linear model. Part of the result is display in Tables 2 and 3. There are only very small differences between the fit presented here and the fit presented in (van der Maas et al. 2003), which are entirely due to differences in optimization algorithm. It should be mentioned that in (van der Maas et al. 2003) the parameter estimates are the coefficients with respect to standardized data. To standardize the data the scale function of R can be used. The column headed "R 2 " in Table 3 lists the squared multiple correlation for the linear regression model and logistic curve model, and the pseudo-R 2 for the cusp catastrophe model.
A visual display of the data fit and diagnostic plots is generated with the command plot(fit); the result is displayed in Figure 3. The figure display the control plane along with the estimated (α i , β i ), i = 1, . . . , n, for each of the observations. Furthermore, it (optionally) displays for each of four regions in the control plane a density estimate of the estimated  Table 3: Model fit statistics for synthetic data example as obtained with summary. The column labeled "R 2 " gives conventional R 2 for the linear and logist model, and the pseudo-R 2 statistic for the cusp model. state variable in that region. These are displayed on the right in the figure: The top two panels reflect the densities in the region left and right of the bifurcation region respectively. These should be positively skewed for the left side and negatively skewed for the right side (compare Figure 2). The next panel display the density estimate for the state estimates below the bifurcation region. Here the density estimate should be approximately uni-modal and less skewed than in the other two regions. The last panel displays the density estimate for estimated state values within the bifurcation region. In this region the density should be bimodal. For the data displayed here the first three densities seem heavily multi-modal. This is due to the fact that the bifurcation factor "political orientation" was measured on a discrete ordinal 10 points scale, thus artifactually introducing modes around these values. The residual versus fitted plot displays the estimated errors as a function of predicted states. As indicated earlier, by default the errors are computed using the delay convention. Although with a good model fit one would expect to see no systematic relationship between the estimated errors and the predicted states, we have observed a consistent negative trend between the two-even when the data were generated from the cusp density in simulations. This may simply result from cases of which the density is strongly skewed. Hence, a slight negative trend should not be taken too seriously as an indication of model misfit.
Note from Figure 3 that the overwhelming majority of cases lie in regions where the cusp density is moderate to strongly skewed. This fact accounts for the negative pseudo-R 2 . Lange et al. (2000) recommend to reject the cusp model entirely when negative pseudo R 2 result. This seems to be an unwarranted recommendation because negative pseudo-R 2 values arebecause of the definition of this measure of fit-entirely possible for the cusp density. In fact, negative pseudo-R 2 statistics are possible for all (strongly) skewed distribution such as for instance a chi-square distribution. Rather than rejecting the cusp catastrophe model on the basis of a negative pseudo-R 2 , one might consider forgetting about this measure of fit entirely, as it clearly is not well-suited for its intended purposes in non-symmetrical distributions. Instead, one can use AIC and BIC to compare the cusp model with competitor models such as the linear regression model and the logistic curve.

Example II
In the second example we use the model specified by Oliva et al. (1987) to demonstrate the use of multiple state variables. The model corresponds to the data in Table 2 of that paper, which displays synthetic data for 50 cases with scores on two (dependent) state variables (Z 1 and Z 2 ), four (independent) bifurcation variables (Y 1 , Y 4 , Y 4 , and Y 4 ), and three (independent) normal variables (X 1 , X 2 , and X 3 ). The data can be loaded from the package with data(oliva).
The "true" model for their synthetic data is For testing purposes, Oliva et al. did not add noise to any of the variables, and hence, the data perfectly comply with their implicit regression equation (3).
Because we are using the stochastic approach of Cobb, we cannot use these deterministic data. We therefore generated data in accordance with the model in equation (14), where X 1 , X 2 , and X 3 were uniformly distributed on the interval (−2, 2), Y 1 , Y 2 and Z 1 were uniformly   Oliva et al. (1987) synthetic data example obtained with summary for the less informed model (Model II; see text for details). Note: "R 2 " value for cusp is pseudo-R 2 .
distributed on (−3, 3), and Y 3 and Y 4 were uniformly distributed on (−5, 5). The states y i were then generated from the cusp density with their respective α and β as normal and splitting factors, and then Z 2 was computed as Z 2i = (y i + 0.52 Z 1i )/(−1.60).
The call for fitting the model of Oliva et al. for the resulting data set is R> oliva.fit <-cusp(y~z1 + z2 -1, + alpha~x1 + x2 + x3 -1, + beta~y1 + y2 + y3 + y4 -1, + data = oliva) Note that in the model there are no intercept coefficients; this is signified in the cusp call by the -1 added to the formula's (see the formula section of the R Manual, R Development Core Team 2009 for details). The data are stored in the data frame oliva in this case. The result from summary(oliva.fit, logist = TRUE) is displayed in the upper part of Table 4.
Some of the estimates differ substantially from the true values. This should however not be too surprising given that there are 9 estimated parameters and only 50 observations, yielding a ratio of 55/9 ≈ 5.5 observations per parameter. Six coefficients differ significantly from zero: a x1 , a x2 , b y1 , b y3 , and both w z1 and w z2 . Using confint to calculate confidence intervals, we obtain for a x1 , a x2 , b y1 , b y3 , w z1 and w z2 95% confidence intervals of respectively (−1.26, −0.3979), (0.2864, 1.0676), (0.1749, 0.7197), (0.5658, 1.07), (0.4074, 0.589) and (1.333, 1.669), which all neatly cover the true values. The same confidence interval for b y4 however, is (−0.4271, 0.134), which does not contain the true value. This may indicate that the estimator is biased. Indeed Hartelman (1997) showed in simulations that the estimators are in fact biased. The control factors (α i 's and β i 's) were however recovered quite decently: Correlations between the actual α i 's and those predicted by the model fit (which can be obtained with the function predict) was .996, and the correlation between actual β i 's and those predicted by the model fit was .924.
In contrast to the model specified in this fit, which is rather informed on the variables and their role in the cusp catastrophe, it is often the case that one doesn't know which variables determine the splitting factor, and which variables determine the normal factor. We therefore also fitted a less informed model, in which both alpha and beta are specified as x1 + x2 + x3 + y1 + y2 + y3 + y4. That is, all independent variables are used to model both the splitting as well as the normal factor. The resulting estimates are given in the second part of Table 4. The same set of estimates differs significantly from zero, demonstrating the usefulness of having significance tests for each estimate separately. The confidence intervals all covered the true parameter value, except for the same parameter as in the informed model.  The fit indices for this less informed model are displayed in Table 5. Note that the pseudo-R 2 statistics indicates that the fit does not differ substantially between the cusp model and the logistic curve model. In fact, it even indicates that the logistic curve model gives a slightly better fit! Thus clearly, the pseudo R 2 is not in all cases a trustworthy guide in selecting a model, and since we are usually unaware of what the "correct" model is, this makes it difficult to rely on it at all. The AIC, AICc, and BIC on the other hand all (correctly) indicate that the cusp is the best model (of the ones compared) for these data. The chi-square Likelihood Ratio test given in the program output indicated that the cusp model fit was significantly better than the linear model with normal errors (χ 2 = 68.71, df = 9, p < .000).
A control plane scatter plot of the synthetic data for the Oliva et al. model fit is presented in Figure 4. A three dimensional display of the model fit as generated with cusp3d is presented in Figure 5. A couple of things may be noted from the scatter plot: First of all, the sizes of the dots differ. In fact the size of the dots, each of which corresponds to a single case, varies according to the observed bivariate density of the control factor values at the location of the point. A second observation is that cases that lie inside the bifurcation set are plotted indiscriminately of whether they lie on the upper surface or on the lower surface. To make it possible to visually distinguish these cases, the color of the points vary according to the value of the state variable; higher values are associate with more intense purple, lower values with more intense green.

Discussion
In this paper we have presented an add-on package for cusp catastrophe modeling in R. The core user interface functions allow the user to easily specify and fit a broad range of models using the maximum likelihood approach of Cobb (1980;Cobb and Watson 1980). It also provides the user a number of tools for the assessment of the cusp catastrophe model fit.
Although the focus of this paper was entirely on quantitative methodology, by no means do we consider the statistically satisfactory fit of a cusp catastrophe model-or any catastrophy model for that matter-definite evidence for the for the presence of dynamical phase transitions. Without concurrent qualitative assessment of the presence of catastrophe flags implied by the model, and without a sound theoretical framework that implicates an underlying gradient system near equilibrium states, any catatrophe model fit remains unconvincing.
Cobb's method is not without it limitations. Deterministic catastrophe classifies systems up to a set of smooth nonlinear scalings of the state variables. Hartelman (1997) points out that this poses a problem for Cobb's approach. The alternative approach offered in Hartelman (1997; see also Wagenmakers et al. 2005a) however, only works for time series. Practical experience with applications to cross-sectional data indicates that Cobb's method is suitable for these cases.
A number of improvement of the package can be thought of. The current package tests for the presence of bifurcation points by fitting a logistic curve that accommodates arbitrarily sudden changes in the state variable as a function of smooth gradual changes in the control variables. More ideally, to test for the presence of bifurcation points one should minimize the negative log-likelihood L of equation (11) subject to the n constraints δ i = α 2 i /4 − β 3 i /27 ≤ 0. The parameter δ i is known as Cardan's Discriminant, named after the 16th century Italian mathematician who first published it (Cobb and Zacks 1985). It is positive, only if the cusp equation (3) has three solutions-i.e., only if there are multiple equilibrium states. One then can use likelihood ratio chi-square testing to compare this model with the unconstrained model. This could be a viable alternative approach to the logistic curve fit, but it would require a different optimization routine. Unfortunately, R currently has no optimization routine that allows for arbitrary nonlinear inequality constraints.
One might further consider a Bayesian approach to estimation. The Bayesian approach however, comes of course at the cost of having to specify a prior belief about the likelihood of parameter values. Given the relations between parameters w 0 , w 1 , . . . , b q−1 , b q and the data in equations (8-10) the latter seems rather laborious in the most general case, and a far from intuitive enterprise. Only (relatively) uninformative priors seem straightforward, hence leaving out the heart of a true Bayesian approach.
These considerations will be explored in future improvements of the package.
parameters. In this regression procedure however, the occurrence of unstable equilibrium states is rewarded just as much as stable states, rather than punished (Hartelman 1997). Alexander et al. (1992) have furthermore criticized this approach because the dependent variable in the regression is a difference score of two variables, one of which is also present as a predictor. As a consequence, the explained variance can be up to 50% for completely random data. Thus the method cannot distinguish between a (cusp) catastrophe model and a linear model. In a reply to Alexander et al. (1992), Guastello (1992) demonstrates that his polynomial regression technique can give an R 2 estimate of 0.55 for purely random data. By constructing the bifurcation and asymmetry factors to be "known", as the author calls it, "a near-perfect [i.e., R 2 = 0.99] cusp model could be obtained" (Guastello 1992, p. 387). On the basis of an, in our view misguided, interpretation of chaos theory, the author argues that from these examples it can be concluded that these random numbers conform to a fold-and even a cusp-catastrophe. Clearly, this cannot be a "cusp" in the sense of Cobb's stochastic version of equation (3). The polynomial regression method of Guastello (1992) estimates the coefficients correctly if the data set is a time series, in which case we can show however, at least for an Ornstein-Uhlenbeck SDE, that theoretically R 2 ranges from 0 for closely spaced samples to 0.5 for widely spaces samples (when samples are nearly uncorrelated). For the cusp SDE a similar result can be demonstrated using simulation. A scenario under which the method would yield correct coefficient estimates and a high R 2 , is when multiple independent systems are observed on two occasions in which the systems are perturbed on the first occasion by at least two standard deviations of the equilibrium noise level.
The multivariate "GEMCAT" methodology of Oliva et al. (Oliva et al. 1987;Lange et al. 2000) assumes that the measured states are at the equilibrium surface and that any discrepancy is due to additive noise (much like Cobb's stochastic extension of catastrophe theory but, as we shall see, not quite the same). Hence, GEMCAT simply minimizes the square of the left hand side of equation (3), summed across all observed states. The qualification "multivariate" refers to the fact that GEMCAT allows for the use of equation (5), whereas the technique of Guastello (1988) presumes that the state variable y is accessible for direct measurement, and is not, as implicated by equation (5), a set of dependent variables that "predict" the state. As is true for the method of Guastello (1988), a problem with this approach is that stable and unstable equilibrium states are treated indiscriminately, and both types of states are "rewarded" for their presence in the model fit. In fact, in the multivariate synthetic data example, for eleven of the cases in Table 2 of (Oliva et al. 1987), the simulated observed states were unstable (i.e., inaccessible) equilibrium states. This is opposite to the predictions from the cusp catastrophe that the unstable equilibrium state is unlikely to be observed, and contrasts with Cobb's conception of stochastic catastrophe theory. Furthermore, as astutely noted in Hartelman (1997), the approach for model evaluation proposed in Oliva et al. (1987) is not valid because their stochastic counterpart of equation (3) as a model, entails an implicit equation to which data should adhere-save for random disturbances. Unlike conventional non-linear regression, implicit equations allow for multiple predictions for each set of values of the predictor variables. Equation (3) and its stochastic counterpart thus do not constitute non-linear regression in the conventional sense. Even more importantly, the conditions under which (asymptotic) statistical inference theory is developed for such implicit models (Seber and Wild 1989) precisely excludes those cases that are central to catastrophe theory (Hartelman 1997). As a consequence, the suggested chi-square statistics initially used in GEMCAT for model comparison and model selection are rendered invalid. In a new version of GEM-