线性混合效应模型(nlme/ lme4)对边缘和条件R^2值的解释

Question

我正在调查细菌多样性，想知道这种多样性是否依赖于pH。我的数据结构如下：

• 我有17份不同研究的数据
• 每项研究都有多样性和pH的价值(研究之间的数值不同)。

现在我正在寻找一种方法来回答这个问题：“在所有的研究中，pH对多样性有影响吗？”

其思想是利用函数lme，将研究设为随机因子。从数据来看，它们似乎更适合于二次项而不是线性回归，所以我尝试用pH的二次项计算模型：

my_model<- lme( fixed = Bacterial_diversity ~ pH +
I(pH^2),
random = ~ pH |Study,
data= my_data)

产出(摘要)是：

> summary(my_model)
Linear mixed-effects model fit by REML
  Data: my_data
       AIC      BIC    logLik
  471.7855 497.7353 -228.8928

Random effects:
 Formula: ~pH | Paper
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 4.4808759 (Intr)
pH          0.4783127 -0.88 
Residual    0.4154606       

Fixed effects:  Bacterial_diversity ~ pH + I(pH^2) 
                 Value Std.Error  DF   t-value p-value
(Intercept)  1.6641091 1.8078372 285  0.920497  0.3581
pH           1.1750097 0.4670426 285  2.515851  0.0124
I(pH^2)     -0.1187954 0.0363455 285 -3.268508  0.0012
 Correlation: 
        (Intr) pH    
pH      -0.895       
I(pH^2)  0.763 -0.959

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-4.2091539 -0.4574020  0.1168270  0.6216446  2.0828655 

Number of Observations: 304
Number of Groups: 17 

我真的不知道如何解释这些结果。当我使用tab_model(my_model)时，我会得到两个R^2值。边缘R^2 (0.025)和条件R^2 (0.974)，给出的p值分别为0.0012 (对于pH的二次项)和0.0124 (对于pH)。是否计算整个模型的R^2值？我如何解释R^2值？我是否可以说，我的数据比线性回归更有可能拟合二次项(因为pH^2的p值更显著)？

甚至是我使用正确的模型，或者应该使用nlme而不是lme4 (我并不真正理解其中的区别！)。

我只想在不同的研究中展示pH与细菌多样性的关系。

如果你有什么想法的话，这会很有帮助的。如果有关于数据、代码或其他方面的问题，请随时询问。

Answer 1

摘要和引文

我建议使用performance包，甚至partR2包，它使用您提到的边缘和条件R2的Nakagawa R2值。我还建议阅读这三篇关于这个问题的论文，说明其价值：

• Nakagawa等人，2013年:原伪R2值论文。
• Nakagawa等人，2017年:伪R2值的扩展。
• Stoffel等人，2021: Part R2值。

拟合模型

我在下面给出了这些价值和解释的工作示例。首先，我在下面加载所需的包和模型。为了简单起见，我在这里删除NA值，但在实际场景中应该小心处理这些值。

#### Libraries ####
library(lmerTest) # for model fitting
library(performance) # for Nakagawa conditional/marginal R2
library(partR2) # for part R2 values
library(tidyverse) # for tidying data

#### Remove NA Values ####
carrots <- carrots %>% 
  drop_na()

#### Fit Model ####
fit <- lmer(Preference 
           ~ sens2 
           + Homesize 
           + (1 | Consumer),
           data=carrots)

边缘和条件R2

然后我们可以运行一个伪R2的简单摘要。

#### Run Pseudo R2 and Part R2 ####
r2_nakagawa(fit)

这给了我们这样的结果：

# R2 for Mixed Models

  Conditional R2: 0.176
     Marginal R2: 0.064

条件R2是整个模型的解释方差量。在这种情况下，固定和随机效应解释了大约17.6%的结果的变化。边际R2解释了这种差异有多大程度上仅归因于固定的影响。这里是一个相当小的数量：.064%。

个人效应:部分R2

然而，这并不能很好地反映个人的影响。为此，partR2可以帮助回答这个问题(在一定程度上)。注意，如果您有丢失的数据或与固定效果匹配的随机斜率，partR2函数将返回一个错误。我已经将引导设置为100，但对于最终模型，它们应该更接近1000 (我这样做只是为了节省时间和计算)。

part <- partR2(fit,
       partvars = c("sens2",
                    "Homesize"),
       nboot = 100)

在这里，您可以在这里得到固定效果大小的摘要：

summary(part)

见下文。您会注意到，它自动提供了一个边缘R2，它与sens2+Homesize的总partR2匹配

R2 (marginal) and 95% CI for the full model: 
 R2     CI_lower CI_upper ndf
 0.0643 0.0407   0.0964   3  

----------

Part (semi-partial) R2:
 Predictor(s)   R2     CI_lower CI_upper ndf
 Model          0.0643 0.0407   0.0964   3  
 sens2          0.0498 0.0262   0.0818   2  
 Homesize       0.0145 0.0000   0.0461   2  
 sens2+Homesize 0.0643 0.0407   0.0964   1  

----------

Inclusive R2 (SC^2 * R2):
 Predictor IR2    CI_lower CI_upper
 sens2     0.0498 0.0307   0.0729  
 Homesize3 0.0145 0.0016   0.0401  

----------

Structure coefficients r(Yhat,x):
 Predictor SC      CI_lower CI_upper
 sens2      0.8804  0.6980   0.9834 
 Homesize3 -0.4749 -0.7164  -0.1814 

----------

Beta weights (standardised estimates)
 Predictor BW      CI_lower CI_upper
 sens2      0.2235  0.1754   0.2703 
 Homesize3 -0.2800 -0.4716  -0.0914 

----------

在这里，我们可以看到，基于更高的解释方差，sens2对模型的贡献大于家庭大小。要获得更多信息，请参考我在上面所包含的文章。