dm "output;clear;log;clear";
/*
* Repeated measures examples from:
* Moser, E.B. 2004. Repeated Measures Modelling with PROC MIXED.
* Proceedings of the 29th SAS Users Group International Conference,
* Montreal, Canada, May 9-12, 2004.
*/
Title1 "Moser - Repeated Measures Examples - SUGI29";
ODS Listing Close;
ODS Results Off;
ODS HTML File=".\SUGI29Repeated.html" GPath=".\SUGI29RepeatedGraphs";
GOptions Device=JPEG NoPrompt /*FText="Helvetica" FTitle="Helvetica" HText=1 HTitle=1*/;
/*
* Clean out the work environment.
*/
Proc DataSets Library=Work MemType=(Data) NoList Kill;
Run;
Quit;
/*
* EXAMPLE 1
*/
/*
* Table 5.1, p 103. Visual acuity and lens strength (taken from Crowder
* and Hand, 1990). Seven subjects had their response times measured when
* a light was flashed into each eye through lenses of powers 6/6, 6/18,
* 6/36, and 6/60 (a lens of power a/b means that the eye will perceive
* as being at a feet an object actually positioned at b feet). Measurements
* were made in milliseconds, and the question of interest was whether
* response time varied with lens strength (Everitt 1995:103).
*
* Everitt, B.S. 1996. Making Sense of Statistics in Psychology:
* A Second-Level Course. Oxford University Press, Oxford, 350pp.
*/
Title2 "Visual Acuity and Lens Strength";
Data Vision; /* From Everitt 1996:103 */
Input Subject @;
Do Eye="Left ","Right";
Do LensStrength="6/06","6/18","6/36","6/60";
Input ResponseTime @;
Output;
End;
End;
DataLines;
1 116 119 116 124 120 117 114 122
2 110 110 114 115 106 112 110 110
3 117 118 120 120 120 120 120 124
4 112 116 115 113 115 116 116 119
5 113 114 114 118 114 117 116 112
6 114 115 94 116 100 99 94 97
7 110 110 105 118 105 105 115 115
;
/*
* Same data set from Lindsey 1993, but the first value
* for subject 6 differs from the corresponding value
* above for subject 6. Pick whichever data set you prefer.
*/
Data Vision; /* From Lindsey 1993:86 */
Input Subject @;
Do Eye="Left ","Right";
Do LensStrength="6/06","6/18","6/36","6/60";
Input ResponseTime @;
Output;
End;
End;
DataLines;
1 116 119 116 124 120 117 114 122
2 110 110 114 115 106 112 110 110
3 117 118 120 120 120 120 120 124
4 112 116 115 113 115 116 116 119
5 113 114 114 118 114 117 116 112
6 119 115 94 116 100 99 94 97
7 110 110 105 118 105 105 115 115
;
Proc Print Data=Vision;
Run;
/*
* Look at boxplots of the data and check for outliers, etc.
*/
Proc Sort Data=Vision;
By Eye LensStrength Subject;
Run;
Proc Boxplot Data=Vision;
Plot ResponseTime*LensStrength (Eye) / BoxStyle=SchematicIdFar Notches;
Id Subject;
Run;
Proc Sort Data=Vision;
By LensStrength Eye Subject;
Run;
Proc Boxplot Data=Vision;
Plot ResponseTime*Eye (LensStrength) / BoxStyle=SchematicIdFar Notches;
Id Subject;
Run;
Proc Sort Data=Vision;
By Eye Subject LensStrength;
Run;
/*
* Construct a profile or spaggetti plot of the repeated measures.
* In this way each subject's data can be studied as a unit.
*/
%AnnoMac; /* bring in the annotate macros */
Data Anno;
%DclAnno;
%System(1,2,4);
Set Vision(Where=(LensStrength In ("6/06","6/60")));
X=(LensStrength="6/60")*100;
Y=ResponseTime;
Labl=Left(Put(Subject,3.));
%Label(X,Y,Labl,BLACK,0,0,1,SWISSXB,6);
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Vision;
By Eye;
Plot ResponseTime*LensStrength=Subject
/ VAxis=Axis1 HAxis=Axis2 NoLegend Annotate=Anno;
Axis1 Label=(A=90 "Response Time")
Order=(90 To 130 By 10);
Axis2 Label=("Lens Strength");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
/*
* These initial examples take the data set and just consider
* parts of it at a time to illustrate the basic underlying
* concepts. The first set of analyses focus on a paired
* experiment.
*/
Proc Sort Data=Vision;
By LensStrength Subject Eye;
Run;
Title3 "Paired Comparison of Eyes at 6/06 Lens Strength";
Title4 "Randomized Block Design (RBD) Model";
Proc Mixed Data=Vision;
Where LensStrength="6/06";
Class Subject Eye;
Model ResponseTime = Eye / Solution;
Random Subject / Solution V VCorr;
Estimate "Mean Difference" Eye 1 -1;
LSMeans Eye / CL;
Run;
Title4 "RBD Model But With Heterogeneous Error Variance";
Proc Mixed Data=Vision;
Where LensStrength="6/06";
Class Subject Eye;
Model ResponseTime = Eye / Solution DDFM=Satterthwaite;
Random Subject / Solution V VCorr;
Repeated / Group=Eye;
Estimate "Mean Difference" Eye 1 -1;
LSMeans Eye;
Run;
Title4 "Multivariate Response Model With Unstructured Error Covariance Matrix";
Proc Mixed Data=Vision;
Where LensStrength="6/06";
Class Subject Eye;
Model ResponseTime = Eye;
Repeated Eye / Subject=Subject Type=UN R RCorr;
LSMeans Eye / PDiff;
Run;
/*
* Now Consider just the Left eye responses.
* Thus, we generalize the paired experiment.
* Fit a randomized complete block design model.
*/
Proc Sort Data=Vision;
By Eye Subject LensStrength;
Run;
Title3 "Comparison of Lens Strengths For Left Eyes";
Title4 "Randomized Block Design (RBD) Model";
Proc Mixed Data=Vision;
Where Eye="Left ";
Class Subject LensStrength;
Model ResponseTime = LensStrength / OutP=Diagnostics;
Random Subject / V VCorr;
LSMeans LensStrength / PDiff;
Run;
/*
* Look at some error residuals diagnostics.
* You should try the Version 9.1+ PROC MIXED diagnostics and ODS Graphics.
*/
Symbol2 C=Black V=Dot H=0.75;
Proc Univariate Data=Diagnostics Normal;
Var Resid;
QQPlot / Normal(Mu=0 Sigma=Est);
Run;
/*
* Construct predicted profile plots.
*/
Proc Sort Data=Diagnostics;
By Subject LensStrength;
Run;
%AnnoMac; /* bring in the annotate macros */
Data Anno;
%DclAnno;
%System(1,2,4);
Set Diagnostics(Where=(LensStrength In ("6/06","6/60")));
X=(LensStrength="6/60")*100;
Y=Pred;
Labl=Left(Put(Subject,3.));
%Label(X,Y,Labl,BLACK,0,0,1,SWISSXB,6);
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Diagnostics;
Plot Pred*LensStrength=Subject
/ VAxis=Axis1 HAxis=Axis2 NoLegend Annotate=Anno;
Axis1 Label=(A=90 "Predicted Response Time")
Order=(90 To 130 By 10);
Axis2 Label=("Lens Strength");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
Title4 "Randomized Block Design (RBD) Heterogeneous Errors Model";
Proc Mixed Data=Vision;
Where Eye="Left "
/* And (Not (Subject=6 And ResponseTime=94)) */ ;
Class Subject LensStrength;
Model ResponseTime = LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Random Subject / V VCorr;
Repeated / Group=LensStrength;
LSMeans LensStrength / PDiff;
Run;
Title4 "Multivariate Response Model With Unstructured Error Covariance Matrix";
Proc Mixed Data=Vision;
Where Eye="Left "
/* And (Not (Subject=6 And ResponseTime=94)) */ ;
Class Subject LensStrength;
Model ResponseTime = LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Repeated LensStrength / Subject=Subject Type=UN HLM HLPS R RCorr;
LSMeans LensStrength / PDiff;
Run;
Title4 "Multivariate Response Model With Compound Symmetry Error Covariance Matrix";
Proc Mixed Data=Vision;
Where Eye="Left "
/* And (Not (Subject=6 And ResponseTime=94)) */ ;
Class Subject LensStrength;
Model ResponseTime = LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Repeated LensStrength / Subject=Subject Type=CS R RCorr;
LSMeans LensStrength / PDiff;
Run;
Title4 "Multivariate Response Model With Heterogeneous Variance Compound Symmetry";
Proc Mixed Data=Vision;
Where Eye="Left "
/* And (Not (Subject=6 And ResponseTime=94)) */ ;
Class Subject LensStrength;
Model ResponseTime = LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Repeated LensStrength / Subject=Subject Type=CSH R RCorr;
LSMeans LensStrength / PDiff;
Run;
/*
* Use the information from both eyes.
* Fit the strip-plot model that permits measurements made
* on same eye to be correlated differently from
* those made on different eyes of the same subject.
*
* This model would correspond to the experiment where both the
* left and right eyes would be measured at a given lens strength
* at the same (or nearly same) time (i.e., left & right are paired).
* The order of the lens strength tests, however, would be randomized
* separately for each subject. Use a random effect to aggregate the
* data according to each unit: Subject, Eye(Subject), LensStrength(Subject).
*/
Proc Sort Data=Vision;
By Subject Eye LensStrength;
Run;
Title3 "Lens Strength And Eyes Comparisons";
Title4 "Strip-plot Model (3 Random Effects)";
Proc Mixed Data=Vision;
* Where Not (Subject=6 And ResponseTime=94);
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Random Subject Eye(Subject) LensStrength(Subject) / V VCorr;
LSMeans Eye*LensStrength / PDiff;
Run;
Symbol2 C=Black V=Dot H=0.75;
Proc Univariate Data=Diagnostics Normal;
Var Resid;
QQPlot / Normal(Mu=0 Sigma=Est);
Run;
/*
* Construct predicted profile plots.
*/
Proc Sort Data=Diagnostics;
By Eye LensStrength;
Run;
%AnnoMac; /* bring in the annotate macros */
Data Anno;
%DclAnno;
%System(1,2,4);
Set Diagnostics(Where=(LensStrength In ("6/06","6/60")));
X=(LensStrength="6/60")*100;
Y=Pred;
Labl=Left(Put(Subject,3.));
%Label(X,Y,Labl,BLACK,0,0,1,SWISSXB,6);
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Diagnostics;
By Eye;
Plot Pred*LensStrength=Subject
/ VAxis=Axis1 HAxis=Axis2 NoLegend Annotate=Anno;
Axis1 Label=(A=90 "Predicted Response Time")
Order=(90 To 130 By 10);
Axis2 Label=("Lens Strength");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
/*
* The above model produces a very large V matrix. By
* use of the SUBJECT effect, the model can be reformulated
* and only a single subject need be printed to study the
* covariance matrix.
*/
Title4 "Re-formulated Strip-plot Model";
Proc Mixed Data=Vision;
* Where Not (Subject=6 And ResponseTime=94);
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Random Intercept Eye LensStrength / Subject=Subject V=1 VCorr=1;
LSMeans Eye*LensStrength / PDiff;
Run;
/*
* Now fit both factors using a correlated errors model.
*/
Proc Sort Data=Vision;
By Subject Eye LensStrength;
Run;
/*
* We might like to use a completely multivariate model.
* Let's see what happens. Notice that the data have been
* ordered by subject, then by the factors that define the
* different levels of the repeated measurement conditions.
*/
Title4 "Unstructured Correlated Errors Model";
Proc Mixed Data=Vision;
* Where (Not (Subject=6 And ResponseTime=94));
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Repeated Eye*LensStrength / Subject=Subject Type=UN R RCorr;
LSMeans Eye*LensStrength / PDiff Slice=(Eye LensStrength);
Run;
/*
* There are not enough subjects to fit the above UN model. Thus,
* we may try a reduced version of the covariance structure model.
*/
Title4 "Unstructured EYE, Compound Symmetry LENSSTRENGTH Covariance Structure";
Proc Mixed Data=Vision;
* Where (Not (Subject=6 And ResponseTime=94));
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Repeated Eye LensStrength / Subject=Subject Type=UN@CS R RCorr;
LSMeans Eye*LensStrength / PDiff Slice=(Eye LensStrength);
Run;
/*
* And a model just to show that the structure can be complicated by
* adding in random effects.
*/
Title4 "Random Subject, Unstructured EYE, Compound Symmetry LENSSTRENGTH Covariance Structure";
Proc Mixed Data=Vision;
* Where (Not (Subject=6 And ResponseTime=94));
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / OutP=Diagnostics DDFM=KenwardRoger;
Random Intercept / Subject=Subject V=1 VCorr=1;
Repeated Eye LensStrength / Subject=Subject Type=UN@CS R RCorr;
LSMeans Eye*LensStrength / PDiff Slice=(Eye LensStrength);
Run;
/*
* An alternative approach would be to treat the lens strength
* set of measurements as a cluster and model them with a
* random effect, but model the measurements on the pair of
* eyes as correlated errors. More subjects are needed for this.
*/
Proc Sort Data=Vision;
By Subject LensStrength Eye;
Run;
Title4 "Random Effects + Covariance Structure on Eyes";
Proc Mixed Data=Vision;
* Where (Not (Subject=6 And ResponseTime=94));
Class Subject Eye LensStrength;
Model ResponseTime = Eye | LensStrength / Solution OutP=Diagnostics DDFM=KenwardRoger;
Random Intercept LensStrength / Subject=Subject Solution G GCorr V=1 VCorr=1;
Repeated Eye / Subject=LensStrength(Subject) Type=CS R RCorr;
LSMeans Eye*LensStrength / PDiff Slice=(Eye LensStrength);
Run;
Proc DataSets Library=Work MemType=(Data) NoList Kill;
Run;
Quit;
/*
* EXAMPLE 2
*/
/*
* Body weights of pigs measured in 9 consecutive weeks of
* followup. Table 3.1, page 35 in Diggle, P. J., K. Liang,
* and S. L. Zeger. 1994. Analysis of Longitudinal Data.
* Oxford Science Publications, 253pp.
*/
Title2 "Pig Weight Growth";
Data Pigs;
Infile Datalines Missover;
Pig+1;
Do Week=1 To 9;
Input Weight @;
Output;
End;
Datalines;
24.0 32.0 39.0 42.5 48.0 54.5 61.0 65.0 72.0
22.5 30.5 40.5 45.0 51.0 58.5 64.0 72.0 78.0
22.5 28.0 36.5 41.0 47.5 55.0 61.0 68.0 76.0
24.0 31.5 39.5 44.5 51.0 56.0 59.5 64.0 67.0
24.5 31.5 37.0 42.5 48.0 54.0 58.0 63.0 65.5
23.0 30.0 35.5 41.0 48.0 51.5 56.5 63.5 69.5
22.5 28.5 36.0 43.5 47.0 53.5 59.5 67.5 73.5
23.5 30.5 38.0 41.0 48.5 55.0 59.5 66.5 73.0
20.0 27.5 33.0 39.0 43.5 49.0 54.5 59.5 65.0
25.5 32.5 39.5 47.0 53.0 58.5 63.0 69.5 76.0
24.5 31.0 40.5 46.0 51.5 57.0 62.5 69.5 76.0
24.0 29.0 39.0 44.0 50.5 57.0 61.5 68.0 73.5
23.5 30.5 36.5 42.0 47.0 55.0 59.0 65.5 73.0
21.5 30.5 37.0 42.5 48.0 52.5 58.5 63.0 69.5
25.0 32.0 38.5 44.0 51.0 59.0 66.0 75.5 86.0
21.5 28.5 34.0 39.5 45.0 51.0 58.0 64.5 72.5
31.0 38.0 48.0 54.0 60.0 62.0 66.5 75.5 84.0
27.5 32.5 36.0 43.0 49.5 52.5 56.0 61.0 64.0
30.0 37.0 45.0 51.0 58.0 63.0 67.5 74.5 81.0
26.0 32.0 40.5 45.5 52.5 55.5 62.5 69.5 74.0
26.0 32.5 39.5 44.0 48.0 54.5 58.0 66.0 73.0
28.5 35.5 41.5 47.5 54.0 59.5 63.5 71.0 78.5
26.5 34.5 42.0 48.5 55.5 62.0 68.0 76.5 85.0
27.5 33.5 41.0 45.0 50.5 56.0 62.5 71.0 78.0
22.5 27.0 33.5 38.5 41.0 49.0 56.0 64.0 68.0
22.0 26.5 32.5 38.5 43.5 50.5 56.5 63.5 68.5
23.5 29.0 35.5 40.0 45.0 50.0 56.5 63.0 67.5
22.5 29.5 36.5 42.0 45.0 55.0 61.0 68.0 72.0
27.5 34.5 42.0 47.5 53.0 63.0 72.0 79.0 85.5
23.5 28.0 33.0 37.0 38.5 48.0 52.5 62.0 64.5
24.5 30.0 38.5 42.0 47.5 54.0 62.5 71.5 77.0
24.5 31.5 40.5 46.5 51.5 61.5 68.5 77.5 84.5
24.5 32.0 39.0 45.0 51.0 55.5 61.5 69.0 75.5
24.0 32.5 40.0 48.0 54.5 61.5 68.0 74.5 81.0
24.0 31.5 38.5 44.0 51.5 57.5 64.0 72.5 79.0
24.5 32.5 39.5 44.5 52.5 56.5 62.0 67.5 72.5
24.5 32.0 38.5 44.0 50.0 56.0 63.5 69.5 76.0
25.5 33.0 41.5 47.0 55.5 60.5 66.5 77.0 82.0
25.5 32.0 39.0 45.5 51.0 57.5 63.5 72.0 78.5
25.0 31.0 36.5 43.0 50.5 55.0 62.5 69.0 75.5
26.5 30.5 33.0 39.0 43.5 49.5 56.5 61.0 65.0
24.0 32.0 39.0 44.5 50.0 56.0 63.0 67.5 74.0
24.5 31.0 37.5 43.5 48.0 56.0 62.5 66.5 70.5
27.0 34.5 42.0 48.5 53.0 60.0 67.0 73.0 76.0
31.0 39.0 47.5 51.0 57.0 64.0 71.0 77.0 80.5
27.0 33.5 40.0 46.5 53.0 60.0 66.5 72.5 80.0
29.5 37.0 46.0 52.5 60.0 67.5 76.0 81.5 88.0
28.5 36.0 42.5 49.0 55.0 63.5 72.0 78.5 85.5
;
Proc Sort Data=Pigs;
By Week;
Run;
Proc Boxplot Data=Pigs;
Plot Weight*Week / BoxStyle=SchematicIdFar Notches;
Id Pig;
Run;
Proc Means Data=Pigs N Mean Std Var;
Class Week;
Var Weight;
Output Out=Profiles(Where=(_TYPE_=1)) Mean=Mean Std=StandardDeviation Var=Variance;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Profiles;
Plot StandardDeviation*Mean=1 / VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Standard Deviation of Weights");
Axis2 Label=("Mean of Weights");
Symbol1 C=Black V=Dot I=RL L=1;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Profiles;
Plot Variance*Mean=1 / VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Variance of Weights");
Axis2 Label=("Mean of Weights");
Symbol1 C=Black V=Dot I=RL L=1;
Run;
Quit;
Proc Sort Data=Pigs;
By Pig Week;
Run;
/*
* Construct profile or spagetti plot
*/
%AnnoMac; /* bring in the annotate macros */
Data Anno;
%DclAnno;
%System(1,2,4);
Set Pigs(Where=(Week In (1,9)));
X=(Week=9)*100;
Y=Weight;
Labl=Left(Put(Pig,3.));
%Label(X,Y,Labl,BLACK,0,0,1,SWISSXB,6);
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Pigs;
Plot Weight*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend /*Annotate=Anno*/;
Axis1 Label=(A=90 "Pig Weight");
Axis2 Label=("Week");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
%Macro Fits(Model);
Data FitStatistics;
Length Structure $10 FitStatistic $4;
Structure="&Model";
Set FitStatistics;
FitStatistic=Compress(Substr(Descr,1,4));
Drop Descr;
Run;
Proc Append Base=PigFits Data=FitStatistics;
Run;
%Mend Fits;
Proc DataSets Library=Work NoList;
Delete PigFits FitStatistics;
Run;
Quit;
/*
* Correlated errors in time model using multivariate model.
*/
Title3 "Multivariate Correlated Errors Model";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week / DDFM=KenwardRoger OutP=UNDiagnostics;
Repeated Week / Subject=Pig Type=UN HLM HLPS R RCorr RI;
LSMeans Week;
ODS Output CovParms=CovarianceParameters FitStatistics=FitStatistics;
Run;
%Fits(UN);
Data Variances;
Set CovarianceParameters;
If Substr(CovParm,4,1) = Substr(CovParm,6,1);
Run;
Proc Print Data=Variances;
Run;
Title3 "Heterogeneous First-order Autoregressive Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week / DDFM=KenwardRoger OutP=ARH1Diagnostics;
Repeated Week / Subject=Pig Type=ARH(1) R RCorr RI;
LSMeans Week;
ODS Output CovParms=CovarianceParameters FitStatistics=FitStatistics;
Run;
%Fits(ARH1);
Title3 "Heterogeneous General-order Autoregressive Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week / DDFM=KenwardRoger OutP=TOEPHDiagnostics;
Repeated Week / Subject=Pig Type=TOEPH R RCorr RI;
LSMeans Week;
ODS Output CovParms=CovarianceParameters FitStatistics=FitStatistics;
Run;
%Fits(TOEPH);
Title3 "Heterogeneous First-order Antedependence";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week / DDFM=KenwardRoger OutP=ANTEDiagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) R RCorr RI;
LSMeans Week;
ODS Output CovParms=CovarianceParameters FitStatistics=FitStatistics;
Run;
%Fits(ANTE1);
Title3;
/*
* Now Print the Table of Model Fits
*/
Proc Transpose Data=PigFits Out=PigTable(Drop=_NAME_);
By Structure NotSorted;
Id FitStatistic;
Var Value;
Run;
Proc Sort Data=PigTable;
By AIC;
Run;
Title3 "Table of Model Fits Ordered By AIC";
Proc Print Data=PigTable NoObs Label;
Label _2R="-2RLogL";
Run;
Title3;
/*
* From the above table, the ANTE(1) structure fits
* best as per AIC. Now, simplify the fixed-effects
* components, but maintain the ANTE(1) structure.
*/
Title3 "ANTE(1) Covariance Structure Model";
/*
* Construct predicted profiles
*/
Proc Sort Data=ANTEDiagnostics;
By Pig Week;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=ANTEDiagnostics;
Plot Pred*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Predicted Pig Weight");
Axis2 Label=("Week");
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
%Macro Fits(Model);
Data FitStatistics;
Length Structure $10 FitStatistic $4;
Structure="&Model";
Set FitStatistics;
FitStatistic=Compress(Substr(Descr,1,4));
Drop Descr;
Run;
Proc Append Base=PigFits Data=FitStatistics;
Run;
%Mend Fits;
Proc DataSets Library=Work NoList;
Delete PigFits FitStatistics;
Run;
Quit;
/*
* Make sure that the data are properly ordered so that
* the covariance matrix can be correctly constructed.
*/
Proc Sort Data=Pigs;
By Pig Week;
Run;
Title4 "Linear Growth Model";
Proc Mixed Data=Pigs Update;
Model Weight = Week / DDFM=KenwardRoger OutP=LinDiagnostics;
Repeated / Subject=Pig Type=ANTE(1) R RCorr;
ODS Output FitStatistics=FitStatistics;
Run;
%Fits(ANTE1);
/*
* Alternatively, a WEEK1 variable could be constructed
* so that WEEK can be included in a class statement to
* obtain the correct covariance matrix.
*/
Data Pigs;
Set Pigs;
Week1=Week;
Week2=Week*Week;
Week3=Week*Week2;
Run;
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week1 / DDFM=KenwardRoger OutP=LinDiagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) R RCorr;
Run;
/*
* Overlay the fitted model on the individual
* pig profiles.
*/
%Annomac;
Data Anno;
%DCLAnno;
%System(2,2,4);
Set LinDiagnostics;
By Pig NotSorted;
If First.Pig Then
Do;
%Move(Week,Weight);
End;
Else
Do;
%Draw(Week,Weight,BLACK,1,1);
End;
Run;
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPlot Data=LinDiagnostics Annotate=Anno;
Where Pig=1;
Plot Pred*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Pig Weight")
Order=(20 To 90 By 10);
Axis2 Label=("Week") Minor=None;
Symbol1 C=RED V=None I=Join L=1 W=5 R=100;
Run;
Quit;
/*
* The fits below show how you might check for Lack Of Fit
* with polynomial models. Note that the covariance parameters
* are being held fixed for the reduced models.
*/
Title4 "Separate Week Means Model";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week / DDFM=KenwardRoger OutP=ANTEDiagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) R RCorr RI;
LSMeans Week;
ODS Output CovParms=CovParms;
Run;
Title4 "Linear Growth Model + LOF";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week1 Week / HTYPE=1 DDFM=KenwardRoger OutP=Diagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) R RCorr;
Parms / ParmsData=CovParms NoIter;
Run;
Title4 "Quadratic Growth Model + LOF";
Proc Mixed Data=Pigs Update;
Class Week;
Model Weight = Week1 Week2 Week / HTYPE=1 DDFM=KenwardRoger OutP=Diagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) R RCorr;
Parms / ParmsData=CovParms NoIter;
Run;
/*
* Looks like the quadratic term is not significant, so we'll stop here.
* You may choose to examine alternative models for growth.
*/
/*
* The profile plots suggests fitting separate regression
* lines to each pigs data. This can be accomplished
* by using a random coefficients regression model.
*/
Title3 "Random Coefficients (Intercept and Slope) Regression Model";
Proc Sort Data=Pigs;
By Pig Week;
Run;
Proc Mixed Data=Pigs Update;
Model Weight = Week / Solution OutPM=PAFit OutP=SSFit;
Random Intercept Week / Solution Subject=Pig Type=UN G GCorr V=1 VCorr=1;
Estimate "Week 1 Ave Weight" Intercept 1 Week 1;
Estimate "Week 2 Ave Weight" Intercept 1 Week 2;
Estimate "Week 3 Ave Weight" Intercept 1 Week 3;
Estimate "Week 4 Ave Weight" Intercept 1 Week 4;
Estimate "Week 5 Ave Weight" Intercept 1 Week 5;
Estimate "Week 6 Ave Weight" Intercept 1 Week 6;
Estimate "Week 7 Ave Weight" Intercept 1 Week 7;
Estimate "Week 8 Ave Weight" Intercept 1 Week 8;
Estimate "Week 9 Ave Weight" Intercept 1 Week 9;
ODS Output FitStatistics=FitStatistics;
Run;
%Fits(RC);
/*
* Construct predicted profiles overlaying the
* Population-averaged Fit with the Subject-specific
* fits.
*/
Proc Sort Data=PAFit;
By Pig Week;
Run;
Proc Sort Data=SSFit;
By Pig Week;
Run;
%Annomac;
Data Anno;
%DCLAnno;
%System(2,2,4);
Set PAFit(Where=(Pig=1));
By Pig NotSorted;
If First.Pig Then
Do;
%Move(Week,Pred);
End;
Else
Do;
%Draw(Week,Pred,RED,1,10);
End;
Run;
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPlot Data=SSFit Annotate=Anno;
Plot Pred*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Predicted Pig Weight")
Order=(20 To 90 By 10);
Axis2 Label=("Week");
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
/*
* Let's go back to the covariance structure models, but now
* consider a method to model repeated measures when time
* is measured on a continuous rather than discrete events
* time scale. Pretend for the moment that each pig could
* be weighed at any given time and for any number of
* measurements. We'll use a spatial structure to accomodate
* continuous time.
*/
Title3 "Continuous Time Covariance Structure -- First-order Autoregressive";
Proc Sort Data=Pigs;
By Pig Week;
Run;
Proc Mixed Data=Pigs Update;
Model Weight = Week / Solution OutPM=PAFit OutP=SSFit;
Random Week / Solution Subject=Pig V VCorr;
Repeated / Subject=Pig Type=SP(POW)(Week) R=1 RCorr=1;
Parms (0.2 To 0.8 By 0.1) (0.5 To 0.9 By 0.1) (4 To 8 By 2);
ODS Output FitStatistics=FitStatistics;
Run;
%Fits(POW+RC);
Title3;
Proc Transpose Data=PigFits Out=PigTable(Drop=_NAME_);
By Structure NotSorted;
Id FitStatistic;
Var Value;
Run;
Proc Sort Data=PigTable;
By AIC;
Run;
Title3 "Table of Model Fits - Ordered by AIC";
Proc Print Data=PigTable NoObs Label;
Label _2R="-2RLogL";
Run;
Title3;
Proc DataSets Library=Work MemType=(Data) NoList Kill;
Run;
Quit;
/*
* Example 2.5
*/
Title2 "Pig Weight Gains";
Data Pigs;
Infile Datalines Missover;
Pig+1;
Array Weights(9) Weight1-Weight9;
Input Weight1-Weight9;
Do Week=1 To 8;
Gain=Weights(Week+1)-Weights(Week);
Output;
End;
Drop Weight1-Weight9;
Datalines;
24.0 32.0 39.0 42.5 48.0 54.5 61.0 65.0 72.0
22.5 30.5 40.5 45.0 51.0 58.5 64.0 72.0 78.0
22.5 28.0 36.5 41.0 47.5 55.0 61.0 68.0 76.0
24.0 31.5 39.5 44.5 51.0 56.0 59.5 64.0 67.0
24.5 31.5 37.0 42.5 48.0 54.0 58.0 63.0 65.5
23.0 30.0 35.5 41.0 48.0 51.5 56.5 63.5 69.5
22.5 28.5 36.0 43.5 47.0 53.5 59.5 67.5 73.5
23.5 30.5 38.0 41.0 48.5 55.0 59.5 66.5 73.0
20.0 27.5 33.0 39.0 43.5 49.0 54.5 59.5 65.0
25.5 32.5 39.5 47.0 53.0 58.5 63.0 69.5 76.0
24.5 31.0 40.5 46.0 51.5 57.0 62.5 69.5 76.0
24.0 29.0 39.0 44.0 50.5 57.0 61.5 68.0 73.5
23.5 30.5 36.5 42.0 47.0 55.0 59.0 65.5 73.0
21.5 30.5 37.0 42.5 48.0 52.5 58.5 63.0 69.5
25.0 32.0 38.5 44.0 51.0 59.0 66.0 75.5 86.0
21.5 28.5 34.0 39.5 45.0 51.0 58.0 64.5 72.5
31.0 38.0 48.0 54.0 60.0 62.0 66.5 75.5 84.0
27.5 32.5 36.0 43.0 49.5 52.5 56.0 61.0 64.0
30.0 37.0 45.0 51.0 58.0 63.0 67.5 74.5 81.0
26.0 32.0 40.5 45.5 52.5 55.5 62.5 69.5 74.0
26.0 32.5 39.5 44.0 48.0 54.5 58.0 66.0 73.0
28.5 35.5 41.5 47.5 54.0 59.5 63.5 71.0 78.5
26.5 34.5 42.0 48.5 55.5 62.0 68.0 76.5 85.0
27.5 33.5 41.0 45.0 50.5 56.0 62.5 71.0 78.0
22.5 27.0 33.5 38.5 41.0 49.0 56.0 64.0 68.0
22.0 26.5 32.5 38.5 43.5 50.5 56.5 63.5 68.5
23.5 29.0 35.5 40.0 45.0 50.0 56.5 63.0 67.5
22.5 29.5 36.5 42.0 45.0 55.0 61.0 68.0 72.0
27.5 34.5 42.0 47.5 53.0 63.0 72.0 79.0 85.5
23.5 28.0 33.0 37.0 38.5 48.0 52.5 62.0 64.5
24.5 30.0 38.5 42.0 47.5 54.0 62.5 71.5 77.0
24.5 31.5 40.5 46.5 51.5 61.5 68.5 77.5 84.5
24.5 32.0 39.0 45.0 51.0 55.5 61.5 69.0 75.5
24.0 32.5 40.0 48.0 54.5 61.5 68.0 74.5 81.0
24.0 31.5 38.5 44.0 51.5 57.5 64.0 72.5 79.0
24.5 32.5 39.5 44.5 52.5 56.5 62.0 67.5 72.5
24.5 32.0 38.5 44.0 50.0 56.0 63.5 69.5 76.0
25.5 33.0 41.5 47.0 55.5 60.5 66.5 77.0 82.0
25.5 32.0 39.0 45.5 51.0 57.5 63.5 72.0 78.5
25.0 31.0 36.5 43.0 50.5 55.0 62.5 69.0 75.5
26.5 30.5 33.0 39.0 43.5 49.5 56.5 61.0 65.0
24.0 32.0 39.0 44.5 50.0 56.0 63.0 67.5 74.0
24.5 31.0 37.5 43.5 48.0 56.0 62.5 66.5 70.5
27.0 34.5 42.0 48.5 53.0 60.0 67.0 73.0 76.0
31.0 39.0 47.5 51.0 57.0 64.0 71.0 77.0 80.5
27.0 33.5 40.0 46.5 53.0 60.0 66.5 72.5 80.0
29.5 37.0 46.0 52.5 60.0 67.5 76.0 81.5 88.0
28.5 36.0 42.5 49.0 55.0 63.5 72.0 78.5 85.5
;
Proc Sort Data=Pigs;
By Week;
Run;
Proc Boxplot Data=Pigs;
Plot Gain*Week / BoxStyle=SchematicIdFar Notches;
Id Pig;
Run;
Proc Means Data=Pigs N Mean Std Var;
Class Week;
Var Gain;
Output Out=Profiles(Where=(_TYPE_=1)) Mean=Mean Std=StandardDeviation Var=Variance;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Profiles;
Plot StandardDeviation*Mean=1 / VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Standard Deviation of Weight Gains");
Axis2 Label=("Mean of Weight Gains");
Symbol1 C=Black V=Dot I=RL L=1;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Profiles;
Plot Variance*Mean=1 / VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Variance of Weight Gains");
Axis2 Label=("Mean of Weight Gains");
Symbol1 C=Black V=Dot I=RL L=1;
Run;
Quit;
Proc Sort Data=Pigs;
By Pig Week;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Pigs;
Plot Gain*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Pig Weight Gain");
Axis2 Label=("Week");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
/*
* Correlated errors in time model using multivariate model.
*/
Title3 "Multivariate Model Error Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Gain = Week / DDFM=KenwardRoger OutP=UNDiagnostics;
Repeated Week / Subject=Pig Type=UN HLM HLPS RCorr RI;
LSMeans Week;
ODS Output CovParms=CovarianceParameters;
Run;
Title3 "Heterogeneous First-order Autoregressive Error Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Gain = Week / DDFM=KenwardRoger OutP=ARH1Diagnostics;
Repeated Week / Subject=Pig Type=ARH(1) RCorr;
LSMeans Week;
Run;
Title3 "Homogeneous First-order Autoregressive Error Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Gain = Week / DDFM=KenwardRoger OutP=AR1Diagnostics;
Repeated Week / Subject=Pig Type=AR(1) RCorr;
LSMeans Week;
Run;
Title3 "Heterogeneous General-order Autoregressive Error Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Gain = Week / DDFM=KenwardRoger OutP=TOEPHDiagnostics;
Repeated Week / Subject=Pig Type=TOEPH RCorr;
LSMeans Week;
Run;
Title3 "Heterogeneous First-order Antedependence Error Structure";
Proc Mixed Data=Pigs Update;
Class Week;
Model Gain = Week / DDFM=KenwardRoger OutP=ANTEDiagnostics;
Repeated Week / Subject=Pig Type=ANTE(1) RCorr;
LSMeans Week;
Run;
Title3 "Random Intercepts Regression Model";
Proc Sort Data=Pigs;
By Pig Week;
Run;
Proc Mixed Data=Pigs Update;
Model Gain = Week / Solution OutPM=RCMeans OutP=RCDiagnostics;
Random Intercept / Solution Subject=Pig Type=UN V VCorr;
Estimate "Week 1 Ave Weight" Intercept 1 Week 1;
Estimate "Week 2 Ave Weight" Intercept 1 Week 2;
Estimate "Week 3 Ave Weight" Intercept 1 Week 3;
Estimate "Week 4 Ave Weight" Intercept 1 Week 4;
Estimate "Week 5 Ave Weight" Intercept 1 Week 5;
Estimate "Week 6 Ave Weight" Intercept 1 Week 6;
Estimate "Week 7 Ave Weight" Intercept 1 Week 7;
Estimate "Week 8 Ave Weight" Intercept 1 Week 8;
Estimate "Week 9 Ave Weight" Intercept 1 Week 9;
Run;
/*
* Construct predicted profiles
*/
Proc Sort Data=RCMeans;
By Pig Week;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=RCDiagnostics;
Plot Pred*Week=Pig / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Predicted Pig Weight Gain")
/*Order=(20 To 90 By 10)*/;
Axis2 Label=("Week");
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
Proc DataSets Library=Work MemType=(Data) NoList Kill;
Run;
Quit;
/*
* Example 3
*/
/*
Table A3: Evolution of AFCR (number of lymphocytes bearing a Fc
receptor in one mm3 of peripheral blood) for three treatments of multiple
sclerosis (Heitjan, 1991b). For each patient, times and AFCR are on pairs
of lines. Negative times are before treatment began.
A.3 Azams data -AFCR
*/
Title2 "Evolution of AFCR in MS Patients";
Data Azams;
Infile DataLines Missover;
Length Treatment $5;
Array Times (*) Time1-Time30;
Array AFCRs (*) AFCR1-AFCR30;
Do s=17,15,16;
Input Treatment;
Do i=1 To s;
Input Subject;
Do j=1 To 30;
Times(j)=.;
AFCRs(j)=.;
End;
Do t=1 To 0 By -1;
j=0;
Input v @;
Do While (v NE .);
j+1;
If t Then Times(j)=v;
Else AFCRs(j)=v;
Input v @;
End;
Input;
End;
Output;
Input;
End;
End;
Stop;
Drop s i j v t;
DataLines;
P+P
1
-27 -13 28 56 84 168 259 331 427 504 672 771 834 945 1008 1092 1289 1306 1351
561 334 157 374 191 465 125 212 232 177 98 207 127 202 143 174 216 245 237
2
-14 -6 58 253 358 508 574 672 855 924
429 587 446 269 131 50 145 634 273 144
3
-14 -7 0 56 84 168 336 420 504 672 756 840 924 1000 1135 1260 1280 1337
231 312 123 127 297 337 225 312 178 111 97 133 239 151 115 297 141 113
4
-14 -7 0 28 56 83 169 252 330 420 504 588 678 756 840 924 1008 1120 1176 1260 1298
351 369 177 107 273 263 266 221 246 226 323 157 260 234 283 211 210 174 335 216 90
5
-22 -8 0 27 83 167 405 504 588 672 756 861 924 1008 1092 1176 1197 1281 1309
347 77 120 379 330 620 798 386 211 176 506 239 312 235 206 430 308 251 151
6
-42 -7 0 28 84 246 343 392 483 644 819 896 1001 1091 1148 1232 1256 1284
23 126 108 350 222 216 312 154 148 132 217 151 102 304 85 462 141 107
7
-14 -7 0 28 56 84 168 261 392 477 561 644 736 812 897 980 1155 1184 1198 1232 1254
273 269 248 345 366 414 341 291 226 331 114 87 100 339 278 133 121 40 116 112 182
8
-7 0 48 85 168 252 335 426 504 588 693 797 923 986 1190
339 431 277 330 602 264 207 165 65 133 202 242 228 147 147
9
-77 -7 0 28 56 83 168 251 342 371 426 518 584 609 672 755 853 924 1008 1091 1176 1214 1317
187 429 465 504 187 199 444 153 489 248 286 120 206 77 43 243 212 171 137 23 283 143 122
10
-14 -7 28 56 112 175 251 336 419 510 595 672 756 924 979 1092 1176 1186 1214 1305
952 572 776 541 622 900 744 940 707 1175 367 353 704 449 459 701 178 644 150 1153
11
-14 -7 20 56 84 168 259 336 489 588 693 784 840 924 1029 1113 1134 1176
90 126 316 245 358 143 148 226 227 118 129 148 298 246 177 156 144 232
12
-14 -7 0 28 56 85 168 335 426 496 595 672 763 840 924 1092 1274
544 120 209 418 262 251 351 282 69 156 228 211 83 236 175 146 262
13
-21 -14 0 28 56 84 167 336 426 505 603 679 762 854 909
94 475 240 139 382 326 176 133 252 137 95 87 212 15 212
14
-80 -73 -31 -3 25 53 137 221 305 389 473 641 725 809 936
288 184 276 483 308 324 173 318 112 196 184 314 217 118 223
15
-7 0 84 252 329 427 497 580 664 756 960 1121 1130 1140 1175
141 149 250 168 189 111 109 162 112 175 89 110 218 127 71
16
-21 -7 0 28 167 252 335 447 504 587 672 840 1092 1098
426 418 404 297 287 60 175 197 289 194 482 315 261 149
17
-175 -168 -162 83 168 252 440 504 713 756 811 917 1056 1064 1137 1151
291 920 144 187 256 91 31 33 144 114 30 113 90 82 146 211
AZ+P
1
-15 -7 0 28 56 84 252 308 405 518 588 671 756 860 945 1008 1112 1196 1259 1281 1309 1333 1358 1400
537 341 328 168 208 228 48 123 19 34 41 18 12 80 14 60 70 11 40 25 138 102 136 88
2
-20 -13 -6 21 50 78 161 246 323 420 497 581 666 750 834 918 1002 1086 1170 1288 1302 1338 1365
272 207 245 151 123 142 12 26 42 46 34 17 10 13 20 77 25 16 32 26 50 238 47
3
-14 -7 0 28 56 84 168 252 323 419 514 588 700 791 868 945 1035 1134 1260 1288 1319 1340
158 114 78 64 65 63 67 44 23 81 76 18 52 32 70 12 10 12 57 42 143 41
4
-14 -7 0 28 56 84 175 252 343 357 504 588 770 847 924 1008 1273 1297 1344
259 213 96 60 48 82 36 58 31 27 54 22 37 80 31 78 58 38 52
5
-14 -7 0 28 84 169 336 420 504 588 673 756 840 924 995 1108 1183 1260 1299
221 130 207 176 380 139 130 33 21 20 115 40 14 11 0 90 18 25 117
6
-14 -7 168 252 336 420 504 658
316 229 145 37 104 16 34 60
7
-14 -7 0 84 168 239 336
442 288 148 201 223 50 95
8
-69 -62 21 84 168 239 420 504 609 778 869 938 1022 1156 1177 1206
748 454 438 230 617 295 195 154 138 87 134 20 60 30 57 34
9
-14 0 28 56 84 253 421 511 609 672 841 924 1008 1100 1176 1331
296 485 165 279 278 163 82 49 88 127 65 32 65 34 41 36
10
-14 -7 0 28 56 84 168 253 336 420 504 588 686 756 840 931 1009 1092 1148 1176 1208 1241
423 518 324 578 391 420 285 96 221 164 130 84 185 132 88 154 285 215 319 111 84 150
11
-20 0 56 86 168 267 336 421 511 672 785 841 980 1093 1176
438 435 204 98 149 78 27 16 17 35 168 23 89 29 64
12
-14 -7 0 56 84 174 258 342 420 504 588 672 756 924
190 250 306 591 263 41 102 23 32 39 20 16 17 26
13
-14 0 29 56 85 169 238 336 420 504 588 672 756 840 945 1029 1092 1120
266 347 316 221 200 40 30 49 22 28 54 46 55 79 48 84 51 57
14
-35 -21 56 84 168 252 343 420 505 603 680 848
224 428 134 171 47 258 100 178 159 203 211 117
15
-14 -7 0 28 84 168 252 336 504
216 388 253 88 139 139 52 78 108
AZ+MP
1
-125 -13 0 28 56 168 239
515 814 469 274 486 230 216
2
-14 28 56 84 168 252 343
267 71 73 80 238 0 0
3
-7 0 28 58 84 175 252 336 504 568 672 770 841 946 1022 1092 1296 1302 1373
408 613 374 273 48 358 84 29 120 293 503 533 205 298 153 131 129 173 304
4
-20 -13 10 28 56 84 168 252 420 504 673 756 840 924 1015 1120
408 263 119 515 177 194 351 112 113 47 40 18 14 60 80 13
5
-14 -7 0 84 156 252 338 428 581 672 840 1128 1191 1260
206 373 262 126 151 162 125 91 99 60 27 34 120 43
6
-14 -7 0 30 56 84 168 252 420 497 588 672 749 854 932 1029 1100 1177 1281 1306 1330
344 468 441 162 485 509 200 66 99 115 69 27 80 20 33 11 70 48 177 84 69
7
-22 -15 -8 28 56 84 162 252 321 419 503 587 671 728 755 853 937 1234 1241 1287 1315
378 349 209 150 470 116 94 37 80 125 56 82 29 65 135 72 18 80 38 59 96
8
-14 -6 0 28 56 84 169 252 316 420 504 581 674 688 875 1017 1094 1212 1239 1267
137 198 219 303 499 182 123 156 114 167 280 406 213 147 178 144 167 389 105 119
9
-14 -7 10 56 78 168 266 420 498 588 673 756 840 924 1043 1157 1352
350 576 364 268 438 296 128 98 17 287 52 74 114 119 123 89 71
10
-14 -7 0 28 56 84 168 259 365 449 533 617 701 1205
372 412 310 457 815 779 597 372 257 228 51 296 40 212
11
-21 -19 0 28 56 84 168 252 336 420 504 594 672 756 840 938 1008 1105 1176 1186 1225
480 626 184 290 462 176 317 109 140 91 62 61 50 106 19 60 13 22 56 133 143
12
-14 0 28 56 86 168 420 504 581 756
436 143 311 277 156 82 40 13 27 25
13
-42 -35 0 28 56 86 340
295 315 176 317 194 108 212
14
-14 -7 0 56 169 240 298 329
88 117 272 305 281 60 189 139
15
-28 0 56 168 252 336 427 518 588 671 756 833 973 1063 1071 1130 1176
228 201 375 143 50 22 31 37 50 34 21 29 37 33 59 105 290
16
-43 0 56 84 168 259 336 434 510 602 672 756 840 945 1091 1101 1119 1151
343 286 159 163 37 45 55 37 59 60 34 33 70 42 26 52 34 121
;
Proc Transpose Data=Azams Out=Times(Rename=(Col1=Time));
By Treatment Subject Notsorted;
Var Time1-Time25;
Run;
Proc Transpose Data=Azams Out=AFCRs(Rename=(Col1=AFCR));
By Treatment Subject Notsorted;
Var AFCR1-AFCR25;
Run;
Data Azams;
Merge Times AFCRs;
If Time=. Then Delete;
Drop _NAME_;
SqrtAFCR=Sqrt(AFCR);
LogAFCR=Log(AFCR+1);
TimeC=Time-501.5;
Run;
Proc Sort Data=Azams;
By Treatment Subject Time;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Azams;
By Treatment;
Plot AFCR*Time=Subject / VAxis=Axis1 HAxis=Axis2 NoLegend HRef=0 LHRef=2;
Axis1 Label=(A=90 "AFCR")
Order=(0 To 1200 By 100);
Axis2 Label=("Days From Start of Treatment")
Order=(-200 To 1400 By 200);
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Azams;
By Treatment;
Plot LogAFCR*Time=Subject / VAxis=Axis1 HAxis=Axis2 NoLegend HRef=0 LHRef=2;
Axis1 Label=(A=90 "Log(AFCR)")
Order=(0 To 8 By 1);
Axis2 Label=("Days From Start of Treatment")
Order=(-200 To 1400 By 200);
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=Azams;
By Treatment;
Plot SqrtAFCR*Time=Subject / VAxis=Axis1 HAxis=Axis2 NoLegend HRef=0 LHRef=2;
Axis1 Label=(A=90 "Sqrt(AFCR)")
/*Order=(0 To 8 By 1)*/;
Axis2 Label=("Days From Start of Treatment")
Order=(-200 To 1400 By 200);
Symbol1 C=Black V=None I=Join L=1 R=100;
Run;
Quit;
Title3 "Square Root Transform - Homogeneous Variance First-order Autoregressive Errors";
Proc Mixed Data=Azams;
Class Treatment Subject;
Model SqrtAFCR = TimeC Treatment Treatment*TimeC TimeC*TimeC Treatment*TimeC*TimeC
/ Solution DDFM=KenwardRoger;
Random Subject(Treatment) / Solution;
Repeated / Subject=Subject(Treatment) Type=SP(POW)(TimeC) R=1,20 RCorr=1,20;
Run;
Quit;
Title3 "Square Root Transform - Homogeneous Variance Independent Errors";
Proc Mixed Data=Azams;
Class Treatment Subject;
Model SqrtAFCR = TimeC Treatment Treatment*TimeC TimeC*TimeC Treatment*TimeC*TimeC
/ Solution DDFM=KenwardRoger;
Random Subject(Treatment) / Solution;
Repeated / Subject=Subject(Treatment) Type=Simple;
Run;
Quit;
Title3 "Log Transform - Homogeneous Variance First-order Autoregressive Errors";
Proc Mixed Data=Azams;
Class Treatment Subject;
Model LogAFCR = TimeC Treatment Treatment*TimeC TimeC*TimeC Treatment*TimeC*TimeC
/ Solution DDFM=KenwardRoger;
Random Subject(Treatment) / Solution;
Repeated / Subject=Subject(Treatment) Type=SP(POW)(TimeC) R=1,20 RCorr=1,20;
Run;
Title3 "Log Transform - Homogeneous Variance Independent Errors";
Proc Mixed Data=Azams;
Class Treatment Subject;
Model LogAFCR = TimeC Treatment Treatment*TimeC TimeC*TimeC Treatment*TimeC*TimeC
/ Solution DDFM=KenwardRoger;
Random Subject(Treatment) / Solution;
Repeated / Subject=Subject(Treatment) Type=Simple;
Run;
Proc DataSets Library=Work MemType=(Data) NoList Kill;
Run;
Quit;
/*
* Example 4
*/
Title2 "Epileptic Fits Data";
Data EpilepticFits;
Subject+1;
Input @17 Treatment Baseline Age @1 @;
Do Period= 1 To 4;
Input NumberOfFits @;
Output;
End;
DataLines;
5 3 3 3 0 11 31
3 5 3 3 0 11 30
2 4 0 5 0 6 25
4 4 1 4 0 8 36
7 18 9 21 0 66 22
5 2 8 7 0 27 29
6 4 0 2 0 12 31
40 20 23 12 0 52 42
5 6 6 5 0 23 37
14 13 6 0 0 10 28
26 12 6 22 0 52 36
12 6 8 5 0 33 24
4 4 6 2 0 18 23
7 9 12 14 0 42 36
16 24 10 9 0 87 26
11 0 0 5 0 50 26
0 0 3 3 0 18 28
37 29 28 29 0 111 31
3 5 2 5 0 18 32
3 0 6 7 0 20 21
3 4 3 4 0 12 29
3 4 3 4 0 9 21
2 3 3 5 0 17 32
8 12 2 8 0 28 25
18 24 76 25 0 55 30
2 1 2 1 0 9 40
3 1 4 2 0 10 19
13 15 13 12 0 47 22
11 14 9 8 1 76 18
8 7 9 4 1 38 32
0 4 3 0 1 19 20
3 6 1 3 1 10 20
2 6 7 4 1 19 18
4 3 1 3 1 24 24
22 17 19 16 1 31 30
5 4 7 4 1 14 35
2 4 0 4 1 11 57
3 7 7 7 1 67 20
4 18 2 5 1 41 22
2 1 1 0 1 7 28
0 2 4 0 1 22 23
5 4 0 3 1 13 40
11 14 25 15 1 46 43
10 5 3 8 1 36 21
19 7 6 7 1 38 35
1 1 2 4 1 7 25
6 10 8 8 1 36 26
2 1 0 0 1 11 25
102 65 72 63 1 151 22
4 3 2 4 1 22 32
8 6 5 7 1 42 25
1 3 1 5 1 32 35
18 11 28 13 1 56 21
6 3 4 0 1 24 41
3 5 4 3 1 16 32
1 23 19 8 1 22 26
2 3 0 1 1 25 21
0 0 0 0 1 13 36
1 4 3 2 1 12 37
;
Proc Sort Data=EpilepticFits;
By Treatment Subject Period;
Run;
GOptions Reset=Symbol Reset=Axis;
Proc GPlot Data=EpilepticFits;
By Treatment;
Plot NumberOfFits*Period=Subject / VAxis=Axis1 HAxis=Axis2 NoLegend;
Axis1 Label=(A=90 "Number of Fits")
Order=(0 To 110 By 10);
Axis2 Label=("Two-week Period");
Symbol1 /*C=Black*/ V=None I=Join L=1 R=100;
Run;
Quit;
/*
* Include the GLIMMIX macro.
* Change the path below, if needed, to point to the location of
* GLIMMIX.sas on your system. If it is not in the SAS/STAT sasmacro
* directory or the SAS/STAT samples directory, then download it from
* the support.sas.com site (do a search there).
*/
%Include "C:\Program Files\SAS\SAS System\9.0\stat\sasmacro\glimmix.sas";
Title3 "Usual Generalized Linear Model W/O Overdispersion";
Proc GenMod Data=EpilepticFits;
Class Treatment Period;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period
/ Link=Log Dist=Poisson Type3;
Run;
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period;
Repeated / Subject=Subject(Treatment) Type=Simple;
Parms (1) / Hold=1;
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Title3 "Usual Generalized Linear Model With Overdispersion";
Proc GenMod Data=EpilepticFits;
Class Treatment Period;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period
/ Link=Log Dist=Poisson Type3 DScale;
Run;
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period / Chisq;
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Title3 "Generalized Linear Model With Clustered Data (GEE Approach)";
Proc GenMod Data=EpilepticFits;
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period
/ Link=Log Dist=Poisson Type3 DScale;
Repeated Subject=Subject(Treatment) / Type=AR(1);
Run;
/*
ProcOpt=Update,
Options=NOTES,
*/
Title4 "Correlated Errors";
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period / DDFM=KenwardRoger;
Repeated Period / Subject=Subject(Treatment) Type=AR(1);
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Title4 "Random Effect";
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period / DDFM=KenwardRoger;
Random Subject(Treatment);
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Title4 "Random Effect and Correlated Errors";
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Period Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period / DDFM=KenwardRoger;
Random Subject(Treatment);
Repeated Period / Subject=Subject(Treatment) Type=AR(1);
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Title4 "Random Coefficient Regression Model";
%Glimmix(Data=EpilepticFits,
Stmts=%Str(
Class Treatment Subject;
Model NumberOfFits = Age BaseLine Treatment Period Treatment*Period / Solution;
Random Intercept Period / Solution Subject=Subject(Treatment) Type=UN;
),
Error=Poisson,
Link=Log,
Out=PoissonResults
)
Run;
Proc DataSets Library=Work Kill NoList;
Run;
Quit;
/*
* Example 5
*/
Title2 "Sitka Tree Growth";
Proc Format;
Value Trt 1="Ozone" 2="Control";
Run;
Data Trees;
Infile DataLines Missover;
Input Time1-Time13;
If Time1=. Then
Do;
Chamber+1;
Delete;
End;
Tree+1;
If Chamber=1 or Chamber=2 Then Treatment=1;
Else Treatment=2;
Format Treatment Trt.;
Label Time1 ="Day 152" Time2 ="Day 174"
Time3 ="Day 201" Time4 ="Day 227"
Time5 ="Day 258" Time6 ="Day 469"
Time7 ="Day 496" Time8 ="Day 528"
Time9 ="Day 556" Time10="Day 579"
Time11="Day 613" Time12="Day 639"
Time13="Day 674";
DataLines;
4.51 4.98 5.41 5.9 6.15 6.16 6.18 6.48 6.65 6.87 6.95 6.99 7.04
4.24 4.2 4.68 4.92 4.96 5.2 5.22 5.39 5.65 5.71 5.78 5.82 5.85
3.98 4.36 4.79 4.99 5.03 5.87 5.88 6.04 6.34 6.49 6.58 6.65 6.61
4.36 4.77 5.1 5.3 5.36 5.53 5.56 5.68 5.93 6.21 6.26 6.2 6.19
4.34 4.95 5.42 5.97 6.28 6.5 6.5 6.79 6.83 7.1 7.17 7.21 7.16
4.59 5.08 5.36 5.76 6.0 6.33 6.34 6.39 6.78 6.91 6.99 7.01 7.05
4.41 4.56 4.95 5.23 5.33 6.13 6.14 6.36 6.57 6.78 6.82 6.81 6.86
4.24 4.64 4.95 5.38 5.48 5.61 5.63 5.82 6.18 6.42 6.48 6.47 6.46
4.82 5.17 5.76 6.12 6.24 6.48 6.5 6.77 7.14 7.26 7.3 6.91 7.28
3.84 4.17 4.67 4.67 4.8 4.94 4.94 5.05 5.33 5.53 5.56 5.57 5.6
4.07 4.31 4.9 5.1 5.1 5.26 5.26 5.38 5.66 5.81 5.84 5.93 5.89
4.28 4.8 5.27 5.55 5.65 5.76 5.77 5.98 6.18 6.39 6.43 6.44 6.41
4.47 4.89 5.23 5.55 5.74 5.99 6.01 6.08 6.39 6.45 6.57 6.57 6.58
4.46 4.84 5.11 5.34 5.46 5.47 5.49 5.7 5.93 6.06 6.15 6.12 6.12
4.6 4.08 4.17 4.35 4.59 4.65 4.69 5.01 5.21 5.38 5.58 5.46 5.5
3.73 4.15 4.61 4.87 4.93 5.24 5.25 5.25 5.45 5.65 5.65 5.76 5.83
4.67 4.88 5.18 5.34 5.49 6.44 6.44 6.61 6.74 7.06 7.11 7.04 7.11
2.96 3.47 3.76 3.89 4.3 4.15 4.15 4.41 4.72 4.76 4.93 4.98 5.07
3.24 3.93 4.76 4.62 4.64 4.63 4.64 4.77 5.08 5.27 5.3 5.43 5.2
4.36 4.77 5.02 5.26 5.45 5.44 5.44 5.49 5.73 5.77 6.01 5.96 5.96
4.04 4.64 4.86 5.09 5.25 5.25 5.27 5.5 5.65 5.69 5.97 5.97 5.89
3.53 4.25 4.68 4.97 5.18 5.64 5.64 5.53 5.74 5.78 5.94 6.18 5.99
4.22 4.69 5.07 5.37 5.58 5.76 5.8 6.11 6.37 6.35 6.58 6.55 6.55
2.79 3.1 3.3 3.38 3.55 3.61 3.65 3.93 4.18 4.13 4.36 4.43 4.39
3.3 3.9 4.34 4.96 5.4 5.46 5.49 5.77 6.03 6.07 6.2 6.26 6.28
3.34 3.81 4.21 4.54 4.86 4.93 4.96 5.15 5.48 5.49 5.7 5.74 5.74
3.76 4.36 4.7 5.44 5.32 5.65 5.67 5.63 6.04 6.02 6.05 6.03 5.91
4.49 4.76 5.15 5.37 5.56 5.73 5.73 5.8 5.97 6.1 6.16 6.22 6.13
4.88 5.14 5.52 6.08 6.17 6.32 6.33 6.37 6.68 6.83 6.94 6.93 6.95
4.88 5.32 5.63 5.75 5.94 6.09 6.09 6.14 6.51 6.61 6.68 6.64 6.74
3.8 4.16 4.45 4.89 5.05 5.06 5.06 5.13 5.32 5.46 5.46 5.42 5.49
4.46 4.62 5.0 5.4 5.49 5.68 5.72 5.95 6.13 6.32 6.33 6.33 6.3
4.29 4.82 5.32 5.46 5.5 5.54 5.54 5.54 5.6 5.57 5.55 5.55 5.55
4.06 4.58 4.81 5.12 5.27 5.48 5.51 5.58 5.93 6.36 6.17 6.07 6.13
5.16 5.43 5.71 6.08 6.21 6.37 6.38 6.41 6.64 6.82 6.89 7.11 7.14
3.81 4.12 4.42 4.62 4.6 4.74 4.76 4.94 5.1 5.21 5.23 5.22 5.23
5.09 5.62 5.9 6.36 6.49 6.72 6.72 6.74 6.87 6.87 6.87 6.83 6.97
4.13 4.71 5.27 5.56 5.72 6.06 6.06 6.21 6.44 6.66 6.71 6.65 6.64
4.85 5.36 5.52 5.96 6.13 6.22 6.24 6.41 6.58 6.78 6.83 6.82 6.8
4.11 4.62 4.95 5.28 5.43 5.8 5.8 5.87 6.2 6.44 6.44 6.4 6.44
4.95 5.39 5.82 6.42 6.48 6.61 6.61 6.66 6.73 6.83 6.9 6.83 6.63
4.36 4.65 5.04 5.38 5.47 5.48 5.48 5.47 5.84 5.97 5.93 5.95 6.01
4.05 4.65 5.09 5.44 5.6 5.79 5.8 6.07 6.14 6.3 6.32 6.34 6.42
3.76 4.27 4.59 5.1 5.25 5.41 5.44 5.48 5.93 5.97 6.08 6.29 6.24
2.84 3.25 3.69 4.16 4.21 4.3 4.3 4.45 4.59 4.74 4.84 4.64 4.64
4.33 4.8 5.09 5.42 5.61 5.85 5.88 6.01 6.22 6.45 6.55 6.55 6.55
3.99 4.55 4.91 5.26 5.3 5.69 5.69 5.9 5.89 5.98 6.05 6.25 6.25
3.5 3.75 3.97 4.71 4.85 5.01 5.02 5.27 5.45 5.59 5.67 5.83 5.86
3.31 3.45 4.16 4.48 4.54 4.72 4.74 4.93 5.07 5.26 5.26 5.35 5.35
3.03 3.55 3.97 4.4 4.58 4.47 4.47 4.66 4.8 5.1 5.08 5.12 5.12
3.27 3.83 4.44 4.8 4.89 5.08 5.09 5.34 5.63 5.81 5.93 5.94 5.94
3.56 4.18 4.7 5.27 5.28 5.5 5.5 5.77 5.98 6.05 6.19 6.14 6.14
3.39 3.73 3.92 4.11 4.15 4.49 4.52 4.82 5.18 5.26 5.32 5.28 5.28
3.72 4.16 4.55 5.03 5.02 5.16 5.16 5.28 5.52 5.7 5.7 5.67 5.67
4.53 5.05 5.18 5.41 5.42 5.71 5.71 5.96 6.17 6.34 6.44 6.46 6.43
4.97 5.32 5.83 6.29 6.45 6.61 6.61 6.79 7.13 7.24 7.32 7.29 7.35
4.37 4.81 5.03 5.19 5.4 5.57 5.59 5.82 6.03 6.17 6.32 6.25 6.29
4.58 4.99 5.37 5.68 5.93 6.14 6.17 6.43 6.56 6.69 6.81 6.82 6.76
4.0 4.5 4.92 5.44 5.87 6.02 6.04 6.11 6.47 6.61 6.66 6.7 6.65
4.73 5.05 5.33 5.92 6.01 6.26 6.26 6.36 6.49 6.63 6.92 6.92 6.92
5.15 5.63 6.11 6.39 6.61 6.82 6.82 6.95 7.11 7.44 7.53 7.46 7.56
4.1 4.46 4.84 5.29 5.48 5.68 5.68 5.85 5.99 6.08 6.25 6.15 6.18
3.22 3.85 4.47 4.85 5.11 5.28 5.28 5.45 5.74 5.95 6.05 6.07 6.02
2.23 2.89 3.16 3.4 3.52 3.89 3.93 4.22 4.51 4.65 4.7 4.73 4.68
3.65 4.36 4.76 5.18 5.44 5.7 5.7 5.89 6.09 6.39 6.57 6.36 6.44
3.4 3.92 4.5 4.97 5.14 5.34 5.35 5.61 5.83 6.03 6.09 5.98 5.95
5.16 5.49 5.74 6.05 6.21 6.37 6.37 6.52 6.65 6.86 6.87 6.88 6.84
4.04 4.52 5.15 5.59 5.87 5.96 5.96 6.17 6.37 6.53 6.6 6.52 6.59
4.52 4.91 5.04 5.71 5.97 6.11 6.12 6.24 6.44 6.54 6.65 6.63 6.64
4.56 5.12 5.4 5.69 5.89 6.16 6.17 6.13 6.44 6.72 6.81 6.87 6.8
4.9 5.35 5.71 6.12 6.25 6.39 6.39 6.52 6.86 7.05 7.09 6.9 6.88
4.83 5.1 5.43 5.59 6.04 6.21 6.21 6.46 6.45 6.59 6.7 6.63 6.66
5.46 5.79 6.12 6.41 6.63 6.73 6.73 6.77 6.68 6.75 6.75 6.62 6.6
4.17 4.67 5.16 5.56 5.75 6.0 6.02 6.14 6.28 6.55 6.66 6.63 6.63
3.35 4.05 4.51 5.22 5.44 5.79 5.82 6.05 6.29 6.22 6.39 6.47 6.42
3.33 3.82 4.38 4.99 5.17 5.4 5.4 5.73 5.85 5.75 5.99 6.1 6.15
3.41 3.68 4.03 4.28 4.54 4.52 4.57 5.01 5.13 5.11 5.3 5.46 5.35
4.5 4.8 5.28 5.83 6.16 6.33 6.34 6.56 6.63 6.75 6.89 6.96 6.94
2.99 3.61 4.48 4.91 5.06 5.23 5.25 5.56 5.95 5.98 6.21 6.28 6.34
;
Proc Sort Data=Trees;
By Tree Chamber Treatment;
Run;
Proc Transpose Data=Trees Out=Sitka;
By Tree Chamber Treatment;
Var Time1-Time13;
Run;
Data Sitka;
Set Sitka;
Time=Input(Substr(_LABEL_,5,3),3.); Drop _NAME_ _LABEL_;
Season=2-(Time <= 496);
Size=Exp(Col1);
Rename Col1=LogSize;
Run;
/*
* Keep only the CONTROL trees for the first season.
* I'll let you write a model to compare treatments.
*/
Proc Sort Data=Sitka Out=Season1;
Where (Season=1) And (Treatment=2);
By Tree Time;
Run;
Title3 "Raw Growth Profiles";
Title4 "First Growing Season Data";
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPLOT Data=Season1;
Plot Size*Time=Tree / NoLegend VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Size")
Order=(0 To 1000 By 200);
Axis2 Label=("Days Since January 1, 1988")
Order=(100 To 500 By 100);
Symbol1 V=None I=Join L=1 /*C=Black*/ W=2 R=1000;
Run;
Quit;
/*
* Smooth spline profiles.
*/
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPLOT Data=Season1;
Plot Size*Time=Tree / NoLegend VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Size")
Order=(0 To 1000 By 200);
Axis2 Label=("Days Since January 1, 1988")
Order=(100 To 500 By 100);
Symbol1 V=None I=Spline L=1 /*C=Black*/ R=1000;
Run;
Quit;
/*
* Fit a von Bertalanfy growth curve to season 1 of tree 70.
* Using just a single subject, we can check to see that the
* model formulation is working and that we have reasonable
* starting values. This can also be a technique for getting
* starting values for the full model.
*/
Title3 "von Bertalanfy Growth Model";
Proc NLIN Data=Season1;
Where (Tree=70);
parameters A=1000 k=1 t0=100;
Model Size = A*(1-exp(-k*(Time-t0)/400));
Output Out=Diagnostics P=PredSize;
Run;
Title4 "Predicted Growth Profiles";
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPLOT Data=Diagnostics;
Where Tree=70;
By Treatment Notsorted;
Plot Size*Time=1 PredSize*Time=2 / Overlay NoLegend VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Size")
Order=(0 To 1000 By 200);
Axis2 Label=("Days Since January 1, 1988")
Order=(100 To 500 By 100);
Symbol1 V=Dot I=None L=1 C=Black;
Symbol2 V=None I=Spline L=1 C=Black;
Run;
Quit;
Title3;
/*
* Include the NLINMIX macro.
* Change the path below, if needed, to point to the location of
* NLINMIX.sas on your system. If it is not in the SAS/STAT sasmacro
* directory or the SAS/STAT samples directory, then download it from
* the support.sas.com site (do a search there).
*/
%Include "D:\Home\Barry\LSU\Projects\MixedModels\SUGI29\nlinmix.sas";
/*
* Based upon the profiles, it appears that a random asymptote model
* is a good place to start.
*/
Title3 "Random Asymptote Model";
%NLINMix(Data=Season1,
Response=Size,
Subject=Tree,
Model=%Str(
expr1=exp(-k*(Time-t0)/400);
expr2=1-expr1;
Pred = (A+b1)*(expr2);
),
Derivs=%Str(
D_A = expr2;
D_k = (A+b1)*(Time-t0)/400*expr1;
D_t0 = -(A+b1)*k/400*expr1;
D_b1 = D_A;
),
Parms=%str(
A =1000
k =3
t0=100),
Random=b1,
Expand=EBLUP
)
Run;
Proc Transpose Data=_Soln Out=Betas(Rename=(D_A=A D_k=k D_t0=t0));
Id Effect;
Var Estimate;
Run;
Data Predictions;
If _N_=1 Then Set Betas;
Set Season1;
By Tree;
If First.Tree Then Set _SolnR(Rename=(Estimate=b1));
expr1=exp(-k*(Time-t0)/400);
expr2=1-expr1;
Pred = (A+b1)*(expr2);
Run;
Title4 "Predicted Growth Profiles";
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPLOT Data=Predictions;
Plot Pred*Time=Tree / NoLegend VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Size")
Order=(0 To 1000 By 200);
Axis2 Label=("Days Since January 1, 1988")
Order=(100 To 500 By 100);
Symbol1 V=None I=Spline L=1 /*C=Black*/ R=1000;
Run;
Quit;
Title3;
/*
* Now also permit the linearized errors to be autocorrelated.
*/
Title3 "Random Asymptote With First-order Autocorrelation Errors Model";
%NLINMix(Data=Season1,
Response=Size,
Subject=Tree,
Model=%Str(
expr1=exp(-k*(Time-t0)/400);
expr2=1-expr1;
Pred = (A+b1)*(expr2);
),
Derivs=%Str(
D_A = expr2;
D_k = (A+b1)*(Time-t0)/400*expr1;
D_t0 = -(A+b1)*k/400*expr1;
D_b1 = D_A;
),
Parms=%str(
A =1000
k =3
t0=100),
Random=b1,
Expand=EBLUP,
RType=SP(POW)(Time),
RSub=Tree
)
Run;
Proc Transpose Data=_Soln Out=Betas(Rename=(D_A=A D_k=k D_t0=t0));
Id Effect;
Var Estimate;
Run;
Data Predictions;
If _N_=1 Then Set Betas;
Set Season1;
By Tree;
If First.Tree Then Set _SolnR(Rename=(Estimate=b1));
expr1=exp(-k*(Time-t0)/400);
expr2=1-expr1;
Pred = (A+b1)*(expr2);
Run;
Title4 "Predicted Growth Profiles";
GOptions Reset=Symbol Reset=Axis Reset=Legend;
Proc GPLOT Data=Predictions;
Plot Pred*Time=Tree / NoLegend VAxis=Axis1 HAxis=Axis2;
Axis1 Label=(A=90 "Size")
Order=(0 To 1000 By 200);
Axis2 Label=("Days Since January 1, 1988")
Order=(100 To 500 By 100);
Symbol1 V=None I=Spline L=1 /*C=Black*/ R=1000;
Run;
Quit;
Title3;
Proc DataSets Library=Work Kill NoList;
Run;
Quit;
ODS HTML Close;
ODS Results On;
ODS Listing;