libname ps "C:\Documents and Settings\tneilands\Desktop\Misc urgent\PS" ; ** --> Compute derived variables <-- ** ; data temp ; set ps.hlp_baseline_n_3818_2_2_10 ; glmadh = 101 - pctadh1 ; glmadh2 = glmadh / 100 ; glmadh3 = glmadh - 1 ; bdisleep = bdi161 ; * Dummy variables for NLMIXED ; bdisleep0 = . ; bdisleep1 = . ; bdisleep2 = . ; bdisleep3 = . ; if bdisleep in (1, 2, 3) then bdisleep0 = 0 ; if bdisleep in (0) then bdisleep0 = 1 ; if bdisleep in (.) then bdisleep0 = . ; if bdisleep in (0, 2, 3) then bdisleep1 = 0 ; if bdisleep in (1) then bdisleep1 = 1 ; if bdisleep in (.) then bdisleep1 = . ; if bdisleep in (0, 1, 3) then bdisleep2 = 0 ; if bdisleep in (2) then bdisleep2 = 1 ; if bdisleep in (.) then bdisleep2 = . ; if bdisleep in (0, 1, 2) then bdisleep3 = 0 ; if bdisleep in (3) then bdisleep3 = 1 ; if bdisleep in (.) then bdisleep3 = . ; * Binary outcome for PROC LOGISTIC ; y_binary = . ; if glmadh = 1 then y_binary = 0 ; if glmadh > 1 then y_binary = 1 ; run ; proc glimmix data = temp ; title "GLIMMIX gamma analysis of HLP data" ; model glmadh = bdisleep1 bdisleep2 bdisleep3 / solution dist = gamma link = log ; estimate 'Exp(bdisleep1)' bdisleep1 1 / exp cl ; estimate 'Exp(bdisleep2)' bdisleep2 1 / exp cl ; estimate 'Exp(bdisleep3)' bdisleep3 1 / exp cl ; run ; ** --> Fit the gamma component of the ZIG <-- ** ; proc glimmix data = temp ; title "GLIMMIX gamma analysis of HLP data for GLMADH > 0" ; where glmadh > 1 ; y = glmadh - 1 ; model y = bdisleep1 bdisleep2 bdisleep3 / solution dist = gamma link = log ; estimate 'Exp(bdisleep1)' bdisleep1 1 / exp cl ; estimate 'Exp(bdisleep2)' bdisleep2 1 / exp cl ; estimate 'Exp(bdisleep3)' bdisleep3 1 / exp cl ; run ; ** --> Fit the logistic component with ZIG <-- ** ; proc logistic data = temp ; title "LOGISTIC analysis of HLP data for GLMADH 0 vs. 1" ; model y_binary(event='0') = bdisleep1 bdisleep2 bdisleep3 / clodds = both ; run ; ** --> Zero-inflated gamma model (ZIG) with no mixing <-- ** ; proc nlmixed data = temp qpoints = 15 lognote = 3 tech = newrap ; parms a0=0.1 a1=0.1 a2=0.1 a3=0.1 b0=0.4 b1=0.2 b2=0.4 b3=0.4 log_theta = 0; title "NLMIXED - Log-Gamma inflated model" ; title2 “Syntax is based on Dale McLerran’s posts on SAS-L” ; y = glmadh-1; linfp = a0 + a1*bdisleep1 + a2*bdisleep2 + a3*bdisleep3 ; infprob = exp(linfp)/(1+exp(linfp)) ; linp = b0 + b1*bdisleep1 + b2*bdisleep2 + b3*bdisleep3 ; mu = exp(linp); theta = exp(log_theta) ; r = mu/theta ; if y = 0 then ll = log(infprob) ; else ll = log(1-infprob) - lgamma(theta) + (theta-1)*log(y) - theta*log(r) - y/r ; model y ~ general(ll); run; data glimmix; length dist $7; set temp; response = (y_binary=1); dist = "Binary"; output; response = glmadh; dist = "Gamma"; output; keep subject bdisleep1 bdisleep2 bdisleep3 response dist; run ; data glimmix2 ; set glimmix ; if dist = "Gamma" and response = 0 then delete ; run ; proc glimmix data=glimmix2 ; class dist ; model response(event='1') = bdisleep1(dist) bdisleep2(dist) bdisleep3(dist) / solution dist=byobs(dist); run; proc glimmix data=glimmix2 method = quad (qpoints = 1) ; class dist ; model response(event='1') = bdisleep1(dist) bdisleep2(dist) bdisleep3(dist) / solution dist=byobs(dist); random int / subject = subject ; nloptions tech = newrap ; run; proc glimmix data=glimmix2 ; class dist ; model response(event='1') = bdisleep1(dist) bdisleep2(dist) bdisleep3(dist) / solution dist=byobs(dist); random _residual_ / type = cs subject = subject ; nloptions tech = nrridg ; run; ** --> End program <-- ** ;