Example 3 for PROC BGLIMM
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: bglmmex3 */
/* TITLE: Example 3 for PROC BGLIMM */
/* DESC: Poisson Regression with Random Effects */
/* REF: STATLIB */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: Generalized linear mixed models */
/* Poisson data */
/* PROCS: BGLIMM */
/* DATA: Scottish lip cancer data */
/* */
/* SUPPORT: Amy Shi */
/* REF: Clayton, D. and Kaldor, J. (1987) */
/* Empirical Bayes Estimates of Age-standardized */
/* Relative Risks for Use in Disease Mapping */
/* Biometrics, 43, 671-681 */
/****************************************************************/
data LipCancer;
input County Observed Expected Employment SMR;
if (Observed > 0) then ExpCount = 100*Observed/SMR;
else ExpCount = Expected;
x = Employment / 10;
LogN = log(ExpCount);
datalines;
1 9 1.4 16 652.2
2 39 8.7 16 450.3
3 11 3.0 10 361.8
4 9 2.5 24 355.7
5 15 4.3 10 352.1
6 8 2.4 24 333.3
7 26 8.1 10 320.6
8 7 2.3 7 304.3
9 6 2.0 7 303.0
10 20 6.6 16 301.7
11 13 4.4 7 295.5
12 5 1.8 16 279.3
13 3 1.1 10 277.8
14 8 3.3 24 241.7
15 17 7.8 7 216.8
16 9 4.6 16 197.8
17 2 1.1 10 186.9
18 7 4.2 7 167.5
19 9 5.5 7 162.7
20 7 4.4 10 157.7
21 16 10.5 7 153.0
22 31 22.7 16 136.7
23 11 8.8 10 125.4
24 7 5.6 7 124.6
25 19 15.5 1 122.8
26 15 12.5 1 120.1
27 7 6.0 7 115.9
28 10 9.0 7 111.6
29 16 14.4 10 111.3
30 11 10.2 10 107.8
31 5 4.8 7 105.3
32 3 2.9 24 104.2
33 7 7.0 10 99.6
34 8 8.5 7 93.8
35 11 12.3 7 89.3
36 9 10.1 0 89.1
37 11 12.7 10 86.8
38 8 9.4 1 85.6
39 6 7.2 16 83.3
40 4 5.3 0 75.9
41 10 18.8 1 53.3
42 8 15.8 16 50.7
43 2 4.3 16 46.3
44 6 14.6 0 41.0
45 19 50.7 1 37.5
46 3 8.2 7 36.6
47 2 5.6 1 35.8
48 3 9.3 1 32.1
49 28 88.7 0 31.6
50 6 19.6 1 30.6
51 1 3.4 1 29.1
52 1 3.6 0 27.6
53 1 5.7 1 17.4
54 1 7.0 1 14.2
55 0 4.2 16 0.0
56 0 1.8 10 0.0
;
proc bglimm data=LipCancer seed=10571042 nmc=10000
outpost=LipCancer_Out;
class County;
model Observed = x / dist=poisson offset=LogN;
random int / sub=County;
run;
data SMR_PRED;
array gamma[56] Intercept__County_1-Intercept__County_56;
array SMR_pred[56];
set LipCancer_Out;
do i = 1 to 56;
set LipCancer(rename=(x=data_x)) point=i;
SMR_pred[i] = 100 * exp(Intercept + x * data_x + gamma[i]);
end;
keep smr_pred:;
run;
%sumint(data=SMR_PRED, var=_numeric_, print=NO, out=SMR_SI)
data combine;
merge LipCancer SMR_SI;
run;
proc sgplot data=combine noautolegend aspect=1;
yaxis label="Predicted SMR" max=700;
xaxis label="Observed SMR" max=700;
text x=SMR y=mean text=employment;
lineparm x=0 y=0 slope=1;
run;
proc bglimm data=LipCancer seed=10571042 nmc=10000
outpost=LipCancer_fOut;
class County;
model Observed = x / dist=poisson offset=LogN;
run;
data SMR_PRED_f;
array SMR_pred[56];
set LipCancer_Out;
do i = 1 to 56;
set LipCancer(rename=(x=data_x)) point=i;
SMR_pred[i] = 100 * exp(Intercept + x * data_x);
end;
keep smr_pred:;
run;
%sumint(data=SMR_PRED_f, var=_numeric_, print=NO, out=SMR_SIf)
data combine_f;
merge LipCancer SMR_SIf;
run;
proc sgplot data=combine_f noautolegend aspect=1;
yaxis label="Predicted SMR" max=700;
xaxis label="Observed SMR" max=700;
text x=SMR y=mean text=employment;
lineparm x=0 y=0 slope=1;
run;
proc bglimm data=LipCancer seed=10571042 nmc=10000 DIC;
class County;
model Observed = x / dist=poisson offset=LogN;
random int / sub=County;
run;
proc bglimm data=LipCancer seed=10571042 nmc=10000 DIC;
class County;
model Observed = x / dist=poisson offset=LogN;
run;
proc bglimm data=LipCancer seed=10571042 nmc=10000;
class County Employment;
model Observed = x / dist=poisson offset=LogN;
random int / sub=Employment;
random int / sub=County;
run;
