2016-06-23 21 views
1

私は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 

答えて

2

モデルの選択プロセスは、どんなかなり標準的ではない、部分的に発明されたデータを使用して

再現性の例あなたがフィットするモデルのクラス、またはあなたが使うパッケージ。

種々の候補共変量を含む任意の回帰モデルについて

、逆方向選択がこの再帰ない:

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-squareモデルが最適です。

+1

ありがとうございます!あなたの答えは単純化プロセスに近づく方法の概要を示していますが、update()やdrop1()のような関数を使用しないと、2つ以上の予測変数を使用するとかなり面倒なことが起こるため、より単純なアプローチ実際のコーディングの – user2390246

+1

もちろん、ループを書くのは解決策ですが、一度やりとりしなければならないことはありません。したがって、私は簡単な方法があるかどうか尋ねます!しかし、はい、それは転倒の可能性があります。もちろん、誰にも私にそれを書くように求めているわけではありません:) – user2390246

+0

glmnetについてはわかりませんが、MCMCglmmのように、系統発生情報を含めることはできません。 – user2390246

関連する問題