Interpreting PROC GLIMMIX output

[sld]

Like @SteveDenham , I played around with code that I thought was pertinent, and I’ll share my new vision with you. The code is attached. Ideally I would take the time to simulate data that would support a full modeling effort—something with random intercepts and random slopes and autocorrelation—but I just don’t have the time for that currently.

The code file includes some comments. I will not go into model comparisons in depth here; you can embark upon your own discovery.

The first three models, which differ only in type—CHOL, UN, or ARH(1)—are re-parameterizations of the exact same statistical model which estimates in different forms (i.e., different parameterizations) the variance among intercepts, the variance among slopes, and the covariance of intercepts and slopes. Clearly none of these models is modeling temporal autocorrelation of residuals.

So my first point is that @Nerdcy's  RANDOM statement is definitely not doing what the OP thinks it is or may want.

My second point is less definitive at the moment. I experimented with models that were intended to model the covariance structure of residuals from the regression of the response on time. My example dataset is not at all sufficient to assess these models adequately. But patterns emerged that are consistent with my mixed model experience: when you have both G-side and R-side random effects, some model formulations that might seem sensible will be over-parameterized. Refer to the Littell et al. reference in the code.

My third point is that this code exploration is based on models that assume a normal distribution for the residuals. Theoretically, you ought to be able to model autocorrelation of those residuals if you have enough data that are well-enough behaved. However, for a generalized linear mixed model, the ability to model R-side random effects (i.e., the covariance structure of residuals) may be limited or (as I understand it) non-existent. In the context of a random coefficients model, I don’t know that that’s a major deal breaker because you are focusing your attention on the variances among intercepts and slopes, which are subject-level statistics.

My fourth point is that these are complicated models. Both Steve and I have experience with mixed models, and we are still discovering things we didn’t know. Just because a model runs, you should not assume it is valid. And it is easy to under-estimate the intricacies of these models.

[Nerdcy]

Thank you. You've given me a lot to think about. I'm going to take another look at my data.

[SteveDenham]

@sld .  Not a bit surprised by any of these models and results.  Overparameterization through inclusion of time as both a continuous and a categorical variable leads to problems, especially if/when time isn't really granular, resulting in a large number of categories for time_c.

I realize that I at first missed, and then continued to miss, that the variable 'time' was continuous.  I have done so many of these with time as a categorical variable that i made yet another mistake in the code.  Here is what I would do:

/*Autoregressive residual:*/

proc glimmix data=wheat2;
    class variety time_c;
    model yield = time_c;
    random int / subject=variety; 
    random time_c / subject=variety type=ar(1) residual;
run;

/*Unstructured residual */
proc glimmix data=wheat2;
    class variety time_c;
    model yield = time_c;
    random time_c / subject=variety type=ar(1) residual;
run;

Comparing AICC values for these two models to all of the previous, I find a substantial reduction(~451 vs ~473), indicating that the models with continuous time were roughly 0.000017 to  0.00012 times as likely to retain the information in the data.when compared to models where time was considered categorical.  I haven't done the PROC AUTOREG work on the residuals, but it is likely that categorization with first-order autoregression results in a fit that is disrupted by also fitting time as a continuous variable.

If any one is looking for source material, see the technical appendix to Littell, Henry and Ammerman, Statistical Analysis of Repeated Measures Data Using SAS Procedures, (1998). J. Anim. Sci. 76:1216–1231. Here is a non-paywall link https://pdfs.semanticscholar.org/7fd4/773cd2d11a0b4c842c1d1bf83b949697a9fa.pdf .

SteveDenham

[sld]

My first point is that the over-parameterization that I noted is due to overlap/redundancy in parameterization of the RANDOM and REPEATED statements in Proc MIXED (and equivalently, in the RANDOM and RANDOM/RESIDUAL statements in Proc GLIMMIX) for certain--actually most--covariance structures. This issue is addressed in the appendix of the link to the Littell, Henry and Ammerman (1998) that @SteveDenham provides. It is also addressed in the Littell, Pendergast & Natajaran (Statistics in Medicine, 2000) paper I referenced in the code I posted previously. In this latter paper, the authors note the distinction between an AR(1) structure and an AR(1)+RE structure, as does Walt Stroup (Generalized Linear Mixed Models (2013)) in Ch 14 (Fig 14.2 does not correctly illustrate the different structures, IMO).

My second point about whether we can model R-side covariance structures is addressed by Stroup in his text on p 435, where he writes

"The primary ambiguity in repeated measures model building for non-Gaussian data occurs when we have a member of the two-parameter exponential family [e.g., gamma]. How do we model repeated measures for beta or negative binomial or other distributions in this family? .... For G-side models, it is not clear how the random ts(a)_ijk [i.e., time x subject(treatment) random effects] effect coexists with f(y|b)'s scale parameters."

In the 10 years or so since Walt wrote his text, he has continued to revisit issues. For example, he and Elizabeth Claassen published a paper in the Journal of Agricultural, Biological, and Environmental Statistics just last month (https://doi.org/10.1007/s13253-020-00402-6) entitled "Pseudo-Likelihood or Quadrature? What We Thought We Knew, What We Think We Know,
and What We Are Still Trying to Figure Out". So he may currently hold a different view about R-side covariance structure modelling with non-Gaussian distributions. Then again, he just retired and could be on to other things 🙂

My third observation is that a random coefficients model (RCM) with time continuous will be quite different from a model with time categorical.

My fourth observation is that although UN, CHOL, and ARH(1) are identical when you are modeling a RCM with random intercepts, random slopes, and one covariance (i.e., a 2 x 2 matrix), you would not want to apply an AR constraint on covariances in a covariance matrix that was larger than 2 x 2 because that would be silly. I'd stick with UN or CHOL, or in my practice VC (or UN(1)) because my datasets are always small, and I have never been able to fit a decent intercept-slope covariance.

Mathematical statistics is not my superpower; my explorations tend to be empirical. My practical experience is consistent with the point that Walt makes in his text, and with the Littell et al. papers.