01-12-2016 05:03 AM
I'm trying to setup a two variables, all with ordinal response scores, doubly repeated measures program where the two variables are evaluated in a number of animals by 3 observers for 3 times (time0=now, time1=after some time form time0, time2=after some time from time1).
It is of interest to state if there is a significant effect of time on scores for the two varables.
I tried the following code (after having read suggestions by Steve Dunham) about similar topic
proc glimmix data=mydata pconv=1e-6; where variable='a'; class animal time evaluator; nloptions maxiter=100; model response = / s link=cumprobit dist=multinomial ddfm=kr2; random int time /s sub=evaluator type=ar(1); run;
Being unsure of correctness of the above code I kindly ask users for useful suggestions.
01-12-2016 11:09 AM
I might suggest something a little different: you need to accommodate animal as a subject as well as specifying time as a fixed effect in the model.
proc glimmix data=mydata pconv=1e-6; where variable='a'; class animal time evaluator; nloptions maxiter=100; model response = time/ s link=cumprobit dist=multinomial ddfm=kr2; random int /s sub=evaluator type=un(1);
random time /s sub=evaluator*animal type=ar(1); run;
This assumes that the time points are equally spaced and that the evaluators are independent. The analysis is a marginal analysis. You may want to use a conditional analysis in this case, where the results are conditional on the random effects, rather than averaged over the random effects. To get that, add METHOD=LAPLACE to the PROC GLIMMIX statement.
01-13-2016 06:55 AM
Dear Steve, thanks for the precious suggestions.
I have applied the syntax you suggested, but convergence was not obtained when covariance structures AR(1) and UN(1) was used.
I set both covariance structures to defaults and convergence was obtained.
What do you think about? Is there a way to get convergence applying covariance structures you suggested ?
I think that a cause of not getting convergence could the small number of animals (N=7).
Is there a significant impaiment of the model if default structures are used?
Here is the program that has worked
proc glimmix data=mydata pconv=1e-6; where variable='a'; class animal time evaluator; nloptions maxiter=100; model response = time/ s link=cumprobit dist=multinomial ddfm=kr2; random int /s sub=evaluator /*type=un(1)*/; random time /s sub=evaluator*animal /*type=ar(1)*/; run;
01-14-2016 02:04 PM
That generally means that you have insufficient data to estimate all the parameters in the G matrix--some come up as zero. This should not deter you from the analysis, as something like 'Final Hessian is not positive definite' would.
For models that aren't converging, one thing you may want to check up on is the behavior of the objective function--if it is behaving relatively smoothly, and you still haven't converged, you might try increasing the maxiter= value, or relaxing the convergence criteria.
01-15-2016 03:36 AM
proc glimmix data=mydata pconv=1e-6; class variable animal time evaluator; nloptions maxiter=100; model response = variable | time / s link=cumprobit dist=multinomial ddfm=kr2; random int /s sub=evaluator type=un(1); random time /s sub=evaluator*animal type=ar(1); run;
Unfortunately the convergence process is not smooth but I think that the information of interest, i.e., difference between times has been obtained (and it is robust between models) with models modified by you suggestions, thanks.
A further improvement could be to analyse all the two response variables (that should be correlated) as a bivariate response; the two variables are correlated and both have multinomial response but one has 4 categories, the other 5.
I have tried a possible analysis for this with responses with each lesser number of categories (I attach the program) but there is no convergence.
Thanks for further help, Fabio
01-15-2016 02:01 PM
The attachment wasn't there, but it is generally a good idea as a last resort, when all else fails, to convert a multinomial (or a pair of multinomials) to a binomial. Have you considered using a different optimizer than the default quasi-Newton? I have had pretty good luck with the ridged Newton-Raphson method for these distributions. To get this add:
to the NLOPTIONS statement.
01-18-2016 08:47 AM
If I convert to binary each of the response variables (class variable =
"variable",having levels "a" and "b") I can get convergence, even with the default
Please I would like to know if the following syntax is correct to analyze
simultaneously (multivariate analysis) the two response variables "a" and "b", each
being a binary response variable, as said above.
Thanks once again, Fabio
title 'bivariate analysis with binary response variables '; proc glimmix data=mydata; nloptions maxiter=100 ; class variable animal time evaluator; model response = variable | time / s dist=binary; random int /s sub=evaluator type=un(1); random time /s sub=evaluator*animal type=ar(1); lsmeans variable*time /diff; run;
01-22-2016 08:34 AM
I agree with this approach as the first attempt. The next step would be to figure out a way to incorporate any correlation between the levels of "variable". For an example, take a look at Joint Modeling of Binary and Count Data located here:
and adapt as needed.
(And then comes the deluge...)