## Documentation Example 6 for ODS Graphics

```/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: ODSGREX6                                            */
/*   TITLE: Documentation Example 6 for ODS Graphics            */
/* PRODUCT: STAT                                                */
/*  SYSTEM: ALL                                                 */
/*    KEYS: graphics, ods                                       */
/*   PROCS:                                                     */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswfk                UPDATE: July 25, 2018         */
/*     REF: ods graphics                                        */
/*    MISC:                                                     */
/*   NOTES: This sample provides DATA step and PROC code        */
/*   from the chapter "Statistical Graphics Using ODS."         */
/****************************************************************/

proc sgplot data=sashelp.heart noautolegend;
reg y=weight x=height / markerattrs=(size=3px) degree=3;
run;

ods graphics on / antialiasmax=6000;
proc sgplot data=sashelp.heart noautolegend;
pbspline y=weight x=height / markerattrs=(size=3px);
run;

ods graphics on / loessmaxobs=6000;
proc sgplot data=sashelp.heart noautolegend;
loess y=weight x=height / markerattrs=(size=3px);
run;

proc sgplot data=sashelp.heart noautolegend;
pbspline y=weight x=height / smooth=0 nknots=5 markerattrs=(size=3px);
run;

%let st = 5;
data a;
s1  = sqrt(&st);
s2  = &st ** (1 / 3);
inc = 0.01 * s2;
do x = -&st to &st by inc;
x1 = x;
x2 = x * x;
x3 = x * x * x;
k1 = -2; k1 = (x > k1) * (x - k1) ** 3;
k2 =  0; k2 = (x > k2) * (x - k2) ** 3;
k3 =  2; k3 = (x > k3) * (x - k3) ** 3;
p0 = -0.5 + 0.01 * x - 0.04 * x2 - 0.01 * x3;
p1 = p0 + 0.1 * k1;
p2 = p1 + - 0.5 * k2;
p3 = p2 + 1.5 * k3;
y  = p3;
g  = 1 + (x > -2) + (x > 0) + (x > 2);
if g le 1 then p1 = .;
if g le 2 then p2 = .;
if g le 3 then p3 = .;
output;
end;
run;

ods path work.templat(update) sashelp.tmplmst(read);
proc template;
define statgraph spplot;
%macro a(c); lineattrs=(thickness=2 color=&c) datatransparency=0 %mend;
begingraph / border=false designwidth=defaultdesignheight;
layout overlayequated / equatetype=square
%let a = offsetmin=0 offsetmax=0;
xaxisopts=(display=(line ticks tickvalues) &a)
yaxisopts=(display=(line ticks tickvalues) &a)
commonaxisopts=(viewmin=-5 viewmax=5
tickvaluelist=(-5 -4 -3 -2 -1 0 1 2 3 4 5));
seriesplot x=x y=y  / lineattrs=(color=yellow thickness=12) group=g;
seriesplot x=x y=p0 / %a(blue);
seriesplot x=x y=p1 / %a(red);
seriesplot x=x y=p2 / %a(green);
seriesplot x=x y=p3 / %a(orange);
referenceline x=-2;
referenceline x=0;
referenceline x=2;
referenceline y=0;
endlayout;
endgraph;
end;
quit;

proc sgrender data=a template=spplot; run;

proc template; delete spplot; quit;

proc sgplot data=sashelp.class(where=(sex='F'));
pbspline y=weight x=height;
run;

proc sgplot data=sashelp.class(where=(sex='F'));
pbspline y=weight x=height / smooth=2e5;
run;

proc sgplot data=sashelp.class(where=(sex='F'));
pbspline y=weight x=height / nknots=20;
run;

proc sgplot data=sashelp.gas;
reg y=nox x=eqratio / degree=3 group=fuel markerattrs=(size=3px) name='a';
keylegend 'a' / location=inside position=topright across=1;
run;

proc sgplot data=sashelp.gas;
pbspline y=nox x=eqratio / group=fuel markerattrs=(size=3px) name='a';
keylegend 'a' / location=inside position=topright across=1;
run;

proc sgplot data=sashelp.gas;
loess y=nox x=eqratio / group=fuel markerattrs=(size=3px) name='a';
keylegend 'a' / location=inside position=topright across=1;
run;

proc sgplot data=sashelp.gas;
pbspline y=nox x=eqratio / group=fuel smooth=0 nknots=5
markerattrs=(size=3px) name='a';
keylegend 'a' / location=inside position=topright across=1;
run;

proc sgplot data=sashelp.gas;
ods output sgplot=sg;
pbspline y=nox x=eqratio / group=fuel smooth=0 nknots=5
markerattrs=(size=3px) name='a';
keylegend 'a' / location=inside position=topright across=1;
run;

data subset(drop=SORT_FUEL_RETAIN_ALL_);
set sg;
Obs = _n_;
by PBSPLINE_EQRATIO_NOX_GROUP_S__GP fuel;
if _N_ gt 169 then do; fuel = '_'; eqratio = ._; nox = ._; end;
if first.fuel or last.fuel or first.PBSPLINE_EQRATIO_NOX_GROUP_S__GP or
last.PBSPLINE_EQRATIO_NOX_GROUP_S__GP or obs = 169 then output;
if lag(first.fuel) or lag(first.PBSPLINE_EQRATIO_NOX_GROUP_S__GP) then do;
call missing(of PBSPLI: Fuel EqRatio NOx obs);
if _n_ gt 169 then do; fuel = '_'; eqratio = ._; nox = ._; end;
output; output; output;
end;
run;

proc print noobs; id obs; run;

proc freq data=sashelp.gas(where=(n(eqratio, nox) eq 2));
tables fuel;
run;

proc transreg data=sashelp.gas nomiss solve ss2 plots=fit(nocli noclm);
ods select anova fitstatistics fitplot;
model identity(nox) = spline(eqratio / nknots=5 evenly after) |
class(fuel / zero=none);
run;

proc means data=sashelp.gas(where=(n(nox, eqratio))) noprint;
class fuel;
var eqratio;
output out=m(where=(_type_ eq 1 and trim(_stat_) in ('MIN', 'MAX')));
run;

proc transpose data=m out=m2(drop=_:);
by fuel;
id _stat_;
var eqratio;
run;

data gas(drop=min max);
if _n_ = 1 then do i = 1 to n;
set m2 nobs=n point=i;
if fuel ne '82rongas' then
do eqratio = min to max by (max - min) / 200; output; end;
end;
set sashelp.gas(where=(n(nox, eqratio)));
output;
run;

proc transreg data=gas solve ss2 plots(interpolate)=fit(nocli noclm);
ods select anova fitstatistics fitplot;
model identity(nox) = class(fuel / zero=none) |
spline(eqratio / nknots=5 after evenly);
run;

data gas(drop=min max);
if _n_ = 1 then do i = 1 to n;
set m2 nobs=n point=i;
do eqratio = min to max by (max - min) / 200; output; end;
end;
set sashelp.gas(where=(n(nox, eqratio)));
output;
run;

proc transreg data=gas solve ss2 plots(interpolate)=fit(nocli noclm);
ods select anova fitstatistics fitplot;
model identity(nox) = class(fuel / zero=none)
spline(eqratio / nknots=5 after evenly);
run;

proc transreg data=gas ss2 plots(interpolate)=fit(nocli noclm);
ods select anova fitstatistics fitplot;
model identity(nox) = class(fuel / zero=none) * pbspline(eqratio / after);
run;

proc transreg data=gas solve ss2 plots(interpolate)=fit(nocli noclm);
ods select anova fitstatistics fitplot;
model identity(nox) = class(fuel / zero=none) *
smooth(eqratio / after sm=60);
run;

data x;
do i = 1 to 100;
g = 1;
x = 10 * uniform(17);
y = x + 2 * sin(x) + normal(17);
output;
g = 2;
x = 10 * uniform(17);
y = 5 - x - 2 * cos(x) + normal(17);
output;
end;
run;

proc sgplot data=x;
title 'Penalized-B-Spline';
pbspline y=y x=x / group=g;
run;
title;

proc transreg data=x ss2 plots=fit(nocli noclm) maxiter=100;
model identity(y) = class(g / zero=none) | mspline(x / nknots=10);
run;

proc transreg data=x ss2 plots=fit(nocli noclm) solve;
model identity(y) = class(g / zero=none) | spline(x / nknots=10 degree=1);
run;

proc transreg data=x ss2 plots=fit(nocli noclm) maxiter=100;
model identity(y) = class(g / zero=none) | mspline(x / nknots=10 degree=1);
run;

proc transreg data=x ss2 plots=fit(nocli noclm) maxiter=100;
model identity(y) = class(g / zero=none) | mspline(x / nknots=10);
output out=msp p replace;
run;

proc sort data=msp; by g x; run;

proc sgplot data=msp;
title 'Transreg Output Data Set';
scatter y=y  x=x / group=g markerattrs=(size=3px);
series  y=py x=x / group=g name='a';
keylegend 'a' / location=inside position=topleft across=1;
run;
title;

proc orthoreg data=sashelp.heart;
effect spl = spline(height / knotmethod=equal(5));
model weight = spl;
effectplot / obs;
run;

proc orthoreg data=sashelp.gas;
effect spl = spline(eqratio / knotmethod=equal(3));
class fuel;
model nox = spl | fuel;
effectplot / obs extend=data;
ods output SliceFitPlot=sp;
run;

proc orthoreg data=sashelp.gas;
effect spl = spline(eqratio / naturalcubic knotmethod=equal(5));
class fuel;
model nox = spl | fuel;
effectplot / obs;
run;

data x2;
do x = -2 to 2 by 0.1;
y = x + sin(x) + normal(17);
output;
end;
run;

proc transreg data=x2 details ss2;
model identity(y) = bspline(x / nknots=3 evenly=3);
run;

proc transreg data=x2 details ss2;
model identity(y) = spline(x / nknots=3 evenly=3);
run;

%let k = -4.000000000001 -3.000000000001 -2.000000000001 -1 0 1
2.000000000001  3.000000000001  4.000000000001;
proc transreg data=x2 details ss2;
model identity(y) = bspline(x / knots=&k);
run;

proc orthoreg data=x2;
effect spl = spline(x / knotmethod=equal(3));
model  y = spl / noint;
run;

proc orthoreg data=x2;
effect spl = spline(x / knotmethod=listwithboundary(&k));
model  y = spl / noint;
run;

proc transreg data=x2 details ss2;
model identity(y) = bspline(x / nknots=3 evenly=3);
output out=b1(keep=x_:) replace;
run;

proc glimmix data=x2 outdesign=b2(drop=x y
rename=(%macro ren; %do i = 1 %to 7; _x&i=x_%eval(&i-1) %end; %mend; %ren));
effect spl = spline(x / knotmethod=equal(3));
model y = spl / noint;
run;

%let k = -4.000000000001 -3.000000000001 -2.000000000001 -1 0 1
2.000000000001  3.000000000001  4.000000000001;

proc iml;
use x2(keep=x); read all into x;
b = bspline(x, 3, {&k});
vname = 'x_0' : rowcatc('x_' || char(ncol(b)-1));
create b3 from b[colname=vname]; append from b;
quit;

data b4(keep=x:);
%let d   = 3;                        /* degree                  */
%let nkn = %eval(&d * 2 + 3);        /* total number of knots   */
%let nb  = %eval(&d + 1 + 3);        /* number of cols in basis */
array k[&nkn] _temporary_ (&k);      /* knots                   */
array b[&nb] x_0 - x_%eval(&nb - 1); /* basis                   */
array w[%eval(2 * &d)];              /* work                    */
set x2;
do i = 1 to &nb; b[i] = 0; end;

* find the index of first knot greater than current data value;
do ki = 1 to &nkn while(k[ki] le x); end;
kki = ki - &d - 1;

* make the basis;
b[1 + kki] = 1;
do j = 1 to &d;
w[&d + j] = k[ki + j - 1] - x;
w[j] = x - k[ki - j];
s = 0;
do i = 1 to j;
t = w[&d + i] + w[j + 1 - i];
if t ne 0.0 then t = b[i + kki] / t;
b[i + kki] = s + w[&d + i] * t;
s = w[j + 1 - i] * t;
end;
b[j + 1 + kki] = s;
end;
run;

options nolabel;
proc compare error note briefsummary criterion=1e-12
data=b1 compare=b2 method=relative(1);
run;

proc compare error note briefsummary criterion=1e-12
data=b1 compare=b3 method=relative(1);
run;

proc compare error note briefsummary criterion=1e-12
data=b1 compare=b4(drop=x) method=relative(1);
run;
options label;

proc sgplot data=b4;
%macro s; %do i = 0 %to &nb-1; series y=x_&i x=x; %end; %mend; %s
yaxis label='B-Spline Basis Functions';
run;

proc transreg data=x2 details ss2;
model identity(y) = pspline(x / name=(p) nknots=3);
output out=p1(keep=p_:) replace;
run;

data both; merge p1 b1; run;

proc print data=both noobs;
format _numeric_ bestd6.5 p_1 p_2 best6.;
run;

proc cancorr data=both;
var x:;
with p:;
run;

data class;
set sashelp.class(rename=(height=Height1)) nobs=n;
output;
if _n_ = n then do;
call missing(name, sex, height1);
do age = 10 to 17; output; end;
end;
run;

proc reg data=class;
model height1 = age;
output out=p1 p=p1y r=r1y;
run;

data class;
f = 1;
set sashelp.class(rename=(height=Height2)) nobs=n;
output;
if _n_ = n then do;
call missing(f, name, sex);
do age = 10 to 17; output; end;
end;
run;

proc reg data=class;
freq f;
model height2 = age;
output out=p2 p=p2y r=r2y;
run;

data all; merge p1 p2; run;

proc print data=all;
var f name sex age height: p: r:;
format p1y p2y r1y r2y 6.3;
run;

proc means min max data=sashelp.gas;
class fuel;
var eqratio;
output out=m(where=(_type_ eq 1 and _stat_ in ('MIN', 'MAX')));
run;

proc transpose data=m out=m2(drop=_:);
var eqratio;
by fuel;
id _stat_;
run;

data score(drop=min max);
set m2;
do eqratio = min to max by (max - min) / 200; output; end;
run;

data gas;
set sashelp.gas(where=(n(nox))) score;
run;

proc glimmix data=gas;
effect spl = spline(eqratio / naturalcubic knotmethod=equal(5));
class fuel;
model nox = spl | fuel;
output out=scored(where=(nmiss(nox))) pred=py;
run;

proc sgplot data=scored;
series y=py x=eqratio / group=fuel;
run;

proc glimmix data=sashelp.gas;
effect spl = spline(eqratio / naturalcubic knotmethod=equal(5));
class fuel;
model nox = spl | fuel;
store SplineModel;
run;

proc plm restore=SplineModel;
score data=score out=scored2 predicted=py;
run;

proc sgplot data=scored2;
series y=py x=eqratio / group=fuel;
run;

data Neuralgia;
input Treatment \$ Sex \$ Age Duration Pain \$ @@;
datalines;
P  F  68   1  No   B  M  74  16  No   P  F  67  30  No   P  M  66  26  Yes
B  F  67  28  No   B  F  77  16  No   A  F  71  12  No   B  F  72  50  No
B  F  76   9  Yes  A  M  71  17  Yes  A  F  63  27  No   A  F  69  18  Yes
B  F  66  12  No   A  M  62  42  No   P  F  64   1  Yes  A  F  64  17  No
P  M  74   4  No   A  F  72  25  No   P  M  70   1  Yes  B  M  66  19  No
B  M  59  29  No   A  F  64  30  No   A  M  70  28  No   A  M  69   1  No
B  F  78   1  No   P  M  83   1  Yes  B  F  69  42  No   B  M  75  30  Yes
P  M  77  29  Yes  P  F  79  20  Yes  A  M  70  12  No   A  F  69  12  No
B  F  65  14  No   B  M  70   1  No   B  M  67  23  No   A  M  76  25  Yes
P  M  78  12  Yes  B  M  77   1  Yes  B  F  69  24  No   P  M  66   4  Yes
P  F  65  29  No   P  M  60  26  Yes  A  M  78  15  Yes  B  M  75  21  Yes
A  F  67  11  No   P  F  72  27  No   P  F  70  13  Yes  A  M  75   6  Yes
B  F  65   7  No   P  F  68  27  Yes  P  M  68  11  Yes  P  M  67  17  Yes
B  M  70  22  No   A  M  65  15  No   P  F  67   1  Yes  A  M  67  10  No
P  F  72  11  Yes  A  F  74   1  No   B  M  80  21  Yes  A  F  69   3  No
;

proc logistic data=Neuralgia outdesign=coded;
class Treatment Sex;
effect Age2      = spline(age / degree=2 knotmethod=equal(0));
effect Duration2 = polynomial(duration / degree=2);
model Pain= Treatment Sex Treatment*Sex Age2 Duration2 / expb;
run;

proc print data=coded(obs=10); run;

proc adaptivereg data=sashelp.gas plots=all details=bases;
class fuel;
model nox = eqratio | fuel;
run;

```