proc import out=ache datafile="/home/grego1/STAT 705/AcheHunting2.txt" replace; run; data ache1; set ache; ltripday=log(tripdays); age=age-45; run; proc glimmix data=ache1; class pid; model kills = age age*age / dist=poisson link=log offset=ltripday solution; run; *Marginal model; *Exchangeable random effect. Marginal model--should be same model in glimmix and genmod, though genmmod uses MOM and glimmix uses MLE; *to estimate variance components. Results are not as significant as for unconstrained case; proc genmod data=ache1; class pid; model kills=age age*age / dist=poisson link=log offset=ltripday; repeated subject=pid / type=exch; ods output GEEEmpPEst=genmod1; run; *Marginal model in GLIMMIX; proc glimmix data=ache1; class pid; model kills = age age*age / dist=poi link=log offset=ltripday solution; *R side random effect; random _residual_ / subject=pid type=cs; *Here is another way to get the R side random effect; *random intercept / subject=pid type=cs residual; ods output ParameterEstimates=glimmix2; run; *Conditional model in PROC NLMIXED; proc sort data=ache1; by pid; run; * need to sort by subject!; proc nlmixed qpoints=100 data=ache1; parms b1=-2.3 b2=0.0251 b3=-0.002 v=1.0; eta=b1+b2*age+b3*age**2+u+ltripday; lambda=exp(eta); model kills ~ poisson(lambda); random u ~ normal(0,v) subject=pid; ods output ParameterEstimates=nlmixed1; run; *Conditional model; *unrestricted covariance matrix for random effect; proc glimmix data=ache1; class pid; model kills = age age*age / dist=poisson link=log offset=ltripday solution; random intercept / subject=pid; *Here is another way to specify random intercept model; *random pid; ods output ParameterEstimates=glimmix1; run; data nlmixed1; set nlmixed1; if Parameter ne 'v'; if Parameter='b1' then Parameter='Intercept'; else if Parameter='b2' then Parameter='Age'; else Parameter='Age*Age'; run; *Conditional model with random class effect--exact same fit as random intercept; proc glimmix data=ache1; class pid; model kills = age age*age / dist=poi link=log offset=ltripday solution; random pid; ods output ParameterEstimates=glimmix4; run; *Let's look at parameter estimates in a tabular format--it takes some work; *There are other approaches (PROC APPEND) that didn't maintain all data features; data parms; set glimmix1 (in=Glimmix1 rename=(tValue=TestStat Probt=PValue) drop=DF) genmod1 (in=genmod1 rename=(z=TestStat ProbZ=PValue Parm=Effect) drop=LowerCL UpperCL) glimmix2 (in=Glimmix2 rename=(tValue=TestStat Probt=PValue) drop=DF) nlmixed1 (in=nlmixed1 rename=(Parameter=Effect tValue=TestStat Probt=PValue StandardError=StdErr) drop=DF Lower Upper alpha gradient) glimmix4 (in=glimmix4 rename=(tValue=TestStat Probt=PValue) drop=DF); if Glimmix1=1 then method="GLIMMIX Conditional Unrestricted"; else if Genmod1=1 then method="GENMOD Marginal Exchangeable"; else if Glimmix2=1 then method="GLIMMIX Marginal Exchangeable"; else if nlmixed1=1 then method="NLMIXED Conditional Exchangeable"; else method="GLIMMIX Conditional Random PID"; run; proc sort data=parms; by Effect; run; proc transpose data=parms out=tparms (drop=_LABEL_); by effect; var Estimate TestStat PValue Stderr; id method; run; proc print data=tparms label noobs; label _NAME_="Summary Measure"; run;