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 | 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 |
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!
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
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
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
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
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?
Use the ilink option.
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.
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.
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.
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
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.
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.
Thank you!
@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).
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).
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.
@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;
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.