Thanks, @SteveDenham, for the kudos.   I'm not seeing the strip-block idea, so I'd be intrigued to see what you might work up.   Greenhouse experiments can be more variable than one might hope or assume, and some greenhouses are better at controlled conditions than are others. So I would not necessarily assume that environmental conditions were similar between the two experiments, and different environmental conditions might differentially affect emergence and growth.   Experiment as a replication could be incorporated in a statistical model as a random effects factor--essentially a super-block of blocks. But then you are trying to estimate a variance among experiments based on only two observations, and blocks within an experiment play the role of subsamples and not true replications.   Experiment could be incorporated as a fixed effect, "randomly" assigned to blocks. This approach could incorporate tests of experiment by treatment interactions; one could argue that those tests are not valid because blocks are not independent replications of a level of "experiment". And experiment as a fixed effect could be difficult to interpret: experiment is just a catch-all for whatever varies among experiments (e.g., light quality, humidity, temperature, day length).   Separate analyses could be done on each experiment. Granted, the comparison of two experiments is subjective, but that seems commensurate with having only two experiments.   This is a nice paper that might indirectly address some of our discussion points: Casler( 2015 ) Fundamentals of Experimental Design: Guidelines for Designing Successful Experiments    ... View more

Would it be incorrect to generalize the analysis to a one-factor rbd despite having an incomplete factorial with two factors, inoculum and fungicide? It would fine. In fact, it would be an appropriate way to deal with the incomplete factorial if you also included appropriate contrasts to test (1) the effects of fungicide (Edit: including control when combined with inoculum), and (2) the effect of inoculum for the control. Edit: These contrasts could be a subset of pairwise comparisons among the 5 means.   If I'm understanding this right, by turning this into a one-factor analysis, my significant indicator (p-value) is the interaction term or close to an interaction term. Your interpretation is not correct. The p-value for treatment in the one-factor analysis addresses the alternative hypothesis that at least one of the five means is different than at least one other of the five means; it is not a test of interaction. With the 5 treatments in this study, you are not actually able to address "interaction" if we are defining interaction as an effect of inoculum that differs among fungicides (or vice versa: effects of fungicide that differ depending on presence or absence of inoculum). You can compare fungicides given inoculum; you can compare inoculum given control. I don't see that you have a pertinent experimental setup for any sort of statistical interaction of fungicide and inoculum.   Would it be improper to use the one-factor analysis for my interaction and the two-factor analysis for my main effects? Yes, it would be incorrect.   To answer your other questions about the bonus question. I did read that paper which I found to be very insightful, but got the impression that pseudoreplication was always a bad or improper form of replication. In my case I don't find this to be true. I would suggest that pseudoreplication of a binary measurement provides greater precision of the treatment. I'm sure there is some reason I am unaware of for why this shouldn't be included in my analysis, which I would be OK with because there are work arounds. I could simply average the binary response variable of the pseudoreplicates and use this average as the response variable. For example 10 of my 18 seedlings emerged in one of the treatments so I'd use 0.55 and ditch the "rep" as one my variables and simplify the analysis whilst providing myself a precise treatment response. Pseudoreplication is always bad. Subsamples are not bad, it's all in how you incorporate them into the statistical model; pseudoreplication means that you are using subsamples as if they are true replications. For your study, you have two options. (1) You can analyze binary data measured at the subsample-level and include a random effect for plot in addition to block, to accommodate the clustering of 18 measurements within a plot; the plot effect can be specified as Block*Fungicide*Inoculum as in the code I suggested previously. Or (2) you can composite the 18 binary responses for a plot into a binomial response (i.e., number emerged out of 18) which is now a plot-level response and so you no longer include a random effect for plot in the statistical model (although you do now need to consider the distinct possibility of overdispersion). This is like your suggestion of an "average" (which is a proportion), but the binomial distribution is theoretically preferable to analyzing proportion data as a continuous variable (e.g., transforming and then assuming normality).   See The arcsine is asinine: the analysis of proportions in ecology  Rethinking the Analysis of Non‐Normal Data in Plant and Soil Science    In response to some of your other questions, my replication would be metareplication where I am simply showing repetition and repeatability. These are greenhouse trials, same location and same year, so more of an exercise of repeatability across time than space perhaps. I understand ideally this would be done in another growing season and another location but I have some limitation. Is the suggestion here that I consider this block 9, 10...to 18 and just add this into the analysis? I was under the impression I would need to include this as another variable to show the "contrast" for lack of better words between the experimental replicate 1 versus 2. Just adding blocks sure would make the analysis easier, just want to make sure I'm following modern/conventional statistical procedures the best I can. I would not recommend treating the second experiment as additional blocks in the first experiment. That approach is not consistent with the experimental design, and it ignores the possibility that the effects of fungicide or inoculum may differ between experiments. (In randomized block designs, we assume that there is no block by treatment interaction.) With only two experiments, in a greenhouse, the first in early season and the second in later season (and where season may matter and is unreplicated), I probably would analyze each experiment separately and subjectively assess whether results are qualitatively similar.   One other question I had was in regard to my error terms. Is there a helpful source on determining these error terms dependent on my design and which factors I choose to run. I know if I do a one-factor rbd or a two factor rbd its just block. But for the binary study as a two-factor with replicates and the binary study if I just do one-factor or average the replicates is where I have some confusion. Generalized linear mixed models are complex. Start with the Stroup paper linked above, and then move to his text Generalized Linear Mixed Models: Modern Concepts, Methods and Applications . This text Analysis of Generalized Linear Mixed Models in the Agricultural and Natural Resources Sciences goes into much less detail than the Stroup text and may be more accessible initially.   I hope this helps move you forward.   ... View more

A few thoughts regarding whether to use a normal distribution for proportions, or a binomial distribution for (number of complications)/(number of cases): (1) I would think that both (number of complications) and (number of cases) would be integers, so I'm puzzled by your statement that the outcome is not an integer. But then I don't fully know your study. (2) The binomial approach is preferable if (number of cases) is variable, or if the proportion of complications is small (or large), or if the (number of cases) is small such that proportion of cases is relatively discrete. (3) Distributional assumptions apply not to the overall (meaning all observations combined) set of response data but rather to the response conditional on the predictors, which is why we often use residual diagnostics to assess distributional assumptions.   Good spotting on the use of time. The solution to that is to have two copies of time, one continuous and one categorical: data have; set have; time_c = time; run; proc mixed data=have; class hospitalid time_c; model complicated_rate=time/solution chisq; random intercept /subject=hospitalid ; repeated time_c / type=ar(1) subject=hospitalid; run; This combination of RANDOM and REPEATED with AR(1) is referred to as "AR(1)+RE"; see Littell et al. (2000) Modelling covariance structure in the analysis of repeated measures data for a nice description of the distinction between this structure and plain AR(1) (which you would implement by keeping the REPEATED statement and omitting the RANDOM statement). In my practice, with small datasets, I usually have estimation trouble with AR(1)+RE.   ... View more

You specify hospitalid in the CLASS statement, but subject = hospitalid_index in the RANDOM statement. What is the difference between hospitalid and hospitalid_index? How many levels of hospitalid are in your dataset?   Do your results change if you use hospitalid as the subject in the RANDOM statement with type=UN?   The RANDOM statement models a 2 x 2 covariance structure, with the intercept variance and slope variance on the diagonal, and the covariance of intercept and slope on the off-diagonal. Using type=ar(1) on the RANDOM statement is nonsensical and does not do what you want. Logically, the off-diagonal elements in the covariance matrix are not autocorrelated (although it's a moot point when there is only one covariance), and more importantly you do not want homogeneous (i.e., equal) variances for intercept and slope. The residuals from the regression of complicated_rate on time might be autocorrelated, but that would be addressed by the REPEATED statement in Proc MIXED.   Have you plotted complicated_rate versus time by hospitalid_index to visually assess whether there is a relationship with time (and whether it appears linear), whether variance among intercepts is appreciable, and whether variance among slopes is appreciable?   You likely would find SAS for Mixed Models: Introduction and Basic Applications by Stroup et al. informative, especially Ch 10 on random coefficients models. You'll notice that in their examples, they do not model autocorrelation of residuals; some people do, some people don't.   Leaping to your questions: Q1: No, using type=ar(1) is wrong. Q2: I would switch to Proc GLIMMIX and use the COVTEST statement: GLIMMIX documentation for COVTEST    I hope this helps move you forward.   ... View more

My first point is that the over-parameterization that I noted is due to overlap/redundancy in parameterization of the RANDOM and REPEATED statements in Proc MIXED (and equivalently, in the RANDOM and RANDOM/RESIDUAL statements in Proc GLIMMIX) for certain--actually most--covariance structures. This issue is addressed in the appendix of the link to the Littell, Henry and Ammerman (1998) that @SteveDenham provides. It is also addressed in the Littell, Pendergast & Natajaran (Statistics in Medicine, 2000) paper I referenced in the code I posted previously. In this latter paper, the authors note the distinction between an AR(1) structure and an AR(1)+RE structure, as does Walt Stroup (Generalized Linear Mixed Models (2013)) in Ch 14 (Fig 14.2 does not correctly illustrate the different structures, IMO).   My second point about whether we can model R-side covariance structures is addressed by Stroup in his text on p 435, where he writes "The primary ambiguity in repeated measures model building for non-Gaussian data occurs when we have a member of the two-parameter exponential family [e.g., gamma]. How do we model repeated measures for beta or negative binomial or other distributions in this family? .... For G-side models, it is not clear how the random ts(a)_ijk [i.e., time x subject(treatment) random effects] effect coexists with f(y|b)'s scale parameters." In the 10 years or so since Walt wrote his text, he has continued to revisit issues. For example, he and Elizabeth Claassen published a paper in the Journal of Agricultural, Biological, and Environmental Statistics just last month (https://doi.org/10.1007/s13253-020-00402-6) entitled "Pseudo-Likelihood or Quadrature? What We Thought We Knew, What We Think We Know, and What We Are Still Trying to Figure Out". So he may currently hold a different view about R-side covariance structure modelling with non-Gaussian distributions. Then again, he just retired and could be on to other things 🙂   My third observation is that a random coefficients model (RCM) with time continuous will be quite different from a model with time categorical.   My fourth observation is that although UN, CHOL, and ARH(1) are identical when you are modeling a RCM with random intercepts, random slopes, and one covariance (i.e., a 2 x 2 matrix), you would not want to apply an AR constraint on covariances in a covariance matrix that was larger than 2 x 2 because that would be silly. I'd stick with UN or CHOL, or in my practice VC (or UN(1)) because my datasets are always small, and I have never been able to fit a decent intercept-slope covariance.   Mathematical statistics is not my superpower; my explorations tend to be empirical. My practical experience is consistent with the point that Walt makes in his text, and with the Littell et al. papers.   ... View more

Like @SteveDenham , I played around with code that I thought was pertinent, and I’ll share my new vision with you. The code is attached. Ideally I would take the time to simulate data that would support a full modeling effort—something with random intercepts and random slopes and autocorrelation—but I just don’t have the time for that currently.   The code file includes some comments. I will not go into model comparisons in depth here; you can embark upon your own discovery.   The first three models, which differ only in type—CHOL, UN, or ARH(1)—are re-parameterizations of the exact same statistical model which estimates in different forms (i.e., different parameterizations) the variance among intercepts, the variance among slopes, and the covariance of intercepts and slopes. Clearly none of these models is modeling temporal autocorrelation of residuals.   So my first point is that @Nerdcy's  RANDOM statement is definitely not doing what the OP thinks it is or may want.   My second point is less definitive at the moment. I experimented with models that were intended to model the covariance structure of residuals from the regression of the response on time. My example dataset is not at all sufficient to assess these models adequately. But patterns emerged that are consistent with my mixed model experience: when you have both G-side and R-side random effects, some model formulations that might seem sensible will be over-parameterized. Refer to the Littell et al. reference in the code.   My third point is that this code exploration is based on models that assume a normal distribution for the residuals. Theoretically, you ought to be able to model autocorrelation of those residuals if you have enough data that are well-enough behaved. However, for a generalized linear mixed model, the ability to model R-side random effects (i.e., the covariance structure of residuals) may be limited or (as I understand it) non-existent. In the context of a random coefficients model, I don’t know that that’s a major deal breaker because you are focusing your attention on the variances among intercepts and slopes, which are subject-level statistics.   My fourth point is that these are complicated models. Both Steve and I have experience with mixed models, and we are still discovering things we didn’t know. Just because a model runs, you should not assume it is valid. And it is easy to under-estimate the intricacies of these models.   ... View more

Temporal autocorrelation is certainly possible. The structure of the autocorrelation might be well captured by AR(1), or other covariance structures might fit the data better (e.g., CS or TOEP, heterogeneous variances versions; or UN).   We might think of the responses as being correlated, but in a design with repeated measures on subjects, we apply these covariance structures to the residuals (R-side random effects)--what's left after the MODEL removes the fixed effect trend of time.   See G-Side and R-Side Random Effects and Covariance Structures  This bit of documentation does suggest that you can do R-side modeling in GLMM; you could try it and see if it worked. ... View more

I am flattered to be confused with @lvm 🙂   ... View more

@SteveDenham I'm good with your first RANDOM statement. I disagree with the second: I do not see why there would be temporal autocorrelation among the subject random effects. If intercepts for subjects are independent (as in the first RANDOM), would not their random effects for linear slopes also be independent?   I'm not sure what the default type is, but for sure you would get covariances among intercepts, linear and quadratic terms with random intercept time time*time / subject = id type=chol; or type=un   Edit: I might agree with your interpretation of the two reported variance components Var(1) and Var(2) if poly (the name of the polynomial effect) replaced time in the RANDOM statement. But I still don't like this use of ARH(1).   ... View more

As you note, you can inverse link estimates of means and predictions. However, you cannot inverse link parameter estimates (I presume that's the table you are referring to). In a GLMM, all the estimates and statistical tests are done on the link scale.   You might find Ben Bolker's response here helpful: https://stats.stackexchange.com/questions/431120/how-to-interpret-parameters-of-glm-output-with-gamma-log-link ; you'll find additional discussions with an internet search.   ... View more