Date:         Tue, 8 Jul 2008 17:35:28 -0400
Reply-To:     Kevin Viel <citam.sasl@GMAIL.COM>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         Kevin Viel <citam.sasl@GMAIL.COM>
Subject:      Simulated ZINB and NLMIXED code

I have been playing with some code to simulate simple zero-inflated
negative binomial (ZINB) data.  I have the following treat:

data ZINB / view = ZINB ;

  call streaminit( 123 ) ;
  seed = 321 ;

  do True_p_logit = 0.0 , 0.05 to 0.95 by 0.05 ;
    do True_k = 1 , 3 , 5, 10 ;
      do True_p_NB = 0.05 to 0.95 by 0.05 ;
        do Repetition = 1 to 2 ;
          do _n_ = 1 to 1000 ;

            /* True_p_logit = 0.0 will be NB not ZINB */
            if True_p_logit ne 0.0 then
              call ranbin( seed , 1 , True_p_logit , zero_inflate ) ;

            if zero_inflate and True_p_logit ne 0.0 then y = 0 ;
             else y = RAND( "NEGBINOMIAL" , True_p_NB , True_k ) ;
            output ;
          end ;
        end ;
      end ;
    end ;
  end ;

run ;


/*****************************************************************/
ods output ParameterEstimates = Viel.PE_NL ;
ods listing close ;

proc nlmixed data = ZINB ;

 /* Initial values */
 parms
       /* logit */
       b_p0               0

       /* NB */
       b0                 0
       k_est              1
       ;


  /* the logit model */
  eta_p =   b_p0
          ;

  p = exp( eta_p ) / ( 1 + exp( eta_p )) ;


  eta =   b0
        ;

  mu = exp( eta ) ;


/*                  gamma( y + 1 / k ) * ( k * mu ) ** y                */
/* Pr( y | mu ) = -----------------------------------------------       */
/*            y! * gamma( 1 / k ) * ( 1 + k * mu ) ** ( y + 1 / k )     */

  /* Negative Binomial */
*  loglike =   lgamma( y + 1 / k_est )
            + y * log( k_est * mu )
            - lgamma( y + 1 )
            - lgamma( 1 / k_est )
            - ( y + 1 / k_est ) * log( 1 + k_est * mu)
            ;

  /* Zero-inflated Negative Binomial */
  if y = 0 then loglike = log( p + ( 1 - p ) *
                              ( 1 + k_est * mu ) ** ( -1 / k_est )
                             ) ;

   else loglike =   log( 1 - p )
                  + lgamma( y + 1 / k_est )
                  + y * log( k_est * mu )
                  - lgamma( y + 1 )
                  - lgamma( 1 / k_est )
                  - ( y + 1 / k_est ) * log( 1 + k_est * mu)
                  ;

  model y ~ general( loglike ) ;

  by True_p_logit True_k True_p_NB Repetition ;

run ;
quit ;

ods output close ;
ods listing ;

data PE_NL_exp ( keep = True_p_logit True_k True_p_NB
                        Repetition true_mean est_: ) ;

  merge Viel.PE_NL ( where  = ( Parameter = "b0"     )
                   )
        Viel.PE_NL ( where  = ( Parameter = "k_est"  )
                     rename = ( estimate = est_k_NL  )
                   )
        Viel.PE_NL ( where  = ( Parameter = "b_p0"   )
                     rename = ( estimate = est_l_p )
                   )
        ;

  by True_p_logit True_k True_p_NB Repetition ;

  est_mean_NL = exp( estimate ) ;
  est_k_NL    = 1 / est_k_NL ;
  est_p_NL    = est_k_NL / ( est_k_NL + est_mean_NL ) ;
  est_log_p   = exp( est_l_p ) / ( 1 + exp( est_l_p )) ;

  true_mean = True_k * ( 1 - True_p_NB ) / True_p_NB ;

run ;

proc print data = PE_NL_exp ;
  var True_p_logit est_log_p True_k est_k_NL True_p_NB
      est_p_NL true_mean est_mean_NL Repetition ;
run ;


As I have not had time to run this code in full (with more repetitions), I
would appreciate any criticism.  So far to me, it looks pretty exciting and
I am looking forward to compare it to the NB model.

Also, how might one select an initial value for k?  I think I read one
suggestion to use data from a GENMOD run (NB regression).

And lastly, how does anyone have experience using NB regression with
rates?  I would appreciate any insight into the use of an offset in either
the logit or NB part.

Many thanks,

Kevin