topic Re: Interpreting PROC GLIMMIX output in Statistical Procedures

Interpreting PROC GLIMMIX output

Nerdcy — Mon, 13 Jul 2020 19:33:13 GMT

Hi, I have conducted a mixed model for longitudinal data using PROC GLIMMIX. The code I used is below:

proc glimmix data=diss method=laplace;
title Model 15: Total Support Final Model;
class carnegie (Ref='0') barrons (Ref='0') flagship(Ref='0');
EFFECT poly = polynomial(time/degree=2);
model totalsupport = 
poly carnegie barrons flagship statelog
poly*carnegie poly*flagship poly*statelog
poly*carnegie*statelog poly*flagship*statelog
/dist=gamma link=log solution;
random intercept time / type=ARH(1) subject = id;
run;

I had two questions about interpreting the output:

As you can see, this is model uses a log link and the one continuous predictor is also transformed onto a log scale (statelog).

Solutions for Fixed Effects
Effect	carnegie	barrons	flagship	Estimate	Standard Error	DF	t Value	Pr > \|t\|
Intercept				1.8144	0.6323	4661	2.87	0.0041
time				1.0161	0.1916	4661	5.30	<.0001
time^2				-0.04453	0.01360	4661	-3.27	0.0011
carnegie	1			-0.6513	0.09370	4661	-6.95	<.0001
carnegie	2			-0.5665	0.1611	4661	-3.52	0.0004
carnegie	0			0	.	.	.	.
barrons		1		0.3896	0.08327	4661	4.68	<.0001
barrons		2		-0.2169	0.09055	4661	-2.40	0.0166
barrons		0		0	.	.	.	.
flagship			1	-0.6884	0.1304	4661	-5.28	<.0001
flagship			0	0	.	.	.	.
statelog				1.5370	0.1524	4661	10.09	<.0001
timecarnegie*	1			1.0983	0.1819	4661	6.04	<.0001
timecarnegie*	2			-0.04387	0.2661	4661	-0.16	0.8691
timecarnegie*	0			0	.	.	.	.
time^2carnegie*	1			-0.07228	0.01423	4661	-5.08	<.0001
time^2carnegie*	2			-0.01247	0.01982	4661	-0.63	0.5293
time^2carnegie*	0			0	.	.	.	.
timeflagship*			1	-1.3081	0.2154	4661	-6.07	<.0001
timeflagship*			0	0	.	.	.	.
time^2flagship*			1	0.08469	0.01663	4661	5.09	<.0001
time^2flagship*			0	0	.	.	.	.
statelogtime*				-0.2541	0.04758	4661	-5.34	<.0001
statelogtime^2*				0.01149	0.003387	4661	3.39	0.0007
statelogtimecarnegie	1			-0.2803	0.04672	4661	-6.00	<.0001
statelogtimecarnegie	2			0.01977	0.06687	4661	0.30	0.7675
statelogtimecarnegie	0			0	.	.	.	.
statelogtime^2carnegie	1			0.01843	0.003647	4661	5.05	<.0001
statelogtime^2carnegie	2			0.002585	0.004956	4661	0.52	0.6019
statelogtime^2carnegie	0			0	.	.	.	.
statelogtimeflagship			1	0.3230	0.05414	4661	5.97	<.0001
statelogtimeflagship			0	0	.	.	.	.
statelogtime^2flagship			1	-0.02112	0.004179	4661	-5.05	<.0001
statelogtime^2flagship			0	0	.	.	.	.

My understanding is that I would interpret this as you would a linear regression model with logged outcome and logged predictor (% change in predictor would lead to a % change in outcome) so for:

statelog

1.5370

0.1524

4661

10.09

<.0001

As statelog increases 1.5370%, logged outcome would increase 1%.

1) Is there a way to convert this data in PROC GLIMMIX so that I can just interpret as increase in statelog = increase in 1 unit of logged outcome. I read about using lsmeans and ilink but that just broke down my categorical predictors.

2) I used an ARH(1) covariance structure for this but I'm unsure how to interpret it. I understand how to interpret the VC or UN but not ARH(1) as I've never had to use it before.

Covariance Parameter Estimates
Cov Parm	Subject	Estimate	Standard Error
Var(1)	id	0.5255	0.04765
Var(2)	id	0.001706	0.000193
ARH(1)	id	-0.4835	0.04940
Residual		0.1098	0.002474

Thank you!

Re: Interpreting PROC GLIMMIX output

SteveDenham — Tue, 14 Jul 2020 15:49:29 GMT

Some things:

You said:

My understanding is that I would interpret this as you would a linear regression model with logged outcome and logged predictor (% change in predictor would lead to a % change in outcome) so for:

statelog

1.5370

0.1524

4661

10.09

<.0001

As statelog increases 1.5370%, logged outcome would increase 1%.

This should be interpreted the other way around, for a unit increase in statelog, you would get a 1.537 unit increase in the log outcome. I wouldn't try to put this on a percentage basis in the log space, because it is something quite different in the probability space. This may answer question 1.

For question 2, a great source in the documentation are the matrix definitions of the various types. For ARH(1), this looks like

ARH(1)

specifies a heterogeneous first-order autoregressive structure,

with cov(zeta_i, zeta_j) = sqrt (sigmasquared_i * sigmasquared_j)* rho ^ abs(i - j). The documentation has this in LaTeX formula that doesn't copy well.

. This covariance structure has the same correlation pattern as the TYPE=AR(1) structure, but the variances are allowed to differ,

So in your example, you have two time points. The first is associated with Var(1) and the second with Var(2), and the correlation is associated with ARH(1). Now, with only two time points, the ARH(1) is equivalent to the unstructured (UN) or the heterogeneous compound symmetry (CSH) structures - it is just a different parameterization. Since these are G-side parameters, they are the variances at the first time point, second time point and correlation between the observed values.

Hope this is relatively understandable.

SteveDenham

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 15:37:37 GMT

Thank you! This was very helpful. I guess I just had one more question. I was told that one of the reasons to choose a glmm with link function instead of transforming the outcome variable and using a lmm is that you can easily back transform the data so that I'm not interpreting the change in log but the change in the original raw form of the variable. How do I do this in glimmmix?

Re: Interpreting PROC GLIMMIX output

sld — Wed, 15 Jul 2020 17:11:06 GMT

To stir the pot a bit:

I do not think that the RANDOM statement is correctly specified

random intercept time / type=ARH(1) subject = id;

This statement estimates random intercepts for subjects (id). We generally consider subjects to be independent; it is implausible that subjects would be temporally autocorrelated.

time is not in the CLASS statement, so it is incorporated in the model as a continuous variable. Consequently, this RANDOM statement is attempting to estimate random slopes of the regression of the response on time. I think it is implausible to assume temporal autocorrelation here either.

Your ability to model covariance structure for residuals in a non-normal GLMM is limited, or even non-existent, as I understand it.

A LMM has a residual variance, but although a GLMM has residuals, it does not have a residual variance because the mean and variance of non-normal distributions are defined by the same parameters so it is not possible to estimate both.

I do not know enough about your study design (for example, are covariates subject-specific or time-varying?) but I'll suggest the following code for a random coefficients model:

random intercept time time*time / subject = id;

If you really wanted to incorporate temporal autocorrelation, you could considering using log(totalsupport) in a LMM assuming a normal distribution. This is not the same model as a GLMM with a gamma distribution with a log link, but it might be "good enough", or it might not:

class ... c_time; /* a categorical version of time */
c_time = time;
...
random intercept time time*time / subject = id;
random c_time / subject=id type=ar(1) residual;

You might find this paper by SAS's Kathleen Kiernan useful: Insights into Using the GLIMMIX Procedure to Model Categorical Outcomes with Random Effects

Re: Interpreting PROC GLIMMIX output

sld — Wed, 15 Jul 2020 17:13:46 GMT

Use the ilink option.

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 17:19:19 GMT

I've been reading about it and I believe you're correct. I guess the question is how do I code ilink. I want the exact output I have just with the estimates using the raw scale not log link. I've read about using estimate or lsmeans but that just seems to partition out my categorical variables instead of providing the same table I put above with the new estimates.

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 17:21:20 GMT

My random statement is correct. I am using a random coefficients model to do a growth curve analysis so the random intercept shows that the subjects vary in starting point and random time shows varying slopes (growth) over time. Because it is longitudinal data, there is expected autocorrelation and heteroscedasticity making ARH(1) the right covariance structure.

Re: Interpreting PROC GLIMMIX output

sld — Wed, 15 Jul 2020 17:41:20 GMT

As you note, you can inverse link estimates of means and predictions. However, you cannot inverse link parameter estimates (I presume that's the table you are referring to). In a GLMM, all the estimates and statistical tests are done on the link scale.

You might find Ben Bolker's response here helpful: https://stats.stackexchange.com/questions/431120/how-to-interpret-parameters-of-glm-output-with-gamma-log-link ; you'll find additional discussions with an internet search.

Re: Interpreting PROC GLIMMIX output

SteveDenham — Wed, 15 Jul 2020 18:01:31 GMT

Actually, if you want a random intercept and a random slopes model, you would have 2 RANDOM statements, like this:

random intercept/subject=id;  /* This gives the random intercept variance component */
random time/type=arh(1) subject=id;  /* And this gives the random slopes variance components * /

If you have sufficient data, this should work (converges and G matrix is positive definite), but if not, consider this:

random time/subject=id type=chol;

In this case, no specific random intercept is modeled, as it is inseparable from the unstructured (using a Cholesky factorization) covariance structure. If you add in that second random statement as above for the random intercept, you will always get a message in the output that the G matrix is not positive definite.

I knew there was something else. Because you used the EFFECT option to get the polynomial with degree 2, the variance components Var(1) and Var(2) refer to the linear and quadratic parts (that's why you can get away with a type= option that is based on categorical values for time). If you use @sld 's suggestion of

random intercept time time*time / subject = id;

you should get similar estimates except that you will not have an estimate of the correlation between linear and quadratic time. At least that is what I think happened in my toy example I was testing things on (stress the part about I think).

SteveDenham.

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 17:49:03 GMT

Thank you! I just thought there was some way to automatically get the exponential but it looks like I'll just have to calculate that manually.

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 17:56:50 GMT

Thank you!

Re: Interpreting PROC GLIMMIX output

sld — Wed, 15 Jul 2020 18:11:22 GMT

@SteveDenham I'm good with your first RANDOM statement. I disagree with the second: I do not see why there would be temporal autocorrelation among the subject random effects. If intercepts for subjects are independent (as in the first RANDOM), would not their random effects for linear slopes also be independent?

I'm not sure what the default type is, but for sure you would get covariances among intercepts, linear and quadratic terms with

random intercept time time*time / subject = id type=chol;

or type=un

Edit: I might agree with your interpretation of the two reported variance components Var(1) and Var(2) if poly (the name of the polynomial effect) replaced time in the RANDOM statement. But I still don't like this use of ARH(1).

Re: Interpreting PROC GLIMMIX output

Nerdcy — Wed, 15 Jul 2020 18:13:30 GMT

Yes, I've been thinking about this. My understanding of the autocorrelation in my data is that responses closer together in time would be more closely related to each other than those farther apart in time (outcome is a dollar amount; for example - the amount received in year 10 is related to the amount received in year 9, but not as related to the amount received in year 2).

Re: Interpreting PROC GLIMMIX output

sld — Wed, 15 Jul 2020 19:03:00 GMT

Temporal autocorrelation is certainly possible. The structure of the autocorrelation might be well captured by AR(1), or other covariance structures might fit the data better (e.g., CS or TOEP, heterogeneous variances versions; or UN).

We might think of the responses as being correlated, but in a design with repeated measures on subjects, we apply these covariance structures to the residuals (R-side random effects)--what's left after the MODEL removes the fixed effect trend of time.

See G-Side and R-Side Random Effects and Covariance Structures This bit of documentation does suggest that you can do R-side modeling in GLMM; you could try it and see if it worked.

Re: Interpreting PROC GLIMMIX output

SteveDenham — Wed, 15 Jul 2020 20:03:12 GMT

@sld : Regarding your concern about this set of random statements

random intercept/subject=id; /* This gives the random intercept variance component */ random time/type=arh(1) subject=id; /* And this gives the random slopes variance components * /

and that the second implies independent slopes by subject, I wonder. I have often used this pair of statements to model autoregressive types of covariance structures, especially R-side (with the residual option in the second statement). A good reference on this would be pages 431 and 432 in Walt Stroup's Generalized Linear Mixed Models, where he presents this particular approach for a GLMM. Theoretical concepts regarding this are covered in section 14.5.1..

I have to admit I am a bit baffled by the behavior.of the type=ARH(1) but by playing a bit with the dataset in Example: Comparing Multiple B-splines, I am now pretty clear as to what was going on.

1. The RANDOM statement: RANDOM int time/type=ARH(1) sub=id; returns two variance components one for id and one for time and a correlation between them, which makes no sense as you pointed out.

2. Two RANDOM statements as I outlined will work, but oddly the default ddfm is Kenward Rogers. Between within returns 0 denominator degrees of freedom for the polynomial terms.

3. You cannot include an effect variable in a RANDOM statement. You get this error: ERROR: Separated or split RANDOM effects are not supported by the GLIMMIX procedure.

So I definitely learned some things this afternoon. Thanks to all of you.
. .

SteveDenham

My code (the dataset spline can be found in the examples):

data test;
set spline;
if _n_ <=25 then grp = 1;
if 26 =< _n_ <=50 then grp=2;
if 51 =< _n_ <=75 then grp=3;
if _n_ >= 76 then grp=4;
run;

proc glimmix data=test;
class group grp;
effect poly=polynomial(x/ degree=3);
effect spl=spline(x);
model y =group poly group*poly/noint;
random int x/type=chol subject=grp;
run;

proc glimmix data=test;
class group grp;
effect poly=polynomial(x/ degree=3);
effect spl=spline(x);
model y =group poly group*poly/noint ;
random int/suject=grp;
random x/type=arh(1) subject=grp;
run;

Re: Interpreting PROC GLIMMIX output

sld — Fri, 17 Jul 2020 04:18:18 GMT

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.

Re: Interpreting PROC GLIMMIX output

Nerdcy — Fri, 17 Jul 2020 13:19:16 GMT

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

Re: Interpreting PROC GLIMMIX output

SteveDenham — Mon, 20 Jul 2020 12:38:05 GMT

@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

Re: Interpreting PROC GLIMMIX output

sld — Tue, 21 Jul 2020 05:06:34 GMT

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.