cusp {cusp}R Documentation

Fit a Cusp Catatrophe Model to Data

Description

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.

Usage

cusp(formula, alpha, beta, data, weights, offset, ..., control = 
    glm.control(), method = "cusp.fit", optim.method = "L-BFGS-B", model = TRUE, 
    contrasts = NULL)

Arguments

formula formula that models the canonical state variable
alpha formula that models the canonical normal/asymmetry factor
beta formula that models the canonical bifurcation/splitting factor
data data.frame that contains all the variables named in the formulas
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 optim
control glm.control object, currently unused
method string, currently unused
optim.method string passed to optim to choose the optimization algorithm
model should the model matrix be returned?
contrasts matrix of contrasts, experimental

Details

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) = (α + β Y(t_ - Y(t)^3)dt + dW(t).

The stationary distribution of the ‘behavioral’, or ‘state’ variable Y, given the control parameters α (‘asymmetry’ or ‘normal’ factor) and β (‘bifurcation’ or ‘splitting’ factor) is

f(y) = Psi exp(α y + β 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[i,1] + cdots + w[p] * Y[i,p],

α[i] = a[0] + a[1] * X[i,1] + cdots + a[p] * X[i,p],

β[i] = b[0] + b[1] * X[i,1] + cdots + b[p] * X[i,p],

in which the a[j]'s, b[j]'s, and w[j]'s are estimated by means of maximum likelihood. Here, the Y[i,j]'s and X[i,j]'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 modelled by specifying a model formula:

y ~ model,

alpha ~ model,

beta ~ model,

in which model can be any valid formula specified in terms of variables that are present in the data.frame.

Value

List with components

coefficients Estimated coefficients
rank rank of Hessian matrix
qr qr decomposition of the Hessian matrix
linear.predictors two column matrix containing the α[i]'s and β[i]'s for each case
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 optim
par parameter estimates for qr standardized data
Hessian Hessian matrix of negative log likelihood function at minimum
hessian.untransformed Hessian matrix of negative log likelihood for qr standardized data
code optim convergence indicator
model list with model design matrices
call function call that created the object
formula list with the formulas
OK logical. TRUE if Hessian matrix is positive definite at the minimum obtained
data original data.frame

Author(s)

Raoul Grasman

References

See cusp-package

See Also

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 α[i] and β[i],

Examples

# 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)

[Package cusp version 2.2 Index]