garma {VGAM} | R Documentation |
Fits GARMA models to time series data.
garma(link = c("identity", "loge", "reciprocal", "logit", "probit", "cloglog", "cauchit"), earg=list(), p.ar.lag = 1, q.lag.ma = 0, coefstart = NULL, step = 1)
link |
Link function applied to the mean response.
By default, the first choice is used, which is suitable for
continuous responses.
The link loge should be chosen if the data are counts.
The links logit , probit ,
cloglog ,
cauchit are suitable for binary responses.
|
earg |
List. Extra argument for the link.
See earg in Links for general information.
In particular, this argument is useful
when the log or logit link is chosen:
for log and logit,
zero values can be replaced by bvalue which
is inputted as earg=list(bvalue = bvalue) .
See loge and logit etc. for specific
information about each link function.
|
p.ar.lag |
A positive integer,
the lag for the autoregressive component.
Called p below.
|
q.lag.ma |
A non-negative integer,
the lag for the moving-average component.
Called q below.
|
coefstart |
Starting values for the coefficients.
For technical reasons, the
argument coefstart in vglm cannot be used.
|
step |
Numeric. Step length, e.g., 0.5 means half-stepsizing.
|
This function draws heavily on Benjamin et al. (1998).
See also Benjamin et al. (2003).
GARMA models extend the ARMA time series model to generalized
responses in the exponential family, e.g., Poisson counts,
binary responses. Currently, this function can handle continuous,
count and binary responses only. The possible link functions
given in the link
argument reflect this, and the user
must choose an appropriate link.
The GARMA(p,q) model is defined by firstly having a response belonging to the exponential family
f(y_t|D_t) = exp [ (y_t theta_t - b(theta_t)) / (phi / A_t) + c(y_t, phi / A_t) ]
where theta_t and phi are the canonical and scale parameters respectively, and A_t are known prior weights. The mean mu_t=E(Y_t|D_t)=b'(theta_t) is related to the linear predictor eta_t by the link function g. Here, D_t={x_t,...,x_1,y_(t-1),...,y_1,mu_(t-1),...,mu_1} is the previous information set. Secondly, the GARMA(p,q) model is defined by
g(mu_t) = eta_t = x_t^T beta + sum_{k=1}^p phi_k (g(y_{t-k}) - x_{t-k}^T beta) + sum_{k=1}^q theta_k (g(y_{t-k}) - eta_{t-k}).
Parameter vectors beta, phi and theta are estimated by maximum likelihood.
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions such as vglm
.
This VGAM family function is `non-standard' in that the model does need some coercing to get it into the VGLM framework. Special code is required to get it running. A consequence is that some methods functions may give wrong results when applied to the fitted object.
This function is unpolished and is requires lots
of improvements. In particular, initialization is quite poor,
and could be improved.
A limited amount of experience has shown that half-stepsizing is
often needed for convergence, therefore choosing crit="coef"
is not recommended.
Overdispersion is not handled.
T. W. Yee
Benjamin, M. A., Rigby, R. A. and Stasinopoulos, M. D. (1998) Fitting Non-Gaussian Time Series Models. Pages 191–196 in: Proceedings in Computational Statistics COMPSTAT 1998 by Payne, R. and P. J. Green. Physica-Verlag.
Benjamin, M. A., Rigby, R. A. and Stasinopoulos, M. D. (2003) Generalized Autoregressive Moving Average Models. Journal of the American Statistical Association, 98: 214–223.
Zeger, S. L. and Qaqish, B. (1988) Markov regression models for time series: a quasi-likelihood approach. Biometrics, 44: 1019–1031.
The site http://www.stat.auckland.ac.nz/~yee contains more documentation about this family function.
# See Zeger and Qaqish (1988) interspike = c(68, 41, 82, 66, 101, 66, 57, 41, 27, 78, 59, 73, 6, 44, 72, 66, 59, 60, 39, 52, 50, 29, 30, 56, 76, 55, 73, 104, 104, 52, 25, 33, 20, 60, 47, 6, 47, 22, 35, 30, 29, 58, 24, 34, 36, 34, 6, 19, 28, 16, 36, 33, 12, 26, 36, 39, 24, 14, 28, 13, 2, 30, 18, 17, 28, 9, 28, 20, 17, 12, 19, 18, 14, 23, 18, 22, 18, 19, 26, 27, 23, 24, 35, 22, 29, 28, 17, 30, 34, 17, 20, 49, 29, 35, 49, 25, 55, 42, 29, 16) spikenum = seq(interspike) bvalue = 0.1 # .Machine$double.xmin # Boundary value fit = vglm(interspike ~ 1, trace=TRUE, garma("loge", earg=list(bvalue=bvalue), p=2, coef=c(4,.3,.4))) summary(fit) coef(fit, matrix=TRUE) Coef(fit) # A bug here ## Not run: plot(interspike, ylim=c(0,120), las=1, font=1, xlab="Spike Number", ylab="Inter-Spike Time (ms)", col="blue") lines(spikenum[-(1:fit@misc$plag)], fitted(fit), col="green") abline(h=mean(interspike), lty=2) ## End(Not run)