leave-one-outからp値を取得するR

Question

私は96の観測（患者）と1098の変数（遺伝子）のデータフレームを持っています。応答はバイナリ（YとN）であり、予測子は数値です。私はleave-one-outの相互検証を実行しようとしていますが、私の興味は標準誤差ではなく、LOOCVから作成された95個のロジスティック回帰モデルのそれぞれの変数のp値です。これらは、これまでの私の試みです：

#Data frame 96 observations 1098 variables 
DF2 

fit <- list() 

for (i in 1:96){ 
    df <- DF2[-i,] 
fit[[i]] <- glm (response ~., data= df, family= "binomial") 
} 
model_pvalues <- data.frame(model = character(), p_value = numeric()) 

これは、16個の要素と30のリストを持つ大規模なリストとしてフィット出力：$係数、$残差、$ fitted.values ....

試み1 ：見積もり、STD： "model_pvalues" 95人の観察（切片及び94個の変数）と4つの変数に

for (i in length(fit)){ 
    model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]]))) 
} 

この出力。エラー、z値、Pr（> | z |）。しかし、私が実際に得ようとしているのは、1097変数すべてのp値です。

試み2：

for (i in length(fit)){ 
    model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]]))[4]) 
} 

私はこれを実行したとき、私は一つの変数のために（ベータ版を想定してどこからかわからない、）1つの番号を取得します。

試み3：推定、STD：私はこれを実行すると

for (i in 1:96){ 
    df <- DF2[-i,] 
    fit[[i]] <- glm (response ~., data= df, family= "binomial") 
    model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]]))) 
} 

が、私は4つの変数の1520回の観測のデータフレームを出します。エラー、z値、Pr（> | z |）。観測は（迎撃）から始まり、その後に82の変数が続きます。その後、このパターンを（Intercept1）と同じ82の変数で、（Intercept15）まで繰り返す。

私の最終目標はLOOCV経由で95個のモデルを作成し、すべてのモデルで使用されているすべての1097個の変数のp値を取得することです。どんな助けも非常に高く評価されるでしょう！

編集：例のデータ（1098個の変数のための本当のDF 96観測）あなたの実際のデータのための（例えば、データでの実際のデータのための96、10）n観測用

Response X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 

P1 N  1 1 1 0 1 0 1 0 2 2 
P2 N  2 1 1 0 2 2 1 2 2 2 
P3 N  2 1 2 1 1 0 1 1 0 1 
P4 Y  1 1 2 0 1 0 0 1 1 1 
P5 N  2 2 1 1 1 0 0 0 1 1 
P6 N  2 1 2 1 1 0 0 0 2 1 
P7 Y  2 1 1 0 2 0 0 0 2 0 
P8 Y  2 1 1 0 2 0 0 1 0 2 
P9 N  1 1 1 0 2 0 0 0 1 0 
P10 N  2 1 2 1 1 0 1 0 0 2 

Answer 1

とp変数（1098、 10の例のデータ）、以下のコードは行をnのp値の列行列で抽出する必要があります。私は、n<<pケースに合うようにしようとすると（パラメタの数に比べてごく少数の観察）統計的性質が極端に悪くなる可能性があり、ペナルティ回帰のようなテクニックを使用しない限り、不可能かもしれないと警告する必要があります。これはおそらく、あなたのパラメータの多くが見積もりに欠けている理由です（つまり、可能な1097変数のうち94個しか得られません）。特に、式のパターンが単純（0,1,2）多数のパラメータが同一直線上にあり、共同して推定することはできません（当初のモデルフィットではNAを多く見たはずです）。私たちはそれらのいくつかはある、なぜならによる共線にp値を抽出し、少し注意する必要があります。この場合、

DF2 <- read.table(row.names=1,header=TRUE,text=" 
Resp. X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 
P1 N 1 1 1 0 1 0 1 0 2 2 
P2 N 2 1 1 0 2 2 1 2 2 2 
P3 N 2 1 2 1 1 0 1 1 0 1 
P4 Y 1 1 2 0 1 0 0 1 1 1 
P5 N 2 2 1 1 1 0 0 0 1 1 
P6 N 2 1 2 1 1 0 0 0 2 1 
P7 Y 2 1 1 0 2 0 0 0 2 0 
P8 Y 2 1 1 0 2 0 0 1 0 2 
P9 N 1 1 1 0 2 0 0 0 1 0 
P10 N 2 1 2 1 1 0 1 0 0 2") 

フィットモデル

n <- nrow(DF2) 
fit <- vector(mode="list",n) ## best to pre-allocate objects 
for (i in 1:n) { 
    df <- DF2[-i,] 
    fit[[i]] <- glm (Resp. ~., data= df, family= "binomial") 
} 

：

例のデータを取得しますmissing - Rは、推定されないパラメータについては係数ベクトル（coef()）にNAを残しますが、同様にサマリの係数テーブルの行を埋め込むことはありません。もちろん

tmpf <- function(x) { 
    ## extract coef vector - has NA values for collinear terms 
    ## [-1] is to drop the intercept 
    r1 <- coef(x)[-1] 
    ## fill in values from p-value vector; leave out intercept with -1, 
    r2 <- coef(summary(x))[-1,"Pr(>|z|)"] 
    r1[names(r2)] <- r2 
    return(r1) 
} 
pvals <- sapply(fit,tmpf) 

、玩具例えば、p値のすべてが1に本質的に等しい...

## round(pvals,4) 
##  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 
## X1 0.9998 0.9998 0.9998 0.9998 0.9998 0.9998 0.9999 0.9998 0.9999 0.9998 
## X2 0.9999 0.9999 0.9999 0.9999  NA 0.9999 0.9999 0.9999 0.9999 0.9999 
## X3 0.9999 0.9999 0.9999 0.9999 0.9999 0.9998 0.9999 0.9999 0.9999 0.9999 
## X4 0.9998 0.9998 0.9998  NA 0.9998 0.9998 0.9998 0.9998 0.9998 0.9998 
## X5  NA 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000  NA 1.0000 
## X6 0.9999  NA 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999 
## X7 1.0000 1.0000 1.0000 1.0000 1.0000  NA 1.0000 1.0000 1.0000 1.0000 
## X8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 
## X9 1.0000 1.0000  NA 1.0000 1.0000 1.0000  NA  NA 1.0000  NA 
## X10  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA