我正在调查细菌多样性，想知道多样性是否依赖于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²的值呢?我能说我的数据更可能适合二次项而不是线性回归(由于pH^2的p值更显著)吗?
甚至是我使用的模型是正确的，或者我应该使用nlme而不是lme4(我真的不明白两者的区别!)
我只是想展示不同研究中pH值和细菌多样性之间的关系
如果你有什么想法，这会很有帮助。如果对数据、代码或其他任何问题有疑问，请随时提出。
摘要及引文
我建议使用performance包甚至partR2包，它们使用您提到的边缘和条件R2的Nakagawa R2值。我还建议阅读这三篇关于这个主题的论文，它们解释了价值是什么:
- Nakagawa et al.， 2013:伪R2值论文原件。
- Nakagawa et al.， 2017: Pseudo R2值的扩展。
- Stoffel et al.， 2021: 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解释了有多少方差单独归因于固定效应。这是一个相当小的数字:0.064%。
个体效果:Part R2
然而，这并不能说明个体的影响。为此，partR2可以帮助回答这个问题(在某种程度上)。请注意，如果缺少与固定效果匹配的数据或随机斜率，则partR2函数将返回一个错误。我已经将靴子设置为100，但对于最终模型，它们应该更接近1000(我这样做只是为了节省时间和计算)。
part <- partR2(fit,
partvars = c("sens2",
"Homesize"),
nboot = 100)
这里是固定效应大小的总结:
summary(part)
所示。您将注意到它自动提供了一个边缘R2，它与sens2+Homesize的总part2相匹配:
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对模型的贡献大于基于更高解释方差的房屋大小。要了解更多信息，请参阅我在上面包含的文章。