Example 38.10 Multiple Trends Correspond to Multiple Extrema in Profile Likelihoods

The following data set contains 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). These data are also used as an example in Chapter 50, The LOESS Procedure, in the SAS/STAT User’s Guide.

   data ENSO;
     input Pressure @@;
     Year = _n_ / 12;
   datalines;
    12.9  11.3  10.6  11.2  10.9   7.5   7.7  11.7
    12.9  14.3  10.9  13.7  17.1  14.0  15.3   8.5
     5.7   5.5   7.6   8.6   7.3   7.6  12.7  11.0
    12.7  12.9  13.0  10.9  10.4  10.2   8.0  10.9
    13.6  10.5   9.2  12.4  12.7  13.3  10.1   7.8
     4.8   3.0   2.5   6.3   9.7  11.6   8.6  12.4
    10.5  13.3  10.4   8.1   3.7  10.7   5.1  10.4
    10.9  11.7  11.4  13.7  14.1  14.0  12.5   6.3
     9.6  11.7   5.0  10.8  12.7  10.8  11.8  12.6
    15.7  12.6  14.8   7.8   7.1  11.2   8.1   6.4
     5.2  12.0  10.2  12.7  10.2  14.7  12.2   7.1
     5.7   6.7   3.9   8.5   8.3  10.8  16.7  12.6
    12.5  12.5   9.8   7.2   4.1  10.6  10.1  10.1
    11.9  13.6  16.3  17.6  15.5  16.0  15.2  11.2
    14.3  14.5   8.5  12.0  12.7  11.3  14.5  15.1
    10.4  11.5  13.4   7.5   0.6   0.3   5.5   5.0
     4.6   8.2   9.9   9.2  12.5  10.9   9.9   8.9
     7.6   9.5   8.4  10.7  13.6  13.7  13.7  16.5
    16.8  17.1  15.4   9.5   6.1  10.1   9.3   5.3
    11.2  16.6  15.6  12.0  11.5   8.6  13.8   8.7
     8.6   8.6   8.7  12.8  13.2  14.0  13.4  14.8
   ;

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=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 restricted log likelihood of the model (, Output 38.10.1). The estimate of the variance of the radial smoother coefficients is .

The "Test of Covariance Parameters" table displays the 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.