One thing to try before you change site to a fixed effect would be to drop the ddfm=kr2 option.  That worked for us just the other day.  You could then adjust for at least some of the bias by using the EMPIRICAL option in the PROC GLIMMIX statement.  If the model converges, be sure to check the denominator degrees of freedom associated with F tests and tests of differences in LSMEANS.  If it looks like the "whole-plot" effects are being tested with the residual degrees of freedom, you may have to specify the denominator degrees of freedom using the ddf= option.  Additionally, don't include 'site' as both a fixed and a random effect, especially if you have a large number of sites (as I did in the code I posted, and for which I apologize).   Another possibility, depending on how many treatments are being examined, is to use an infinite number for the denominator degrees of freedom by specifying the CHISQ option rather than the ddfm= option in the MODEL statement.  This will give both F tests and chi-squared tests for the effects in the model. With 52 weekly timepoints, the whole-plot df is large enough to consider this option.  If you go to monthly averages, it still may be a very good approximation.   SteveDenham ... View more

I would try the following code:   proc glimmix data=<weekly_summarized data> order=data; method=reml; class site id trt week; model meanweekly=site trt week trt*week / ddfm=kr2; random int / subject=site; random week/ residual subject=id(site) type=ar(1); run; this is close to your Option 2, but rearranged a bit. @StatsMan makes a good point about fitting weekly averages, but even that may require more memory than is available to you.  The V matrix including the random and residual effects may be too large, which might lead to using monthly averages, or using month as a by variable and doing 12 separate analyses on a monthly basis and combining what you get with the monthly average analysis.   I also switched the ddfm=kr to ddfm=kr2, as the latter is better "behaved".   SteveDenham   ... View more

If there is a reason for using a log10 transform for this field, then the backtransformed (10**estimate, 10**stderr) values are what you should probably report to be consistent with the literature.  If you truly need the expected value and the standard deviation of the expected value (standard error), there are two methods that we have used.  These are the delta method and the omega method.  The latter is spelled out in the PROC GLIMMIX documentation, referring to a lognormal distribution.  You can use the formulas there, after converting the log10 values to natural logs (multiply by a constant).  For the delta method, it turns out that Wikipedia has good formulas for E(X) and Var(X).   SteveDenham ... View more

It does.  That model is one we have used with pretty good success for safety studies with opthalmic drug administration.   SteveDenham ... View more

One thing that jumps out at me about the summing to get an exposure variable in a longitudinal study is that the shape of the exposure likely will have a profound effect on the response - suppose all the mass for exposure is at time=0 and no other exposure compared to equal exposure at the various time points such that the sum is identical in both cases.  That would almost certainly lead to different responses over time.   I don't know if this would work or not, but you might consider a multimember EFFECT, where the inputs in the EFFECT statement are the current exposure level and lagged values of the exposure level.  Run something that calculates the autoregression to get the number of lags you might need, stopping when the autocorrelation levels out (a sill in a semivariogram).  If there are a lot of time points, you might need to thin them for inputs.  I would start with patients that completed the trial, and then work on developing methods to deal with missing measurements, either at random or due to drop out.   SteveDenham ... View more

Start out with   repeated /subject=day(subject); and see what you get.  I know that seems overly simple, but at this point, you probably need something that should work where you can look at results to see if you have adequately captured what is going on.   SteveDenham   ... View more

That looks like it should work.  Alternatively, you could fit a Kronecker product (UN@CS, UN@AR(1), UN@UN) with eye and week as the factors and animal as the  subject.  For repeated in time, I would start with UN@AR(1) as the UN@UN may exceed your memory limits.   SteveDenham ... View more

Be aware that the same issue that @ballardw outlined can happen even with similar operating systems and the same version of SAS/STAT.  Check to make sure that your CLASS and MODEL statements align with each other and with the sort in the dataset you are analyzing.  We just ran into this issue and a sort, followed by a CLASS statement with the by variables from the sort in the same order, and a MODEL statement that reflects the design solved our differences.   SteveDenham ... View more

You need to use either the CONTROLL or CONTROLU difftype to get one-sided tests and confidence intervals.  Check the LSMEANS documentation in the Shared Concepts section of the documentation.   SteveDenham ... View more

One of the difficulties in finding effect size calculations in mixed models is that there is very little agreement on how to calculate them, and if they are meaningful at all when there are multiple variance components (i.e. it is difficult to calculate a pooled standard deviation to use, when that may depend on which class, school and location are utilized.).   Confidence bounds on parameters are easy enough, just add CL to the MODEL statement.  Confidence bounds on the difference between conditional means is a bit harder.  You'll need to use an LSMESTIMATE statement with a CL option, and correctly specify the means you are interested in.  That would be easier in GLIMMIX where there is a SLICEDIFF option for the LSMEANS statement, thus enabling to compare between levels of a factor of interest at each level of some other factor or combination of factors.   SteveDenham ... View more

Well, I went digging through the documentation and came across two things that might affect this.  The first is that the standard error of the mean is the square root of the variance.  Note that this is different to the square root of the variance divided by n, which is what I think you were expecting.  The second thing is how the variance is estimated.  The default is a Taylor series expansion, and it looks like the formulas do involve various values for n by cluster.  However, later on when the documentation talks about degrees of freedom it says: Taylor Series Variance Estimation For the Taylor series method, PROC SURVEYMEANS calculates the degrees of freedom for the t test as the number of clusters minus the number of strata. If there are no clusters, then the degrees of freedom equals the number of observations minus the number of strata. If the design is not stratified, then the degrees of freedom equals the number of PSUs minus one.   To me, that says that the divisor when calculating the variance is the number of clusters minus the number of strata.  If not stratified, then it is the number of primary sampling units minus one.  There is nothing about the sample size within clusters entering into this calculation.  So to match up the t test calculation to the variance calculation, I would expect you would have to use the same logic. I am really willing to believe I might be in error on this, but based on your example, that is what it looks like is happening.   SteveDenham   ... View more

I like @jiltao 's suggestion.  Another one would be to try a different technique as an optimizer. i believe that the default method is being used, which is probably QUANEW. Check the log - if it is not converging the error message could say something with regards to QUANEW, and the optimization information table where the technique should be listed.  Failure to converge might be the source of the missing estimates.  Another would be something pathological with the data, so that the Hessian runs into issues.  In either case, try:   nloptions tech=congra; The conjugate gradient method requires only that the first derivatives exist (same as QUANEW) but is appropriate when the objective function and gradient are faster to evaluate than the Hessian.  Watch out for local optima with this method - you may have to search a parameter grid to check that it doesn't get stuck in a local stationary point.   SteveDenham   ... View more

Yes, those plots look exactly like the ones I used to illustrate an interaction years ago when I taught an introductory stats class.   SteveDenham ... View more