Fits Aster Models.
aster(x, ...)# S3 method for default
aster(x, root, pred, fam, modmat, parm,
type = c("unconditional", "conditional"), famlist = fam.default(),
origin, origin.type = c("model.type", "unconditional", "conditional"),
method = c("trust", "nlm", "CG", "L-BFGS-B"), fscale, maxiter = 1000,
nowarn = TRUE, newton = TRUE, optout = FALSE, coef.names, ...)
# S3 method for formula
aster(formula, pred, fam, varvar, idvar, root,
data, parm, type = c("unconditional", "conditional"),
famlist = fam.default(),
origin, origin.type = c("model.type", "unconditional", "conditional"),
method = c("trust", "nlm", "CG", "L-BFGS-B"), fscale, maxiter = 1000,
nowarn = TRUE, newton = TRUE, optout = FALSE, ...)
aster returns an object of class inheriting from "aster".
aster.formula, returns an object of class "aster" and
subclass "aster.formula".
The function summary (i.e., summary.aster) can
be used to obtain or print a summary of the results, the function
anova (i.e., anova.aster)
to produce an analysis of deviance table, and the function
predict (i.e., predict.aster)
to produce predicted values and standard errors.
An object of class "aster" is a list containing at least the
following components:
a named vector of coefficients.
the numeric rank of the fitted generalized linear model
part of the aster model (i.e., the rank of modmat).
up to a constant, minus twice the maximized log-likelihood. (Note the minus. This is somewhat counterintuitive, but cannot be changed for reasons of backward compatibility.)
the number of iterations used by the optimization method.
logical. Was the optimization algorithm judged to have converged?
integer. The convergence code returned by the optimization method.
The gradient vector of minus the log likelihood at the
fitted coefficients vector.
The Hessian matrix of minus the log likelihood
(i.e., the observed Fisher information) at the
fitted coefficients vector.
This is also the expected Fisher information when
type == "unconditional".
Expected Fisher information at the fitted coefficients
vector.
The object returned by the optimization routine
(trust, nlm, or optim).
Only returned when the argument optout is TRUE.
Calls to aster.formula return a list also containing:
the matched call.
the formula supplied.
the terms object used.
the data argument.
an nind by nnode matrix, the data for an
aster model. The rows are independent and identically modeled
random vectors. See details below for further requirements.
aster.formula constructs such an x from the response
in its formula. Hence data for aster.formula must be a
vector of length nind * nnode.
an object of the same shape as x, the root data.
For aster.default an nind by nnode matrix,
For aster.formula an nind * nnode vector.
an integer vector of length nnode determining
the dependence
graph of the aster model. pred[j] is
the index of the predecessor of
the node with index j unless the predecessor is a root
node, in which case pred[j] == 0. See details below for
further requirements.
an integer vector of length nnode determining
the exponential family structure of the aster model. Each element
is an index into the vector of family specifications given by
the argument famlist.
an nind by nnode by ncoef
three-dimensional array, the model matrix.
aster.formula constructs such a modmat from
its formula, the data frame data, and the variables
in the environment of the formula.
usually missing. Otherwise a vector of length ncoef
giving a starting point for the optimization.
type of model. The value of this argument can be abbreviated.
a list of family specifications (see families).
Distinguished point in parameter space. May be missing,
in which case an unspecified default is provided. See details below
for further explanation. This is what lm,
glm and other functions that do regression call
“offset” but we don't change our name for reasons of backward
compatibility.
Parameter space in which specified distinguished point
is located. If "conditional" then argument "origin" is
a conditional canonical parameter value.
If "unconditional" then argument "origin" is
an unconditional canonical parameter value.
If "model.type" then the type is taken from argument "type".
The value of this argument can be abbreviated.
optimization method. If "trust" then the
trust function is used. If "nlm" then the
nlm function is used. Otherwise the
optim function is used with the specified method
supplied to it.
The value of this argument can be abbreviated.
an estimate of the size of the log likelihood at the maximum.
Defaults to nind.
maximum number of iterations. Defaults to '1000'.
if TRUE (the default), suppress warnings from
the optimization routine.
if TRUE (the default), do one Newton iteration
on the result produced by the optimization routine, except when
method == "trust" when no such Newton iteration is done,
regardless of the value of newton, because trust
always terminates with a Newton iteration when it converges.
if TRUE, save the entire result of the optimization
routine (trust, nlm, or optim,
as the case may be).
names of the regression coefficients. If missing,
dimnames(modmat)[[3]] is used. In aster.formula these
are produced automatically by the R formula machinery.
other arguments passed to the optimization method.
a symbolic description of the model to be fit. See
lm, glm, and
formula for discussions of the R formula mini-language.
a variable of the same length as the response in
the formula that is a factor whose levels are character strings
treated as variable names. The number of variable names is nnode.
Must be of the form rep(vars, each = nind) where vars is
a vector of variable names. Usually found in the data frame data
when this is produced by the reshape function.
a variable of the same length as the response in
the formula that indexes individuals. The number
of individuals is nind.
Must be of the form rep(inds, times = nnode) where inds is
a vector of labels for individuals. Usually found in the data frame
data when this is produced by the reshape function.
an optional data frame containing the variables
in the model. If not found in data, the variables are taken
from environment(formula), typically the environment from
which aster is called. Usually produced by
the reshape function.
It was almost always wrong for aster model data to have NA values.
Although theoretically possible for the R formula mini-language to do the
right thing for an aster model with NA values in the data, usually
it does some wrong thing. Thus, since version 0.8-20, this function and
the reaster function give errors when used with data having
NA values. Users must remove all NA values (or replace them
with what they should be, perhaps zero values) “by hand”.
Even if the result of this function
has component converges equal to TRUE, the result will be
nonsense if there are one or more directions of recession. These are
not detected by this function, but rather by the summary function
applied to the result of this function,
for which see summary.aster.
The vector pred must satisfy all(pred < seq(along = pred)),
that is, each predecessor must precede in the order given in pred.
The vector pred defines a function \(p\).
The joint distribution of the data matrix x is a product of conditionals
$$\prod_{i \in I} \prod_{j \in J} \Pr \{ X_{i j} | X_{i p(j)} \}$$
When \(p(j) = 0\), the notation \(X_{i p(j)}\)
means root[i, j]. Other elements of the matrix root are
not used.
The conditional distribution
\(\Pr \{ X_{i j} | X_{i p(j)} \}\)
is the \(X_{i p(j)}\)-fold convolution of the \(j\)-th family
in the vector fam, a one-parameter exponential family
(i.e., the sum of \(X_{i p(j)}\) i.i.d. terms having
this one-parameter exponential family distribution).
For type == "conditional" the canonical parameter vector
\(\theta_{i j}\) is modeled in GLM fashion as
\(\theta = a + M \beta\) where \(M\) is the model
matrix modmat and \(a\) is the distinguished point origin.
Since the “vector” \(\theta\) is
actually a matrix, the “matrix” \(M\) must correspondingly
be a three-dimensional array. So \(\theta = a + M \beta\)
written out in full is
$$\theta_{i j} = a_{i j} + \sum_{k \in K} m_{i j k} \beta_k$$
This specifies the log likelihood.
For type == "unconditional" the canonical parameter vector
for an unconditional model is modeled in GLM fashion as
\(\varphi = a + M \beta\) (where the notation is as above).
The unconditional canonical parameters are then specified in terms of
the conditional ones by
$$\varphi_{i j} = \theta_{i j} - \sum_{k \in S(j)} \psi_k(\theta_{i k})$$
where \(S(j)\) denotes the set of successors of \(j\),
the \(k\) such that \(p(k) = j\), and \(\psi_k\) is the
cumulant function for the \(k\)-th exponential family.
This rather crazy looking formulation is an invertible change of parameter
and makes \(\varphi\)
the canonical parameter and \(x\) the canonical statistic of a full
flat unconditional exponential family.
Again, this specifies the log likelihood.
In versions of aster prior to version 0.6 there was no \(a\) in the model specification, which is the same as specifying \(a = 0\) in the current specification. If \(a\) is in the column space of the model matrix, that is, if there exists an \(\alpha\) such that \(a = M \alpha\), then there is no difference in the model specified with \(a\) and the one with \(a = 0\). The maximum likelihood regression coefficients \(\hat{\beta}\) will be different, but the maximum likelihood estimates of all other parameters (conditional and unconditional, canonical and mean value) will be the same. This is the usual case and explains why “linear” models (with \(a = 0\)) as opposed to “affine” models (with general \(a\)) are popular. In the unusual case where \(a\) is not in the column space of the design matrix, then affine models are a generalization of linear models: the two are not equivalent, their maximum likelihood estimates are not the same in any parameterization.
In order to use the R model formula mini-language we must flatten
the dimensionality, making the model matrix modmat two-dimensional
(a true matrix). This must be done as if by
matrix(modmat, ncol = ncoef),
which imposes the requirements on varvar and idvar
given in the arguments section: they must look like row(x) and
col(x) modulo relabeling.
Then x and root
become one-dimensional, done as if by as.numeric(x)
and as.numeric(root).
The standard way to do this in R is to use the reshape
function on a data frame in which the columns of the x matrix
are variables in the data frame. reshape automatically puts
things in the right order and creates varvar and idvar.
This is shown in the examples section below and in the package vignette
titled "Aster Package Tutorial".
Geyer, C. J., Wagenius, S., and Shaw, R. G. (2007) Aster Models for Life History Analysis. Biometrika, 94, 415--426. tools:::Rd_expr_doi("10.1093/biomet/asm030").
Shaw, R. G., Geyer, C. J., Wagenius, S., Hangelbroek, H. H., and Etterson, J. R. (2008) Unifying Life History Analysis for Inference of Fitness and population growth. American Naturalist, 172, E35--E47. tools:::Rd_expr_doi("10.1086/588063").
anova.aster,
summary.aster,
and
predict.aster
### see package vignette for explanation ###
library(aster)
data(echinacea)
vars <- c("ld02", "ld03", "ld04", "fl02", "fl03", "fl04",
"hdct02", "hdct03", "hdct04")
redata <- reshape(echinacea, varying = list(vars), direction = "long",
timevar = "varb", times = as.factor(vars), v.names = "resp")
redata <- data.frame(redata, root = 1)
pred <- c(0, 1, 2, 1, 2, 3, 4, 5, 6)
fam <- c(1, 1, 1, 1, 1, 1, 3, 3, 3)
hdct <- grepl("hdct", as.character(redata$varb))
redata <- data.frame(redata, hdct = as.integer(hdct))
level <- gsub("[0-9]", "", as.character(redata$varb))
redata <- data.frame(redata, level = as.factor(level))
aout <- aster(resp ~ varb + level : (nsloc + ewloc) + hdct : pop,
pred, fam, varb, id, root, data = redata)
summary(aout, show.graph = TRUE)
Run the code above in your browser using DataLab