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
- how the PhiPS2 value varies with habitat and months(sampling_no). And since I have different species, I think I should consider looking at
- 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