cusp {cusp}

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 = (\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.

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 \alpha_i's and \beta_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 \alpha[i] and \beta[i],

Examples

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)