Fit Cumulative Link Models
clm.fit.RdA direct fitter of cumulative link models.
Usage
clm.fit(y, ...)
# Default S3 method
clm.fit(y, ...)
# S3 method for class 'factor'
clm.fit(y, X, S, N, weights = rep(1, nrow(X)),
offset = rep(0, nrow(X)), S.offset = rep(0, nrow(X)),
control = list(), start, doFit=TRUE,
link = c("logit", "probit", "cloglog", "loglog", "cauchit",
"Aranda-Ordaz", "log-gamma"),
threshold = c("flexible", "symmetric", "symmetric2", "equidistant"),
...)Arguments
- y
for the default method a list of model components. For the factor method the response variable; a factor, preferably and ordered factor.
- X, S, N
optional design matrices for the regression parameters, scale parameters and nominal parameters respectively.
- weights
optional case weights.
- offset
an optional offset.
- S.offset
an optional offset for the scale part of the model.
- control
a list of control parameters, optionally a call to
clm.control.- start
an optional list of starting values of the form
c(alpha, beta, zeta)for the thresholds and nominal effects (alpha), regression parameters (beta) and scale parameters (zeta).- doFit
logical for whether the model should be fit or the model environment should be returned.
- link
the link function.
- threshold
the threshold structure, see further at
clm.- ...
currently not used.
Details
This function does almost the same thing that clm does:
it fits a cumulative link model. The main differences are that
clm.fit does not setup design matrices from formulae and only
does minimal post processing after parameter estimation.
Compared to clm, clm.fit does little to warn the
user of any problems with data or model. However, clm.fit will
attempt to identify column rank defecient designs. Any unidentified
parameters are indicated in the aliased component of the fit.
clm.fit.factor is not able to check if all thresholds are
increasing when nominal effects are specified since it needs access to
the terms object for the nominal model. If the terms object for the
nominal model (nom.terms) is included in y, the default
method is able to chech if all thresholds are increasing.
Value
A list with the following components:
aliased, alpha, coefficients, cond.H, convergence, df.residual,
edf, fitted.values, gradient, Hessian, logLik, maxGradient, message,
n, niter, nobs, tJac, vcov
and optionally
beta, zeta
These components are documented in clm.
Examples
## A simple example:
fm1 <- clm(rating ~ contact + temp, data=wine)
summary(fm1)
#> formula: rating ~ contact + temp
#> data: wine
#>
#> link threshold nobs logLik AIC niter max.grad cond.H
#> logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> contactyes 1.5278 0.4766 3.205 0.00135 **
#> tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Threshold coefficients:
#> Estimate Std. Error z value
#> 1|2 -1.3444 0.5171 -2.600
#> 2|3 1.2508 0.4379 2.857
#> 3|4 3.4669 0.5978 5.800
#> 4|5 5.0064 0.7309 6.850
## get the model frame containing y and X:
mf1 <- update(fm1, method="design")
names(mf1)
#> [1] "y" "y.levels" "X" "offset" "terms" "contrasts"
#> [7] "xlevels" "weights" "doFit" "control" "link" "threshold"
#> [13] "start" "formulas"
res <- clm.fit(mf1$y, mf1$X) ## invoking the factor method
stopifnot(all.equal(coef(res), coef(fm1)))
names(res)
#> [1] "Hessian" "Theta" "aliased" "alpha"
#> [5] "beta" "coefficients" "cond.H" "convergence"
#> [9] "df.residual" "edf" "fitted.values" "gradient"
#> [13] "logLik" "maxGradient" "message" "n"
#> [17] "niter" "nobs" "tJac" "vcov"
## Fitting with the default method:
mf1$control$method <- "Newton"
res2 <- clm.fit(mf1)
stopifnot(all.equal(coef(res2), coef(fm1)))