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, conserving only one replicate for Exclosure situation in Summer time. 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

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

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 (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...

Thanks in advance!!!

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...