The GLIMMIX procedure has a specific COVTEST for homogeneity of variance components.  Check the documentation for specifics on how to implement this, as it requires a group= option in the RANDOM statement.   Note that the test is a log likelihood ratio test, which in a fully fixed effect model is not as powerful as Levene's test.  But since GLIMMIX allows you to easily fit a heterogeneous variance situation, it's not an issue (at least in my field).   I also caution you about fitting a change from baseline response with baseline as a covariate. The model is algebraically identical to fitting the response with baseline as a covariate save that the coefficient for the covariate in the change from baseline is equal to the coefficient minus one for the straight response variable.   SteveDenham ... View more

Elapsed times often result in a gamma distribution.   SteveDenham ... View more

This paper served as the source for our current analysis of ocular data. Consider a design of 4 treatment levels, 6 time points and measures taken on each eye of an animal. Given that, you might consider the following approach using PROC GLIMMIX:   proc glimmix data=eye_data; class treatment time eye subject_id; model response = treatment|time/ddfm=kr2; random eye/subject=subject_id type=un; random time/residual subject=eye(subject_id) type=ar(1); lsmeans <means of interest>; lsmestimate <comparisons of interest>/adjust=simulate(seed=1) adjdfe=row; run; There is a lot to unpack here, but the critical assumption is that eye is a repeated measure on the subject, but only adds variability to the fixed effects. By specifying type=un, you will get a measure of the variance due to each eye and the covariance between them. GLIMMIX really does not care for doubly repeated measures on R side variables, so this has been our approach.    Using the Kronecker product covariance structure in PROC MIXED is another possibility. This code may be useful:   proc mixed data=eye_data; class treatment time eye subject_id; model response = treatment|teye|time/ddfm=kr2; repeated eye time/ subject=subject_id type=un@ar(1); lsmeans <means of interest>; lsmestimate <comparisons of interest>/adjust=simulate(seed=1) adjdfe=row; run; The final thing I wanted to thrrow in there is that Dunnett's adjustment is really directed to comparison to a single control group mean, and that is why we tend to use Edwards and Berry's simulate method.  If you ever have the Dunnett-Hsu method fail to converge (which may be the case for repeated measures designs), there is a message in the log to consider the adjust=simulate method.   SteveDenham       ... View more

I would think that PROC POWER for comparison of two binomials and some additional options ought to work.  Check the examples in the documentation.   SteveDenham ... View more

Noting that you "keep getting errors" without giving us what is found in the log is like cutting off one leg before running a race.  Please post the log, being sure to including any DATA step messages prior to the errors.   SteveDenham ... View more

You could run a cross-tab in PROC FREQ to identify avistn and grpn combinations where the resp values are all in one table of trt01pn*sitegr1.   For instance:   ods output crosstabfreqs=crosstabsfreqs; proc freq data = ana2 descending; by avisitn grpn; tables resp*trt01pn*sitegr1; run; Post processing the crosstabfreqs dataset should identify the separated combinations, which in turn could be removed from the PROC LOGISTIC analysis. SteveDenham   ... View more

Note my soapbox rant above regarding the WARNINGS. Truly the only way I know to suppress them is as follows:   WARNING: Stopped because of infinite likelihood.     There are several things that might cause this, but the most common is multiple observations for the subject for a given cell defined by the fixed effects, e.g. two identical timepoints for a subject in a repeated measures analysis.   WARNING: Unable to make Hessian positive definite.   This almost always is due to a model being more complex than the data will support. Unless you can collect more data, you should probably not include a model that results in this, no matter the convergence status. Odd things like infinite or zero degrees of freedom often turn up with this error.   WARNING: Did not converge.   This may because the iteration has stopped at the default maximum. Try adding a MAXITER= option to your PROC MIXED statement. If that does not solve the problem, there are some things to try, but they require shifting to PROC GLIMMIX to get at other optimization techniques or rescaling variables.   For a great discussion of this, see this paper by Kiernan, Tao and Gibbs: Tips and Strategies for Mixed Modeling with SAS/STAT® Procedures    SteveDenham ... View more

If the data come from a survey, you may want to change to one of the appropriate survey procedures such as SURVEYMEANS or SURVEYREG. SURVEYMEANS can provide geometric means directly, while SURVEYREG would require transformation/back-transformation. Shifting to the SURVEY procedures requires you to have knowledge of the complex survey design used to collect the data.   There are several ways to convert continuous variables such as age to categorical variables in a DATA step, but beware of loss of power/increase in standard error when doing so.  There is a lot of discussion in the literature about the dangers of categorization of continuous variables. From the short example you present, you would have to convince me that the difference between a 12 year old and a 13 year old (different categories) was markedly greater than between a 12 year old and a 6 year old (same category).   SteveDenham ... View more

If you want an adjusted geometric mean, I would suggest the following code:   proc glimmix data=NH.demo0118_phth_fv2; ods output lsmeans=geomeans; class timephth; model value=timephth creatinine_phth/dist=lognormal; lsmeans timephth/cl; run; data geomeans2; set geomeans; mean=exp(estimate); run; This will give you the geometric mean at the ARITHMETIC mean of creatinine_phth. You can also get confidence bounds by exponentiating the variables lower and upper in the geomeans dataset.   SteveDenham   ... View more

With regards to this question: Is there a rationale why Laplace integration is inappropriate here? It does allow model selection in this fashion, but with your new code Laplace provides F-stats/Pvalues for all terms except those with Time, and the error “Estimated G matrix is not positive definite” The missing F stats and the Estimated G matrix statements point to an issue that the model is too complex to fit with this amount of data and a single point quadrature. My normal reaction there is to drop back to RSPL methods, and sacrifice IC calculations. One thing you might try is method=quad(fastquad). I am not up to speed on using that option yet, so it didn't occur to me that it might help. When I try, I get an ERROR: Insufficient resources to perform adaptive quadrature with 3 quadrature points. METHOD=LAPLACE, corresponding to a single point, may provide a computationally less intensive possibility If you have a cloud account, you might be able to increase your MEMSIZE to use this method.   Given that, you could try selecting a covariance structure based on the generalized chi-square/DF ratio. Look for the structure that minimizes the absolute value of ln(generalized chi-squared/DF). That would give the structure that results in the smallest amount of over- or under-dispersion (and no, I don't have a reference for that method, so it is pretty much an ad hoc kind of thing. Any structure selected that way is crucially dependent on the distribution/link used and conditional on the fixed effect solution.   On to the covariate shrub as a mediator that serves as a surrogate for year - here I believe you are correct. I would proceed with caution though, as I doubt the number of shrubs is a complete surrogate for the year effect. Right now I am trying to fit a spline to the variable covar as maybe something better, but it is slow going.   SteveDenham ... View more

I think we can address most of your questions by rearranging some of the code. It appears to me that you have 3 fixed effect design factors implemented for each plot within each site. If that is the case, then the following would be my first attempt:   proc glimmix /*method = laplace*/ data = exper1; class SITE Herb veg time plot year; nloptions maxiter=5000; model count = herb|veg|time / dist=negbin link=log ddfm=kr2; random intercept plot/ subject=site; random year / subject= plot(site) type=cs; run; This nests plots within each year. The issue here that needs some clarification is whether the plots are identical within site from year to year. The code above assumes that they are. If not, the easiest thing to do would be to recode the plots so that they are unique across site and year. Since not all plots within a site. were measured on all years. Some changes made were to use the default pseudo-likelihood method of optimization rather than the Laplace integration, and to specify the type of covariance structure for year as compound symmetric as at most 2 years were available for estimation on any plot.   SteveDenham     ... View more

There is a bit of a right tail, as shown in the histogram and at the upper right part of the QQ plot, but definitely not enough to make the assumption of normality of the residuals untrue. One thing I do see is that it looks like one of your model factors has three levels, rather than the two you mentioned in the original post.   So here is a way to analyze this using time as a repeated factor at the residual level, using PROC GLIMMIX:   proc glimmix data=ild.all; class intervention visit study_id; model mean_steps =intervention visit intervention*visit; random visit /residual type=un subject=study_id; run; Adding a covariate is then easy to do by adding the variable to the MODEL statement. If the covariate is categorical, you will also need to add it to the CLASS statement.   SteveDenham   ... View more

Things have changed since the Lambert paper. GENMOD does not employ an EM algorithm any longer.   SteveDenham. ... View more

I can't find the original post, but I don't see any reason that a cumulative problit link would not be valid with DIST=MULTINOMIAL. Note that for either cumulative link, the first category is assumed to be the reference in this parameterization. If you want to compare somehow to a different level as a reference and still use a cumulative link, you will probably have to use an LSMEANS statement with the diff option, or (much better) use an LSMESTIMATE statement. Comparisons will be on the linked scale. If you want comparisons on the original scale, you ought to look into the NLmeans and NLest macros.  @StatDave has posted many times about this.   However, since I don't see the original post, and @PaigeMiller 's response refers to the use of a generalized logit, I think that maybe your data does not really reflect an ordinal scale, which would then require the use of a reference category.   SteveDenham ... View more

First off, you may want to check the parameterization of the model that you are using. The example at the link below uses a more "interpretable" parameterization that you might find also has better properties for your goal.   If you don't have a good idea from prior knowledge as to what the best starting value(s) might be, you might try a grid search over a range of reasonable values. You don't need the exact starting values this way, and if several combinations of starting values lead to the same optimum solution, you get what you are looking for without knowing a precise set of values.  This also can help you avoid converging to a saddle point rather than an extremum.  See the documentation for the PARAMETER statement in NLIN  https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.4/statug/statug_nlin_syntax12.htm and this example which goes into various parameterizations for an ED50 model  https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.4/statug/statug_nlin_examples04.htm    Good luck.    SteveDenham ... View more