```Date: Tue, 8 Jul 2008 17:35:28 -0400 Reply-To: Kevin Viel Sender: "SAS(r) Discussion" From: Kevin Viel 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 ```

Back to: Top of message | Previous page | Main SAS-L page