data igg;
  input IgG Age;
  Age2=age*age;
  datalines;
      1.5         .5
      2.7         .5
      1.9         .5
        4         .5
      1.9         .5
      4.4         .5
      1.5         .5
      2.2         .5
      1.6         .5
      4.7         .5
      1.6         .5
      1.4         .5
        3         .5
      2.5         .5
        1         .5
      4.3         .5
      4.7         .5
      1.7         .5
      1.9   .5833333
       .9   .5833333
      4.1   .5833333
      2.8   .5833333
      2.2   .5833333
      5.4   .6666667
      8.4   .6666667
        2   .6666667
      5.1   .6666667
      1.5   .6666667
      3.2   .6666667
      7.7        .75
      4.5        .75
      6.6        .75
      4.2        .75
      2.1        .75
      4.2        .75
      3.8        .75
      5.7   .8333333
        3   .8333333
      3.2   .9166667
      5.1   .9166667
      2.1   .9166667
      2.3   .9166667
      3.4   .9166667
      3.9   .9166667
      4.3   .9166667
      5.3   .9166667
      7.2   .9166667
      3.8   .9166667
      5.6   .9166667
      1.5          1
        7          1
      4.6          1
      3.7          1
      4.5          1
      4.5          1
        5          1
      5.5          1
      5.5          1
      3.2          1
      3.2          1
      2.2          1
      2.3          1
      3.8   1.083333
      3.5   1.083333
      5.8   1.083333
        4   1.083333
      2.6   1.083333
      5.1   1.166667
      4.4   1.166667
      3.1   1.166667
        5   1.166667
      1.4   1.166667
      6.7   1.166667
      5.3   1.166667
      1.7   1.166667
      6.6   1.166667
      2.9   1.166667
      6.1       1.25
        4       1.25
      5.5       1.25
      4.7       1.25
      6.1       1.25
        4       1.25
      7.4   1.333333
      4.7   1.333333
      3.9   1.333333
      4.5   1.333333
      5.1   1.416667
      3.4   1.416667
      3.5   1.416667
      3.7        1.5
      5.8        1.5
      4.1        1.5
      9.8        1.5
      2.8        1.5
      5.8   1.583333
        7   1.583333
      3.1   1.583333
      4.2   1.666667
      5.4   1.666667
      5.7   1.666667
      4.4   1.666667
      5.8       1.75
      4.1       1.75
        4       1.75
      5.3       1.75
        5   1.833333
        6   1.833333
      5.7   1.833333
        7   1.833333
      2.5   1.833333
      6.8   1.833333
      6.3   1.916667
      5.3   1.916667
      4.7          2
        7          2
      4.2          2
      5.7          2
      3.4          2
      6.7          2
      4.6          2
      5.6          2
      1.8          2
      3.5          2
      4.3          2
      4.3          2
      5.4          2
      4.9          2
      5.4          2
      5.6          2
      3.8   2.083333
      7.3   2.083333
      3.4   2.083333
      4.4   2.083333
      3.7   2.083333
      3.3   2.083333
      3.5   2.083333
        5   2.083333
      2.7   2.083333
      4.4   2.166667
        8   2.166667
      6.2       2.25
      3.3   2.333333
      5.8   2.333333
      7.5   2.333333
      5.5   2.333333
      5.5        2.5
      6.1        2.5
      4.9   2.583333
      7.2   2.583333
      3.5   2.583333
      5.9   2.666667
      3.2       2.75
      6.1       2.75
      3.7       2.75
      7.3   2.833333
      3.3   2.916667
      1.8   2.916667
     10.4   2.916667
      4.2   2.916667
      4.2          3
      6.1          3
      7.8          3
      4.4          3
        5          3
      5.3          3
      5.6          3
      4.5          3
        6          3
      7.3   3.083333
      4.5   3.083333
      4.7   3.083333
      3.9   3.083333
        4   3.083333
      4.8   3.083333
      3.3   3.083333
      5.7   3.083333
      4.3   3.083333
        7   3.083333
     13.4   3.083333
        4   3.083333
      5.8   3.166667
      6.3   3.166667
      8.8   3.166667
      4.8   3.166667
      5.3       3.25
      4.6       3.25
      6.9       3.25
      5.7   3.333333
      6.5   3.416667
      6.3   3.416667
      6.8   3.416667
      3.9        3.5
      7.8   3.583333
        8   3.583333
      5.4   3.666667
      6.2   3.666667
      6.1   3.833333
      3.9   3.833333
        6   3.833333
      3.5   3.833333
      4.2   3.833333
      3.6   3.833333
      4.3   3.833333
      5.4   3.916667
      5.8          4
      7.5          4
      7.1          4
        6          4
      3.2   4.083333
      6.9   4.083333
      7.9   4.083333
      3.4   4.083333
      9.5   4.166667
      3.8   4.166667
      8.3   4.166667
      6.7   4.166667
      7.1       4.25
      7.8       4.25
      7.2       4.25
      6.6   4.333333
      2.5   4.333333
      2.1   4.333333
        4   4.333333
      3.7   4.333333
      5.6   4.333333
      5.6   4.333333
      4.5   4.416667
      5.9   4.416667
      6.3   4.416667
      4.7   4.416667
        8   4.416667
      4.6   4.416667
        4   4.416667
      1.9   4.416667
        4        4.5
      8.6        4.5
      2.6        4.5
      3.9        4.5
      6.4   4.583333
      7.8   4.583333
      3.8   4.583333
      5.5   4.583333
      7.1   4.583333
     10.2   4.666667
        7   4.666667
      7.4   4.666667
      9.4   4.666667
      6.8   4.666667
      9.1   4.666667
      5.2       4.75
      4.3       4.75
      2.7       4.75
       11       4.75
      9.6   4.833333
     12.6   4.833333
      8.9   4.833333
      3.8   4.916667
      6.1   5.083333
      7.5   5.083333
      7.4   5.083333
     10.9   5.083333
      8.3   5.166667
      8.2   5.166667
      4.7   5.166667
      8.1   5.166667
      9.3   5.166667
      4.4   5.166667
      8.7   5.166667
      9.8       5.25
      7.1       5.25
      8.1       5.25
      7.9       5.25
      8.4   5.333333
      8.2   5.333333
     10.4   5.333333
      9.7   5.333333
      8.1   5.416667
      4.8   5.416667
      4.9        5.5
     12.5        5.5
      3.8        5.5
      8.8        5.5
     10.4   5.583333
      4.7   5.583333
      3.3   5.583333
      5.6   5.583333
      4.6   5.583333
     14.4   5.583333
      9.1   5.666667
      6.3   5.666667
      6.1       5.75
      5.6       5.75
      5.4       5.75
        7   5.833333
      7.6   5.916667
      3.1          6
      6.8          6
;
run;


ods graphics on;
proc sgplot data=igg;
  histogram igg / nbins=50;
run;


proc sgplot data=igg;
  scatter x=Age y=IgG;
run;


proc mcmc data=igg
          seed=5263
          propcov=congra
          ntu=1000
          mintune=10
          nmc=30000;
  begincnst;
    p=0.95;
  endcnst;
  parms (b0-b2) 0;
  prior b: ~ general(0);
  mu= b0 + b1*age + b2*age2;
  u = igg - mu;
  ll = log(p)+log(1-p) - 0.5*(abs(u)+(2*p-1)*u);
  model igg ~ general(ll);
run;


data stackloss;
  input stackloss airflow watertemp acidconc;
  datalines;
42       80         27         89
37       80         27         88
37       75         25         90
28       62         24         87
18       62         22         87
18       62         23         87
19       62         24         93
20       62         24         93
15       58         23         87
14       58         18         80
14       58         18         89
13       58         17         88
11       58         18         82
12       58         19         93
 8       50         18         89
 7       50         18         86
 8       50         19         72
 8       50         19         79
 9       50         20         80
15       56         20         82
15       70         20         91
;
run;


proc sgplot data=stackloss;
  histogram stackloss / nbins=10;
  xaxis label="Stack Loss";
run;


proc sgplot data=stackloss;
  scatter y=stackloss x=airflow;
run;


proc sgplot data=stackloss;
  scatter y=stackloss x=watertemp;
run;


proc sgplot data=stackloss;
  scatter y=stackloss x=acidconc;
run;


data by_stackloss;
  set stackloss;
  do p = .05 to .95 by .05;
    output;
  end;
run;


proc sort data=by_stackloss;
  by p;
run;


ods output postsummaries=by_ps postintervals=by_pi;
proc mcmc data=by_stackloss
          seed=73625
          propcov=congra
          ntu=1000
          nmc=30000
          mintune=17;
  by p;
  parms (b0-b3) 0;
  prior b: ~ general(0);
  mu= b0 + b1*airflow + b2*watertemp + b3*acidconc;
  u = stackloss - mu;
  ll = log(p)+log(1-p) - 0.5*(abs(u)+(2*p-1)*u);
  model stackloss ~ general(ll);
run;


data process;
  merge by_ps by_pi;
run;


proc sort data=process out=process;
  by parameter p;
run;


proc sgplot data=process(where=(parameter="b0"));
  title "Estimated Parameter by Quantile";
  title2 "With 95% HPD Interval";
  series x=p y=mean / markers legendlabel="Intercept (b0)";
  band x=p lower=hpdlower upper=hpdupper /
       transparency=.5 legendlabel="HPD Interval";
  yaxis label="Intercept (b0)";
  xaxis label="Quantile";
  refline 0 / axis=y;
run;


proc print data=process(where=(parameter="b0")) noobs;
  title "Estimated Parameter b0 by Quantile";
  title2 "with 95% HPD Interval";
  var p mean hpdlower hpdupper;
run;


proc sgplot data=process(where=(parameter="b1"));
  title "Estimated Parameter by Quantile";
  title2 "With 95% HPD Interval";
  series x=p y=mean / markers legendlabel="Airflow (b1)";
  band x=p lower=hpdlower upper=hpdupper /
       transparency=.5 legendlabel="HPD Interval";
  yaxis label="Airflow (b1)";
  xaxis label="Quantile";
  refline 0 / axis=y;
run;


proc print data=process(where=(parameter="b1")) noobs;
  title "Estimated Parameter b1 by Quantile";
  title2 "with 95% HPD Interval";
  var p mean hpdlower hpdupper;
run;


proc sgplot data=process(where=(parameter="b2"));
  title "Estimated Parameter by Quantile";
  title2 "With 95% HPD Interval";
  series x=p y=mean / markers legendlabel="Water Temperature (b2)";
  band x=p lower=hpdlower upper=hpdupper /
       transparency=.5 legendlabel="HPD Interval";
  yaxis label="Water Temperature (b2)";
  xaxis label="Quantile";
  refline 0 / axis=y;
run;


proc print data=process(where=(parameter="b2")) noobs;
  title "Estimated Parameter b2 by Quantile";
  title2 "with 95% HPD Interval";
  var p mean hpdlower hpdupper;
run;


proc sgplot data=process(where=(parameter="b3"));
  title "Estimated Parameter by Quantile";
  title2 "With 95% HPD Interval";
  series x=p y=mean / markers legendlabel="Acid Concentration (b3)";
  band x=p lower=hpdlower upper=hpdupper /
       transparency=.5 legendlabel="HPD Interval";
  yaxis label="Acid Concentration (b3)";
  xaxis label="Quantile";
  refline 0 / axis=y;
run;


proc print data=process(where=(parameter="b3")) noobs;
  title "Estimated Parameter b3 by Quantile";
  title2 "with 95% HPD Interval";
  var p mean hpdlower hpdupper;
run;