LISTSERV at the University of Georgia
Menubar Imagemap
Home Browse Manage Request Manuals Register
Previous messageNext messagePrevious in topicNext in topicPrevious by same authorNext by same authorPrevious page (July 2008, week 2)Back to main SAS-L pageJoin or leave SAS-L (or change settings)ReplyPost a new messageSearchProportional fontNon-proportional font
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


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