hi Kyle @Laser_Taco_ ,   Just a few things that occurred to me over the weekend. In GLIMMIX, for ratios you should probably specify a distribution, and change the method to LAPLACE with tech=nrridg. So perhaps something like this:   Proc GLIMMIX data=Long_version method=laplace; nloptions maxiter=1000 tech=nrridg; class Pen Trt diet Block Feed instance; model y=Trt Feed instance Feed*instance/dist=beta ddfm=KR; random instance/subject=pen*diet type=ar(1); random feed/subject=pen type=chol; random block; lsmeans Feed/pdiff lines adjust=tukey; run; The first thing to do is to re-arrange the data into a long format, where there is one dependent value in each observation (here it is y). The other variables in each line would be pen, trt, diet (although this looks completely confounded with feed), block, feed and instance, which is indexed from 1 to 5. So instead of: PEN DIET Treatment Rep Feed pref1 pref2 pref3 pref4 pref5 1 1 3 1 800 68.7500 57.6471 67.8571 57.9545 78.9157 1 3 3 1 100 31.2500 42.3529 32.1429 42.0455 21.0843   You would have (for the first pen, and all other pens would follow this pattern: PEN DIET Treatment Rep Feed instance y 1 1 3 1 800 1 0.6875 1 1 3 1 800 2 0.576471 1 1 3 1 800 3 0.678571 1 1 3 1 800 4 0.579545 1 1 3 1 800 5 0.789157 1 3 3 1 100 1 0.3125 1 3 3 1 100 2 0.423529 1 3 3 1 100 3 0.321429 1 3 3 1 100 4 0.420455 1 3 3 1 100 5 0.210843  And so on. Note that all of the y values have been rescaled to the (0, 1) interval. So now there are two levels of feed for each pen and five levels of instance for each pen by diet combination.   If I were better informed about how to use it, I might recommend PROC COPULA for the analysis, but that is only based on what the examples look like. I thought I had a way to address this, but so long as each pref within a treatment is such that the values for the two diets within a treatment sum to 1, all of the above will have some major issues due to lack of independence. If the preference response was converted to a ratio or some other type of value, so that a single value is obtained per pen at each preference measurement (time), you should be able to improve the analysis.   SteveDenham     ... View more

I hope I understand this correctly. Each pen is randomly assigned to one of 3 treatments (A, B, C). Thus, a pen will have one value. I don't see the advantage of calculating a value for each as a proportion. Have you considered fitting the raw consumption values (which should be independent), and then calculating the ratios from the least squares means?   The other point to consider is that ratios probably do not yield normally distributed residuals, as the cut-offs at 0 and 1 may be a problem. If you are going to analyze the proportions as the dependent variable, a beta distribution for the values may be more appropriate. That would mean shifting to GLIMMIX or GEE.   And last, failure to converge. That very likely can be fixed by adding MAXITER= and MAXFUNC= options to the PROC MIXED statement, provided there are no other issues. Check both the .log and .lst file for notes, warnings and notifications that indicate that the algorithm has difficulties other than slow convergence.   SteveDenham ... View more

Just a quick comment on covariance structure selection. If you use any of the pseudo-likelihood methods, the information criteria probably should not be used for selection, as the pseudo-likelihood estimates aren't the same under various structures. Thus @StatsMan 's comments re LAPLACE or QUADRATURE. If you truly want to use pseudo-likelihood methods, then probably the best you can do for covariance structure selection is look at the Gener. Chi-Square / DF value, and pick the structure that has the least over- or under-dispersion. You should note that this measure will get closer to 1 the more variables are estimated, and there is no penalization for this as there is for the information criteria, so "Caveat emptor" - let the user (buyer) beware.   SteveDenham ... View more

Hmm, that does depend on what your objective is, but I will take a swing at it without knowing what your data is (site means, SDs and Ns, or raw data).   The point I want to address is whether to treat site as a fixed effect (inference is only to the sites included) or a random effect (inference is to all possible sites). All of this assumes using PROC GLIMMIX, as I am fairly certain you are not looking at data with normally distributed residuals if it is primary endpoint clinical data. So, considering SITE as a fixed effect, here is code that will yield a variance estimate for each site:   proc glimmix data=your_data; class site treatment; model response = treatment/dist= <insert appropriate distribution here>; random _residual_/group=site; run; For SITE as a random effect, you could try:   proc glimmix data=your_data; class site treatment; model response = treatment/dist= <insert appropriate distribution here>; random intercept/subject=site solution; run;  Using ODS output statements, you can get the variance estimates for each site from code like this. Weights can then be calculated as 1/variance, merged back against the original dataset, and the data refit using a WEIGHT statement. The Type 3 F test for this analysis reflects any differences in treatment effect after weighting by site.   I am not a clinical biostatistician, but the meta-analyses I have seen generally treat site as a fixed effect. You might want to talk with others about the choice.   SteveDenham     ... View more

Non-estimable main effect lsmeans in situations like this are usually due to "missing" data/cells. Run a PROC FREQ or PROC MEANS to see which combinations of innocyear*harvestyear*TreeSp are totally missing observations.   How to get around this issue? I come from the Milliken and Johnson school of Analysis of Messy Data, and what I learned is to fit a one-way model, with only the highest order interaction in the MODEL statement. ESTIMATE statements or even better LSMESTIMATE statements can be used to get marginal means (equivalent to lsmeans) for the effects currently listed as non-estimable.   On to the question about the zeroes. That is a result of the use of non-full rank (GLM) parameterization by PROC GLIMMIX. It is also why the solution vector is NOT what you want to examine for significance testing. Instead, the Type 3 F tests (not shown) and differences between LSMEANS/LSMESTIMATES provide that information.   Last thing, be sure you want to use the shrinkage estimates obtained through the Kenward-Rogers ddfm option. Looking at the degrees of freedom in the fixed effect solution vector, I believe you should be in good shape with its use.   SteveDenham ... View more

As @ballardw says, it comes from the t distribution part. I get 12.71 for a two-tailed t value of 0.05. With a standard error of 25, then a half interval is 12.71*25 =317.25, so your confidence bounds around a mean of 25 would look like -292.25 to 342.25. A more precise t value will yield what you are seeing.   SteveDenham ... View more

@PaigeMiller  That is an excellent point. So we now have two pretty good reasons to fit all of the data, rather than parts of it, and trying to make sense of the parts.   SteveDenham ... View more

If you fit a model with a covariate equal to the response variable, odd things can happen, mostly because you are fitting some sort of adjusted response - essentially the response being fit is now equal to response_variable - response_variable * beta_hat, where beta_hat reflects a slope estimate. Since your means are relatively close, I suspect beta_hat is near, but not equal to, zero.   SteveDenham ... View more

One thing I didn't see in @sbxkoenk 's response deals with the default number of iterations (20) that GLIMMIX uses when fitting. Add an NLOPTIONS statement setting the number of iterations at some value greater than 20, which is the default for GLIMMIX. A good value, in my opinion, is at least 500; if that doesn't solve the issue, there are many other adjustments to the algorithm that may help, almost all of which are listed in @sbxkoenk 's post..   Oh yeah, the statement to add is  NLOPTIONS maxiter=550; SteveDenham     ... View more

@Ksharp 's post hits the mark: intraclass correlation is a good measure of intra-rater reliability. It also points out that the kappa statistic is for categorical variables.  An example can be found at  https://documentation.sas.com/doc/en/statug/15.2/statug_mixed_examples02.htm    About two-thirds of the way through the example, it covers how to calculate an intraclass correlation.  And near the end, it fits a heterogeneous variances model, which can be used to calculate intraclass correlations for each level used in the GROUP= statement.   SteveDenham     ... View more

Have you looked at Stroup & Claassen's 2020 paper?   <PAYWALL> Stroup, W., Claassen, E. Pseudo-Likelihood or Quadrature? What We Thought We Knew, What We Think We Know, and What We Are Still Trying to Figure Out. JABES 25, 639–656 (2020). https://doi.org/10.1007/s13253-020-00402-6.   There is a tremendous amount of simulation involved and much discussion of bias and how it might be overcome. I may be misinterpreting your aims, but if you haven't read this paper you are missing out.   SteveDenham ... View more

Binomial distribution: bounded below by zero, above by 1 (inclusive at both bounds). It is the proportion (or rate) of subjects (or trials) having a given characteristic divided by the number of subjects (or trials) observed (N). The variance is a function of the proportion (=N*p*(1-p)). Poisson distribution: bounded below by zero (inclusive). It is the count of subjects having a given characteristic. The population variance is equal to the population mean.   SteveDenham     ... View more