Model validity and specific contrasts in mixed model

Question

I have a design where mice are distributed in a two-way anova setup. With genotype (WT and KO) and treatment (Ctrl and Treat). From each mouse 5 different tissues have been extracted and the number of cells positive for a mark have been counted.

I am interested in testing if the ratio of positive to negative cells are affected by genotype and treatment and to find the effects of treatment within each genotype. All of these tests I want to do for each tissue.

The table used has this format:

Animal Genotype Treatment Tissue Positive Count
M1     WT       Ctrl      TA     Yes      25
M1     WT       Ctrl      TA     No       10
...
M12    KO       Treat     EDL    No       5

I have fitted the following model, using glmer from the lme4 package:

ratioMod <- glmer(Positive ~ 0 + Tissue + (Treatment * Genotype):Tissue + (Tissue | Animal), data = myData, weights = Count, family = binomial)

To find the effects of treatment and genotype I used the emmeans package:

emRes <- emmeans(ratioMod, specs = c("Treatment", "Genotype"), by = "Tissue")

This notices that there is nesting of Treatment and Genotype in Tissue, but that seems correct to me. I found the effects (and significance) of the main effects by:

emmeans(emRes, specs = c("Treatment"), by = "Tissue", contr = "pairwise")

emmeans(emRes, specs = c("Genotype"), by = "Tissue", contr = "pairwise")

This warns that NOTE: Results may be misleading due to involvement in interactions.

Finally I wrote my own little hack to find only the pairwise comparisons I was interested in (by subsetting the matrix generated by pairwise.emmc here called testFn), and found the interactions I was interested in by:

contrast(emRes, interaction = c("Treatment", "Genotype"), method = "testFn")

Is this the correct way to analyse this experiment? Is there anything I should check or worry about (such as the warning about interactions).

Answer 1

Welcome to CrossValidated!

The nesting you mentioned is suggested by the model, which does not contain a main effect for either genotype or treatment. However, I want you to consider whether you have the right model. My understanding from the beginning of the question is that treatment and genotype are the factors of primary interest, so it would be curious to have them be subclasses of tissue. It appears (contrary to an earlier draft of this answer) that tissue is also a fixed factor of primary interest, but I don't think there is any nesting relationship among these factors; they are all crossed.

I'd suggest the model

mod <- glmer(Positive ~ Treatment * Genotype * Tissue + 
    (1 | Animal), data = myData, 
    weights = Count, family = binomial)

This model includes the intercept, and allows for random intercepts for each mouse. The fixed effects are for Treatment, Genotype, Tissue, and their interaction. This would be a much more parsimonious and conventional model for this kind of experiment. Having random Tissue effects for each animal is a much more complex model with tons of parameters (the tissue variances and the pairwise correlations of all these tissue effects) and seems unnecessary. The predictions from this model (and also yours) are on the logit scale, which is defined as the log of the odds, where the odds is Pr(pos)/Pr(neg).

Post hoc analysis would depend on whether the Treatment:Genotype interaction is significant. If it's not significant, then you can do something like

emmeans(mod, pairwise ~ Treatment | Tissue)
emmeans(mod, pairwise ~ Genotype | Tissue)

If you add type = "response" to these, you will get the results on the odds and odds-ratio scales instead of the logit and difference-of-logits scales.

If the interaction is significant, just add by variables:

emmeans(mod, pairwise ~ Treatment | Genotype * Tissue)
emmeans(mod, pairwise ~ Genotype | Treatment * Tissue)

Also, I wonder if you actually tried that last contrast() call. It isn't right. I tested a similar example, and I got an error message that says

Error ... contrast 'Treatment.emmc' not found

This is because the interaction argument needs to be name(s) of contrast functions. You don't need to re-specify the factors, because they are already contained in the object emRes that you passed in the first argument.

I hope these comments help you find a way through this. I'm happy to follow-up if I misunderstand something.