Are mixed model results valid when several missing replicates?

Question

I want to know if it is correct to take as valid the results of a mixed model (lme) test for a triple factor experiment with several missing replicates in only one level factor situation.

My objective is to test the effect of grazing exclosure over a soil property in a silvopastoral system (along different seasons and tree cover situations). I am mainly interested in the existence of Exclosure effects; the estimates are secondary for me.

The 3 factors of my proposed experiment are (with 3 replications):

• Grazing (2 levels): Exclosure / continuous (main interest factor)
• Tree_cover (2 levels): Below trees / Between trees
• Season (2 levels): Summer / Winter

The problem is that I've lost some samples during my summer trip, conserving only one replicate (of 3) for Exclosure situation in Summer time.

    Season Random_site    Tree_cover    Grazing soil_property
1   Summer           1   Below trees Contiunuos         7.396
2   Summer           1   Below trees  Exclosure            NA
3   Summer           1 Between trees Contiunuos         8.612
4   Summer           1 Between trees  Exclosure            NA
5   Summer           2   Below trees Contiunuos         6.942
6   Summer           2   Below trees  Exclosure         8.661
7   Summer           2 Between trees Contiunuos        13.795
8   Summer           2 Between trees  Exclosure        15.768
9   Summer           3   Below trees Contiunuos         5.702
10  Summer           3   Below trees  Exclosure            NA
11  Summer           3 Between trees Contiunuos         7.393
12  Summer           3 Between trees  Exclosure            NA
13  Winter           1   Below trees Contiunuos         6.702
14  Winter           1   Below trees  Exclosure         7.421
15  Winter           1 Between trees Contiunuos         5.058
16  Winter           1 Between trees  Exclosure         5.886
17  Winter           2   Below trees Contiunuos         8.596
18  Winter           2   Below trees  Exclosure         9.714
19  Winter           2 Between trees Contiunuos         5.657
20  Winter           2 Between trees  Exclosure        14.918
21  Winter           3   Below trees Contiunuos         7.722
22  Winter           3   Below trees  Exclosure         6.941
23  Winter           3 Between trees Contiunuos         5.436
24  Winter           3 Between trees  Exclosure         7.897

In the figure can clearly be noted the problem (asterisks represents the measurements and missing data situation is rounded in red).

Despite this problem I tried to fit a mixed model (lme) with heteroskedasticity along the "Tree cover" factor (observed in residuals and levene test).

The mixed model (lme) is: response ~ Grazing * Tree_cover * Season

model <- nlme::lme(soil_property ~ Grazing*Tree_cover*Season,
               random = ~1|Random_site,
               weights = varIdent(form = ~1|Tree_cover),
               data = df,
               na.action = na.omit,
               method = "REML")

Anova test gives me a not significant (but very close to be) Grazing:Tree_cover interaction:

car::Anova(model)

Analysis of Deviance Table (Type II tests)

Response: soil_property
                            Chisq Df Pr(>Chisq)   
Grazing                    2.6069  1   0.106399   
Season                     0.2685  1   0.604328   
Tree_cover                 2.4448  1   0.117914   
Grazing:Season             0.3711  1   0.542394   
Grazing:Tree_cover         3.7907  1   0.051537 . 
Season:Tree_cover          7.9251  1   0.004875 **
Grazing:Season:Tree_cover  0.0000  1   0.995744   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

But if I run the "post hoc" test, it tells me that there is actually a grazing effect under the "Between trees" situation:

emmeans::emmeans (model, pairwise ~ Grazing | Tree_cover)

Tree_cover = Below trees:
 contrast               estimate    SE df lower.CL upper.CL t.ratio p.value
 Exclosure - Continuous    0.722 0.635 10  -0.6928     2.14 1.137   0.2819 

Tree_cover = Between trees:
 contrast               estimate    SE df lower.CL upper.CL t.ratio p.value
 Exclosure - Continuous    4.565 2.007 10   0.0927     9.04 2.274   0.0462 

Results are averaged over the levels of: Season 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

NOW THE QUESTION: Should I disregard this result (because of the several missing data, which seems not to be random but in fact actually is) or could it be considered as valid?

My doubt arose when I saw the graph of the adjusted mixed model (lme), which assigned greater dispersion to the situation where I lost repetitions, giving me the impression that the model considers this problem (as it can be seen in the next figure, red arrows points out the wider error bars under missing data situations).

How does lme package manage this kind of problems? It seems that somehow it considers missing data... I revised Pinhero & Bates (2000) book, but did not found anything about this issue ...

PD: I am aware that this issue have some problems and it seems me to be forcing it to have results. I just want to receive some advice about it whether to decide to discard the experiment or to report this results ...

Answer 1

Mixed models provide valid results when you have missing data in your outcome that are of the missing completely at random and missing at random type. The former assumption means that the reasons why you have missing data is completely unrelated to the outcome you are measuring. The latter assumption means that the reasons why you have missing data you have in a site are related to the observed outcome data in this site.