Interpreting PROC GLIMMIX output

[Nerdcy]

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
time*carnegie	1			1.0983	0.1819	4661	6.04	<.0001
time*carnegie	2			-0.04387	0.2661	4661	-0.16	0.8691
time*carnegie	0			0	.	.	.	.
time^2*carnegie	1			-0.07228	0.01423	4661	-5.08	<.0001
time^2*carnegie	2			-0.01247	0.01982	4661	-0.63	0.5293
time^2*carnegie	0			0	.	.	.	.
time*flagship			1	-1.3081	0.2154	4661	-6.07	<.0001
time*flagship			0	0	.	.	.	.
time^2*flagship			1	0.08469	0.01663	4661	5.09	<.0001
time^2*flagship			0	0	.	.	.	.
statelog*time				-0.2541	0.04758	4661	-5.34	<.0001
statelog*time^2				0.01149	0.003387	4661	3.39	0.0007
statelog*time*carnegie	1			-0.2803	0.04672	4661	-6.00	<.0001
statelog*time*carnegie	2			0.01977	0.06687	4661	0.30	0.7675
statelog*time*carnegie	0			0	.	.	.	.
statelog*time^2*carnegie	1			0.01843	0.003647	4661	5.05	<.0001
statelog*time^2*carnegie	2			0.002585	0.004956	4661	0.52	0.6019
statelog*time^2*carnegie	0			0	.	.	.	.
statelog*time*flagship			1	0.3230	0.05414	4661	5.97	<.0001
statelog*time*flagship			0	0	.	.	.	.
statelog*time^2*flagship			1	-0.02112	0.004179	4661	-5.05	<.0001
statelog*time^2*flagship			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!

[SteveDenham · Accepted solution]

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

[Nerdcy]

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?

[sld]

Use the ilink option.

[Nerdcy]

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.

[sld]

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-gamm... ; you'll find additional discussions with an internet search.

[Nerdcy]

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. 

[sld]

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 

[Nerdcy]

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.

[SteveDenham]

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.

[Nerdcy]

Thank you!

[sld]

@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).

[Nerdcy]

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).

[sld]

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.

[SteveDenham]

@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;