R 混合模型：lme、lmer 或两者兼而有之？哪一项与我的数据相关，为什么？

机器算法验证 r lme4-nlme 混合模式
2022-03-23 20:15:35

我有一个数据是通过每年（1994-2015）在一个西非国家进行的森林清查获得的。从未管理的天然林中选择 10 个大小相等的地块（每个 1 公顷），然后识别和计数树木和灌木的种类。计算了丰度、香农、辛普森等生物多样性指数。我只选择了在所有 10 个地块中收集数据的 9 年，并且我丢弃了不完整的年份并将“年份”作为因素。

数据结构如下：

str(BIData)
'data.frame':   90 obs. of  9 variables:
 $ Year          : Factor w/ 9 levels "1994","1995",..: 1 1 1 1 1 1 1 1 1 1 ...
     $ Plot          : Factor w/ 10 levels "Bas Kolel","Bougou",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Richness      : int  8 21 13 14 8 10 6 10 8 20 ...
     $ Abundance     : int  286 1471 1121 466 242 97 250 790 208 2015 ...
 $ Shannon       : num  1.33 1.79 1.55 1.68 1.44 1.71 1.35 1.27 1.27 1.86 ...
     $ Simpson       : num  0.656 0.71 0.682 0.694 0.665 0.714 0.66 0.647 0.649 0.718 ...
 $ InverseSimpson: num  2.91 3.45 3.14 3.28 2.99 3.52 2.95 2.83 2.86 3.54 ...
     $ Topography    : Factor w/ 3 levels "Plateau","Slope",..: 3 1 1 3 3 2 2 2 3 1 ...
 $ Land_use      : Factor w/ 2 levels "Cultivated","Pasture": 1 2 2 2 1 1 2 2 1 2 ...

此外，地块位于不同的地形（坡地、山谷、高原）和土地利用（耕地、牧场）。

我在 lmer 和 lme 中有以下两个模型：

model=lmer(Abundance~Year+Topography+Land_use+(1|Plot), method="ML", data=BIData)
model=lme(Abundance~Year+Topography+Land_use, random=~1|Plot, method="ML", data=BIData)

我得到了完全不同的结果：我的问题？

我不是专家，但我发现它lme提供了一种带有 p 值的“美丽”结果。我可以看到许多重要因素，例如年份、地形和土地利用，而lmer只有 t 值没有 p 值。我不知道哪一个对我的数据是正确的。在这两种情况下，它都显示了良好且可接受的残差图。

请帮助我了解哪一个对我的数据是正确的。

谢谢@fcoppens。不，我没有尝试该参数。这是 和 的lme输出lmer。

lmer

model=lmer(Abundance~Year+Topography+Land_use+(1|Plot), method="ML", data=BIData)
summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: Abundance ~ Year + Topography + Land_use + (1 | Plot)
   Data: BIData

REML criterion at convergence: 1106.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5754 -0.5024 -0.0186  0.4015  3.4341 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept) 51753    227.5   
 Residual             48592    220.4   
Number of obs: 90, groups:  Plot, 10

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       1073.15     252.41   4.252
Year1995             0.40      98.58   0.004
Year1996           -32.70      98.58  -0.332
Year1998          -198.10      98.58  -2.010
Year1999          -341.90      98.58  -3.468
Year2002          -295.80      98.58  -3.001
Year2004          -324.90      98.58  -3.296
Year2010          -291.60      98.58  -2.958
Year2015          -371.00      98.58  -3.763
TopographySlope   -756.87     206.36  -3.668
TopographyValley  -645.82     236.71  -2.728
Land_usePasture    178.07     200.85   0.887

lme

model=lme(Abundance~Year+Topography+Land_use, random=~1|Plot, method="ML", data=BIData)
summary(model)
Linear mixed-effects model fit by maximum likelihood
 Data: BIData 
       AIC      BIC    logLik
  1264.675 1299.673 -618.3377

Random effects:
 Formula: ~1 | Plot
        (Intercept) Residual
StdDev:    171.5578 209.1232

Fixed effects: Abundance ~ Year + Topography + Land_use 
                     Value Std.Error DF   t-value p-value
(Intercept)      1073.1495  213.5506 72  5.025271  0.0000
Year1995            0.4000  100.4595 72  0.003982  0.9968
Year1996          -32.7000  100.4595 72 -0.325504  0.7457
Year1998         -198.1000  100.4595 72 -1.971938  0.0525
Year1999         -341.9000  100.4595 72 -3.403360  0.0011
Year2002         -295.8000  100.4595 72 -2.944469  0.0044
Year2004         -324.9000  100.4595 72 -3.234138  0.0018
Year2010         -291.6000  100.4595 72 -2.902661  0.0049
Year2015         -371.0000  100.4595 72 -3.693029  0.0004
TopographySlope  -756.8671  171.7008  6 -4.408058  0.0045
TopographyValley -645.8214  196.9543  6 -3.279041  0.0168
Land_usePasture   178.0654  167.1213  6  1.065486  0.3276
Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.6851599 -0.5159528 -0.0222693  0.4401886  3.6493837 

Number of Observations: 90
Number of Groups: 10
1个回答

看起来这些不是“完全不同的结果”；两个输出之间的固定效应系数相同，标准误差和随机效应的估计值略有不同。这并不奇怪，因为lmer拟合使用了受限最大似然 (REML)，而lme拟合使用了最大似然 (ML)。有关这些方法之间差异的介绍，请参阅此交叉验证页面。

输出中缺少plmer值是包的作者有意识的选择，正如包的文档和此交叉验证页面中所讨论的那样。论据是，在此类分析中通常为系数计算的p值具有误导性。lme4加载包后，该命令提供help("pvalues")了继续操作的指南，如当前小插图的第 35 页所述。

后续查询：

感谢您的回答。我可以看到，在两个输出中，固定效应系数是相同的。t 值几乎相同。虽然我对绘图（随机）效果不感兴趣，但结果并不相同（残差）。我的问题：（1）在 lme4 包（lmer）中，作者故意省略了一些警告的 p 值，那么 nlme pacakge 呢？(2) p 值是否来自 lme 可靠？这些 p 值符合我数据的趋势。(3) 我的数据是否有任何理由表明 lmer 比 lme 更相关？请？

回复：

这里有2个问题需要考虑。首先是模型拟合的 REML 和 ML 方法之间的选择。（您给出的“方法”参数lmer无法识别；要适应 ML，lmer您必须指定“REML = FALSE”，如@f_coppen 所述。您可以lme使用参数 'method="REML "' 或者只是省略其中的参数，lme因为 REML 是 ) 的默认值lme。REMLlmer拟合与 MLlme拟合几乎可以肯定地解释了估计随机效应的差异、估计的系数误差的差异以及由此产生的t 值差异。要在 REML 和 ML 之间做出选择，请参阅此交叉验证页面. 对于像这样的相对较小的研究（至少相对于您尝试估计的参数数量），您可能需要 REML。

第二个问题是如何计算p值。一个主要的潜在问题是在考虑了估计参数的数量之后，如何估计模型中剩余的自由度；您需要知道将t 值转换为p值的自由度。我想报告提供的p值没有任何问题，lme只要您指定您曾经lme获得它们。这些p值可能太低了，但读者至少能够知道如何思考它们及其局限性。help("pvalues")但是，使用包下推荐的方法lme4或附件包提供的方法可能会更好lmerTest. 无论哪种方式，重要的是要考虑潜在的问题，而不是简单地接受统计软件的标准输出。

后续查询：

太感谢了。你所有的建议都奏效了。我在没有指定方法的情况下以默认形式运行 lme 和 lmer 模型（REML 是默认形式），我发现结果几乎相同，除了与 DF 差异相关的 p 值略有不同。如果您有兴趣，我会发布结果吗？

model1=lme(Abundance~Year+Topography+land_use, random=~1|Plot, data=BIData) model2=lmer(Abundance~Year+Topography+land_use+(1|Plot), data=BIData)

其它你可能感兴趣的问题
上一篇时间序列的最佳预测窗口 下一篇R：检验与协变量的相关性？