観測データとシミュレーションデータのカイ2乗の不一致を計算し、ベイズ推定を使用してモデルの適合度を評価しようとしています。観測されたデータセットには、欠落している（ "NA"）値が含まれています。ただし、シミュレートされた値には欠落値はありません。したがって、それらの間の不一致の統計を比較することはできません。
以下に示すコードは例であり、私の仕事に似ています。
p <- array(runif(3000*195*6, 0, 1), c(3000, 195, 6))
N <- array(rpois(3000*195, 10), c(3000, 195))
y <- array(0, c(195, 6))
for(j in 1:195){
for(k in 1:6){
y[j,k] <- (rbinom(1, N[j], p[1,j,k]))
}
}
foo <- runif(50, 1, 195)
bar <- runif(50, 1, 6)
for(i in 1:50){
y[foo[i], bar[i]] <- NA
}
コードは、いくつかの欠落値（ "NA"）を含む応答変数yを導出します。次に、データ「y」とシミュレートされた「理想的な」データセット「y.new」のカイ二乗を計算しました。それどころか、y.newには欠落している値はありません。したがって、EとE.newの合計を比較しようとすると、y.newではなくyの欠落データを除外すると、E.newは常に大きくなるはずです。
eval <- array(NA, c(3000, 195, 6))
E <- array(NA, c(3000, 195, 6))
E.new <- array(NA, c(3000, 195, 6))
y.new <- array(NA, c(195, 6))
for(i in 1:3000){
for(j in 1:195){
for(k in 1:6){
eval[i,j,k] <- p[i,j,k]*N[i,j]
E[i,j,k] <- ((y[j,k] - eval[i,j,k])^2) / (eval[i,j,k] + 0.5)
y.new[i,j,k] <- rbinom(1, N[i,j], p[i,j,k]) # Create new "ideal" dataset
E.new[i,j,k] <- ((y.new[i,j,k] - eval[i,j,k])^2) / (eval[i,j,k] + 0.5)
}
}
} # very slow! think about how to vectorize instead of nested for loops
fit <- sum(E)
fit.new <- sum(E.new)
さて、私の質問は、不足している値をどのように処理するかです。現在、上記のコードでは、値が欠落しているため、yからevalを減算できません。たとえそれができたとしても、fitとfit.newは比較できません。私の考えは、yで欠落している値の場所を見つけて、使用している他のすべての配列から同じ[j、k]値を削除することです。これを行うための最善の方法に関する提案はありますか？
編集：私は非常に奇妙な結果を得ています。上記のようにコードを実行する場合でも、以下のように（スイープを使用して）実行する場合でも、E [1,,]はE[>,,]よりもはるかに小さくなります。特に奇妙なのは、eval [1 ,、]とeval [> 1 ,、]が同じように見えることです。問題となったのが異なるサイズの行列の処理であるかどうかを確認するために、y [j、k]を複製して各y [i ,,]が等しいy[i、j、k]にすることも試みました。なぜこれが当てはまるのか誰かが知っていますか？理論的には、このシミュレートされたデータを使用すると、E[i,,]とE.new[i,,]のすべての反復はある程度類似しているはずです。以下は、私が話していることを示すためのいくつかの要約情報です。これは新しい質問のようですが、私の元の質問に関連しています。問題を引き起こしているのはNAであるに違いないと思ったのですが、それだけではないようです。
> summary(eval[1,,])
V1 V2 V3 V4
Min. : 0.01167 Min. : 0.01476 Min. : 0.0293 Min. : 0.01953
1st Qu.: 2.60909 1st Qu.: 2.35093 1st Qu.: 2.5239 1st Qu.: 1.85789
Median : 4.85460 Median : 5.12719 Median : 5.2480 Median : 4.35639
Mean : 5.09371 Mean : 5.39451 Mean : 5.3891 Mean : 4.72061
3rd Qu.: 6.91273 3rd Qu.: 7.44676 3rd Qu.: 7.5431 3rd Qu.: 7.06119
Max. :15.81298 Max. :14.94309 Max. :14.9851 Max. :16.25751
> summary(eval1[2,,])
V1 V2 V3 V4
Min. : 0.06346 Min. : 0.06468 Min. : 0.2092 Min. : 0.006769
1st Qu.: 2.44825 1st Qu.: 1.93702 1st Qu.: 2.4226 1st Qu.: 2.426689
Median : 4.16865 Median : 4.01536 Median : 5.0771 Median : 4.833679
Mean : 4.85646 Mean : 4.64887 Mean : 5.3450 Mean : 5.169656
3rd Qu.: 6.64691 3rd Qu.: 6.96278 3rd Qu.: 7.7034 3rd Qu.: 7.229125
Max. :13.00335 Max. :13.79093 Max. :17.2673 Max. :17.915080
> summary(E[1,,])
V1 V2 V3 V4
Min. :0.00001 Min. :0.00000 Min. :0.000003 Min. :0.000008
1st Qu.:0.02744 1st Qu.:0.02723 1st Qu.:0.023008 1st Qu.:0.035854
Median :0.11750 Median :0.11889 Median :0.109138 Median :0.146706
Mean :0.39880 Mean :0.41636 Mean :0.353876 Mean :0.479533
3rd Qu.:0.46435 3rd Qu.:0.40993 3rd Qu.:0.390625 3rd Qu.:0.604021
Max. :4.43466 Max. :4.83871 Max. :6.254577 Max. :5.231650
NA's :10 NA's :8 NA's :8 NA's :10
> summary(E[2,,])
V1 V2 V3
Min. : 0.0000 Min. : 0.00003 Min. : 0.00002
1st Qu.: 0.8213 1st Qu.: 0.42091 1st Qu.: 0.36853
Median : 2.0454 Median : 2.31697 Median : 2.39892
Mean : 8.0619 Mean : 9.40838 Mean : 6.38919
3rd Qu.: 5.6755 3rd Qu.: 6.34782 3rd Qu.: 4.89749
Max. :395.9499 Max. :172.83324 Max. :120.93648
NA's :10 NA's :8 NA's :8
ありがとう、ダン