Title: | Cusp-Catastrophe Model Fitting Using Maximum Likelihood |
---|---|
Description: | Cobb's maximum likelihood method for cusp-catastrophe modeling (Grasman, van der Maas, and Wagenmakers (2009) <doi:10.18637/jss.v032.i08>; Cobb (1981), Behavioral Science, 26(1), 75-78). Includes a cusp() function for model fitting, and several utility functions for plotting, and for comparing the model to linear regression and logistic curve models. |
Authors: | Raoul P. P. P. Grasman [aut, cre, cph]
|
Maintainer: | Raoul P. P. P. Grasman <[email protected]> |
License: | GPL-2 |
Version: | 2.3.8 |
Built: | 2025-02-14 03:40:32 UTC |
Source: | https://github.com/cran/cusp |
Fits cusp catastrophe to data using Cobb's maximum likelihood method with a different algorithm. The package contains utility functions for plotting, and for comparing the model to linear regression and logistic curve models. The package allows for multivariate response subspace modeling in the sense of the GEMCAT software of Oliva et al.
Package: | cusp |
Type: | Package |
Version: | 2.0 |
Date: | 2008-02-14 |
License: | GNU GPL v2 (or higher) |
This package helps fitting Cusp catastrophe models to data, as advanced in Cobb et al. (1985). The main functions are
cusp |
Fit Cobb's Cusp catastrophe model; see example below. |
summary.cusp |
Summary statistics of cusp model fit. |
confint.cusp |
Confidence intervals for parameter estimates |
plot.cusp |
Diagnostic plots for cusp model fit |
cusp3d |
3D graphical display of cusp model fit (experimental). |
dcusp |
Density of Cobb's cusp distribution |
pcusp |
Cumulative probability function of Cobb's cusp distribution |
qcusp |
Quantile function of Cobb's cusp distribution |
rcusp |
Sample from Cobb's cusp distribution. |
cusp.logist |
Fit logistic model for bifurcation testing (experimental) |
Raoul Grasman <[email protected]>
L. Cobb and S. Zacks (1985) Applications of Catastrophe Theory for Statistical Modeling in the Biosciences (article), Journal of the American Statistical Association, 392:793–802.
P. Hartelman (1996). Stochastic Catastrophe Theory. Unpublished PhD-thesis.
H. L. J. van der Maas, R. Kolstein, and J van der Pligt (2003). Sudden Transitions in Attitudes, Sociological Methods and Research, 32:125-152.
Oliva, DeSarbo, Day, and Jedidi. (1987) GEMCAT : A General Multivariate Methodology for Estimating Catastrophe Models, Behavioral Science, 32:121-137.
R. P. P. P. Grasman, H. L. J. van der Maas, and E-J. Wagenmakers (2009). Fitting the Cusp Catastrophe in R: A cusp Package Primer. Journal of Statistical Software 32(8), 1-28. URL https://www.jstatsoft.org/v32/i08/.
set.seed(123) # fitting cusp to cusp data x <- rcusp(100, alpha=0, beta=1) fit <- cusp(y ~ x, alpha ~ 1, beta ~ 1) print(fit) # example with regressors ## Not run: x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit) plot(fit) cusp3d(fit) ## End(Not run) # use of OK npar <- length(fit$par) ## Not run: while(!fit$OK) # refit if necessary until convergence is OK fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data, start=rnorm(npar)) ## End(Not run) ## Not run: # example 1 from paper data(attitudes) data(attitudeStartingValues) fit.attitudes <- cusp(y ~ Attitude, alpha ~ Orient + Involv, beta ~ Involv, data = attitudes, start=attitudeStartingValues) summary(fit.attitudes) plot(fit.attitudes) cusp3d(fit.attitudes, B = 0.75, Y = 1.35, theta = 170, phi = 30, Yfloor = -9) ## End(Not run)
set.seed(123) # fitting cusp to cusp data x <- rcusp(100, alpha=0, beta=1) fit <- cusp(y ~ x, alpha ~ 1, beta ~ 1) print(fit) # example with regressors ## Not run: x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit) plot(fit) cusp3d(fit) ## End(Not run) # use of OK npar <- length(fit$par) ## Not run: while(!fit$OK) # refit if necessary until convergence is OK fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data, start=rnorm(npar)) ## End(Not run) ## Not run: # example 1 from paper data(attitudes) data(attitudeStartingValues) fit.attitudes <- cusp(y ~ Attitude, alpha ~ Orient + Involv, beta ~ Involv, data = attitudes, start=attitudeStartingValues) summary(fit.attitudes) plot(fit.attitudes) cusp3d(fit.attitudes, B = 0.75, Y = 1.35, theta = 170, phi = 30, Yfloor = -9) ## End(Not run)
Data set reflecting bistability in political attitudes
data(attitudes) data(attitudeStartingValues)
data(attitudes) data(attitudeStartingValues)
A data frame with 1387 observations on the following 3 variables.
Orient
a numeric vector
Involv
a numeric vector
Attitude
a numeric vector
The format of attitudeStartingValues is: num [1:7] 0.153 -0.453 -0.097 -0.124 -0.227 ...
The data set was taken from (van der Maas, Kolstein, & van der Pligt, 2003). It 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”. Responses of some 1387 Dutch respondents are included who indicated their level of agreement with this statement on a 5 point scale (1 = total ly agree, 5 = total ly 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 starting values provided here for a cusp analysis of the attitude
data set give proper convergence in one run. They were found after many trial starting values that yielded improper convergence.
van der Maas HLJ, Kolstein R, van der Pligt J (2003). Sudden Transitions in Attitudes. Sociological Methods & Research, 23(2), 125152.
van der Maas HLJ, Kolstein R, van der Pligt J (2003). Sudden Transitions in Attitudes. Sociological Methods & Research, 23(2), 125152.
data(attitudes) data(attitudeStartingValues) ## Not run: fit <- cusp(y ~ Attitude, alpha ~ Orient + Involv, beta ~ Involv, data = attitudes, start=attitudeStartingValues) ## End(Not run) ## maybe str(attitudeStartingValues) ; plot(attitudeStartingValues) ...
data(attitudes) data(attitudeStartingValues) ## Not run: fit <- cusp(y ~ Attitude, alpha ~ Orient + Involv, beta ~ Involv, data = attitudes, start=attitudeStartingValues) ## End(Not run) ## maybe str(attitudeStartingValues) ; plot(attitudeStartingValues) ...
This function fits a cusp catastrophe model to data using the maximum likelihood method of Cobb. Both the state variable may be modeled by a linear combination of variables and design factors, as well as the normal/asymmetry factor alpha
and bifurction/splitting factor beta
.
cusp(formula, alpha, beta, data, weights, offset, ..., control = glm.control(), method = "cusp.fit", optim.method = "L-BFGS-B", model = TRUE, contrasts = NULL)
cusp(formula, alpha, beta, data, weights, offset, ..., control = glm.control(), method = "cusp.fit", optim.method = "L-BFGS-B", model = TRUE, contrasts = NULL)
formula |
|
alpha |
|
beta |
|
data |
|
weights |
vector of weights by which each data point is weighted (experimental) |
offset |
vector of offsets for the data (experimental) |
... |
named arguments that are passed to |
control |
|
method |
string, currently unused |
optim.method |
string passed to |
model |
should the model matrix be returned? |
contrasts |
matrix of |
cusp
fits a cusp catastrophe model to data. Cobb's definition for the canonical form of the stochastic cusp catastrophe is the stochastic differential equation
The stationary distribution of the ‘behavioral’, or ‘state’ variable , given the control parameters
(‘asymmetry’ or ‘normal’ factor) and
(‘bifurcation’ or ‘splitting’ factor) is
where is a normalizing constant.
The behavioral variable and the asymmetry and bifurcation factors are usually not directly related to the dependent and independent variables in the data set. These are therefore used to predict the state variable and control parameters:
in which the 's,
's, and
's are estimated by means of maximum likelihood.
Here, the
's and
's are variables constructed from variables in the data set. Variables predicting the
's and
's need not be the same.
The state variable and control parameters can be modeled by specifying a model formula
:
in which model
can be any valid formula
specified in terms of variables that are present in the data.frame
.
List with components
coefficients |
Estimated coefficients |
rank |
rank of Hessian matrix |
qr |
|
linear.predictors |
two column matrix containing the |
deviance |
sum of squared errors using Delay convention |
aic |
AIC |
null.deviance |
variance of canonical state variable |
iter |
number of optimization iterations |
weights |
weights provided through weights argument |
df.residual |
residual degrees of freedom |
df.null |
degrees of freedom of constant model for state variable |
y |
predicted values of state variable |
converged |
convergence status given by |
par |
parameter estimates for |
Hessian |
Hessian matrix of negative log likelihood function at minimum |
hessian.untransformed |
Hessian matrix of negative log likelihood for |
code |
|
model |
list with model design matrices |
call |
function call that created the object |
formula |
list with the formulas |
OK |
logical. |
data |
original data.frame |
Raoul Grasman
See cusp-package
summary.cusp
for summaries and model assessment.
The generic functions coef
, effects
, residuals
, fitted
, vcov
.
predict
for estimated values of the control parameters and
,
set.seed(123) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit) ## Not run: plot(fit) cusp3d(fit) ## End(Not run) # useful use of OK ## Not run: while(!fit$OK) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data, start=rnorm(fit$par)) # use different starting values ## End(Not run)
set.seed(123) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit) ## Not run: plot(fit) cusp3d(fit) ## End(Not run) # useful use of OK ## Not run: while(!fit$OK) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data, start=rnorm(fit$par)) # use different starting values ## End(Not run)
Given bifurcation/splitting factor values this function computes the border values of the normal/symmetry factor for the bifurcation set of the cusp catastrophe.
cusp.bifset(beta)
cusp.bifset(beta)
beta |
values of the bifurcation/splitting factor at which the border values of the normal/symmetry factor is computed |
Matrix with columns named beta
, alpha.l
, alpha.u
. The latter two columns give respectively the lower and upper border values of the normal/symmetry factor. Negative values of beta
give NaN
values for the normal factor.
Raoul Grasman
See cusp-package
cusp.bifset(-3:3)
cusp.bifset(-3:3)
This function computes the locations of the extrema of the cusp catastrophe potential function.
cusp.extrema(alpha, beta)
cusp.extrema(alpha, beta)
alpha |
(single) value of normal/symmetry factor |
beta |
(single) value of bifurcation/splitting factor |
The locations are determined by computing the solutions to the equation
Ordered vector with locations of extremes.
Use Vectorize
to allow for array input.
Raoul Grasman
http://www.scholarpedia.org/article/Cusp_bifurcation
# simple use cusp.extrema(2,3) # using vectorize to allow for array input; # returns a matrix with locations in each column Vectorize(cusp.extrema)(-3:3, 2)
# simple use cusp.extrema(2,3) # using vectorize to allow for array input; # returns a matrix with locations in each column Vectorize(cusp.extrema)(-3:3, 2)
This function fits a logistic curve model to data using maximum likelihood under the assumption of normal errors (i.e., nonlinear least squares). Both the response variable may be modeled by a linear combination of variables and design factors, as well as the normal/asymmetry factor alpha
and bifurcation/splitting factor beta
.
cusp.logist(formula, alpha, beta, data, ..., model = TRUE, x = FALSE, y = TRUE)
cusp.logist(formula, alpha, beta, data, ..., model = TRUE, x = FALSE, y = TRUE)
formula , alpha , beta
|
|
data |
|
... |
named arguments that are passed to |
model , x , y
|
logicals. If |
A nonlinear regression is carried out of the model
for ,
where
in which the 's, and
's, are estimated. The
's are variables in the data set
and specified by
formula
; the 's are variables in the data set and are specified in
alpha
and beta
. Variables in alpha
and beta
need not be the same. The 's are estimated implicitly
using concentrated likelihood methods, and are not returned explicitly.
List with components
minimum |
Objective function value at minimum |
estimate |
Coordinates of objective function minimum |
gradient |
Gradient of objective function at minimum. |
code |
Convergence |
iterations |
Number of iterations used by |
coefficients |
A named vector of estimates of |
linear.predictors |
Estimates of |
fitted.values |
Predicted values of |
residuals |
Residuals |
rank |
Numerical rank of matrix of predictors for |
deviance |
Residual sum of squares. |
logLik |
Log of the likelihood at the minimum. |
aic |
Akaike's information criterion |
rsq |
R Squared (proportion of explained variance) |
df.residual |
Degrees of freedom for the residual |
df.null |
Degrees of freedom for the Null residual |
rss |
Residual sum of squares |
hessian |
Hessian matrix of objective function at the minimum if |
Hessian |
Hessian matrix of log-likelihood function at the minimum (currently unavailable) |
qr |
QR decomposition of the |
converged |
Boolean indicating if optimization convergence is proper (based on exit code |
weights |
|
call |
the matched call |
y |
If requested (the default), the matrix of response variables used. |
x |
If requested, the model matrix used. |
null.deviance |
The sum of squared deviations from the mean of the estimated |
Raoul Grasman
Hartelman PAI (1997). Stochastic Catastrophe Theory. Amsterdam: University of Amsterdam, PhD thesis.
A family of functions that return the normalization constant for the cusp density given the values of the bifurcation and asymmetry parameters (default), or returns the moment of a specified order (cusp.nc
).
cusp.nc(alpha, beta, mom.order = 0, ...) cusp.nc.c(alpha, beta, ..., keep.order = TRUE) cusp.nc.C(alpha, beta, subdivisions = 100, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, aux = NULL, keep.order = TRUE) cusp.nc.vec(alpha, beta, ..., keep.order = FALSE)
cusp.nc(alpha, beta, mom.order = 0, ...) cusp.nc.c(alpha, beta, ..., keep.order = TRUE) cusp.nc.C(alpha, beta, subdivisions = 100, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, aux = NULL, keep.order = TRUE) cusp.nc.vec(alpha, beta, ..., keep.order = FALSE)
alpha |
the asymmetry parameter in Cobb's cusp density (see |
beta |
the bifurcation parameter in Cobb's cusp density (see |
mom.order |
the moment order to be computed (see details below) |
subdivisions , rel.tol , abs.tol , stop.on.error , aux
|
Arguments used by the internal integration routine of R (see |
keep.order |
Logical, that indicates whether the order of the output should be the same as the order of the input. |
... |
Extra arguments in |
The function cusp.nc
returns if
mom.order = 0
and times the moment of order
mom.order
otherwise.
The function cusp.nc
is internally used if the C-routine symbol "cuspnc"
is not loaded.
The functions cusp.nc.c
and cusp.nc.C
call this C routine, which is considerably faster than
cusp.nc
.
These functions are not intended to be called directly by the user.
cusp.nc, cusp.nc.c, cusp.nc.vec
return a numeric vector of the same length as alpha
and beta
with normalizing constants, or the indicated moments times the normalization constant (cusp.nc
only).
cusp.nc.C
returns a list with vectors with the results obtained from integrate
.
cusp.nc.c
first sorts the input in such a way that the numerical integrals can be evaluated more quickly than in arbitrary order
Raoul Grasman
(Negative) log-likelihood for Cobb's cusp probability density function used by cusp
. This function is not to be called by the user. See help(cusp)
.
cusp.nlogLike(p, y, X.alpha, X.beta = X.alpha, ..., verbose = FALSE) cusp.nlogLike.c(p, y, X.alpha, X.beta = X.alpha, ..., verbose = FALSE) cusp.logLike(p, x, verbose = FALSE)
cusp.nlogLike(p, y, X.alpha, X.beta = X.alpha, ..., verbose = FALSE) cusp.nlogLike.c(p, y, X.alpha, X.beta = X.alpha, ..., verbose = FALSE) cusp.logLike(p, x, verbose = FALSE)
p |
parameter vector |
x |
vector of observed values for the state variable in the cusp ( |
y |
design matrix predicting state values at which the likelihood is evaluated |
X.alpha |
design matrix predicting alpha in the model |
X.beta |
design matrix predicting beta in the model |
... |
unused extra arguments |
verbose |
logical, if |
cusp.nlogLike
is the R version of the corresponding C function wrapped by cusp.nlogLike.c
These functions are not intended to be called directly by the user.
The value of the negative log-likelihood function (cusp.nlogLike
, cusp.nlogLike.c
),
the value of the log-likelihood function (cusp.logLike
).
The functions are not to be called by the user directly.
Raoul Grasman
See cusp-package
This function generates a 3D display of the fit (object) of a cusp model.
cusp3d(y, alpha = if (!missing(y) && is.list(y)) y$lin[, "alpha"], beta = if (!missing(y) && is.list(y)) y$lin[, "beta"], w = 0.03, theta = 170, phi = 35, B = 4, Y = 3, Yfloor = -15, np = 180, n.surface = 30, surface.plot = TRUE, surf.alpha = 0.75, surf.gamma = 1.5, surf.chroma = 35, surf.hue = 240, surf.ltheta = 0, surf.lphi = 45, ...)
cusp3d(y, alpha = if (!missing(y) && is.list(y)) y$lin[, "alpha"], beta = if (!missing(y) && is.list(y)) y$lin[, "beta"], w = 0.03, theta = 170, phi = 35, B = 4, Y = 3, Yfloor = -15, np = 180, n.surface = 30, surface.plot = TRUE, surf.alpha = 0.75, surf.gamma = 1.5, surf.chroma = 35, surf.hue = 240, surf.ltheta = 0, surf.lphi = 45, ...)
y |
object returned by |
alpha |
vector of normal/symmetry factor values corresponding to the state values in |
beta |
vector of bifurcation/splitting factor values corresponding to the state values in |
w |
number that specifies the size of the data points plotted on the cusp surface |
theta , phi
|
Angles defining the viewing direction: |
B |
range of the splitting factor axis |
Y |
range of the state variable axis |
Yfloor |
location on state variable axis where the control surface is plotted |
np |
factor that determines the fineness of the drawing |
n.surface |
factor that determines the fineness of the rendered surface |
surface.plot |
plot the surface? |
surf.alpha |
transparency level of rendered surface |
surf.gamma |
factor that determines the shading of surface facets ( |
surf.chroma , surf.hue
|
chroma and hue of surface color (see |
surf.ltheta , surf.lphi
|
the surface is shaded as though it was being illuminated from the direction specified by azimuth |
... |
named parameters that are passed to |
This function is experimental.
cusp3d
returns the viewing transformation matrix, say VT
, a 4 x 4 matrix suitable for projecting 3D coordinates (x,y,z) into the 2D plane using homogeneous 4D coordinates (x,y,z,t). It can be used to superimpose additional graphical elements on the 3D plot, by lines() or points(), using the simple function trans3d().
Currently still somewhat buggy.
Raoul Grasman
See cusp-package
persp
, plot.cusp
, cusp3d.surface
set.seed(123) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) cusp3d(fit)
set.seed(123) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) cusp3d(fit)
This function generates a 3D display of the cusp equilibrium surface.
cusp3d.surface(alpha = c(-5, 5), beta = c(-3, 3), y = 41, xlim = range(alpha), ylim = range(beta), zlim = c(-5, 4), xlab = expression(alpha), ylab = expression(beta), zlab = "equilibrium states", main = NULL, sub = NULL, phi = 20, theta = 160, r = sqrt(3), d = 1, scale = TRUE, expand = 1, hue = 240, chroma = 35, surf.alpha = 0.75, gamma = 1.5, bcol = NA, lcol = "gray", ltheta = 90, lphi = 70, box = TRUE, axes = FALSE, nticks = 5, ticktype = "simple", floor.lines = TRUE, ...)
cusp3d.surface(alpha = c(-5, 5), beta = c(-3, 3), y = 41, xlim = range(alpha), ylim = range(beta), zlim = c(-5, 4), xlab = expression(alpha), ylab = expression(beta), zlab = "equilibrium states", main = NULL, sub = NULL, phi = 20, theta = 160, r = sqrt(3), d = 1, scale = TRUE, expand = 1, hue = 240, chroma = 35, surf.alpha = 0.75, gamma = 1.5, bcol = NA, lcol = "gray", ltheta = 90, lphi = 70, box = TRUE, axes = FALSE, nticks = 5, ticktype = "simple", floor.lines = TRUE, ...)
alpha |
numeric 2-vector specifying the normal/symmetry factor axis range |
beta |
numeric 2-vector specifying the bifurcation/splitting factor axis range |
y |
numeric specifying the iso contours used to render the surface (see details below) |
xlim , ylim , zlim
|
numeric 2-vectors (see |
xlab , ylab , zlab , main , sub
|
strings (see |
phi , theta
|
numeric, determine viewing direction (see |
r |
numeric, distance to center of the plotting box (see |
d |
numeric, strength of perspective transformation (see |
scale , expand
|
logical, see |
hue , chroma , surf.alpha
|
hue, chroma and alpha (transparency) of the surface segments (see |
gamma |
gamma for shading of surface (see |
bcol |
color, |
lcol |
color of the lines on the floor of the plotting cube |
ltheta , lphi
|
numeric, direction of illumination of the surface (similar to |
box , axes , nticks , ticktype
|
(see |
floor.lines |
logical, if |
... |
If y
has length 1, it is interpreted as the number of contours. Otherwise it is interpreted as a vector of contour levels from which the surface must be determined. If y
is a number, the exact range of y
is determined by the ranges of alpha
and beta
through the cusp equilibrium equation below.
The surface is constructed from the iso-contours of the cusp equilibrium surface that makes up the solutions to
as a (multi-)function of the asymmetry variable and bifurcation variable
. For each possible solution
the iso-contours are given by the equation
which are linear in . For each value of
the values of
are determined for the end points of the
range specified by
beta
. The two 3D coordinates (,
,
) are projected onto the 2D canvas using the
persp
transformation matrix and used for drawing the lines and polygons.
cusp3d.surface
returns the viewing transformation matrix, say VT
, a 4 x 4 matrix suitable for projecting 3D coordinates (x,y,z) into the 2D plane using homogeneous 4D coordinates (x,y,z,t). It can be used to superimpose additional graphical elements on the 3D plot, by lines() or points(), using the simple function trans3d().
This function is an alternative to cusp3d
which uses a different method of rendering and also plots fitted points on the surface.
Raoul Grasman
See cusp-package
, cusp3d
## Not run: p = cusp3d.surface(chroma=40,lcol=1,surf.alpha=.95,phi=30,theta=150, bcol="surface",axes=TRUE,main="Cusp Equilibrium Surface") lines(trans3d(c(5,5), c(3,3), c(-5,4), p), lty=3) # replot some of the box outlines lines(trans3d(c(-5,5), c(3,3), c(4,4), p), lty=3) ## End(Not run)
## Not run: p = cusp3d.surface(chroma=40,lcol=1,surf.alpha=.95,phi=30,theta=150, bcol="surface",axes=TRUE,main="Cusp Equilibrium Surface") lines(trans3d(c(5,5), c(3,3), c(-5,4), p), lty=3) # replot some of the box outlines lines(trans3d(c(-5,5), c(3,3), c(4,4), p), lty=3) ## End(Not run)
Functions for the cusp distribution.
dcusp(y, alpha, beta) pcusp(y, alpha, beta, subdivisions = 100, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, aux = NULL, keep.order = TRUE) qcusp(p, alpha, beta) rcusp(n, alpha, beta)
dcusp(y, alpha, beta) pcusp(y, alpha, beta, subdivisions = 100, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, aux = NULL, keep.order = TRUE) qcusp(p, alpha, beta) rcusp(n, alpha, beta)
y |
vector of quantiles |
p |
vector of probabilities |
n |
number of observations. |
alpha |
normal/asymmetry factor value of cusp density |
beta |
bifurcation/splitting factor value of cusp density |
subdivisions |
See |
rel.tol |
See |
abs.tol |
See |
stop.on.error |
See |
aux |
See |
keep.order |
logical. If true the order of the output values is the same as those of the input values |
The cusp distribution is defined by
where is the normalizing constant.
rcusp
uses rejection sampling to generate samples.
qcusp
implements binary search and is rather slow.
dcusp
gives the density function, pcusp
gives the distribution function, qcusp
gives the quantile function, and rcusp
generates observations.
Raoul Grasman
See cusp-package
, integrate
# evaluate density and distribution dcusp(0,2,3) pcusp(0,2,3) pcusp(qcusp(0.125,2,3),2,3) # = 0.125 # generate cusp variates rcusp(100, 2, 3) # generate cusp variates for random normal and splitting factor values alpha = runif(20, -3, 3) beta = runif(20, -3, 3) Vectorize(rcusp)(1, alpha, beta)
# evaluate density and distribution dcusp(0,2,3) pcusp(0,2,3) pcusp(qcusp(0.125,2,3),2,3) # = 0.125 # generate cusp variates rcusp(100, 2, 3) # generate cusp variates for random normal and splitting factor values alpha = runif(20, -3, 3) beta = runif(20, -3, 3) Vectorize(rcusp)(1, alpha, beta)
Add a miniature bifurcation set for the cusp catastrophe to an existing plot.
draw.cusp.bifset(rx = par("usr")[1:2], ry = par("usr")[3:4], xpos = min(rx) + 0.01 * diff(rx)[1], ypos = max(ry) - 0.01 * diff(ry)[1], xscale = 0.1 * diff(rx), yscale = 0.1 * diff(ry) / xscale, aspect = 1, mark = 1, col = hsv(0.7, s = 0.8, alpha = 0.5), border = NA, density = NA, bifurcation.set.fill = gray(0.8), background = hsv(0.1, s = 0.1, alpha = 0.5), ..., X)
draw.cusp.bifset(rx = par("usr")[1:2], ry = par("usr")[3:4], xpos = min(rx) + 0.01 * diff(rx)[1], ypos = max(ry) - 0.01 * diff(ry)[1], xscale = 0.1 * diff(rx), yscale = 0.1 * diff(ry) / xscale, aspect = 1, mark = 1, col = hsv(0.7, s = 0.8, alpha = 0.5), border = NA, density = NA, bifurcation.set.fill = gray(0.8), background = hsv(0.1, s = 0.1, alpha = 0.5), ..., X)
rx |
x-axis range of the plot window |
ry |
y-axis range of the plot window |
xpos |
x-axis position of drawing |
ypos |
y-axis position of drawing |
xscale |
scaling applied to drawing along x-axis |
yscale |
scaling applied to drawing along y-axis |
aspect |
aspect ratio |
mark |
0, 1, 2, 3, or 4; indicates which part of the cusp surface should be marked |
col |
color used for marking a part of the cusp surface |
border |
color used for the marked part of the cusp surface. See |
density |
the density of shading lines of the marked part of the cusp surface, in lines per inch. The default value of NULL means that no shading lines are drawn. See |
bifurcation.set.fill |
color for marking the bifurcation set |
background |
background color of the cusp surface |
... |
arguments passed to |
X |
|
This function is mainly intended for internal use by cusp.plot
.
No return value. Called for its side effect.
Raoul Grasman
http://www.scholarpedia.org/article/Cusp_bifurcation
## Not run: plot(1:10) draw.cusp.bifset(mark=0) # no marking ## End(Not run)
## Not run: plot(1:10) draw.cusp.bifset(mark=0) # no marking ## End(Not run)
Synthetic ‘multivariate’ data from the cusp catastrophe as generated from the equations specified by Oliva et al. (1987).
data(oliva)
data(oliva)
A data frame with 50 observations on the following 12 variables.
x1
splitting factor predictor
x2
splitting factor predictor
x3
splitting factor predictor
y1
the bifurcation factor predictor
y2
the bifurcation factor predictor
y3
the bifurcation factor predictor
y4
the bifurcation factor predictor
z1
the state factor predictor
z2
the state factor predictor
alpha
the true 's
beta
the true 's
y
the true state variable values
The data in Oliva et al. (1987) are obtained from the equations
Here the 's are uniformly distributed on (-2,2), and the
's and
are
uniform on (-3,3).
The states
were then generated from the cusp density, using
rcusp
, with their respective
's and
's as normal and splitting factors, and then
was computed as
Oliva T, DeSarbo W, Day D, Jedidi K (1987). GEMCAT: A general multivariate methodology for estimating catastrophe models. Behavioral Science, 32(2), 121137.
Oliva T, Desarbo W, Day D, Jedidi K (1987). GEMCAT: A general multivariate methodology for estimating catastrophe models. Behavioral Science, 32(2), 121137.
data(oliva) set.seed(121) fit <- cusp(y ~ z1 + z2 - 1, alpha ~ x1 + x2 + x3 - 1, ~ y1 + y2 + y3 + y4 - 1, data = oliva, start = rnorm(9)) summary(fit) ## Not run: cusp3d(fit, B=5.25, n.surf=50, theta=150) # B modifies the range of beta (is set here to 5.25 to make # sure all points lie on the surface) ## End(Not run)
data(oliva) set.seed(121) fit <- cusp(y ~ z1 + z2 - 1, alpha ~ x1 + x2 + x3 - 1, ~ y1 + y2 + y3 + y4 - 1, data = oliva, start = rnorm(9)) summary(fit) ## Not run: cusp3d(fit, B=5.25, n.surf=50, theta=150) # B modifies the range of beta (is set here to 5.25 to make # sure all points lie on the surface) ## End(Not run)
This function generates diagnostic graphical displays of fits of a cusp catastrophe model to data obtained with cusp
## S3 method for class 'cusp' plot(x, what = c("all", "bifurcation", "residual", "densities"), ...)
## S3 method for class 'cusp' plot(x, what = c("all", "bifurcation", "residual", "densities"), ...)
x |
Object returned by |
what |
1-character string giving the type of plot desired. The following values are possible: |
... |
named arguments that are passed to lower level plotting function |
These diagnostic plots help to identify problems with the fitted model. In optimal cases the fitted locations in the parameter plane are dispersed over regions of qualitatively different behavior. Within each region the fitted dependent values have a density of the appropriate shape (e.g., bimodal in the bifurcation set).
No return value. Called for its side effect.
Raoul Grasman
See cusp-package
plotCuspBifurcation
, plotCuspResidfitted
, plotCuspDensities
set.seed(20) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) ## Not run: plot(fit) # just densities layout(matrix(1:4,2)) plot(fit, what="densities") ## End(Not run)
set.seed(20) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) ## Not run: plot(fit) # just densities layout(matrix(1:4,2)) plot(fit, what="densities") ## End(Not run)
Displays fitted data points on the control plane of cusp catastrophe. The function takes a fit object obtained with cusp
and generates a plot. Different diagnostic plots may be chosen, or all can be combined in a single plot (the default).
plotCuspBifurcation(object, xlim = a + c(-0.3, 0.3), ylim = b + c(-0.1, 0.1), xlab = expression(alpha), ylab = expression(beta), hue = 0.5 + 0.25 * tanh(object$y), col = hsv(h = hue, s = 1, alpha = 0.4), cex.xlab = 1.55, cex.ylab = cex.xlab, axes = TRUE, box = TRUE, add = FALSE, bifurcation.set.fill = gray(0.8), cex.scale = 15, cex = (cex.scale/log(NROW(ab))) * dens/max(dens), pch = 20)
plotCuspBifurcation(object, xlim = a + c(-0.3, 0.3), ylim = b + c(-0.1, 0.1), xlab = expression(alpha), ylab = expression(beta), hue = 0.5 + 0.25 * tanh(object$y), col = hsv(h = hue, s = 1, alpha = 0.4), cex.xlab = 1.55, cex.ylab = cex.xlab, axes = TRUE, box = TRUE, add = FALSE, bifurcation.set.fill = gray(0.8), cex.scale = 15, cex = (cex.scale/log(NROW(ab))) * dens/max(dens), pch = 20)
object |
object returned by |
xlim |
the x limits (x1, x2) of the plot. |
ylim |
the y limits of the plot. |
xlab |
a label for the x axis. |
ylab |
a label for the x axis. |
hue |
hue of points (see |
col |
color used in plots |
cex.xlab , cex.ylab
|
see |
axes |
logical. Should the axes be displayed? |
box |
logical. Should a box be drawn around the plot? |
add |
logical. Add to current plot? |
bifurcation.set.fill |
1-character string. Color used to fill the bifurcation set (see |
cex.scale , cex , pch
|
see |
The default hue of each dot is a function of the height of the cusp surface to which it is closest. This is especially useful in the bifurcation set. Purple dots are higher than green dots.
The size of the dots depends on the density of dots at its location. The higher the density the larger the dot.
No return value. Called for its side effect.
Raoul Grasman
See cusp-package
set.seed(20) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) ## Not run: plot(fit, what='bifurcation', box=TRUE, axes=FALSE) ## End(Not run)
set.seed(20) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) ## Not run: plot(fit, what='bifurcation', box=TRUE, axes=FALSE) ## End(Not run)
Plot density of state variables conditioned on their location on the cusp control surface.
plotCuspDensities(object, main = "Conditional density", ...)
plotCuspDensities(object, main = "Conditional density", ...)
object |
cusp fit object returned by |
main |
title of plot |
... |
named arguments that are passed to |
This function is mainly intended for internal use by plot.cusp
.
No return value. Called for its side effect.
Raoul Grasman
Plot Residuals against Fitted Values for a Cusp Model Fit.
plotCuspResidfitted(object, caption = "Residual vs Fitted", xlab = paste("Fitted (", colnames(fitted(object))[1], " convention)", sep = ""), ylab = "Residual", ...)
plotCuspResidfitted(object, caption = "Residual vs Fitted", xlab = paste("Fitted (", colnames(fitted(object))[1], " convention)", sep = ""), ylab = "Residual", ...)
object |
cusp fit object returned by |
caption |
plot caption |
xlab |
label for x-axis |
ylab |
label for y-axis |
... |
named arguments that are passed to |
This function is mainly intended for internal use by plot.cusp
.
No return value. Called for its side effect.
Raoul Grasman
Predicted values based on a cusp model object.
## S3 method for class 'cusp' predict(object, newdata, se.fit = FALSE, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, method = c("delay", "maxwell", "expected"), keep.linear.predictors = FALSE, ...)
## S3 method for class 'cusp' predict(object, newdata, se.fit = FALSE, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, method = c("delay", "maxwell", "expected"), keep.linear.predictors = FALSE, ...)
object |
Object of class " |
newdata |
An optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used. |
se.fit |
See |
interval |
See |
level |
See |
type |
See |
terms |
See |
na.action |
See |
pred.var |
See |
weights |
See |
method |
Type of prediction convention to use. Can be abbreviated. ( |
keep.linear.predictors |
Logical. Should the linear predictors (alpha, beta, and y) be returned? |
... |
further arguments passed to or from other methods. |
predict.cusp
produces predicted values, obtained by evaluating the regression functions from the
cusp object
in the frame newdata
using predict.lm
. This results in linear
predictors for the cusp control variables alpha
, and beta
, and, if method = "delay"
,
for the behavioral cusp variable y
. These are then used to compute predicted values: If
method = "delay"
these are the points on the cusp surface defined by
that are closest to y
. If method = "maxwell"
they are
the points on the cusp surface corresponding to the minimum of the associated potential function
.
A vector of predictions. If keep.linear.predictors
the return value has a "data"
attribute
which links to newdata
augmented with the linear predictors alpha
, beta
, and, if
method = "delay"
, y
. If method = "expected"
, the expected value from the equilibrium distribution of the stochastic process
where is
a Wiener proces (aka Brownian motion) is returned. (This distribution is implemented in
dcusp
.)
Currently method = "expected"
should not be trusted.
Raoul Grasman
See cusp-package
.
set.seed(123) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) newdata = data.frame(x1 = runif(10), x2 = runif(10), z = 0) predict(fit, newdata)
set.seed(123) # example with regressors x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) newdata = data.frame(x1 = runif(10), x2 = runif(10), z = 0) predict(fit, newdata)
summary
method for class “cusp”
## S3 method for class 'cusp' summary(object, correlation = FALSE, symbolic.cor = FALSE, logist = FALSE, ...) ## S3 method for class 'summary.cusp' print(x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"), ...)
## S3 method for class 'cusp' summary(object, correlation = FALSE, symbolic.cor = FALSE, logist = FALSE, ...) ## S3 method for class 'summary.cusp' print(x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"), ...)
object |
Object returned by |
x |
‘ |
correlation |
logical; if |
symbolic.cor |
logical; currently unused |
logist |
logical. If |
digits |
numeric; the number of significant digits to use when printing. |
signif.stars |
logical. If |
... |
further arguments passed to or from other methods. |
print.summary.cusp
tries to be smart about formatting the coefficients, standard errors, etc. and additionally gives significance stars if signif.stars
is TRUE
.
Correlations are printed to two decimal places (or symbolically): to see the actual correlations print summary(object)$correlation
directly.
The function summary.cusp
computes and returns a list of summary statistics of the fitted linear model given in object, using the components (list elements) “call
” and “terms
” from its argument, plus
call |
the matched call |
terms |
the |
deviance |
sum of squared residuals of cusp model fit |
aic |
Akaike Information Criterion for cusp model fit |
contrasts |
contrasts used |
df.residual |
degrees of freedom for the residuals of the cusp model fit |
null.deviance |
variance of canonical state variable |
df.null |
degrees of freedom of constant model for state variable |
iter |
number of optimization iterations |
deviance.resid |
residuals computed by |
coefficients |
a |
aliased |
named logical vector showing if the original coefficients are aliased. |
dispersion |
always 1 |
df |
3-vector containing the rank of the model matrix, residual degrees of freedom, and model degrees of freedom. |
resid.name |
string specifying the convention used in determining the residuals (i.e., "Delay" or "Maxwell"). |
cov.unscaled |
the unscaled (dispersion = 1) estimated covariance matrix of the estimated coefficients. |
r2lin.r.squared |
where |
r2lin.dev |
residual sums of squares of the linear model |
r2lin.df |
degrees of freedom for the linear model |
r2lin.logLik |
value of the log-likelihood for the linear model assuming normal errors |
r2lin.npar |
number of parameters in the linear model |
r2lin.aic |
AIC for the linear model |
r2lin.aicc |
corrected AIC for the linear model |
r2lin.bic |
BIC for the linear model |
r2log.r.squared |
|
r2log.dev |
if |
r2log.df |
ditto, degrees of freedom for the logistic model |
r2log.logLik |
ditto, value of log-likelihood function for the logistic model assuming normal errors. |
r2log.npar |
ditto, number of parameters for the logistic model |
r2log.aic |
ditto, AIC for logistic model |
r2log.aicc |
ditto, corrected AIC for logistic model |
r2log.bic |
ditto, BIC for logistic model |
r2cusp.r.squared |
pseudo-
This value can be negative. |
r2cusp.dev |
residual sums of squares for cusp model |
r2cusp.df |
residual degrees of freedom for cusp model |
r2cusp.logLik |
value of the log-likelihood function for the cusp model |
r2cusp.npar |
number of parameters in the cusp model |
r2cusp.aic |
AIC for cusp model fit |
r2cusp.aicc |
corrected AIC for cusp model fit |
r2cusp.bic |
BIC for cusp model fit. |
Raoul Grasman
Cobb L, Zacks S (1985). Applications of Catastrophe Theory for Statistical Modeling in the Biosciences. Journal of the American Statistical Association, 80(392), 793–802.
Hartelman PAI (1997). Stochastic Catastrophe Theory. Amsterdam: University of Amsterdam, PhD thesis.
Cobb L (1998). An Introduction to Cusp Surface Analysis.
https://www.aetheling.com/models/cusp/Intro.htm.
set.seed(97) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit, logist=FALSE) # set logist to TRUE to compare to logistic fit
set.seed(97) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) print(fit) summary(fit, logist=FALSE) # set logist to TRUE to compare to logistic fit
Returns an estimate of the variance-covariance matrix of the main parameters of a fitted cusp model object.
## S3 method for class 'cusp' vcov(object, ...) ## S3 method for class 'cusp' confint(object, parm, level = 0.95, ...)
## S3 method for class 'cusp' vcov(object, ...) ## S3 method for class 'cusp' confint(object, parm, level = 0.95, ...)
object |
a fitted cusp model object. |
parm |
a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
the confidence level required. |
... |
additional arguments for method functions. |
The variance-covariance matrix is estimated by the inverse of the Hessian matrix of the log-likelihood at the maximum likelihood estimate (vcov
).
Normal theory confidence intervals are computed for all parameters in the cusp model object using vcov
to obtain the standard errors (confint
).
The variance-covariance matrix (vcov
).
A matrix (or vector) with columns giving lower and upper confidence limits for each parameter. These will be labeled as (1-level)/2 and 1 - (1-level)/2 in
Raoul Grasman
Seber, Wild (2005) Nonlinear regression. New York: Wiley
set.seed(123) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) vcov(fit)
set.seed(123) x1 = runif(150) x2 = runif(150) z = Vectorize(rcusp)(1, 4*x1-2, 4*x2-1) data <- data.frame(x1, x2, z) fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data) vcov(fit)
Data sets with measurements from different physical instances of Zeeman's Catastrophe Machine
data(zeeman1) data(zeeman2) data(zeeman3)
data(zeeman1) data(zeeman2) data(zeeman3)
A data frame with 150/198/282 observations on the following 3 variables.
x
A control plane variable that is manipulable by the experimentalist.
y
A control plane variable that is manipulable by the experimentalist.
z
The state variable of the machine: the shortest distance to the longitudinal axis of the machine.
The behavior Zeeman's catastrophe machine is archetypal for the Cusp catastrophe. This device consists of a wheel is tethered by an elastic chord to a fixed point. Another elastic, also attached to the wheel is moved about in the ‘control plane’ area opposite to the fixed point. The shortest distance between the strap point on the wheel and the axis defined by the fixed point and the control plane is recorded as a function of the position in the control plane. (In the original machine the angle between this axis and the line through the wheel center and the strap point is used.) See https://www.math.stonybrook.edu/~tony/whatsnew/column/catastrophe-0600/cusp4.html for a vivid demonstration. These data sets were obtained from 3 different physical instances of this machine, made by different people.
Measurements were made by systematically sampling different points in the control plane.
See vignette for example analysis with all three data sets.
For pictures of the machines, see
zeeman1
is due to Noemi Schuurman
zeeman2
is due to Karin Visser
zeeman3
is due to Mats Nagel & Joris ?
See https://sites.google.com/site/zeemanmachine/data-repository
Zeeman (1976).
data(zeeman1) data(zeeman2) data(zeeman3) ## Not run: fit <- cusp(y~z, alpha~x+y, beta~x+y, data=zeeman1) plot(fit) cusp3d(fit, surf.hue = 40, theta=215, phi=37.5, B=5.25) ## End(Not run)
data(zeeman1) data(zeeman2) data(zeeman3) ## Not run: fit <- cusp(y~z, alpha~x+y, beta~x+y, data=zeeman1) plot(fit) cusp3d(fit, surf.hue = 40, theta=215, phi=37.5, B=5.25) ## End(Not run)