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), nsimEIM=100, zero = 1)

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

    lphi: Link function for the parameter phi. 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'.

 nsimEIM: See 'CommonVGAMffArguments'.

    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. See 'CommonVGAMffArguments' for more information.

_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:

     This 'VGAM' family function requires care for a successful
     application. In particular, good initial values are required
     because of the presence of local solutions. Therefore running this
     function with several different combinations of arguments such as
     'iphi', 'il1', 'il2', 'qmu' is highly recommended. Graphical
     methods such as 'hist' can be used as an aid.

     With grouped data (i.e., using the 'weights' argument) one has to
     use a large value of 'nsimEIM'; see the example below.

_N_o_t_e:

     The response must be integer-valued since 'dpois' is invoked.

     Fitting this model successfully to data can be difficult due to
     local solutions 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. The default control argument 'trace=TRUE' is to
     encourage monitoring convergence.

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

     T. W. Yee

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

     'rpois', 'poissonff', 'mix2normal1'.

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

     # Example 1: simulated data
     n = 1000
     mu1 = exp(2.5) # also known as lambda1
     mu2 = exp(3)
     (phi = logit(-0.5, inverse=TRUE))
     y = ifelse(runif(n) < phi, rpois(n, mu1), rpois(n, mu2))
     fit = vglm(y ~ 1, mix2poisson)
     coef(fit, matrix=TRUE)

     # Compare the results with the truth
     round(rbind('Estimated'=Coef(fit), 'Truth'=c(phi, mu1, mu2)), dig=2)

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

     # Example 2: London Times data (Lange, 1997, p.31)
     deaths = 0:9
     freq = c(162, 267, 271, 185, 111, 61, 27, 8, 3, 1)
     y = rep(deaths, freq)

     # Usually this does not work well unless nsimEIM is large
     fit = vglm(deaths ~ 1, weight = freq,
                mix2poisson(iphi=0.3, il1=1, il2=2.5, nsimEIM=5000))

     # This works better in general
     fit = vglm(y ~ 1, mix2poisson(iphi=0.3, il1=1, il2=2.5))

     coef(fit, matrix=TRUE)
     Coef(fit)

