This function fits a cusp catatrophe model to data using the maximum likelihood method of Cobb. Both the state variable may be modelled 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)
List with components
Estimated coefficients
rank of Hessian matrix
qr decomposition of the Hessian matrix
two column matrix containing the \(\alpha_i\)'s and \(\beta_i\)'s for each case
sum of squared errors using Delay convention
AIC
variance of canonical state variable
number of optimization iterations
weights provided through weights argument
residual degrees of freedom
degrees of freedom of constant model for state variable
predicted values of state variable
convergence status given by optim
parameter estimates for qr standardized data
Hessian matrix of negative log likelihood function at minimum
Hessian matrix of negative log likelihood for qr standardized data
optim convergence indicator
list with model design matrices
function call that created the object
list with the formulas
logical. TRUE if Hessian matrix is positive definite at the minimum obtained
original data.frame
formula that models the canonical state variable
formula that models the canonical normal/asymmetry factor
formula that models the canonical bifurcation/splitting factor
data.frame that contains all the variables named in the formulas
vector of weights by which each data point is weighted (experimental)
vector of offsets for the data (experimental)
named arguments that are passed to optim
glm.control object, currently unused
string, currently unused
string passed to optim to choose the optimization algorithm
should the model matrix be returned?
matrix of contrasts, experimental
Raoul Grasman
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 $$dY_t = (\alpha + \beta Y_t - Y_t^3)dt + dW_t.$$ The stationary distribution of the ‘behavioral’, or ‘state’ variable \(Y\), given the control parameters \(\alpha\) (‘asymmetry’ or ‘normal’ factor) and \(\beta\) (‘bifurcation’ or ‘splitting’ factor) is $$ f(y) = \Psi \exp(\alpha y + \beta y^2/2 - y^4/4), $$ where \(\Psi\) 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: $$y_i = w_0 + w_1 Y_{i1} + \cdots + w_p Y_{ip}$$ $$\alpha_i = a_0 + a_1 X_{i1} + \cdots + a_p X_{ip}$$ $$\beta_i = b_0 + b_1 X_{i1} + \cdots + b_q X_{iq}$$ in which the \(a_j\)'s, \(b_j\)'s, and \(w_j\)'s are estimated by means of maximum likelihood. Here, the \(Y_{ij}\)'s and \(X_{ij}\)'s are variables constructed from variables in the data set. Variables predicting the \(\alpha\)'s and \(\beta\)'s need not be the same.
The state variable and control parameters can be modelled by specifying a model formula: $$\code{y ~ model},$$ $$\code{alpha ~ model},$$ $$\code{beta ~ model},$$ in which model can be any valid formula specified in terms of variables that are present in the data.frame.
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)
print(fit)
summary(fit)
if (FALSE) {
plot(fit)
cusp3d(fit)
}
# useful use of OK
if (FALSE) {
while(!fit$OK)
fit <- cusp(y ~ z, alpha ~ x1+x2, beta ~ x1+x2, data,
start=rnorm(fit$par)) # use different starting values
}
Run the code above in your browser using DataLab