This example conducts an analysis similar to Example 15 in Examples: GLIMMIX Procedure: Examples: GLIMMIX Procedure. It uses simulated data to perform multiple comparisons among predicted values in a model with group-specific trends that
are modeled through regression splines. The estimable functions are formed using nonpositional syntax with constructed effects.
Consider the data in the following DATA step. Each of the 100 observations for the continuous response variable y
is associated with one of two groups.
data spline; input group y @@; x = _n_; datalines; 1 -.020 1 0.199 2 -1.36 1 -.026 2 -.397 1 0.065 2 -.861 1 0.251 1 0.253 2 -.460 2 0.195 2 -.108 1 0.379 1 0.971 1 0.712 2 0.811 2 0.574 2 0.755 1 0.316 2 0.961 2 1.088 2 0.607 2 0.959 1 0.653 1 0.629 2 1.237 2 0.734 2 0.299 2 1.002 2 1.201 1 1.520 1 1.105 1 1.329 1 1.580 2 1.098 1 1.613 2 1.052 2 1.108 2 1.257 2 2.005 2 1.726 2 1.179 2 1.338 1 1.707 2 2.105 2 1.828 2 1.368 1 2.252 1 1.984 2 1.867 1 2.771 1 2.052 2 1.522 2 2.200 1 2.562 1 2.517 1 2.769 1 2.534 2 1.969 1 2.460 1 2.873 1 2.678 1 3.135 2 1.705 1 2.893 1 3.023 1 3.050 2 2.273 2 2.549 1 2.836 2 2.375 2 1.841 1 3.727 1 3.806 1 3.269 1 3.533 1 2.948 2 1.954 2 2.326 2 2.017 1 3.744 2 2.431 2 2.040 1 3.995 2 1.996 2 2.028 2 2.321 2 2.479 2 2.337 1 4.516 2 2.326 2 2.144 2 2.474 2 2.221 1 4.867 2 2.453 1 5.253 2 3.024 2 2.403 1 5.498 ;
The following statements fit a model with separate trends for the two groups; the trends are modeled as B-splines.
proc orthoreg data=spline; class group; effect spl = spline(x); model y = group spl*group / noint; store ortho_spline; title 'B-splines Comparisons'; run;
Results from this analysis are shown in Output 87.6.1. The "Parameter Estimates" table shows the estimates for the spline coefficients in the two groups.
Output 87.6.1: Results for Group-Specific Spline Model
Parameter | DF | Parameter Estimate | Standard Error |
t Value | Pr > |t| |
---|---|---|---|---|---|
(group='1') | 1 | 9.70265463962039 | 3.1341899987 | 3.10 | 0.0026 |
(group='2') | 1 | 6.30619220563569 | 2.6299147768 | 2.40 | 0.0187 |
spl*group 1 1 | 1 | -11.1786451718041 | 3.7008097395 | -3.02 | 0.0033 |
spl*group 1 2 | 1 | -20.1946092746139 | 3.9765046236 | -5.08 | <.0001 |
spl*group 2 1 | 1 | -9.53273697995301 | 3.2575832048 | -2.93 | 0.0044 |
spl*group 2 2 | 1 | -5.85652496534967 | 2.7906116773 | -2.10 | 0.0388 |
spl*group 3 1 | 1 | -8.96118371893294 | 3.0717508806 | -2.92 | 0.0045 |
spl*group 3 2 | 1 | -5.55671605245205 | 2.5716715573 | -2.16 | 0.0335 |
spl*group 4 1 | 1 | -7.26153231478755 | 3.243690314 | -2.24 | 0.0278 |
spl*group 4 2 | 1 | -4.36778889738236 | 2.7246809593 | -1.60 | 0.1126 |
spl*group 5 1 | 1 | -6.44615256510896 | 2.9616955361 | -2.18 | 0.0323 |
spl*group 5 2 | 1 | -4.03801618914902 | 2.4588839125 | -1.64 | 0.1042 |
spl*group 6 1 | 1 | -4.63816959094139 | 3.7094636319 | -1.25 | 0.2146 |
spl*group 6 2 | 1 | -4.30290104395061 | 3.0478540171 | -1.41 | 0.1616 |
spl*group 7 1 | 0 | 0 | . | . | . |
spl*group 7 2 | 0 | 0 | . | . | . |
By default, the ORTHOREG procedure constructs B-splines with seven knots. Since B-spline coefficients satisfy a sum-to-one constraint and since the model contains group-specific intercepts, the last spline coefficient for each group is redundant and is set to 0.
The following statements make a prediction for the input data set by using the SCORE statement with PROC PLM and graph the observed and predicted values in the two groups:
proc plm restore=ortho_spline; score data=spline out=ortho_pred predicted=p; run; proc sgplot data=ortho_pred; series y=p x=x / group=group name="fit"; scatter y=y x=x / group=group; keylegend "fit" / title="Group"; run;
The prediction plot in Output 87.6.2 suggests that there is some separation of the group trends for small values of x
and for values that exceed about x = 40.
Output 87.6.2: Observed Data and Predicted Values by Group
In order to determine the range on which the trends separate significantly, the PLM procedure is executed in the following
statements with an ESTIMATE
statement that applies group comparisons at a number of values for the spline variable x
:
%macro GroupDiff; %do x=0 %to 75 %by 5; "Diff at x=&x" group 1 -1 group*spl [1,1 &x] [-1,2 &x], %end; 'Diff at x=80' group 1 -1 group*spl [1,1 80] [-1,2 80] %mend; proc plm restore=ortho_spline; show effects; estimate %GroupDiff / adjust=simulate seed=1 stepdown; run;
For example, the following ESTIMATE statement compares the trends between the two groups at x = 25:
estimate 'Diff at x=25' group 1 -1 group*spl [1,1 25] [-1,2 25];
The nonpositional syntax is used for the group
*spl
effect. For example, the specification requests that the spline be computed at x = 25 for the second level of variable group
. The resulting coefficients are added to the vector for the estimate after being multiplied with –1.
Because comparisons are made at a large number of values for x
, a multiplicity correction is in order to adjust the p-values to reflect familywise error control. Simulated p-values with step-down adjustment are used here.
Output 87.6.3 displays the "Store Information" for the item store and information about the spline effect (the result of the SHOW statement).
Output 87.6.3: Spline Details
Output 87.6.4 displays the results from the ESTIMATE statement.
Output 87.6.4: Estimate Results with Multiplicity Correction
Estimates Adjustment for Multiplicity: Holm-Simulated |
||||||
---|---|---|---|---|---|---|
Label | Estimate | Standard Error | DF | t Value | Pr > |t| | Adj P |
Diff at x=0 | 12.4124 | 4.2130 | 86 | 2.95 | 0.0041 | 0.0206 |
Diff at x=5 | 1.0376 | 0.1759 | 86 | 5.90 | <.0001 | <.0001 |
Diff at x=10 | 0.3778 | 0.1540 | 86 | 2.45 | 0.0162 | 0.0545 |
Diff at x=15 | 0.05822 | 0.1481 | 86 | 0.39 | 0.6952 | 0.9101 |
Diff at x=20 | -0.02602 | 0.1243 | 86 | -0.21 | 0.8346 | 0.9565 |
Diff at x=25 | 0.02014 | 0.1312 | 86 | 0.15 | 0.8783 | 0.9565 |
Diff at x=30 | 0.1023 | 0.1378 | 86 | 0.74 | 0.4600 | 0.7418 |
Diff at x=35 | 0.1924 | 0.1236 | 86 | 1.56 | 0.1231 | 0.2925 |
Diff at x=40 | 0.2883 | 0.1114 | 86 | 2.59 | 0.0113 | 0.0450 |
Diff at x=45 | 0.3877 | 0.1195 | 86 | 3.24 | 0.0017 | 0.0096 |
Diff at x=50 | 0.4885 | 0.1308 | 86 | 3.74 | 0.0003 | 0.0024 |
Diff at x=55 | 0.5903 | 0.1231 | 86 | 4.79 | <.0001 | <.0001 |
Diff at x=60 | 0.7031 | 0.1125 | 86 | 6.25 | <.0001 | <.0001 |
Diff at x=65 | 0.8401 | 0.1203 | 86 | 6.99 | <.0001 | <.0001 |
Diff at x=70 | 1.0147 | 0.1348 | 86 | 7.52 | <.0001 | <.0001 |
Diff at x=75 | 1.2400 | 0.1326 | 86 | 9.35 | <.0001 | <.0001 |
Diff at x=80 | 1.5237 | 0.1281 | 86 | 11.89 | <.0001 | <.0001 |
Notice that the "Store Information" in Output 87.6.3 displays the classification variables (from the CLASS statement in PROC ORTHOREG), the constructed effects (from the EFFECT statement in PROC ORTHOREG), and the model effects (from the MODEL statement in PROC ORTHOREG). Output 87.6.4 shows that at the 5% significance level the trends are significantly different for and for . Between those values you cannot reject the hypothesis of trend congruity.
To see this effect more clearly, you can filter the results by adding the a filtering statement to the previous PROC PLM step:
proc plm restore=ortho_spline; estimate %GroupDiff / adjust=simulate seed=1 stepdown; filter adjp > 0.05; run;
This produces Output 87.6.5, which displays the subset of the results in Output 87.6.4 that meets the condition in the FILTER expression.
Output 87.6.5: Filtered Estimate Results
B-splines Comparisons |
Estimates Adjustment for Multiplicity: Holm-Simulated |
||||||
---|---|---|---|---|---|---|
Label | Estimate | Standard Error | DF | t Value | Pr > |t| | Adj P |
Diff at x=10 | 0.3778 | 0.1540 | 86 | 2.45 | 0.0162 | 0.0545 |
Diff at x=15 | 0.05822 | 0.1481 | 86 | 0.39 | 0.6952 | 0.9101 |
Diff at x=20 | -0.02602 | 0.1243 | 86 | -0.21 | 0.8346 | 0.9565 |
Diff at x=25 | 0.02014 | 0.1312 | 86 | 0.15 | 0.8783 | 0.9565 |
Diff at x=30 | 0.1023 | 0.1378 | 86 | 0.74 | 0.4600 | 0.7418 |
Diff at x=35 | 0.1924 | 0.1236 | 86 | 1.56 | 0.1231 | 0.2925 |