#linear-regression #lme4 #mixed-models #emmeans #posthoc
Вопрос:
У меня есть набор данных измерений «Y» в разных местах, и я пытаюсь определить, как переменная Y зависит от переменных A, B и D, запустив lmer()
модель и проанализировав результаты. Однако, когда я достигаю этапа post hoc, я получаю ошибку при попытке анализа.
Вот пример моих данных:
table <- " ID location A B C D Y
1 1 AA 0 0.6181587 -29.67 14.14 168.041
2 2 AA 1 0.5816176 -29.42 14.21 200.991
3 3 AA 2 0.4289670 -28.57 13.55 200.343
4 4 AA 3 0.4158891 -28.59 12.68 215.638
5 5 AA 4 0.3172721 -28.74 12.28 173.299
6 6 AA 5 0.1540603 -27.86 14.01 104.246
7 7 AA 6 0.1219355 -27.18 14.43 128.141
8 8 AA 7 0.1016643 -26.86 13.75 179.330
9 9 BB 0 0.6831649 -28.93 17.03 210.066
10 10 BB 1 0.6796935 -28.54 18.31 280.249
11 11 BB 2 0.5497743 -27.88 17.33 134.023
12 12 BB 3 0.3631052 -27.48 16.79 142.383
13 13 BB 4 0.3875498 -26.98 17.81 136.647
14 14 BB 5 0.3883785 -26.71 17.56 142.179
15 15 BB 6 0.4058061 -26.72 17.71 109.826
16 16 CC 0 0.8647298 -28.53 11.93 220.464
17 17 CC 1 0.8664036 -28.39 11.59 326.868
18 18 CC 2 0.7480748 -27.61 11.75 322.745
19 19 CC 3 0.5959143 -26.81 13.27 170.064
20 20 CC 4 0.4849077 -26.77 14.68 118.092
21 21 CC 5 0.3584687 -26.65 15.65 95.512
22 22 CC 6 0.3018285 -26.33 16.11 71.717
23 23 CC 7 0.2629121 -26.39 16.16 60.052
24 24 DD 0 0.8673077 -27.93 12.09 234.244
25 25 DD 1 0.8226558 -27.96 12.13 244.903
26 26 DD 2 0.7826429 -27.44 12.38 252.485
27 27 DD 3 0.6620447 -27.23 13.84 150.886
28 28 DD 4 0.4453213 -27.03 15.73 102.787
29 29 DD 5 0.3720257 -27.13 16.27 109.201
30 30 DD 6 0.6040217 -27.79 16.41 101.509
31 31 EE 0 0.8770987 -28.62 12.72 239.036
32 32 EE 1 0.8504547 -28.47 12.92 220.600
33 33 EE 2 0.8329484 -28.45 12.94 174.979
34 34 EE 3 0.8181102 -28.37 13.17 138.412
35 35 EE 4 0.7942685 -28.32 13.69 121.330
36 36 EE 5 0.7319724 -28.22 14.62 111.851
37 37 EE 6 0.7014828 -28.24 15.04 110.447
38 38 EE 7 0.7286984 -28.15 15.18 121.831"
#Create a dataframe with the above table
df <- read.table(text=table, header = TRUE)
df
# Make sure location is a factor
df$location<-as.factor(df$location)
Вот моя модель:
# Load libraries
library(ggplot2)
library(pscl)
library(lmtest)
library(lme4)
library(car)
mod = lmer(Y ~ A * B * poly(D, 2) * (1|location), data = df)
summary(mod)
plot(mod)
Теперь мне нужно определить, какие переменные существенно влияют на Y. Таким образом, я бежал Anova()
из пакета car
(вывод вставлен здесь).
Anova(mod)
# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Y
# Chisq Df Pr(>Chisq)
# A 8.2754 1 0.004019 **
# B 0.0053 1 0.941974
# poly(D, 2) 40.4618 2 1.636e-09 ***
# A:B 0.1709 1 0.679348
# A:poly(D, 2) 1.6460 2 0.439117
# B:poly(D, 2) 5.2601 2 0.072076 .
# A:B:poly(D, 2) 0.6372 2 0.727175
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Это говорит о том, что:
A существенно влияет на Y
B не оказывает существенного влияния на Y
D существенно влияет на Y
Поэтому далее я бы провел специальный тест для каждой из этих переменных, но именно здесь я столкнулся с проблемами. Я пробовал использовать оба emmeans
multcomp
пакета и ниже:
library(emmeans)
emmeans(mod, list(pairwise ~ A), adjust = "tukey")
# NOTE: Results may be misleading due to involvement in interactions
# Error in if ((misc$estType == "pairs") amp;amp; (paste(c("", by), collapse = ",") != :
# missing value where TRUE/FALSE needed
pairs(emmeans(mod, "A"))
# NOTE: Results may be misleading due to involvement in interactions
# Error in if ((misc$estType == "pairs") amp;amp; (paste(c("", by), collapse = ",") != :
# missing value where TRUE/FALSE needed
library(multcomp)
summary(glht(mod, linfct = mcp(A = "Tukey")), test = adjusted("fdr"))
# Error in h(simpleError(msg, call)) :
# error in evaluating the argument 'object' in selecting a method for function 'summary': Variable(s) ‘depth’ of class ‘integer’ is/are not contained as a factor in ‘model’.
Это первый раз, когда я запускаю тест ANOVA/post hoc для lmer()
модели, и хотя я прочитал несколько вводных сайтов для этой модели, я не уверен, что тестирую ее правильно. Любая помощь будет признательна.
Ответ №1:
Если я правильно смотрю на данные, A-это переменная, имеющая значения 0, 1,…, 7. Теперь посмотрите на свою таблицу anova, где вы видите, что A имеет только 1 d.f., а не 7, как должно быть для фактора, имеющего 8 уровней. Это означает, что ваша модель принимает А за числовой предиктор, что довольно бессмысленно. Превратите А в фактор и перестройте его модель. Тебе повезет больше.
Я также думаю, что вы имели в виду, чтобы (1|location)
в конце формулы модели были случайные эффекты, взаимодействующие с некоторыми полиномиальными эффектами.
Комментарии:
1. Привет, я пытался внести эти коррективы, но я все еще получаю ошибки при попытке запустить специальные тесты.
2. Вам нужно заново подогнать модель. Затем снова запустите Anova (), потому что тесты изменятся, и подтвердите, что у A все еще нет 1 d.f.