garma                  package:VGAM                  R Documentation

_G_A_R_M_A (_G_e_n_e_r_a_l_i_z_e_d _A_u_t_o_r_e_g_r_e_s_s_i_v_e _M_o_v_i_n_g-_A_v_e_r_a_g_e) _M_o_d_e_l_s

_D_e_s_c_r_i_p_t_i_o_n:

     Fits GARMA models to time series data.

_U_s_a_g_e:

     garma(link = c("identity", "loge", "reciprocal",
                    "logit", "probit", "cloglog", "cauchit"),
           earg=list(),
           p.ar.lag = 1, q.lag.ma = 0,
           coefstart = NULL, step = 1)

_A_r_g_u_m_e_n_t_s:

    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.

_D_e_t_a_i_l_s:

     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.

_V_a_l_u_e:

     An object of class '"vglmff"' (see 'vglmff-class'). The object is
     used by modelling functions such as 'vglm'.

_W_a_r_n_i_n_g:

     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.

_N_o_t_e:

     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.

_A_u_t_h_o_r(_s):

     T. W. Yee

_R_e_f_e_r_e_n_c_e_s:

     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.

_S_e_e _A_l_s_o:

     The site <URL: http://www.stat.auckland.ac.nz/~yee> contains more
     documentation about this family function.

_E_x_a_m_p_l_e_s:

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

