使用 lmer 进行预测：错误和缺失 p 值答案

【问题标题】：Predicting with lmer: errors and missing p-values使用 lmer 进行预测：错误和缺失 p 值
【发布时间】：2020-05-27 10:54:49
【问题描述】：

在一项实验中，雌性鱼被暴露在两个级别的光周期（环境和压缩）和两个级别的温度（4 和 7）中。它们在四个罐中（每个光周期两个罐，光周期内每个温度一个罐）。数据中有九个用 time_date 表示的样本。在其他响应中是“k”。我的兴趣是光周期、温度和时间日期对“k”的影响。面临的挑战：不平衡设计（在采样期间未采样一个光周期或温度水平），伪复制（每个罐是一个处理（光周期内的温度掩蔽））。通过一些阅读，我遇到了混合模型。我已经尝试过 lmer （更重要的是：我不确定是否正确）并且陷入了没有 p 值的警告和输出。我感谢您的帮助。提前谢谢你。

这是样本数据

fem.fish <- structure(list(time_date = structure(c(8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 
12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 
14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 
14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 
15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L), .Label = c("30-Jan-18", 
"11-Apr-18", "13-Jun-18", "07-Aug-18", "19-Sep-18", "30-Oct-18", 
"28-Nov-18", "03-Jan-19", "17-Jan-19", "31-Jan-19", "14-Feb-19", 
"28-Feb-19", "14-Mar-19", "27-Mar-19", "10-Apr-19", "24-Apr-19"
), class = "factor"), photo = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Ambient", 
"Compress"), class = "factor"), temp = structure(c(2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("4", 
"7"), class = "factor"), tank = structure(c(2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("T1", 
"T2", "T3", "T4"), class = "factor"), k = c(5.041791145, 5.408503999, 
5.535282299, 5.346402317, 5.376649977, 5.072021484, 6.097412109, 
4.390658006, 5.13676712, 4.472827193, 5.381892125, 4.882544582, 
4.655393586, 5.435528121, 4.985185185, 4.548431822, 5.041791145, 
5.408503999, 5.535282299, 5.346402317, 5.376649977, 5.072021484, 
6.097412109, 4.390658006, 5.13676712, 4.472827193, 5.381892125, 
4.882544582, 4.655393586, 5.435528121, 4.985185185, 4.548431822, 
5.517125816, 4.772205603, 5.928149807, 4.152323266, 4.666037968, 
4.638984928, 4.044444444, 4.720296599, 5.315500686, 4.967790359, 
3.520804755, 4.722326417, 5.051895044, 4.807450844, 5.096461818, 
5.28703008, 5.653368614, 6.357164944, 3.979492188, 3.928861374, 
5.632685221, 5.264668498, 5.281464786, 5.387205387, 4.332381668, 
5.250388878, 4.580237638, 4.650926114, 5.65951009, 4.401587625, 
5.194587481, 4.184813255, 4.44738449, 5.829977261, 4.331985587, 
4.827988338, 4.022222222, 3.672891297, 5.148148148, 4.068381688, 
5.71922963, 4.566763848, 5.330442907, 2.422536369, 5.346580575, 
4.971865289, 5.018922289, 5.513702624, 4.432146456, 5.692296224, 
4.738120151, 4.896057489, 5.50365439, 5.249023438, 5.737818961, 
4.260276996, 5.242507722, 4.580758017, 5.021888504, 5.013662642, 
4.308286338, 5.50840192, 4.732342764, 4.672289386, 5.715557782, 
3.827088497, 4.632069971, 4.935541824, 4.008746356, 4.963859809, 
4.836806618, 4.46244856, 4.839677641, 4.498269896, 4.88357943, 
4.984069185, 4.596844478, 5.196200195, 5.165529005, 14.74622771, 
5.397084548, 7.983198678, 5.691090246, 5.707491082, 5.187172012, 
6.297376093, 4.647178889, 4.282407407, 4.333496094, 4.773656052, 
4.770999725, 4.092207407, 3.917638484, 5.193905817, 3.704833984, 
5.571239611, 4.226680384, 3.65230095, 4.78515625, 5.603027344, 
4.159218067, 4.719370009, 4.437016946, 4.407713499, 4.284050303, 
4.676783265, 4.311689337, 4.540625, 4.864470022, 4.668176455, 
5.221193416, 4.997084123, 4.112752873, 5.587217586, 6.045051626, 
4.605417744, 4.35030714, 5.185252617, 4.752696927, 4.446670562, 
4.268256569, 4.30372087, 4.025205761, 5.696474074, 4.068342788, 
3.5212701, 4.544646911, 5.212620027, 5.31978738, 4.879910442, 
4.606482493, 4.33502906, 5.294067215, 5.770262391, 4.264308136, 
4.501028807, 2.944958848, 4.180638577, 4.120435057, 3.833076111, 
4.496793003, 4.232167131, 3.783896334, 5.070553936, 4.825776352, 
4.643534043, 6.318587106, 5.66205358, 5.194631597, 4.72557037, 
4.195096521, 4.956238551, 3.503093444, 5.24857851, 4.792524005, 
4.44229595, 5.285131195, 4.335878892, 4.170953361, 4.045779268
)), row.names = c(NA, -192L), class = "data.frame")

我尝试了什么和第一个警告

    fit1 <- lmer(k ~ 0 + photo*temp*time_date + (1|tank), data = fem.fish, REML = FALSE)

fixed-effect model matrix is rank deficient so dropping 12 columns / coefficients
boundary (singular) fit: see ?isSingular

我的总结和对相关矩阵的另一个警告

summary(fit1)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: k ~ 0 + photo * temp * time_date + (1 | tank)
   Data: fem.fish

     AIC      BIC   logLik deviance df.resid 
   551.2    635.9   -249.6    499.2      166 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7467 -0.4380 -0.0447  0.3663  9.7226 

Random effects:
 Groups   Name        Variance Std.Dev.
 tank     (Intercept) 0.0000   0.0000  
 Residual             0.7883   0.8879  
Number of obs: 192, groups:  tank, 4

Fixed effects:
                                         Estimate Std. Error t value
photoAmbient                            5.284e+00  3.139e-01  16.832
photoCompress                           4.937e+00  3.139e-01  15.728
temp7                                  -1.218e-14  4.439e-01   0.000
time_date17-Jan-19                     -9.116e-02  4.439e-01  -0.205
time_date31-Jan-19                     -9.798e-02  4.439e-01  -0.221
time_date14-Feb-19                      1.264e-01  4.439e-01   0.285
time_date28-Feb-19                     -3.986e-01  4.439e-01  -0.898
time_date14-Mar-19                      3.655e-01  4.439e-01   0.823
time_date27-Mar-19                     -3.979e-01  4.439e-01  -0.896
time_date10-Apr-19                     -4.122e-01  4.439e-01  -0.929
time_date24-Apr-19                     -2.184e-01  4.439e-01  -0.492
photoCompress:temp7                     8.874e-15  6.278e-01   0.000
photoCompress:time_date31-Jan-19       -2.957e-01  6.278e-01  -0.471
photoCompress:time_date28-Feb-19        1.575e+00  6.278e-01   2.509
photoCompress:time_date14-Mar-19       -6.073e-01  6.278e-01  -0.967
temp7:time_date17-Jan-19               -4.121e-02  6.278e-01  -0.066
temp7:time_date31-Jan-19                2.382e-01  6.278e-01   0.379
temp7:time_date14-Feb-19               -2.024e-01  6.278e-01  -0.322
temp7:time_date28-Feb-19               -1.441e+00  6.278e-01  -2.295
temp7:time_date14-Mar-19               -1.104e+00  6.278e-01  -1.759
temp7:time_date27-Mar-19               -4.306e-01  6.278e-01  -0.686
temp7:time_date10-Apr-19               -7.885e-01  6.278e-01  -1.256
temp7:time_date24-Apr-19               -5.872e-01  6.278e-01  -0.935
photoCompress:temp7:time_date14-Mar-19  9.077e-01  8.879e-01   1.022

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

fit warnings:
fixed-effect model matrix is rank deficient so dropping 12 columns / coefficients
convergence code: 0
boundary (singular) fit: see ?isSingular

我对 t 值的理解一点都不好，所以我无法确定是否存在显着影响，甚至无法确定相互作用是否显着。

我会感谢您对建模的建议（拟合正确的模型？）以及您认为有用的更多建议

非常感谢大家。

【问题讨论】：

你可以使用一些包来计算pvalue......但是......正如你从错误中看到的那样，你没有足够的观察来适应那个复杂的模型......你确定要做什么这？喜欢这个模型lmer(k ~ 0 + photo*temp*time_date + (1|tank) 合理吗？
您消除拦截是否有原因（在开始时指定 0+？）。您遇到了一些问题。一是你有一定水平的预测因子完全一致，导致排名不足警告。这意味着您的预测变量的某些级别没有进入模型。混合模型默认没有p值，一般建议用似然比检验比较两个嵌套模型
非常感谢您的建议。关于模型是否合理，我只是在阅读后尝试过（所以我对这些模型完全陌生-如果有更好的建议-我也可以考虑一下）。我已经尝试过 lmertest，现在我可以看到 p 值了。我已经删除了零，实际上我没有理由放它（非常感谢）。现在我正在尝试通过似然比测试——希望我能通过它得到一个解决方案——但是，欢迎提出可能更好的方法的建议。非常感谢您的宝贵时间。

标签： r lme4 mixed-models

【解决方案1】：

尝试导入“lmerTtest”包。

在拟合你的模型之前导入这个包，这样你就会看到p值和意义的“*”：

library("lme4")
library("lmerTest")

【讨论】：

非常感谢！我已经做到了，现在我可以看到它们了。
很好！所以请您将答案标记为有用。谢谢

【解决方案2】：

我已将您的数据用于以下示例。我认为由于你所有的术语都是分类的，你得到了排名不足的模型。我建议您使用时间作为连续预测器，从而摆脱排名不足的警告。

library(lme4)
library(parameters)
library(performance)

levels(fem.fish$time_date) <- 1:nlevels(fem.fish$time_date)
fem.fish$time_date <- as.numeric(fem.fish$time_date)

fit1 <- lmer(
  k ~ 1 + photo * temp * time_date + (1 | tank),
  data = fem.fish,
  REML = FALSE
)
#> boundary (singular) fit: see ?isSingular

关于奇异拟合的第二个警告（现在是第一个也是唯一的警告）是因为您在不同组中的结果实际上没有可变性（由tank 表示）。这意味着这里的随机效应模型不会比简单的线性模型带来更多好处。

ranef(fit1)
#> $tank
#>    (Intercept)
#> T1           0
#> T2           0
#> T3           0
#> T4           0
#> 
#> with conditional variances for "tank"

最后，您可以使用包 parameters 和 performance 获得全面的模型摘要（包括不同的 p 值近似值，如 Satterthwaite 或 Kenward-Roger、标准化参数或（集群）稳健标准误差）或模型拟合指数（如 r2）。

parameters::model_parameters(fit1)
#> Parameter                                 | Coefficient |   SE |        95% CI |     t |  df |      p
#> -----------------------------------------------------------------------------------------------------
#> (Intercept)                               |        5.56 | 0.61 | [ 4.36, 6.76] |  9.06 | 182 | < .001
#> photo [Compress]                          |       -1.46 | 1.04 | [-3.50, 0.58] | -1.41 | 182 | 0.160 
#> temp [7]                                  |        0.69 | 0.94 | [-1.16, 2.53] |  0.73 | 182 | 0.467 
#> time_date                                 |       -0.04 | 0.05 | [-0.13, 0.06] | -0.74 | 182 | 0.461 
#> photo [Compress] * temp [7]               |        0.73 | 1.52 | [-2.24, 3.70] |  0.48 | 182 | 0.631 
#> photo [Compress] * time_date              |        0.12 | 0.09 | [-0.06, 0.31] |  1.35 | 182 | 0.178 
#> temp [7] * time_date                      |       -0.09 | 0.07 | [-0.23, 0.05] | -1.29 | 182 | 0.198 
#> (photo [Compress] * temp [7]) * time_date |       -0.07 | 0.13 | [-0.33, 0.19] | -0.52 | 182 | 0.604

performance::r2(fit1)
#> Warning: Can't compute random effect variances. Some variance components equal zero.
#>   Solution: Respecify random structure!
#> Random effect variances not available. Returned R2 does not account for random effects.
#> # R2 for Mixed Models
#> 
#>   Conditional R2: NA
#>      Marginal R2: 0.088

既然您的时间变量是连续的，您可以考虑时间趋势的非线性关系。您可以使用 spline 包对此进行建模，并使用ggeffects 获取效果图。当然，这适用于具有线性时间趋势以及其他曲线时间趋势的模型。

library(ggeffects)
pr <- ggpredict(fit1, c("time_date", "photo", "temp"))
plot(pr)

library(splines)
fit2 <- lmer(
  k ~ 1 + photo * temp * bs(time_date) + (1 | tank),
  data = fem.fish,
  REML = FALSE
)
#> boundary (singular) fit: see ?isSingular

pr <- ggpredict(fit2, c("time_date [all]", "photo", "temp"))
plot(pr)

希望有帮助！

【讨论】：

非常感谢@Daniel?，当然可以！