gnm | R Documentation |
gnm
fits generalised nonlinear models using an
over-parameterized representation. Nonlinear terms are specified by
calls to functions of class "nonlin"
.
gnm(formula, eliminate = NULL, ofInterest = NULL, constrain = numeric(0),
constrainTo = numeric(length(constrain)), family = gaussian,
data = NULL, subset, weights, na.action, method = "gnmFit",
checkLinear = TRUE, offset, start = NULL, etastart = NULL,
mustart = NULL, tolerance = 1e-06, iterStart = 2, iterMax = 500,
trace = FALSE, verbose = TRUE, model = TRUE, x = TRUE,
termPredictors = FALSE, ridge = 1e-08, ...)
formula |
a symbolic description of the nonlinear predictor. |
eliminate |
a factor to be included as the first term in the
model. |
ofInterest |
optional coefficients of interest, specified by a
regular expression, a numeric vector of indices, a character vector of
names, or "[?]" to select from a Tk dialog. If |
constrain |
(non-eliminated) coefficients to constrain, specified by a regular expression, a numeric vector of indices, a logical vector, a character vector of names, or "[?]" to select from a Tk dialog. |
constrainTo |
a numeric vector of the same length as
|
family |
a specification of the error distribution and link function
to be used in the model. This can be a character string naming
a family function; a family function, or the result of a call
to a family function. See |
data |
an optional data frame containing the variables in the model.
If not found in |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used in the fitting process. |
na.action |
a function which indicates what should happen when the data
contain |
method |
the method to be used: either |
checkLinear |
logical: if |
offset |
this can be used to specify an a priori known component to
be added to the predictor during fitting. |
start |
a vector of starting values for the parameters in the
model; if a starting value is |
etastart |
starting values for the linear predictor. |
mustart |
starting values for the vector of means. |
tolerance |
a positive numeric value specifying the tolerance level for convergence. |
iterStart |
a positive integer specifying the number of start-up iterations to perform. |
iterMax |
a positive integer specifying the maximum number of main iterations to perform. |
trace |
a logical value indicating whether the deviance should be printed after each iteration. |
verbose |
logical: if |
model |
logical: if |
x |
logical: if |
termPredictors |
logical: if |
ridge |
numeric, a positive value for the ridge constant to be used in the fitting algorithm |
... |
further arguments passed to fitting function. |
Models for gnm
are specified by giving a symbolic description
of the nonlinear predictor, of the form response ~ terms
. The
response
is typically a numeric vector, see later in this
section for alternatives. The usual symbolic language may be used to
specify any linear terms, see formula
for details.
Nonlinear terms may be specified by calls to functions of class
"nonlin". There are several "nonlin" functions in the gnm
package. Some of these specify simple
mathematical functions of predictors: Exp
, Mult
, and
Inv
. Others specify more specialised nonlinear terms, in
particular MultHomog
specifies homogeneous multiplicative
interactions and Dref
specifies diagonal reference terms. Users
may also define their own "nonlin" functions, see
nonlin.function
for details.
The eliminate
argument may be used to specify a factor that
is to be included as the first term in the model (since an intercept
is then redundant, none is fitted). The structure of the factor is
exploited to improve computational efficiency — substantially so if
the eliminate
d factor has a large number of levels. Use of
eliminate
is designed for factors that are required in the
model but are not of direct interest (e.g., terms needed to fit
multinomial-response models as conditional Poisson models). See
backPain
for an example.
The ofInterest
argument may be used to specify coefficients of
interest, the indices of which are returned in the ofInterest
component of the model object. print()
displays of the model
object or its components obtained using accessor functions such as
coef()
etc, will only show these coefficients. In addition
methods for "gnm"
objects which may be applied to a subset of
the parameters are by default applied to the coefficients of interest.
See ofInterest
for accessor and replacement functions.
For contingency tables, the data may be provided as an object of class
"table"
from which the frequencies will be extracted to use
as the response. In this case, the response should be specified as
Freq
in the model formula. The "predictors"
,
"fitted.values"
, "residuals"
, "prior.weights"
,
"weights"
, "y"
and "offset"
components of
the returned gnm
fit will be tables with the same format as the
data, completed with NA
s where necessary.
For binomial models, the response
may be specified as a factor
in which the first level denotes failure and all other levels denote
success, as a two-column matrix with the columns giving the numbers
of successes and failures, or as a vector of the proportions of
successes.
The gnm
fitting algorithm consists of two stages. In the start-up
iterations, any nonlinear parameters that are not specified by either the
start
argument of gnm
or a plug-in function are
updated one parameter at a time, then the linear parameters are
jointly updated before the next iteration. In the main iterations, all
the parameters are jointly updated, until convergence is reached or
the number or iterations reaches iterMax
. To solve the
(typically rank-deficient) least squares problem at the heart of the
gnm
fitting algorithm, the design matrix is standardized and
regularized (in the Levenberg-Marquardt sense) prior to solving; the
ridge
argument provides a degree of control over the
regularization performed (smaller values may sometimes give faster
convergence but can lead to numerical instability).
Convergence is judged by comparing the squared components of the score vector
with corresponding elements of the diagonal of the Fisher information
matrix. If, for all components of the score vector, the ratio is less
than tolerance^2
, or the corresponding diagonal element of the
Fisher information matrix is less than 1e-20, iterations cease. If the
algorithm has not converged by iterMax
iterations,
exitInfo
can be used to print information on the
parameters which failed the convergence criteria at the last iteration.
By default, gnm
uses an over-parameterized representation of
the model that is being fitted. Only minimal identifiability constraints
are imposed, so that in general a random parameterization is obtained.
The parameter estimates are ordered so that those for any linear terms
appear first.
getContrasts
may be used to obtain estimates of
specified scaled contrasts, if these contrasts are identifiable. For
example, getContrasts
may be used to estimate the contrasts
between the first level of a factor and the rest, and obtain standard
errors.
If appropriate constraints are known in advance, or have been
determined from a gnm
fit, the model may be (re-)fitted using
the constrain
argument to specify coefficients which should be
set to values specified by constrainTo
. Constraints should only
be specified for non-eliminated parameters. update
provides a convenient way of re-fitting a gnm
model with new
constraints.
If method = "gnmFit"
, gnm
returns NULL
if the
algorithm has failed and an object of class "gnm"
otherwise. A
"gnm"
object inherits first from "glm"
then "lm"
and is a list containing the following components:
call |
the matched call. |
formula |
the formula supplied. |
constrain |
a numeric vector specifying any coefficients that were constrained in the fitting process. |
constrainTo |
a numeric vector of the same length as
|
family |
the |
prior.weights |
the case weights initially supplied. |
terms |
the |
data |
the |
na.action |
the |
xlevels |
a record of the levels of the factors used in fitting. |
y |
the response used. |
offset |
the offset vector used. |
coefficients |
a named vector of non-eliminated coefficients,
with an attribute |
eliminate |
the |
ofInterest |
a named numeric vector of indices corresponding
to non-eliminated coefficients, or |
predictors |
the fitted values on the link scale. |
fitted.values |
the fitted mean values, obtained by transforming the predictors by the inverse of the link function. |
deviance |
up to a constant, minus twice the maximised log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero. |
aic |
Akaike's An Information Criterion, minus twice the maximized log-likelihood plus twice the number of parameters (so assuming that the dispersion is known). |
iter |
the number of main iterations. |
conv |
logical indicating whether the main iterations
converged, with an attribute for use by |
weights |
the working weights, that is, the weights used in the last iteration. |
residuals |
the working residuals, that is, the residuals from the last iteration. |
df.residual |
the residual degrees of freedom. |
rank |
the numeric rank of the fitted model. |
The list may also contain the components model
, x
,
or termPredictors
if requested in the arguments to gnm
.
If a table was passed to data
and the default for
na.action
was not overridden, the list will also contain a
table.attr
component, for use by the extractor functions.
If a binomial gnm
model is specified by giving a two-column
response, the weights returned by prior.weights
are the total
numbers of cases (factored by the supplied case weights) and the
component y
of the result is the proportion of successes.
The function summary.gnm
may be used to obtain and print
a summary of the results, whilst plot.gnm
may be used
for model diagnostics.
The generic functions formula
, family
,
terms
, coefficients
,
fitted.values
, deviance
,
extractAIC
, weights
,
residuals
, df.residual
,
model.frame
, model.matrix
,
vcov
and termPredictors
maybe used to
extract components from the object returned by gnm
or to
construct the relevant objects where necessary.
Note that the generic functions weights
and
residuals
do not act as straight-forward accessor
functions for gnm
objects, but return the prior weights and
deviance residuals respectively, as for glm
objects.
Regular expression matching is performed using grep
with
default settings.
Heather Turner and David Firth
Cautres, B, Heath, A F and Firth, D (1998). Class, religion and vote in Britain and France. La Lettre de la Maison Francaise 8.
formula
for the symbolic language used to specify
formulae.
Diag
and Symm
for specifying special types
of interaction.
Exp
, Mult
, Inv
, MultHomog
,
Dref
and nonlin.function
for incorporating
nonlinear terms in the formula
argument to gnm
.
residuals.glm
and the generic functions
coef
, fitted
, etc. for extracting
components from gnm
objects.
exitInfo
to print more information on last iteration
when gnm
has not converged.
getContrasts
to estimate (identifiable) scaled contrasts
from a gnm
model.
### Analysis of a 4-way contingency table
set.seed(1)
print(cautres)
## Fit a "double UNIDIFF" model with the religion-vote and class-vote
## interactions both modulated by nonnegative election-specific
## multipliers.
doubleUnidiff <- gnm(Freq ~ election:vote + election:class:religion
+ Mult(Exp(election), religion:vote) +
Mult(Exp(election), class:vote), family = poisson,
data = cautres)
## Examine the multipliers of the class-vote log odds ratios
ofInterest(doubleUnidiff) <- pickCoef(doubleUnidiff, "class:vote[).]")
coef(doubleUnidiff)
## Coefficients of interest:
## Mult(Exp(.), class:vote).election1
## -0.38357138
## Mult(Exp(.), class:vote).election2
## 0.29816599
## Mult(Exp(.), class:vote).election3
## 0.06580307
## Mult(Exp(.), class:vote).election4
## -0.02174104
## Re-parameterize by setting first multiplier to zero
getContrasts(doubleUnidiff, ofInterest(doubleUnidiff))
## estimate SE
## Mult(Exp(.), class:vote).election1 0.0000000 0.0000000
## Mult(Exp(.), class:vote).election2 0.6817374 0.2401644
## Mult(Exp(.), class:vote).election3 0.4493745 0.2473521
## Mult(Exp(.), class:vote).election4 0.3618301 0.2534754
## quasiSE quasiVar
## Mult(Exp(.), class:vote).election1 0.22854401 0.052232363
## Mult(Exp(.), class:vote).election2 0.07395886 0.005469913
## Mult(Exp(.), class:vote).election3 0.09475938 0.008979340
## Mult(Exp(.), class:vote).election4 0.10934798 0.011956981
## Same thing but with last multiplier as reference category:
getContrasts(doubleUnidiff, rev(ofInterest(doubleUnidiff)))
## estimate SE
## Mult(Exp(.), class:vote).election4 0.00000000 0.0000000
## Mult(Exp(.), class:vote).election3 0.08754436 0.1446833
## Mult(Exp(.), class:vote).election2 0.31990727 0.1320022
## Mult(Exp(.), class:vote).election1 -0.36183013 0.2534754
## quasiSE quasiVar
## Mult(Exp(.), class:vote).election4 0.10934798 0.011956981
## Mult(Exp(.), class:vote).election3 0.09475938 0.008979340
## Mult(Exp(.), class:vote).election2 0.07395886 0.005469913
## Mult(Exp(.), class:vote).election1 0.22854401 0.052232363
## Re-fit model with first multiplier set to zero
doubleUnidiffConstrained <-
update(doubleUnidiff, constrain = ofInterest(doubleUnidiff)[1])
## Examine the multipliers of the class-vote log odds ratios
coef(doubleUnidiffConstrained)[ofInterest(doubleUnidiff)]
## ...as using 'getContrasts' (to 4 d.p.).
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.