negbin                  package:aod                  R Documentation

_N_e_g_a_t_i_v_e-_B_i_n_o_m_i_a_l _M_o_d_e_l _f_o_r _C_o_u_n_t_s

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

     The function fits a negative-binomial log linear model accounting
     for overdispersion in counts y.

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

       negbin(formula, random, data, phi.ini = NULL, warnings = FALSE, 
              na.action = na.omit, fixpar = list(),
              hessian = TRUE, ...)
       

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

 formula: A formula for the fixed effects. The left-hand side of the
          formula must be the counts 'y' i.e.,  positive integers ('y
          >= 0'). The right-hand side can involve an offset term.

  random: A right-hand formula for the overdispersion parameter(s) phi.

    data: A data frame containing the response ('y') and explanatory
          variable(s).

 phi.ini: Initial values for the overdispersion parameter(s) phi.
          Default to 0.1.

warnings: Logical to control printing of warnings occurring during
          log-likelihood maximization. Default to  FALSE (no printing).

na.action: A function name. Indicates which action should be taken in
          the case of missing value(s).

  fixpar: A list with 2 components (scalars or vectors) of the same
          size, indicating which parameters are  fixed (i.e., not
          optimized) in the global parameter vector (b, phi) and the
          corresponding fixed values.
           For example, 'fixpar = list(c(4, 5), c(0, 0))' means that
          4th and 5th parameters of the model are set to 0.

 hessian: A logical. When set to 'FALSE', the hessian and the
          variances-covariances matrices of the  parameters are not
          computed.

     ...: Further arguments passed to 'optim'.

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

     For a given count y, the model is:

                     y | lambda ~ Poisson(lambda)

     with lambda following a Gamma distribution Gamma(r, theta).
      If G denote the gamma function, then:

 P(lambda) = r^{-theta} * lambda^{theta - 1} * exp(-lambda / r) / G(theta)


                        E[lambda] = theta * r


                      Var[lambda] = theta * r^2

     The marginal negative-binomial distribution is:

 P(y) = G(y + theta) * (1 / (1 + r))^theta * (r / (1 + r))^y / (y! * G(theta))

     The function uses the recrystallization mu = theta * r = exp(X b)
     = exp(eta) and phi = 1 / theta,  where X is a design-matrix, b is
     a vector of fixed effects, eta = X b is the linear predictor and 
     phi the overdispersion parameter.
      The marginal mean and variance are:

                              E[y] = mu


                       Var[y] = mu + phi * mu^2

     The parameters b and phi are estimated by maximizing the
     log-likelihood of the marginal model (using the  function
     'optim()'). Several explanatory variables are allowed in b. Only
     one is allowed in phi.
      An offset can be specified in the 'formula' argument to model
     rates y/T. The offset and the marginal mean are log(T) and mu =
     exp(log(T) + eta), respectively.

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

     An object of formal class "glimML": see 'glimML-class' for
     details.

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

     Matthieu Lesnoff matthieu.lesnoff@cirad.fr, Renaud Lancelot
     renaud.lancelot@cirad.fr

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

     Lindsey, J.K., 1999, 2nd ed. _Models for Repeated Measurements_.
     Oxford, London.
      Venables, W.N., Ripley, B.D., 2002, 4th ed. _Modern applied
     statistics with S_. Springer-Verlag, New York.

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

     'glimML-class', 'glm' and 'optim',
      'glm.nb' in the recommended package 'MASS',
      'gnlr' in package 'gnlm' available at <URL:
     www.luc.ac.be/~jlindsey/rcode.html>.

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

       # without offset
       data(salmonella)
       negbin(y ~ log(dose + 10) + dose, ~ 1, salmonella)
       library(MASS) # function glm.nb in MASS
       fm.nb <- glm.nb(y ~ log(dose + 10) + dose,
                       link = log, data = salmonella)
       coef(fm.nb)
       1 / fm.nb$theta # theta = 1 / phi
       c(logLik(fm.nb), AIC(fm.nb))
       # with offset
       data(dja)
       negbin(y ~ group + offset(log(trisk)), ~ group, dja)
       # phi fixed to zero in group TREAT
       negbin(y ~ group + offset(log(trisk)), ~ group, dja,
         fixpar = list(4, 0))
       # glim without overdispersion
       summary(glm(y ~ group + offset(log(trisk)),
         family = poisson, data = dja))
       # phi fixed to zero in both groups
       negbin(y ~ group + offset(log(trisk)), ~ group, dja,
         fixpar = list(c(3, 4), c(0, 0))) 
       

