Каков наиболее подходящий метод для контрастов в lme (nlme)?
У меня есть вопрос, касающийся проведения специальных контрастов в смешанной модели glmm. По сути, мои данные представляют собой серию переменных ответа с фиксированными коэффициентами "Тип объекта" (три типа субъектов), "Обработка" (четыре курса) и "Год" (два года), с "Номер объекта" и " Местоположение 'как случайные эффекты.
Как я буду обсуждать оценки ниже, вот средства и стандартные отклонения по типу субъекта и лечение:
Моя модель, которая работает нормально, выглядит следующим образом:
library(nlme)
ctrl <- lmeControl(opt='optim')
myModel <- lme(y ~ Treatment * SubjectType * Year, random=list(SubjectNumber=~1, Location=~1), control=ctrl, method="REML", data = dframe1, na.action="na.exclude")
Затем я получаю следующий вывод, который выглядит нормально:
anova.lme(myModel)
numDF denDF F-value p-value
(Intercept) 1 168 3199.405 <.0001
Treatment 3 168 5.837 0.0008
SubjectType 2 168 7.502 0.0008
Year 1 168 9.410 0.0025
Treatment:SubjectType 6 168 0.933 0.4725
Treatment:Year 3 168 3.858 0.0106
SubjectType:Year 2 168 6.322 0.0023
Treatment:SubjectType:Year 6 168 1.383 0.2239
Тем не менее, теперь я хочу провести контрастные контрасты, чтобы увидеть, где находятся различия в основных эффектах "Лечение" и "Тип объекта" (т. Е. "Является ли тип объекта 1 отличным от обоих типов объектов 2 и 3, или просто 2/3"?".
Я провел немало исследований и попробовал несколько методов, но я, кажется, немного растерялся, и все дает мне разные результаты, поэтому некоторые разъяснения будут с благодарностью.
Сначала я использовал Tukey со следующим кодом, который меня порадовал, но мне посоветовали, что это не учитывает различия в других основных эффектах. Вот код и выходные данные, которые я использовал изначально, и lsmeans предоставлены для ссылки, когда я говорю об оценках различий позже в этом вопросе:
library(lsmeans)
lsmeans(myModel, pairwise~SubjectType, adjust="tukey")
NOTE: Results may be misleading due to involvement in interactions
$lsmeans
SubjectType lsmean SE df lower.CL upper.CL
SubjectType1 0.5531300 0.01853118 168 0.5165460 0.5897140
SubjectType2 0.6078394 0.01853118 168 0.5712554 0.6444234
SubjectType3 0.6545389 0.01853118 168 0.6179549 0.6911229
Results are averaged over the levels of: Treatment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
SubjectType1 - SubjectType2 -0.05470938 0.02620704 168 -2.088 0.0955
SubjectType1 - SubjectType3 -0.10140890 0.02620704 168 -3.870 0.0005
SubjectType2 - SubjectType3 -0.04669951 0.02620704 168 -1.782 0.1788
Results are averaged over the levels of: Treatment
P value adjustment: tukey method for comparing a family of 3 estimates
lsmeans(myModel, pairwise~Treatment, adjust="tukey")
NOTE: Results may be misleading due to involvement in interactions
$lsmeans
Treatment lsmean SE df lower.CL upper.CL
Treatment1 0.6379318 0.02139796 168 0.5956883 0.6801754
Treatment2 0.5299945 0.02139796 168 0.4877510 0.5722380
Treatment3 0.6123516 0.02139796 168 0.5701081 0.6545952
Treatment4 0.6403998 0.02139796 168 0.5981563 0.6826434
Results are averaged over the levels of: SubjectType
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Treatment1 - Treatment2 0.107937332 0.03026128 168 3.567 0.0026
Treatment1 - Treatment3 0.025580191 0.03026128 168 0.845 0.8327
Treatment1 - Treatment4 -0.002468028 0.03026128 168 -0.082 0.9998
Treatment2 - Treatment3 -0.082357141 0.03026128 168 -2.722 0.0358
Treatment2 - Treatment4 -0.110405360 0.03026128 168 -3.648 0.0020
Treatment3 - Treatment4 -0.028048219 0.03026128 168 -0.927 0.7905
Results are averaged over the levels of: SubjectType
P value adjustment: tukey method for comparing a family of 4 estimates
Оценки в этом действительно выглядят правильными на основе данных, и те, которые отмечены как значимые, на графике отличаются друг от друга. Я также провел аналогичное сравнение с пакетом 'multcomp', который дает разные результаты, и предполагаемые различия не выглядят корректными:
library("multcomp")
Contrasts2 <- glht(myModel, linfct = mcp(SubjectType = "Tukey"))
Warning message:
In mcp2matrix(model, linfct = linfct) :
covariate interactions found -- default contrast might be inappropriate
summary(Contrasts2)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
SubjectType2 - SubjectType1 == 0 0.26395 0.16575 1.593 0.249
SubjectType3 - SubjectType1 == 0 0.31642 0.16575 1.909 0.136
SubjectType3 - SubjectType2 == 0 0.05247 0.16575 0.317 0.946
(Adjusted p values reported -- single-step method)
ContrastsTreatment <- glht(myModel, linfct = mcp(Treatment = "Tukey"))
Warning message:
In mcp2matrix(model, linfct = linfct) :
covariate interactions found -- default contrast might be inappropriate
summary(ContrastsTreatment)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Treatment2 - Treatment1 == 0 -0.49321 0.16575 -2.976 0.01542 *
Treatment3 - Treatment1 == 0 0.08369 0.16575 0.505 0.95794
Treatment4 - Treatment1 == 0 -0.02263 0.16575 -0.137 0.99909
Treatment3 - Treatment2 == 0 0.57690 0.16575 3.481 0.00294 **
Treatment4 - Treatment2 == 0 0.47059 0.16575 2.839 0.02373 *
Treatment4 - Treatment3 == 0 -0.10631 0.16575 -0.641 0.91856
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Затем я попытался создать контрастную матрицу с glht, которая снова дает довольно разные результаты, и снова оценки не выглядят правильными на основе данных:
contrast.matrix.SubjectType <- rbind(
`SubjectType1 vs. Soanta` = c(1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`SubjectType1 vs. SubjectType3` = c(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`SubjectType2 vs. SubjectType3` = c(0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
compsSubjectType <- glht(myModel, contrast.matrix.SubjectType)
summary(compsSubjectType)
Simultaneous Tests for General Linear Hypotheses
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
SubjectType1 vs. Soanta == 0 0.33426 0.26207 1.275 0.401
SubjectType1 vs. SubjectType3 == 0 0.28180 0.26207 1.075 0.522
SubjectType2 vs. SubjectType3 == 0 -0.05247 0.16575 -0.317 0.945
(Adjusted p values reported -- single-step method)
contrast.matrix.Treatment <- rbind(
`Treatment1 vs. Treatment 1` = c(1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`Treatment1 vs. Treatment 2` = c(1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`Treatment1 vs. Treatment4` = c(1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`Treatment 1 vs. Treatment 2` = c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`Treatment 1 vs. Treatment4` = c(0, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
`Treatment 2 vs. Treatment4` = c(0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
compsTreatment <- glht(myModel, contrast.matrix.Treatment)
summary(compsTreatment)
Simultaneous Tests for General Linear Hypotheses
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Treatment1 vs. Treatment 1 == 0 1.0914 0.2621 4.165 < 0.001 ***
Treatment1 vs. Treatment 2 == 0 0.5145 0.2621 1.963 0.19465
Treatment1 vs. Treatment4 == 0 0.6208 0.2621 2.369 0.07935 .
Treatment 1 vs. Treatment 2 == 0 -0.5769 0.1657 -3.481 0.00268 **
Treatment 1 vs. Treatment4 == 0 -0.4706 0.1657 -2.839 0.02195 *
Treatment 2 vs. Treatment4 == 0 0.1063 0.1657 0.641 0.91584
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Изменение коррекции, очевидно, влияет на результаты, но, опять же, оценки выглядят неверно:
summary(glht(myModel, contrast.matrix.SubjectType), test = adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
SubjectType1 vs. Soanta == 0 0.33426 0.26207 1.275 0.202
SubjectType1 vs. SubjectType3 == 0 0.28180 0.26207 1.075 0.282
SubjectType2 vs. SubjectType3 == 0 -0.05247 0.16575 -0.317 0.752
(Adjusted p values reported -- none method)
summary(glht(myModel, contrast.matrix.Treatment), test = adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Fit: lme.formula(fixed = y ~ Treatment * SubjectType * Year, data = dframe1,
random = list(SubjectNumber = ~1, SubjectLocation = ~1), method = "REML",
na.action = "na.exclude", control = ctrl)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Treatment1 vs. Treatment 1 == 0 1.0914 0.2621 4.165 3.12e-05 ***
Treatment1 vs. Treatment 2 == 0 0.5145 0.2621 1.963 0.04961 *
Treatment1 vs. Treatment4 == 0 0.6208 0.2621 2.369 0.01784 *
Treatment 1 vs. Treatment 2 == 0 -0.5769 0.1657 -3.481 0.00050 ***
Treatment 1 vs. Treatment4 == 0 -0.4706 0.1657 -2.839 0.00452 **
Treatment 2 vs. Treatment4 == 0 0.1063 0.1657 0.641 0.52125
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- none method)
И, наконец, сравнение моделей в сводке модели, кажется, соответствует тому, что показывают данные, но опять же предоставляет другие значения p, и оценки, кажется, не соответствуют данным:
summary(myModel)
Linear mixed-effects model fit by REML
Data: dframe1
AIC BIC logLik
0.9785978 104.0406 26.5107
Random effects:
Formula: ~1 | SubjectNumber
(Intercept)
StdDev: 0.0001623009
Formula: ~1 | SubjectLocation %in% SubjectNumber
(Intercept) Residual
StdDev: 0.0001623009 0.2029986
Fixed effects: y ~ Treatment * SubjectType * Year
Value Std.Error DF t-value p-value
(Intercept) 0.5982163 0.11720132 168 5.104177 0.0000
TreatmentTreatment2 -0.4932126 0.16574769 168 -2.975683 0.0034
TreatmentTreatment3 0.0836861 0.16574769 168 0.504900 0.6143
TreatmentTreatment4 -0.0226273 0.16574769 168 -0.136517 0.8916
SubjectTypeSubjectType2 0.2639539 0.16574769 168 1.592504 0.1132
SubjectTypeSubjectType3 0.3164198 0.16574769 168 1.909045 0.0580
Year 0.0082773 0.07412461 168 0.111667 0.9112
TreatmentTreatment2:SubjectTypeSubjectType2 0.0807378 0.23440263 168 0.344441 0.7309
TreatmentTreatment3:SubjectTypeSubjectType2 -0.0983580 0.23440263 168 -0.419611 0.6753
TreatmentTreatment4:SubjectTypeSubjectType2 0.1477391 0.23440263 168 0.630279 0.5294
TreatmentTreatment2:SubjectTypeSubjectType3 0.3533741 0.23440263 168 1.507552 0.1335
TreatmentTreatment3:SubjectTypeSubjectType3 -0.2917698 0.23440263 168 -1.244738 0.2150
TreatmentTreatment4:SubjectTypeSubjectType3 0.0481905 0.23440263 168 0.205589 0.8374
TreatmentTreatment2:Year 0.2465353 0.10482803 168 2.351807 0.0198
TreatmentTreatment3:Year -0.0791173 0.10482803 168 -0.754734 0.4515
TreatmentTreatment4:Year -0.0326548 0.10482803 168 -0.311508 0.7558
SubjectTypeSubjectType2:Year -0.1651260 0.10482803 168 -1.575208 0.1171
SubjectTypeSubjectType3:Year -0.1671907 0.10482803 168 -1.594904 0.1126
TreatmentTreatment2:SubjectTypeSubjectType2:Year -0.0516816 0.14824922 168 -0.348613 0.7278
TreatmentTreatment3:SubjectTypeSubjectType2:Year 0.0718031 0.14824922 168 0.484341 0.6288
TreatmentTreatment4:SubjectTypeSubjectType2:Year -0.0043490 0.14824922 168 -0.029336 0.9766
TreatmentTreatment2:SubjectTypeSubjectType3:Year -0.2067818 0.14824922 168 -1.394826 0.1649
TreatmentTreatment3:SubjectTypeSubjectType3:Year 0.2071016 0.14824922 168 1.396983 0.1643
TreatmentTreatment4:SubjectTypeSubjectType3:Year 0.0218841 0.14824922 168 0.147617 0.8828
Итак, мой общий вопрос: что лучше всего использовать для определения сравнений между видами лечения / типами предметов? Почему два метода Тьюки дают разные результаты? Извините за столь многословный вопрос, и большое спасибо за ваше время и советы. Если бы ответы были даны в относительно простой терминологии, я был бы признателен! Спасибо!
1 ответ
Во-первых, я бы посоветовал вам взглянуть на emmeans Vignettes (ранее lsmeans). Есть несколько хороших примеров, которым вы можете следовать, чтобы лучше понять, что делает пакет и что означают результаты. Я считаю, что эти виньетки являются одним из лучших ресурсов, чтобы понять, как работает пакет. Что касается того, почему ваши результаты варьируются от метода к методу. У меня нет достаточных статистических знаний, чтобы попытаться ответить на этот вопрос.
Во-вторых, для моих моделей я следовал указаниям, приведенным в этом вопросе от StatsExchange, чтобы выполнить конкретные парные сравнения.