quasibin                 package:aod                 R Documentation

_Q_u_a_s_i-_L_i_k_e_l_i_h_o_o_d _M_o_d_e_l _f_o_r _P_r_o_p_o_r_t_i_o_n_s

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

     The function fits the generalized linear model "II" proposed by
     Williams (1982) accounting  for overdispersion in clustered
     binomial data (n, y).

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

     quasibin(formula, data, link = "logit", phi = NULL, tol =  0.001)

_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 of the form  'cbind(y, n - y)' where the
          modelled probability is 'y/n'.

    link: The link function for the mean p: "logit" or "cloglog".

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

     phi: When 'phi' is NULL (the default), the overdispersion
          parameter phi is estimated from the data.  Otherwise, its
          value is considered as fixed.

     tol: A positive scalar (default to 0.001). The algorithm stops at
          iteration r + 1 when the condition chi{^2}[r+1] - chi{^2}[r]
          <= tol is met by the chi-squared statistics .

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

     For a given cluster (n, y), the model is:

                   y | lambda ~ Binomial(n, lambda)

     with lambda a random variable of mean E[lambda] = p and variance
     Var[lambda] = phi * p * (1 - p).
      The marginal mean and variance are:

                               E[y] = p


              Var[y] = p * (1 - p) * [1 + (n - 1) * phi]

     The overdispersion parameter phi corresponds to the intra-cluster
     correlation coefficient, which is here restricted to be positive.
      The function uses the function 'glm' and the parameterization: p
     = h(X b) = h(eta), where h is the inverse of a given link
     function, X is a design-matrix, b is a vector of fixed effects and
     eta = X b is the linear predictor. 
      The estimate of b maximizes the quasi log-likelihood of the
     marginal model. The parameter phi is estimated with the moment
     method or can be set to a constant (a regular _glim_ is fitted
     when phi is set to zero). The literature recommends to estimate
     phi from the saturated model. Several explanatory variables are
     allowed in b. None is allowed in phi.

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

     An object of formal class "glimQL": see 'glimQL-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:

     Moore, D.F., 1987, _Modelling the extraneous variance in the
     presence of extra-binomial variation_.  Appl. Statist. 36, 8-14.
      Williams, D.A., 1982, _Extra-binomial variation in logistic
     linear models_. Appl. Statist. 31, 144-148.

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

     'glm', 'geese' in the contributed package 'geepack', 
     'glm.binomial.disp' in the contributed package 'dispmod'.

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

       data(orob2) 
       fm1 <- glm(cbind(y, n - y) ~ seed * root,
                  family = binomial, data = orob2)
       fm2 <- quasibin(cbind(y, n - y) ~ seed * root,
                       data = orob2, phi = 0)
       fm3 <- quasibin(cbind(y, n - y) ~ seed * root,
                       data = orob2)
       rbind(fm1 = coef(fm1), fm2 = coef(fm2), fm3 = coef(fm3))
       # show the model
       fm3
       # dispersion parameter and goodness-of-fit statistic
       c(phi = fm3@phi, 
         X2 = sum(residuals(fm3, type = "pearson")^2))
       # model predictions
       predfm1 <- predict(fm1, type = "response", se = TRUE)
       predfm3 <- predict(fm3, type = "response", se = TRUE)
       New <- expand.grid(seed = levels(orob2$seed),
                          root = levels(orob2$root))
       predict(fm3, New, se = TRUE, type = "response")
       data.frame(orob2, p1 = predfm1$fit, 
                         se.p1 = predfm1$se.fit,
                         p3 = predfm3$fit,
                         se.p3 = predfm3$se.fit)
       fm4 <- quasibin(cbind(y, n - y) ~ seed + root,
                       data = orob2, phi = fm3@phi)
       # Pearson's chi-squared goodness-of-fit statistic
       # compare with fm3's X2
       sum(residuals(fm4, type = "pearson")^2)
       

