Hi Katy,   Thank you for the update.    There are rules of thumb for the ratio of number of observations to number of estimated parameters that range from 10:1 to 30:1, or more or less. (It's a "rule of thumb", so that's what we can expect.) Regardless, with only 67 observations spread over an unknown but smaller number of horses, your current model is attempting too much (i.e., it is overfitted). I suggest that you try to find a more parsimonious model with fewer predictor variables.    You do want to assess whether data may be overdispersed with respect to the binomial distribution. But first you should consider your other RANDOM statement:    random int/ subject=horse type=ar(1); This statement specifies that horses are not independent--it specifies that horses are correlated according to an AR(1) structure, which makes no sense to me. What are you attempting to do with this RANDOM statement?   ... View more

Before diving into graphics for predicted values, you'll want to identify a correct model that works for your data. I have some concerns about both data and model.   Does your actual dataset have more than 24 horses and 33 observations? Your model is trying to estimate 13 fixed effects parameters, and you would need a lot more data for the model to do that well.   In the MadeUp dataset, there are two observations with Freq_PMN > TNCC. These observations are rejected by the binomial distribution because Freq_PMN must be <= TNCC. Is this a problem that appears in your actual dataset?   The model is regressing Freq_PMN/TNCC on resp (plus the other continuous predictors). Are you content to fit a random intercepts model and consequently to assume that the slopes of the regressions do not vary by horse (i.e., to not incorporate random slopes)? In the MadeUp dataset, only 5 of the 24 horses have multiple measures, ranging from n= 2 to 4, which is not enough information to fit random slopes. If the actual dataset has little by way of repeated measures on horses, then random slopes is a moot point.   I doubt that your RANDOM statements are actually doing what you intended them to do. The best resource for this topic is Walt Stroup's text https://www.crcpress.com/Generalized-Linear-Mixed-Models-Modern-Concepts-Methods-and-Applications/Stroup/p/book/9781439815120   One way to acquire predicted values for graphing is to score a new data set using the PLM procedure. See https://support.sas.com/resources/papers/proceedings10/258-2010.pdf and http://support.sas.com/kb/33/307.html and https://support.sas.com/resources/papers/proceedings15/SAS1919-2015.pdf Pay attention to the information about scoring a mixed model in this last link.   Hope this helps some.   ... View more

I would think of each subplot being divided into two sub-subplots, with each level of depth being "randomly" assigned to a sub-subplot. The model now is a split-split plot design:   proc mixed; class tillage CR N replication depth ; model SOC= tillage | CR | N | depth ; random replication(tillage) cr*n*replication(tillage); lsmeans tillage CR N / alpha=0.01 diff adjust=simulate ; lsmeans tillage | cr | n | depth; run; Note that I removed "replication" from RANDOM because based on your description, the whole plots are not clustered in blocks. "replication(tillage)" is the whole plot variance. "cr*n*replication(tillage)" is the subplot variance. The sub-subplot variance is residual.   I do not recommend applying Type I error adjustment (like the simulate method) to interaction means. (I am fine with it for main effects.) The adjustment will account for ALL pairwise comparisons among the interaction means, most of which you have no interest in. Adjusting for Type I error increases Type II error (and so, reduces power) and to control for more comparisons than necessary is an unnecessary reduction in power--unless you care only about the risk of Type I error and you care nothing about the risk of Type II error. For most of us, I seriously doubt that is the case.   In addition, interaction generally is more complicated than mere pairwise comparisons. I think you have to take a more thoughtful and context-specific approach to interpretation of interaction.  My usual approach is to identify post-hoc contrasts that are pertinent in the context of the study and apply Type I error control to those. But that's my opinion, and yours may differ.   Hope this helps.   Edit: It could be that "location" is a block, and that "location" is the same as "replication", in which would require a different RANDOM statement. The design is not clear from your description.   ... View more

The Poisson distribution is bounded by zero and has no upper bound. So if your count can be no greater than 5, the Poisson quite likely is not appropriate.   @Rick_SAS recommended the binomial distribution in a previous response in this thread. I agree with his suggestion. The binomial will produce an analysis of the proportion of activities engaged in: 0/5, 1/5, 2/5, 3/5, 4/5, 5/5.   ... View more

Cool! Let me know how it sorts out. ... View more

I'll play devil's advocate: Are you even the tiniest bit skeptical that Well 29, which was similar to the other two wells at the other 3 sampling dates, should be non-detectable at N16? What biological/hydrological mechanisms would explain that? (I don't need to know, but you do.)   Are there procedures in your protocol, from sampling to inserting into the thermocycler, where something might have gone awry? Are the duplicate samples processed entirely separately, or are they two aliquots of the same sample? Did you count other genes, and do any of them show a similar pattern?   Once you've resolved the data validity issues to your satisfaction, you can return to the statistical analysis. Making the zeros be anything "small" will produce the same thing you've seen so far and will trigger significant results for well, date and interaction. If you want to see the impact of the detectibility-replacement, run the analysis without N16. The heterogeneous variances model probably is about as good as is feasible; you won't see any difference in the residual plots, but you'll note differences in SEs for different dates and fit statistics.     ... View more

Hmm. Do you think that the gene was truly not there or there at very low levels, or do you think you might have had a problem with the process by which counts are acquired? What is the detectability limit for your process? You can clearly see that how you deal with non-detectable results has a big impact on the statistical analysis. This is not a context that I am not very familiar with, so I cannot offer you much advice about an analysis approach.   I'm struck by the low variability among the duplicate samples at N16, and I wonder whether your process was running well for these samples. ... View more

The lsmean for well 29 in N16 is oddly low, and I am suspicious about the validity of those data. In fact, I'm a bit suspicious of the rest of the N16 data. What do you know about that?     It is a CRD because each duplicate sample is "randomly assigned" to a level of well and a level of date. Obviously, we are taking liberties with "randomly assigned", and you have to decide whether treating this study as a quasi-experiment is appropriate for your purposes.  ... View more

I concur, the residuals look awful.   With the range of values taken by cDNA, the distributional characteristics probably are going to be a challenge. Try   proc glimmix data=first plots=(studentpanel boxplot(fixed)); class well date; model cDNA = well date well*date / dist=lognormal; lsmeans well date / diff adjust=tukey lines; lsmeans well*date / plot=meanplot(sliceby=well join cl) slice=(well date) slicediff=(well date) adjust=tukey; random _residual_ / group=date; run;   The RANDOM will fit a separate variance for each date.   Please post the entire output.   Generally (with some exceptions), if the interaction is significant, then main effects are nonsensical. Interaction implies that the story for main effect A is not the same at all levels of main effect B; a main effect for factor A applies the same story to all levels of factor B. So depending on the nature of the interaction, you might be looking at only the interaction. But first, you need a model that adequately meets assumptions.   This is the nature of analysis of real data 🙂   ... View more

Using duplicate samples as replicates, with well and date as categorical fixed effects factors, the statistical model is a two-way factorial (well x date) in a completely randomized design. Try something like   proc glimmix data=first plots=(studentpanel boxplot(fixed)); class well date; model cDNA = well date well*date / dist=lognormal; lsmeans well date / diff adjust=tukey lines; lsmeans well*date / plot=meanplot(sliceby=well join cl) slice=(well date) slicediff=(well date) adjust=tukey; run;   and see how the residual plots look with respect to the assumption of lognormal distribution. You could also try a Poisson or negative binomial using data on the count scale.  ... View more