私はRのパッケージMCMCglmmで実行している混合効果モデルを持っています。私が調査に興味を持っている可能性のある予測子はいくつかあります。 モデルから簡潔な情報を除外して予測モデルをモデルから除外するにはどうすればよいですか?他のコンテキストでMCMCglmmのモデル単純化はどのように実行できますか?
、私はつまり、私が興味を持って予測因子のすべてが含まれている「フル」モデルで始まり、そして無益変数を見つけることがdrop1()
を使用してからそれらを削除するにはupdate()
続い回帰後方段階的に使用していますモデル。
これらの機能のいずれも、MCMCglmmオブジェクトでは動作しません。私の最善の選択肢は何ですか?私は、drop1()
が動作しない1つの理由は、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
ありがとうございます!あなたの答えは単純化プロセスに近づく方法の概要を示していますが、update()やdrop1()のような関数を使用しないと、2つ以上の予測変数を使用するとかなり面倒なことが起こるため、より単純なアプローチ実際のコーディングの – user2390246
もちろん、ループを書くのは解決策ですが、一度やりとりしなければならないことはありません。したがって、私は簡単な方法があるかどうか尋ねます!しかし、はい、それは転倒の可能性があります。もちろん、誰にも私にそれを書くように求めているわけではありません:) – user2390246
glmnetについてはわかりませんが、MCMCglmmのように、系統発生情報を含めることはできません。 – user2390246