Как я могу выполнить упрощение модели для MCMCglmm?
У меня есть модель со смешанными эффектами, которую я использую в пакете MCMCglmm в R. Есть несколько возможных предикторов, которые я заинтересован в исследовании. Как я могу выполнить упрощение модели, чтобы исключить неинформативные предикторы из модели?
В других контекстах я использовал пошаговую регрессию назад, то есть, начиная с "полной" модели, которая содержит все интересующие меня предикторы, и используя drop1() найти неинформативные переменные, а затем update() удалить их из модели.
Ни одна из этих функций не работает с объектами MCMCglmm. Какая моя лучшая альтернатива? Я думаю, что одна из причин drop1() не работает, так как объект MCMCglmm не имеет значения AIC. Может ли значение DIC использоваться аналогичным образом? Я не знаю много о значениях DIC, но в приведенном ниже примере они, похоже, не ведут себя так, как AIC: после удаления переменной из модели DIC уменьшается более чем на 100.
За update()Ясно, что можно задать модель заново с нуля, но она включает в себя копирование / вставку довольно длинного вызова функции и становится еще более сложной с более длинными формулами модели, которые включают взаимодействия. Я задавался вопросом, есть ли более простая альтернатива.
Мне также интересно, есть ли причина, по которой не был написан метод для этих функций для работы с объектами MCMCglmm - является ли то, что я пытаюсь сделать, каким-то образом недействительным?
Воспроизводимый пример с использованием частично изобретенных данных:
set.seed(1234)
library(MCMCglmm)
data(bird.families)
n <- Ntip(bird.families)
# Create some dummy variables
d <- data.frame(taxon = bird.families$tip.label,
X1 = rnorm(n),
X2 = rnorm(n),
X3 = sample(c("A", "B", "C"), n, replace = T),
X4 = sample(c("A", "B", "C"), n, replace = T))
# Simulate a phenotype composed of phylogenetic, fixed and residual effects
d$phenotype <- rbv(bird.families, 1, nodes="TIPS") +
d$X1*0.7 +
ifelse(d$X3 == "B", 0.5, 0) +
ifelse(d$X3 == "C", 0.8, 0) +
rnorm(n, 0, 1)
# Inverse matrix of shared phyloegnetic history
Ainv <- inverseA(bird.families)$Ainv
# Set priors
prior <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = 1, nu = 0.002)))
# Run model
model <- MCMCglmm(phenotype ~ X1 + X2 + X3 + X4,
random = ~taxon,
ginverse = list(taxon=Ainv),
data = d,
prior = prior,
verbose = FALSE)
summary(model)
# drop1(model, test = "Chisq") # doesn't work
# update(model, ~.- X2) # doesn't work
model2 <- MCMCglmm(phenotype ~ X1 + X3 + X4,
random = ~taxon,
ginverse = list(taxon=Ainv),
data = d,
prior = prior,
verbose = FALSE)
summary(model2)
Полный вывод модели:
Iterations = 3001:12991
Thinning interval = 10
Sample size = 1000
DIC: 145.2228
G-structure: ~taxon
post.mean l-95% CI u-95% CI eff.samp
taxon 1.808 0.2167 3.347 17.88
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.6712 0.0003382 1.585 20.98
Location effects: phenotype ~ X1 + X2 + X3 + X4
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -0.6005279 -1.5753606 0.2689091 1000.0 0.198
X1 0.7506249 0.5220491 1.0033785 506.8 <0.001 ***
X2 -0.0339347 -0.2452923 0.2085228 900.7 0.712
X3B 0.6276852 0.1238884 1.1727409 1000.0 0.022 *
X3C 1.0533919 0.4244092 1.5095804 1000.0 <0.001 ***
X4B -0.0002833 -0.5099753 0.5103896 478.3 0.978
X4C 0.0723440 -0.4435520 0.6193011 1000.0 0.776
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Уменьшенный выход модели:
Iterations = 3001:12991
Thinning interval = 10
Sample size = 1000
DIC: 296.4755
G-structure: ~taxon
post.mean l-95% CI u-95% CI eff.samp
taxon 1.542 0.07175 3.1 26.33
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.8155 0.0004526 1.661 21.89
Location effects: phenotype ~ X1 + X3 + X4
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -0.606920 -1.380267 0.335225 1000.0 0.172
X1 0.750929 0.540667 1.011041 1000.0 <0.001 ***
X3B 0.633123 0.149282 1.243790 1000.0 0.032 *
X3C 1.047085 0.518664 1.570641 1000.0 <0.001 ***
X4B 0.007019 -0.528080 0.512161 509.2 0.998
X4C 0.063093 -0.490103 0.590965 1000.0 0.818
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
0 ответов
Процесс выбора модели достаточно стандартный, независимо от того, какой класс моделей вам подходит или какой пакет вы используете.
Для любой регрессионной модели с участием различных ковариат-кандидатов обратный отбор делает это рекурсивно:
1. fit a model with `p` covariates;
2. drop the least insignificant covariate based on p-value (in your case, the pMCMC value);
3. `p = p - 1` and go back to 1.
The process goes on until there is no insignificant covariates to drop.
Точно так же существует прямой выбор, который начинается с перехвата и непрерывно добавляет ковариаты. Я не буду объяснять это здесь.
Выбор модели также может быть основан на информационных критериях, таких как AIC / BIC. Для вывода MCMC используется DIC. Лучшая модель - та, с наименьшим DIC.
Что я предлагаю, так это то, что вы выполняете обратный выбор. Пока вы устанавливаете различные модели, запишите их значения DIC. В итоге вы получите таблицу, подобную этой:
Model_1 DIC_1
Model_2 DIC_2
. .
. .
В лучшем случае обратный выбор и выбор DIC дают одну и ту же модель. Это означает, что эта модель имеет убедительные доказательства по сравнению с остальными. Ну, ты золотой!
Вы также можете самостоятельно вычислить (откорректировать) R-квадрат, исходя из остатков, исходных данных и предполагаемой стандартной ошибки. Модель с самым большим R-квадратом - лучшая.