A User-Defined Trend Model
/*--------------------------------------------------------------
SAS Sample Library
Name: ssmex05.sas
Description: Example program from SAS/ETS User's Guide,
The SSM Procedure
Title: A User-Defined Trend Model
Product: SAS/ETS Software
Keys: State space models
PROC: SSM
Notes:
--------------------------------------------------------------*/
title;
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
;
data Cows;
set Cows;
ironInf = "No Iron and No Infection";
if iron=1 and infection=1 then ironInf = "Iron and Infection";
else if iron=1 and infection=0 then ironInf = "Iron and No Infection";
else if iron=0 and infection=1 then ironInf = "No Iron and Infection";
else ironInf = "No Iron and No Infection";
run;
proc sort data=Cows;
by tpoint ;
run;
proc ssm data=Cows;
where iron=1 and infection=1;
id tpoint;
parms var1 var2 / lower=(1.e-8 1.e-8);
array tMat{2,2};
tMat[1,1] = 1;
tMat[2,2] = 1;
tMat[1,2] = _ID_DELTA_;
array covMat{2,2};
covMat[1,1] = var1*_ID_DELTA_ + var2*_ID_DELTA_**3/3;
covMat[1,2] = var2*_ID_DELTA_**2/2;
covMat[2,1] = covMat[1,2] ;
covMat[2,2] = var2*_ID_DELTA_;
state harveyLL(2) T(g)=(tMat) cov(g)=(covMat) a1(2);
component trend = harveyLL[1];
component slope = harveyLL[2];
irregular wn;
model weight = trend wn;
output out=for;
run;
proc sgplot data=For;
title "Model Fit: Two-Parameter Polynomial Spline of Order 2";
series x=tpoint y=smoothed_trend;
scatter x=tpoint y=weight;
run;
proc sgplot data=For;
title "Smoothed Estimate of Growth-Rate";
series x=tpoint y=smoothed_slope;
run;
title;