Description Usage Arguments Details Value Author(s) References See Also Examples
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
.
1 2 3 
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
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) = Ψ \exp(α y + β y^2/2  y^4/4),
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:
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
:
\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
.
List with components
coefficients 
Estimated coefficients 
rank 
rank of Hessian matrix 
qr 

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 
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 cusppackage
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],
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22  set.seed(123)
# example with regressors
x1 = runif(150)
x2 = runif(150)
z = Vectorize(rcusp)(1, 4*x12, 4*x21)
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)

Call: cusp(formula = y ~ z, alpha = alpha ~ x1 + x2, beta = beta ~ x1 + x2, data = data)
Coefficients:
a[(Intercept)] a[x1] a[x2] b[(Intercept)] b[x1]
2.93474 6.15488 0.16132 2.23399 0.59292
b[x2] w[(Intercept)] w[z]
3.81285 0.07288 0.95201
Degrees of Freedom: 149 Total (i.e. Null); 142 Residual
Null Deviance: 192.7
Delay Deviance: 67.49 AIC: 242.4
Call:
cusp(formula = y ~ z, alpha = alpha ~ x1 + x2, beta = beta ~
x1 + x2, data = data)
Deviance Residuals:
Min 1Q Median 3Q Max
2.47861 0.42193 0.01327 0.37423 1.64024
Coefficients:
Estimate Std. Error z value Pr(>z)
a[(Intercept)] 2.93474 0.61408 4.779 1.76e06 ***
a[x1] 6.15488 0.95874 6.420 1.36e10 ***
a[x2] 0.16132 0.50064 0.322 0.7473
b[(Intercept)] 2.23399 0.89596 2.493 0.0127 *
b[x1] 0.59292 1.16289 0.510 0.6101
b[x2] 3.81285 0.66732 5.714 1.11e08 ***
w[(Intercept)] 0.07288 0.11371 0.641 0.5216
w[z] 0.95201 0.05280 18.032 < 2e16 ***

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Null deviance: 192.732 on 149 degrees of freedom
Linear deviance: 82.579 on 146 degrees of freedom
Logist deviance: NA on NA degrees of freedom
Delay deviance: 67.486 on 142 degrees of freedom
R.Squared logLik npar AIC AICc BIC
Linear model 0.6116778 168.0745 4 344.1491 344.4250 356.1916
Cusp model 0.6521824 113.1932 8 242.3865 243.4078 266.4716

Note: R.Squared for cusp model is Cobb's pseudoR^2. This value
can become negative.
Chisquare test of linear vs. cusp model
Xsquared = 109.8, df = 4, pvalue = 0
Number of optimization iterations: 51
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.