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 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],
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*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)
|
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.76e-06 ***
a[x1] 6.15488 0.95874 6.420 1.36e-10 ***
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.11e-08 ***
w[(Intercept)] 0.07288 0.11371 0.641 0.5216
w[z] 0.95201 0.05280 18.032 < 2e-16 ***
---
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 pseudo-R^2. This value
can become negative.
Chi-square test of linear vs. cusp model
X-squared = 109.8, df = 4, p-value = 0
Number of optimization iterations: 51
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.