zipebcom                package:VGAM                R Documentation

_E_x_c_h_a_n_g_e_a_b_l_e _b_i_v_a_r_i_a_t_e _c_l_o_g_l_o_g _o_d_d_s-_r_a_t_i_o _m_o_d_e_l _f_r_o_m _a
_z_e_r_o-_i_n_f_l_a_t_e_d _P_o_i_s_s_o_n _d_i_s_t_r_i_b_u_t_i_o_n

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

     Fits an exchangeable bivariate odds-ratio model to two binary
     responses with a complementary log-log link. The data are assumed
     to come from a zero-inflated Poisson distribution that has been
     converted to presence/absence.

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

     zipebcom(lmu12="cloglog", lphi12="logit", loratio="loge",
              emu12=list(), ephi12=list(), eoratio=list(),
              imu12=NULL, iphi12=NULL, ioratio = NULL, 
              zero=2:3, tol=0.001, addRidge=0.001)

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

lmu12, emu12, imu12: Link function, extra argument and optional initial
          values for the first (and second) marginal probabilities.
          Arguments 'lmu12' and 'emu12' should be left alone. Argument
          'imu12' may be of length 2 (one element for each response).

  lphi12: Link function applied to the phi parameter of the
          zero-inflated Poisson distribution (see 'zipoisson'). See
          'Links' for more choices.

 loratio: Link function applied to the odds ratio. See 'Links' for more
          choices.

iphi12, ioratio: Optional initial values for phi and the odds ratio.
          See 'CommonVGAMffArguments' for more details. In general,
          good initial values (especially for 'iphi12') are often
          required, therefore use these arguments if convergence
          failure occurs. If inputted, the value of 'iphi12' cannot be
          more than the sample proportions of zeros in either response.

ephi12, eoratio: List. Extra argument for each of the links. See 'earg'
          in 'Links' for general information.

    zero: Which linear/additive predictor is modelled as an intercept
          only? A 'NULL' means none. The default has both phi and the
          odds ratio as not being modelled as a function of the
          explanatory variables (apart from an intercept).

     tol: Tolerance for testing independence. Should be some small
          positive numerical value.

addRidge: Some small positive numerical value. The first two diagonal
          elements of the working weight matrices are  multiplied by
          '1+addRidge' to make it diagonally dominant, therefore
          positive-definite.

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

     This 'VGAM' family function fits an exchangeable bivariate odds
     ratio model ('binom2.or') with a 'cloglog' link. The data are
     assumed to come from a zero-inflated Poisson (ZIP) distribution
     that has been converted to presence/absence. Explicitly, the
     default model is

             cloglog[P(Y_j=1)/(1-phi)] =  eta_1,   j=1,2

     for the (exchangeable) marginals, and

                         logit[phi] =  eta_2,

     for the mixing parameter, and

 log[P(Y_{00}=1) P(Y_{11}=1) / (P(Y_{01}=1) P(Y_{10}=1))] =  eta_3,

     specifies the dependency between the two responses. Here, the
     responses equal 1 for a success and a 0 for a failure, and the
     odds ratio is often written psi=p00 p11 / (p10 p01). We have p10 =
     p01 because of the exchangeability.

     The second linear/additive predictor models the phi parameter (see
     'zipoisson'). The third linear/additive predictor is the same as
     'binom2.or', viz., the log odds ratio.

     Suppose a dataset1 comes from a Poisson distribution that has been
     converted to presence/absence, and that both marginal
     probabilities are the same (exchangeable). Then
     'binom2.or("cloglog", exch=TRUE)' is appropriate. Now suppose a
     dataset2 comes from a _zero-inflated_ Poisson distribution. The
     first linear/additive predictor of 'zipebcom()' applied to
     dataset2 is the same as that of 'binom2.or("cloglog", exch=TRUE)'
     applied to dataset1. That is, the phi has been taken care of by
     'zipebcom()' so that it is just like the simpler 'binom2.or'.

     Note that, for eta_1, 'mu12 = prob12 / (1-phi12)' where 'prob12'
     is the probability of a 1 under the ZIP model. Here, 'mu12'
     correspond to 'mu1' and 'mu2' in the 'binom2.or'-Poisson model.

     If phi=0 then 'zipebcom()' should be equivalent to
     'binom2.or("cloglog", exch=TRUE)'. Full details are given in Yee
     and Dirnbock (2009).

     The leading 2 x 2 submatrix of the expected information matrix
     (EIM) is of rank-1, not 2! This is due to the fact that the
     parameters corresponding to the first two linear/additive
     predictors are unidentifiable. The quick fix around this problem
     is to use the 'addRidge' adjustment. The model is fitted by
     maximum likelihood estimation since the full likelihood is
     specified. Fisher scoring is implemented.

     The default models eta2 and eta3 as single parameters only, but
     this can be circumvented by setting 'zero=NULL' in order to model
     the  phi and odds ratio as a function of all the explanatory
     variables.

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

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

     When fitted, the 'fitted.values' slot of the object contains the
     four joint probabilities, labelled as (Y1,Y2) = (0,0), (0,1),
     (1,0), (1,1), respectively. These estimated probabilities should
     be extracted with the 'fitted' generic function.

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

     The fact that the EIM is not of full rank may mean the model is
     naturally ill-conditioned. Not sure whether there are any negative
     consequences wrt theory. For now it is certainly safer to fit
     'binom2.or' to bivariate binary responses.

_N_o_t_e:

     The '"12"' in the argument names reinforce the user about the
     exchangeability assumption. The name of this 'VGAM' family
     function stands for _zero-inflated Poisson exchangeable bivariate
     complementary log-log odds-ratio model_ or ZIP-EBCOM.

     See 'binom2.or' for details that are pertinent to this 'VGAM'
     family function too. Even better initial values are usually needed
     here.

     The 'xij' (see 'vglm.control') argument enables environmental
     variables with different values at the two time points to be
     entered into an exchangeable 'binom2.or' model. See the author's
     webpage for sample code.

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

     Thomas W. Yee

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

     Yee, T. W. and Dirnbock, T. (2009) Models for analysing species'
     presence/absence data at two time points. Journal of Theoretical
     Biology, *259*. In press.

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

     'binom2.or', 'zipoisson', 'cloglog', 'CommonVGAMffArguments'.

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

     mydat = data.frame(x = seq(0, 1, len=(nsites <- 2000)))
     mydat = transform(mydat, eta1 =  -3 + 5 * x,
                              phi1 = logit(-1, inverse=TRUE),
                              oratio = exp(2))
     mydat = transform(mydat, mu12 = cloglog(eta1, inverse=TRUE) * (1-phi1))
     tmat =  with(mydat, rbinom2.or(nsites, mu1=mu12, oratio=oratio, exch=TRUE))
     mydat = transform(mydat, ybin1 = tmat[,1], ybin2 = tmat[,2])

     with(mydat, table(ybin1,ybin2)) / nsites  # For interest only
     ## Not run: 
     # Various plots of the data, for interest only
     par(mfrow=c(2,2))
     with(mydat, plot(jitter(ybin1) ~ x, col="blue"))

     with(mydat, plot(jitter(ybin2) ~ jitter(ybin1), col="blue"))

     with(mydat, plot(x, mu12, col="blue", type="l", ylim=0:1,
          ylab="Probability", main="Marginal probability and phi"))
     with(mydat, abline(h=phi1[1], col="red", lty="dashed"))

     tmat2 = with(mydat, dbinom2.or(mu1=mu12, oratio=oratio, exch=TRUE))
     with(mydat, matplot(x, tmat2, col=1:4, type="l", ylim=0:1,
          ylab="Probability", main="Joint probabilities"))
     ## End(Not run)

     # Now fit the model to the data.
     fit = vglm(cbind(ybin1,ybin2) ~ x, fam=zipebcom, dat=mydat, trace=TRUE)

     coef(fit, matrix=TRUE)
     summary(fit)
     vcov(fit)

