Documentation Example 6 for PROC GLIMMIX
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: gmxex06 */
/* TITLE: Documentation Example 6 for PROC GLIMMIX */
/* Radial Smoothing of Repeated Measures Data */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: Generalized linear mixed models */
/* Mixed model smoothing */
/* Smooth temporal trends */
/* PROCS: GLIMMIX, SORT, SGPANEL */
/* DATA: Diggle, P.J., Liang, K.-Y., and Zeger, S.L. (1994) */
/* Analysis of Longitudinal Data */
/* Oxford, UK: Oxford University Press */
/* */
/* SUPPORT: Oliver Schabenberger */
/* REF */
/* MISC: */
/****************************************************************/
data times;
input time1-time23;
datalines;
122 150 166 179 219 247 276 296 324 354 380 445
478 508 536 569 599 627 655 668 723 751 781
;
data cows;
if _n_ = 1 then merge times;
array t{23} time1 - time23;
array w{23} weight1 - weight23;
input cow iron infection weight1-weight23 @@;
do i=1 to 23;
weight = w{i};
tpoint = (t{i}-t{1})/10;
output;
end;
keep cow iron infection tpoint weight;
datalines;
1 0 0 4.7 4.905 5.011 5.075 5.136 5.165 5.298 5.323
5.416 5.438 5.541 5.652 5.687 5.737 5.814 5.799
5.784 5.844 5.886 5.914 5.979 5.927 5.94
2 0 0 4.868 5.075 5.193 5.22 5.298 5.416 5.481 5.521
5.617 5.635 5.687 5.768 5.799 5.872 5.886 5.872
5.914 5.966 5.991 6.016 6.087 6.098 6.153
3 0 0 4.868 5.011 5.136 5.193 5.273 5.323 5.416 5.46
5.521 5.58 5.617 5.687 5.72 5.753 5.784 5.784
5.784 5.814 5.829 5.872 5.927 5.9 5.991
4 0 0 4.828 5.011 5.136 5.193 5.273 5.347 5.438 5.561
5.541 5.598 5.67 . 5.737 5.844 5.858 5.872
5.886 5.927 5.94 5.979 6.052 6.028 6.12
5 1 0 4.787 4.977 5.043 5.136 5.106 5.298 5.298 5.371
5.438 5.501 5.561 5.652 5.67 5.737 5.784 5.768
5.784 5.784 5.829 5.858 5.914 5.9 5.94
6 1 0 4.745 4.868 5.043 5.106 5.22 5.298 5.347 5.347
5.416 5.501 5.561 5.58 5.687 5.72 5.737 5.72
5.737 5.753 5.768 5.784 5.844 5.844 5.9
7 1 0 4.745 4.905 5.011 5.106 5.165 5.273 5.371 5.416
5.416 5.521 5.541 5.635 5.687 5.704 5.784 5.768
5.768 5.814 5.829 5.858 5.94 5.94 6.004
8 0 1 4.942 5.106 5.136 5.193 5.298 5.347 5.46 5.521
5.561 5.58 5.635 5.704 5.784 5.823 5.858 5.9
5.94 5.991 6.016 6.064 6.052 6.016 5.979
9 0 1 4.605 4.745 4.868 4.905 4.977 5.22 5.165 5.22
5.22 5.247 5.298 5.416 5.501 5.521 5.58 5.58
5.635 5.67 5.72 5.753 5.799 5.829 5.858
10 0 1 4.7 4.868 4.905 4.977 5.011 5.106 5.165 5.22
5.22 5.22 5.273 5.384 5.438 5.438 5.501 5.501
5.541 5.598 5.58 5.635 5.687 5.72 5.704
11 0 1 4.828 5.011 5.075 5.165 5.247 5.323 5.394 5.46
5.46 5.501 5.541 5.609 5.687 5.704 5.72 5.704
5.704 5.72 5.737 5.768 5.858 5.9 5.94
12 0 1 4.7 4.828 4.905 5.011 5.075 5.165 5.247 5.298
5.298 5.323 5.416 5.505 5.561 5.58 5.561 5.635
5.687 5.72 5.72 5.737 5.784 5.814 5.799
13 0 1 4.828 5.011 5.075 5.136 5.22 5.273 5.347 5.416
5.438 5.416 5.521 5.628 5.67 5.687 5.72 5.72
5.799 5.858 5.872 5.914 5.94 5.991 6.016
14 0 1 4.828 4.942 5.011 5.075 5.075 5.22 5.273 5.298
5.323 5.298 5.394 5.489 5.541 5.58 5.617 5.67
5.704 5.753 5.768 5.814 5.872 5.927 5.927
15 0 1 4.745 4.905 4.977 5.075 5.193 5.22 5.298 5.323
5.394 5.394 5.438 5.583 5.617 5.652 5.687 5.72
5.753 5.768 5.814 5.844 5.886 5.886 5.886
16 0 1 4.7 4.868 5.011 5.043 5.106 5.165 5.247 5.298
5.347 5.371 5.438 5.455 5.617 5.635 5.704 5.737
5.784 5.768 5.814 5.844 5.886 5.94 5.927
17 1 1 4.605 4.787 4.828 4.942 5.011 5.136 5.22 5.247
5.273 5.247 5.347 5.366 5.416 5.46 5.541 5.481
5.501 5.635 5.652 5.598 5.635 5.635 5.598
18 1 1 4.828 4.977 5.011 5.136 5.273 5.298 5.371 5.46
5.416 5.416 5.438 5.557 5.617 5.67 5.72 5.72
5.799 5.858 5.886 5.914 5.979 6.004 6.028
19 1 1 4.7 4.905 4.942 5.011 5.043 5.136 5.193 5.193
5.247 5.22 5.323 5.338 5.371 5.394 5.438 5.416
5.501 5.561 5.541 5.58 5.652 5.67 5.704
20 1 1 4.745 4.905 4.977 5.043 5.136 5.273 5.347 5.394
5.416 5.394 5.521 5.617 5.617 5.617 5.67 5.635
5.652 5.687 5.652 5.617 5.687 5.768 5.814
21 1 1 4.787 4.942 4.977 5.106 5.165 5.247 5.323 5.416
5.394 5.371 5.438 5.521 5.521 5.561 5.635 5.617
5.687 5.72 5.737 5.737 5.768 5.768 5.704
22 1 1 4.605 4.828 4.828 4.977 5.043 5.165 5.22 5.273
5.247 5.22 5.298 5.375 5.371 5.416 5.501 5.501
5.521 5.561 5.617 5.635 5.72 5.737 5.768
23 1 1 4.7 4.905 5.011 5.075 5.106 5.22 5.22 5.298
5.323 5.347 5.416 5.472 5.501 5.541 5.598 5.598
5.598 5.652 5.67 5.704 5.737 5.768 5.784
24 1 1 4.745 4.942 5.011 5.075 5.106 5.247 5.273 5.323
5.347 5.371 5.416 5.481 5.501 5.541 5.598 5.598
5.635 5.687 5.704 5.72 5.829 5.844 5.9
25 1 1 4.654 4.828 4.828 4.977 4.977 5.043 5.136 5.165
5.165 5.165 5.193 5.204 5.22 5.273 5.371 5.347
5.46 5.58 5.635 5.67 5.753 5.799 5.844
26 1 1 4.828 4.977 5.011 5.106 5.165 5.22 5.273 5.323
5.371 5.394 5.46 5.576 5.652 5.617 5.687 5.67
5.72 5.784 5.784 5.784 5.829 5.814 5.844
;
proc glimmix data=cows;
t2 = tpoint / 100;
class cow iron infection;
model weight = iron infection iron*infection tpoint;
random t2 / type=rsmooth subject=cow
knotmethod=kdtree(bucket=100 knotinfo);
output out=gmxout pred(blup)=pred;
nloptions tech=newrap;
run;
data plot;
set gmxout;
length group $26;
if (iron=0) and (infection=0) then group='Control Group (n=4)';
else if (iron=1) and (infection=0) then group='Iron - No Infection (n=3)';
else if (iron=0) and (infection=1) then group='No Iron - Infection (n=9)';
else group = 'Iron - Infection (n=10)';
run;
proc sort data=plot; by group cow;
run;
proc sgpanel data=plot noautolegend;
title 'Radial Smoothing With Cow-Specific Trends';
label tpoint='Time' weight='log(Weight)';
panelby group / columns=2 rows=2;
scatter x=tpoint y=weight;
series x=tpoint y=pred / group=cow lineattrs=GraphFit;
run;
ods select OptInfo FitStatistics CovParms;
proc glimmix data=cows;
t2 = tpoint / 100;
class cow iron infection;
model weight = iron infection iron*infection tpoint;
random t2 / type=rsmooth
subject=cow
group=iron*infection
knotmethod=kdtree(bucket=100);
nloptions tech=newrap;
run;
ods select FitStatistics CovParms;
proc glimmix data=cows;
t2 = tpoint / 100;
class cow iron infection;
model weight = iron infection iron*infection tpoint;
random t2 / type=rsmooth
subject=iron*infection
knotmethod=kdtree(bucket=100);
random intercept / subject=cow;
output out=gmxout pred(blup)=pred;
nloptions tech=newrap;
run;
data plot;
set gmxout;
length group $26;
if (iron=0) and (infection=0) then group='Control Group (n=4)';
else if (iron=1) and (infection=0) then group='Iron - No Infection (n=3)';
else if (iron=0) and (infection=1) then group='No Iron - Infection (n=9)';
else group = 'Iron - Infection (n=10)';
run;
proc sort data=plot; by group cow;
run;
proc sgpanel data=plot noautolegend;
title 'Group-Specific Smoothing and Random Cow Intercepts';
label tpoint='Time' weight='log(Weight)';
panelby group / columns=2 rows=2;
scatter x=tpoint y=weight;
series x=tpoint y=pred / group=cow lineattrs=GraphFit;
run;