Прогнозирование с помощью lmer: ошибки и отсутствующие p-значения - PullRequest
Новая
       137

Прогнозирование с помощью lmer: ошибки и отсутствующие p-значения

0 голосов
/ 27 мая 2020

В эксперименте самки fi sh подвергались воздействию двух уровней фотопериода (окружающей и сжатой) и двух уровней температуры (4 и 7). Они находились в четырех резервуарах (два резервуара для каждого фотопериода, по одному резервуару для каждой температуры в течение фотопериода). Всего было девять выборок, обозначенных в данных time_date. Среди других ответов стоит «k». Меня интересуют эффекты фотопериода, температуры и time_date на "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-значений совсем не очень хорошее, поэтому я не могу установить sh, есть ли существенные эффекты или даже взаимодействия значимы или нет.

Я буду ценю ваши предложения по моделированию (Подобрать правильную модель?) и многое другое из того, что вы считаете полезным

Большое спасибо всем.

Ответы [ 2 ]

0 голосов
/ 29 мая 2020

Я использовал ваши данные для следующего примера. Я думаю, что из-за того, что все ваши термины категоричны, вы получите модель с дефицитом ранга. Я бы посоветовал вам использовать время в качестве непрерывного предсказателя, тем самым вы избавитесь от предупреждения о недостаточном ранге. 

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"

Наконец, вы можете использовать пакеты параметры и производительность чтобы получить исчерпывающие сводные данные модели (включая различные аппроксимации 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)

Надеюсь, что это поможет!

0 голосов
/ 27 мая 2020

Попробуйте импортировать пакет «lmerTtest».

Перед тем, как подгонять вашу модель, импортируйте этот пакет, таким образом вы увидите p-значение и значимость «*»: 

library("lme4")
library("lmerTest")
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...