It might be that be setting group 4 to be the reference, that the indexing gets shifted such that your ESTIMATE statements are now messed up somehow.  Try adding the 'E' option to print out the L matrix.  The other thing to recall is that the estimate statement combines multiple entries from the solution vector to calculate the estimates, so there may be some small differences in the estimated value, and certainly some differences in the standard error, which in the estimate is a standard error of a linear combination of the parameters, and so should be a bit larger. This will affect the t values and thus the associated p values.   SteveDenham ... View more

Your approach to homogeneity is exactly what I wish was done by validated programs, but it seems that there is a great attraction somehow to keep doing this and using it to decide whether a parametric or a nonparametric testwill be implemented.  The same goes for normality tests.  At least for linear mixed models (not so much for generalized linear mixed models, where it is just, well, stupid), the problem with tests for normality being so dependent on sample size.   So what I hear through all of this is that what the mixed model procedures need is some decent interval estimates for the variance components.  I know that I just ran something yesterday with the "cl" option in the RANDOM statement in an attempt to get the limits you are talking about, and it was seemingly ignored by PROC GLIMMIX, with no error or warning in the log or message in the output.   SteveDenham ... View more

@WillTheKiwi - very interesting.  Just out of curiosity, I ran the data through PROC GLIMMIX with no separation of the two groups as in your method to get a likelihood ratio test of variance homogeneity for the two groups using this code:   proc glimmix data=dat1 chol; class subjectid group; model y=group/ddfm=sat solution noint; random _residual_/group=group solution; lsmeans group/diff; covtest homogeneity; run; The test had a chisquared value of 2.46 and an associate p value of 0.1165 for the analysis with 10 per group.  Residual VC's for the two groups were control = 0.9428 and exptal = 2.7497.   Then I cranked this up to 100 per treatment group.  Now the heterogeneity was obvious and significant (p<0.0001).  The least squares means were much closer to the simulated base values of 0 and 5, while the variance components went to ~1 and ~5.  This approach always yielded positive VC's.  So my feeling here is that I have been approaching this issue in a substantially different way, where 0 or negative values lead to a nonpositive definite G matrix, at which point I consider re-parameterizing the random effects.  I will keep this work in my "I learned some more about mixed models" folder for future reference.  Thanks for the effort and the code!   SteveDenham   ... View more

This is a great question, and I'm sorry the datamethods group wasn't really helpful (need to remember that the group is R oriented for a lot of things). Without sounding cruel, could you please tell me how a negative value for a variance has meaning in the context of real data? That question, by itself, led to the development of REML methods, as the Henderson type I methods (which can give negative values) were severely biased in estimating breeding values (animal breeding was the first big adopter of mixed model methods, so far as I know). We know ML estimators are biased, and REML estimators less biased, and are asymptotically unbiased. So this is a learning experience for me - I can't wrap my head around how to interpret a negative variance. Either the true value is a boundary condition (=0) or the model is misspecified in some sense. OK, I'll put the soapbox away for a while and try answering the first part of your post.   A copy of the output (at least for the random effects and the model) would be really helpful in addressing your questions.  In return, I refer you to the following papers in hopes that it will help:   Stroup&Claassen Pseudolikelihood vs Quadrature    Kiernan SGF GLIMMIX categorical outcomes with random effects    Note that there is a clear contradiction between Kiernan and Stroup&Claassen regarding bias of estimates, when comparing REML (pseudo-likelihood) to ML (.Laplace or quadrature methods).  Personally, if my interest was primarily in fixed effects, I would use REML methods, and if primarily interested in random effects, I would use ML methods, and accept the inherent bias of MLE's as compared to REML estimates.  I suppose that is where the negative variances arise.  You might consider using some of the empirical shrinkage methods in that case, to see if it "improves" those estimates.   SteveDenham ... View more

Well, perhaps this will help.  I calculated the means of the numerator and denominator that you are fitting for the class variable Age.   I got this:   While this ignores the other factors in the model, this demonstrates that simply removing Age from the model (_type_ = 0) results in very different segmentation of the data.  Without Age in the model, count=1682.625, n = 10074.54, so the estimated probabilities, and the consequent differences, are based on the averaged Age response = 0.167.  With Age in the model: Age=0 bases the estimates on an average of 1976.25/10125.3333 = 0.195, while Age=1 is based on 1389/10023.75 = 0.139.  If everything were linear, the removal would make no difference.  However, the logistic model works on the logit scale, and logit(0.167) does not equal the mean of logit(0.195) and logit(0.139).  So it should not be surprising that the difference in differences end up slightly different, as the two Age categories do not receive equal weight.   This is, perhaps, an oversimplified way of looking at this.  Actually all of the logit estimates for the other terms in the model are changed to accommodate the difference due to Age, which is ignored when it is excluded from the model. SteveDenham ... View more

In the code presented, Age is still in the model.  It only has been removed from the lsmeans statement.  Since it is still in the model statement, the reduced lsmeans give equal weight to each of the two Age groups, so that the results are different.   If the code presented is not what you ran, then there are other issues to address.   SteveDenham ... View more

I am going to assume you have no background in numerical analysis, and in particular,numeric minimization.. The Marquardt method can be much more sensitive to the initial conditions than the Newton method.  So what, you say?  It means that if there are multiple minima for the function in question, that your initial parameter guesses had better be close to the values at the global minimum.  If not, then you could converge to a set of values that result in a larger error term.  The other issue that Newton's method is a bit better at handling highly correlated parameters, such as those you have here.  Does that help at all?   SteveDenham   PS.  There are some reasons that other optimizers are available while Marquardt's is not for other optimization-based procedures, but that is material I have not dealt with since graduate school, so I am liable to end up with completely wrong responses. You might get better answers if you posted the question in the Mathematical Optimization forum. ... View more

Maybe ridging will help, but to do it requires switching to NLMIXED.  Since I don't have access to the exp dataset, I haven't gone through this for any errors or odd behavior:   proc nlmixed data=exp technique=nrridg; parms RLTE=0.00132 B=0.5 C=0.8 A=0.08 T=3 TU=20 s2e=0.1; R=RLTE*predlcm**C*EXP(A*temp)*(1-EXP(T*(temp-TU))); s=sow**(1-B)-R*(1-B)*time; if s>0 then stx=s**(1/(1-B)); else stx=0; model stw ~ normal(stx,s2e); run; I assume the variables 'sow' and 'predlcm' are in the 'exp' dataset.  What concerns me a bit about it is the units for all of these. R should have units of RLTE*predclm, so 'sow' should have units like RLTE*predclm*time.  The issue is that sow is exponentiated, and so has to accommodate that.  It looks like the definition of 'stx' handles that with the exponentiation there equal to the reciprocal of the first term in the definition of 's'.  So now we have to deal with the second term in the definition of s after exponentiation.  I don't see how to make that rather complicated expression unitless, so that the units of 'stx' aren't some horrible mixture of things.  I presume this equation is standard in the literature for the digestive physiology of fish, or you wouldn't have chosen it.  I just think some rearranging of terms or something should help with the correlation which is responsible for the extended confidence interval.   SteveDenham   ... View more

The examples contain sample programs, data and output.   SteveDenham ... View more