mix2poisson               package:VGAM               R Documentation

_M_i_x_t_u_r_e _o_f _T_w_o _P_o_i_s_s_o_n _D_i_s_t_r_i_b_u_t_i_o_n_s

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

     Estimates the three parameters of a mixture of two Poisson
     distributions by maximum likelihood estimation.

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

     mix2poisson(lphi = "logit", llambda = "loge",
                 ephi=list(), el1=list(), el2=list(),
                 iphi = 0.5, il1 = NULL, il2 = NULL,
                 qmu = c(0.2, 0.8), zero = 1)

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

    lphi: Link function for the parameter phi. See below for more
          details. See 'Links' for more choices.

 llambda: Link function applied to each lambda parameter. See 'Links'
          for more choices.

ephi, el1, el2: List. Extra argument for each of the links. See 'earg'
          in 'Links' for general information.

    iphi: Initial value for phi, whose value must lie between 0 and 1.

il1, il2: Optional initial value for lambda1 and lambda2. These values
          must be positive. The default is to compute initial values
          internally using the argument 'qmu'.

     qmu: Vector with two values giving the probabilities relating to
          the sample quantiles for obtaining initial values for lambda1
          and lambda2. The two values are fed in as the 'probs'
          argument into 'quantile'.

    zero: An integer specifying which linear/additive predictor is
          modelled as intercepts only.  If given, the value must be
          either 1 and/or 2 and/or 3, and the default is the first one
          only, meaning phi is a single parameter even when there are
          explanatory variables. Set 'zero=NULL' to model all
          linear/additive predictors as functions of the explanatory
          variables.

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

     The probability function can be loosely written as 

     P(Y=y) = phi * Poisson(lambda1) + (1-phi) * Poisson(lambda2)

     where phi is the probability an observation belongs to the first
     group, and y=0,1,2,.... The parameter phi satisfies 0 < phi < 1.
     The mean of Y is phi*lambda1 + (1-phi)*lambda2 and this is
     returned as the fitted values. By default, the three
     linear/additive predictors are (logit(phi), log(lambda1),
     log(lambda2))^T.

_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'.

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

     Numerical problems can occur. Half-stepping is not uncommon. If
     failure to converge occurs, try obtaining better initial values,
     e.g., by using 'iphi' and 'qmu' etc.

     This function uses a quasi-Newton update for the working weight
     matrices (BFGS variant). It builds up approximations to the weight
     matrices, and currently the code is not fully tested. In
     particular, results based on the weight matrices (e.g., from
     'vcov' and 'summary') may be quite incorrect, especially when the
     arguments 'weights' is used to input prior weights.

     This 'VGAM' family function should be used with caution.

_N_o_t_e:

     Fitting this model successfully to data can be difficult due to
     numerical problems and ill-conditioned data.  It pays to fit the
     model several times with different initial values, and check that
     the best fit looks reasonable. Plotting the results is
     recommended.  This function works better as lambda1 and lambda2
     become more different.

     Convergence is often slow, especially when the two component
     distributions are not well separated. The control argument 'maxit'
     should be set to a higher value, e.g., 200, and use 'trace=TRUE'
     to monitor convergence.

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

     T. W. Yee

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

     'rpois', 'mix2normal1'.

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

     n = 3000
     mu1 = exp(2.4) # also known as lambda1
     mu2 = exp(3.1)
     phi = 0.3
     y = ifelse(runif(n) < phi, rpois(n, mu1), rpois(n, mu2))

     fit = vglm(y ~ 1, mix2poisson, maxit=200) # good idea to have trace=TRUE
     coef(fit, matrix=TRUE)
     Coef(fit) # the estimates
     c(phi, mu1, mu2) # the truth

     ## Not run: 
     # Plot the results
     ty = table(y)
     plot(names(ty), ty, type="h", main="Red=estimate, blue=truth")
     abline(v=Coef(fit)[-1], lty=2, col="red")
     abline(v=c(mu1, mu2), lty=2, col="blue")
     ## End(Not run)

