### Example 41.10 Multiple Trends Correspond to Multiple Extrema in Profile Likelihoods

Observations for a period of 168 months for the Southern Oscillation Index, measurements of monthly averaged atmospheric pressure differences between Easter Island and Darwin, Australia (Kahaner, Moler, and Nash 1989, Ch. 11.9; National Institute of Standards and Technology 1998) is available in the data set ENSO in the Sashelp library. These data are also used as an example in Chapter 53: The LOESS Procedure. The following statements print the first 10 observations of this data set in Output 41.10.1.

proc print data=Sashelp.enso(obs=10);
run;


Output 41.10.1: El Niño Southern Oscillation Data

Obs Month Year Pressure
1 1 0.08333 12.9
2 2 0.16667 11.3
3 3 0.25000 10.6
4 4 0.33333 11.2
5 5 0.41667 10.9
6 6 0.50000 7.5
7 7 0.58333 7.7
8 8 0.66667 11.7
9 9 0.75000 12.9
10 10 0.83333 14.3

Differences in atmospheric pressure create wind, and the differences recorded in the data set ENSO drive the trade winds in the southern hemisphere. Such time series often do not consist of a single trend or cycle. In this particular case, there are at least two known cycles that reflect the annual weather pattern and a longer cycle that represents the periodic warming of the Pacific Ocean (El Niño).

To estimate the trend in these data by using mixed model technology, you can apply a mixed model smoothing technique such as TYPE=RSMOOTH or TYPE=PSPLINE. The following statements fit a radial smoother to the ENSO data and obtain profile likelihoods for a series of values for the variance of the random spline coefficients:

data tdata;
do covp1=0,0.0005,0.05,0.1,0.2,0.5,
1,2,3,4,5,6,8,10,15,20,50,
75,100,125,140,150,160,175,
200,225,250,275,300,350;
output;
end;
run;

ods select FitStatistics CovParms CovTests;
proc glimmix data=sashelp.enso noprofile;
model pressure = year;
random year / type=rsmooth knotmethod=equal(50);
parms (2) (10);
covtest tdata=tdata / parms;
ods output covtests=ct;
run;


The tdata data set contains value for the variance of the radial smoother variance for which the profile likelihood of the model is to be computed. The profile likelihood is obtained by setting the radial smoother variance at the specified value and estimating all other parameters subject to that constraint.

Because the model contains a residual variance and you need to specify nonzero values for the first covariance parameter, the NOPROFILE is added to the PROC GLIMMIX statements. If the residual variance is profiled from the estimation, you cannot fix covariance parameters at a given value, because they would be reexpressed during model fitting in terms of ratios with the profiled (and changing) variance.

The PARMS statement determines starting values for the covariance parameters for fitting the (full) model. The PARMS option in the COVTEST statement requests that the input parameters be added to the output and the output data set. This is useful for subsequent plotting of the profile likelihood function.

The Fit Statistics table displays the –2 restricted log likelihood of the model (897.76, Output 41.10.2). The estimate of the variance of the radial smoother coefficients is 3.5719.

The Test of Covariance Parameters table displays the –2 restricted log likelihood for each observation in the tdata set. Because the tdata data set specifies values for only the first covariance parameter, the second covariance parameter is free to vary and the values for –2 Res Log Like are profile likelihoods. Notice that for a number of values of CovP1 the chi-square statistic is missing in this table. For these values the –2 Res Log Like is smaller than that of the full model. The model did not converge to a global minimum of the negative restricted log likelihood.

Output 41.10.2: REML and Profile Likelihood Analysis

The GLIMMIX Procedure

Fit Statistics
-2 Res Log Likelihood 897.76
AIC (smaller is better) 901.76
AICC (smaller is better) 901.83
BIC (smaller is better) 897.76
CAIC (smaller is better) 899.76
HQIC (smaller is better) 897.76
Generalized Chi-Square 1554.38
Gener. Chi-Square / DF 9.36

Covariance Parameter Estimates
Cov Parm Estimate Standard Error
Var[RSmooth(Year)] 3.5719 3.7672
Residual 9.3638 1.3014

Tests of Covariance Parameters
Based on the Restricted Likelihood
Label DF -2 Res Log Like ChiSq Pr > ChiSq Input Parameters Note
CovP1 CovP2
WORK.TDATA 1 893.01 . 1.0000 0 9.3638 MI
WORK.TDATA 1 892.76 . 1.0000 0.000500 9.3638 DF
WORK.TDATA 1 897.34 . 1.0000 0.05000 9.3638 DF
WORK.TDATA 1 898.53 0.77 0.3816 0.1000 9.3638 DF
WORK.TDATA 1 899.38 1.62 0.2038 0.2000 9.3638 DF
WORK.TDATA 1 899.49 1.73 0.1888 0.5000 9.3638 DF
WORK.TDATA 1 898.83 1.07 0.3016 1.0000 9.3638 DF
WORK.TDATA 1 898.04 0.28 0.5967 2.0000 9.3638 DF
WORK.TDATA 1 897.79 0.03 0.8693 3.0000 9.3638 DF
WORK.TDATA 1 897.77 0.01 0.9145 4.0000 9.3638 DF
WORK.TDATA 1 897.86 0.10 0.7517 5.0000 9.3638 DF
WORK.TDATA 1 897.99 0.23 0.6311 6.0000 9.3638 DF
WORK.TDATA 1 898.27 0.51 0.4761 8.0000 9.3638 DF
WORK.TDATA 1 898.49 0.73 0.3919 10.0000 9.3638 DF
WORK.TDATA 1 898.70 0.94 0.3318 15.0000 9.3638 DF
WORK.TDATA 1 898.45 0.69 0.4068 20.0000 9.3638 DF
WORK.TDATA 1 892.63 . 1.0000 50.0000 9.3638 DF
WORK.TDATA 1 887.44 . 1.0000 75.0000 9.3638 DF
WORK.TDATA 1 883.79 . 1.0000 100.00 9.3638 DF
WORK.TDATA 1 881.55 . 1.0000 125.00 9.3638 DF
WORK.TDATA 1 880.72 . 1.0000 140.00 9.3638 DF
WORK.TDATA 1 . . . 150.00 9.3638
WORK.TDATA 1 880.07 . 1.0000 160.00 9.3638 DF
WORK.TDATA 1 879.85 . 1.0000 175.00 9.3638 DF
WORK.TDATA 1 . . . 200.00 9.3638
WORK.TDATA 1 880.21 . 1.0000 225.00 9.3638 DF
WORK.TDATA 1 880.80 . 1.0000 250.00 9.3638 DF
WORK.TDATA 1 881.56 . 1.0000 275.00 9.3638 DF
WORK.TDATA 1 882.44 . 1.0000 300.00 9.3638 DF
WORK.TDATA 1 884.41 . 1.0000 350.00 9.3638 DF

DF: P-value based on a chi-square with DF degrees of freedom.
MI: P-value based on a mixture of chi-squares.

The following statements plot the –2 restricted profile log likelihood (Output 41.10.3):

proc sgplot data=ct;
series y=objective x=covp1;
run;


Output 41.10.3: –2 Restricted Profile Log Likelihood for Smoothing Variance

The local minimum at which the optimization stopped is clearly visible, as are a second local minimum near zero and the global minimum near 180.

The observed and predicted pressure differences that correspond to the three minima are shown in Output 41.10.4. These results were produced with the following statements:

proc glimmix data=sashelp.enso;
model pressure = year;
random year / type=rsmooth knotmethod=equal(50);
parms (0) (10);
output out=gmxout1 pred=pred1;
run;
proc glimmix data=sashelp.enso;
model pressure = year;
random year / type=rsmooth knotmethod=equal(50);
output out=gmxout2 pred=pred2;
parms (2) (10);
run;
proc glimmix data=sashelp.enso;
model pressure = year;
random year / type=rsmooth knotmethod=equal(50);
output out=gmxout3 pred=pred3;
parms (200) (10);
run;
data plotthis; merge gmxout1 gmxout2 gmxout3;
run;
proc sgplot data=plotthis;
scatter x=year y=Pressure;
series  x=year y=pred1 /
lineattrs   = (pattern=solid thickness=2)
legendlabel = "Var[RSmooth] = 0.0005"
name        = "pred1";
series  x=year y=pred2 /
lineattrs   = (pattern=dot thickness=2)
legendlabel = "Var[RSmooth] = 3.5719"
name        = "pred2";
series  x=year y=pred3 /
lineattrs   = (pattern=dash thickness=2)
legendlabel = "Var[RSmooth] = 186.71"
name        = "pred3";
keylegend "pred1" "pred2" "pred3" / across=2;
run;


Output 41.10.4: Observed and Predicted Pressure Differences

The one-year cycle () and the El Niño cycle () are clearly visible. Notice that a larger smoother variance results in larger BLUPs and hence larger adjustments to the fixed-effects model. A large smoother variance thus results in a more wiggly fit. The third local minimum at applies only very small adjustments to the linear regression between pressure and time, creating slight curvature.