How to interpret and visualise output from lmer model r?

Question

I have data from 310 individual plants, from various habitats (5 type), species (10) occuring in different habitats, measured photosynthesis performance (PhiPS2) in different months (sampling_no_).

I aim to look at

1. how the PhiPS2 value varies with habitat and months(sampling_no). And since I have different species, I think I should consider looking at
2. how it changes with different species in different habitats and months (sampling_no). I have PhiPS2 as a numeric continuous variable and habitat and sampling_no as categorical values.

I tried using a linear mixed-effect model: M5 <- lmer(PhiPS2~habitat*sampling_no + (1+ species|habitat), REML=F, data =physio2)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PhiPS2 ~ habitat * sampling_no + (1 + species | habitat)
   Data: physio2

     AIC      BIC   logLik deviance df.resid 
 -1338.5  -1058.2    744.2  -1488.5      235 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0414 -0.5610 -0.0908  0.5200  3.2800 

Random effects:
 Groups   Name                              Variance  Std.Dev.  Corr                                                 
 habitat  (Intercept)                       0.000e+00 0.000e+00                                                      
          speciesArtemesia brevifolia       1.136e-03 3.371e-02   NaN                                                
          speciesAster flaccidus            9.169e-04 3.028e-02   NaN  1.00                                          
          speciesLactuca tatarica           2.085e-12 1.444e-06   NaN -0.62 -0.62                                    
          speciesLeontopodium ochroleucum   3.035e-03 5.509e-02   NaN  1.00  1.00 -0.62                              
          speciesPotentilla pamerica        2.755e-03 5.249e-02   NaN -1.00 -1.00  0.62 -1.00                        
          speciesPrimula macrophylla        4.410e-04 2.100e-02   NaN -1.00 -1.00  0.62 -1.00  1.00                  
          speciesPsychrogeton andryaloides  8.150e-04 2.855e-02   NaN  1.00  1.00 -0.62  1.00 -1.00 -1.00            
          speciesSchistophyllidium bifurcum 3.639e-03 6.032e-02   NaN  1.00  1.00 -0.62  1.00 -1.00 -1.00  1.00      
          speciesWaldhemia tridactylites    6.230e-04 2.496e-02   NaN  1.00  1.00 -0.62  1.00 -1.00 -1.00  1.00  1.00
 Residual                                   4.576e-04 2.139e-02                                                      
Number of obs: 310, groups:  habitat, 5

Fixed effects:
                                  Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                      1.320e-01  3.627e-03  2.900e+02  36.398  < 2e-16 ***
habitatRuderal                  -1.361e-02  1.023e-02  3.068e+02  -1.330  0.18457    
habitatSemi-desert              -1.281e-02  1.023e-02  3.068e+02  -1.252  0.21166    
habitatSteppe                    6.537e-03  6.061e-03  2.920e+02   1.079  0.28165    
habitatSubnival                 -1.599e-02  7.683e-03  2.578e+02  -2.082  0.03836 *  
sampling_no2                    -1.571e-02  5.114e-03  3.070e+02  -3.073  0.00231 ** 
sampling_no3                    -2.060e-02  5.114e-03  3.070e+02  -4.028 7.09e-05 ***
sampling_no4                    -3.249e-02  5.114e-03  3.070e+02  -6.353 7.63e-10 ***
habitatRuderal:sampling_no2      2.431e-02  1.446e-02  3.070e+02   1.681  0.09378 .  
habitatSemi-desert:sampling_no2  2.791e-02  1.446e-02  3.070e+02   1.930  0.05454 .  
habitatSteppe:sampling_no2       3.954e-03  7.922e-03  3.070e+02   0.499  0.61804    
habitatSubnival:sampling_no2     2.763e-02  1.085e-02  3.070e+02   2.547  0.01136 *  
habitatRuderal:sampling_no3     -1.660e-02  1.446e-02  3.070e+02  -1.148  0.25200    
habitatSemi-desert:sampling_no3 -1.800e-02  1.446e-02  3.070e+02  -1.244  0.21428    
habitatSteppe:sampling_no3      -1.712e-02  7.922e-03  3.070e+02  -2.161  0.03147 *  
habitatSubnival:sampling_no3    -7.386e-03  1.085e-02  3.070e+02  -0.681  0.49649    
habitatRuderal:sampling_no4      3.049e-02  1.446e-02  3.070e+02   2.108  0.03587 *  
habitatSemi-desert:sampling_no4 -1.771e-02  1.446e-02  3.070e+02  -1.225  0.22162    
habitatSteppe:sampling_no4      -9.943e-04  7.922e-03  3.070e+02  -0.126  0.90021    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 19 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')````

````anova(M5)
Missing cells for: habitatSubnival:sampling_no1.  
Interpret type III hypotheses with care.
Type III Analysis of Variance Table with Satterthwaite's method
                      Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
habitat             0.006643 0.0016608     4 217.20  3.6290 0.0069352 ** 
sampling_no         0.042669 0.0142230     3 306.97 31.0786 < 2.2e-16 ***
habitat:sampling_no 0.016647 0.0015134    11 306.97  3.3069 0.0002599 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The habitat, sampling_no, and their interaction is significant. I have the following doubts:

Is the model I used correct for my purpose? if yes, how to interpret the result, and graphically represnt it? How do I model/graphically represnt the species level differences? I have not done any transformation of the response variable, is it needed? i am not sure.

i also tried other models such as

M1 <- lmer(PhiPS2~habitat+(1|species), REML = F, data =physio2) ##(not significant)
M2 <- lmer(PhiPS2~habitat+(1+habitat|species), REML=F, data =physio2) ##(not significant)
M3 <- lmer(PhiPS2~habitat +sampling_no + (1|species), REML=F, data =physio2) (habitat = slightly significant, sampling_no = significant)
M4 <- lmer(PhiPS2~habitat*sampling_no + (1|species), REML=F, data =physio2) ##(habitat = slightly significant, sampling_no = significant, interaction = significant)
M5 <- lmer(PhiPS2~habitat*sampling_no +(1+ species|habitat), REML=F, data =physio2) ##(habitat = significant, sampling_no = significant, interaction = significant)

Please suggest. Thanks

Answer 1

With data such as yours, it is quite common to use mixed effects models, so you are on the right track. However, the way you treat species in your model(s) is not something I have seen before. Please note that I work in psychology and education but I see a fair number of ecology posts on this forum. Many people would treat your data as hierarchical with repeated measures nested within species that are nested within habitats. This implies the following base structure for a model:

M0 <- lmer(PhiPS2 ~ 1 + (1|habitat) + (1|species), REML=T, data =physio2)

Note that I have included two random (or varying) intercepts - one for habitat and one for species. This is because one would expect that PhiPS2 measurements should be correlated among observations from the same species and from the same habitats. Consequently, variance in PhiPS2 is now partitioned into three sources - 1) among habitats, 2) among species within habitats, and 3) among observations within species and habitats (the residual). You can look at the proportion of variation in PhiPS2 at each of these levels - performance::icc(M0, by_group=TRUE). The other thing I've changed is that I turned on REML estimation. This is because you have a small number of habitats and REML is preferred when the number of clusters/groupings is small.

Now that you've correctly partitioned outcome variance, you can go about examining how PhiPS2 varies as a function of time (sampling_no). Your currently treat sampling_no as categorical - is that necessary? Might you be able to treat it as continuous and then model its relation w/ PhiPS2 as linear or quadratic? You will want to look at plots of PhiPS2 against time for different species and habitats to see if that is a reasonable approach. This simplifies the modeling relative to treating sampling_no as categorical - you get one (linear) or two (linear + quadratic) slopes vs. three. It makes interpretation a bit easier as well as estimation.

Either way, your next step is similar. You want to allow for a fixed (non-varying) population slope(s) for sampling_no and then random (varying) slopes for sampling_no by species and habitat:

M1 <- lmer(PhiPS2 ~ 1 + sampling_no + (1+sampling_no|habitat) +
      (1+sampling_no|species), REML=T, data =physio2)

This will lead to new parameter estimates of fixed slopes for sampling_no and then random effect variances (and covariances) for sampling_no slopes. As I said above, if you treat sampling_no as categorical, it could lead to estimation problems due to the large number of random effects that need to be estimated. You will need to be on the lookout for this. You might wonder how you can figure out whether allowing the slope(s) for sampling_no to vary across habitats and species leads to non-trivial improvements in model fit. For that, you can use the anova function to compare nested models. I would run the following set of models:

M0 <- lmer(PhiPS2 ~ 1 + (1|habitat) + (1|species), REML=T, data =physio2)
performance::icc(M0, by_group=TRUE)
M1 <- lmer(PhiPS2 ~ 1 + sampling_no + (1|habitat) +
      (1|species), REML=T, data =physio2)
M2 <- lmer(PhiPS2 ~ 1 + sampling_no + (1+ sampling_no|habitat) +
      (1|species), REML=T, data =physio2)
# test of whether sampling_no slope(s) vary by habitats ignoring species variation
anova(M2,M1) 
M3 <- lmer(PhiPS2 ~ 1 + sampling_no + (1|habitat) +
          (1+sampling_no|species), REML=T, data =physio2)
# test of whether sampling_no slope(s) vary by species ignoring habitat variation
anova(M3,M2)    
M4 <- lmer(PhiPS2 ~ 1 + sampling_no + (1+sampling_no|habitat) +
          (1+sampling_no|species), REML=T, data =physio2)
# test of whether sampling_no slope(s) vary by species, after accounting for habitat variation 
anova(M4,M3)