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

[![enter image description here][1]][1]
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,
                   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 

<br />
**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). 

[![enter image description here][2]][2]

How does lme package manage this kind of problems? It seems that somehow it considers missing data... I revised [Pinhero & Bates (2000)][3] 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 ...

  [1]: https://i.sstatic.net/wA00h.png
  [2]: https://i.sstatic.net/JF7uV.png
  [3]: https://link.springer.com/book/10.1007/b98882