aster {aster}

Aster Models

Description

Fits Aster Models.

Usage

aster(x, ...)

## Default S3 method:
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 class '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, ...)

Arguments

x

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.

root

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.

pred

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.

fam

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.

modmat

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.

parm

usually missing. Otherwise a vector of length ncoef giving a starting point for the optimization.

type

type of model. The value of this argument can be abbreviated.

famlist

a list of family specifications (see families).

origin

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.

origin.type

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.

method

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.

fscale

an estimate of the size of the log likelihood at the maximum. Defaults to nind.

maxiter

maximum number of iterations. Defaults to '1000'.

nowarn

if TRUE (the default), suppress warnings from the optimization routine.

newton

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.

optout

if TRUE, save the entire result of the optimization routine (trust, nlm, or optim, as the case may be).

coef.names

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.

formula

a symbolic description of the model to be fit. See lm, glm, and formula for discussions of the R formula mini-language.

varvar

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.

idvar

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.

data

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.

Details

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".

Value

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:

coefficients

a named vector of coefficients.

rank

the numeric rank of the fitted generalized linear model part of the aster model (i.e., the rank of modmat).

deviance

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.)

iter

the number of iterations used by the optimization method.

converged

logical. Was the optimization algorithm judged to have converged?

code

integer. The convergence code returned by the optimization method.

gradient

The gradient vector of minus the log likelihood at the fitted coefficients vector.

hessian

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".

fisher

Expected Fisher information at the fitted coefficients vector.

optout

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:

call

the matched call.

formula

the formula supplied.

terms

the terms object used.

data

the data argument.

NA Values

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”.

Directions of Recession

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.

References

Geyer, C. J., Wagenius, S., and Shaw, R. G. (2007) Aster Models for Life History Analysis. Biometrika, 94, 415–426. 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. doi:10.1086/588063.

See Also

anova.aster, summary.aster, and predict.aster

Examples

### 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)