我正在尝试计算观察数据和模拟数据的卡方差异，并使用贝叶斯推理评估模型拟合度。观察到的数据集包含缺失(“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
}
代码派生响应变量 y，其中包括一些缺失值(“NA”)。然后，我计算了数据“y”和模拟“理想”数据集“y.new”的卡方。相反，y.new 没有任何缺失值。因此，当我尝试比较 E 和 E.new 的总和时，如果我忽略 y 中的缺失数据而不是 y.new，则 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[>1,,] 小得多。特别奇怪的是 eval[1,,] 和 eval[>1,,] 看起来是相同的。我什至尝试复制 y[j,k] 使其成为 y[i,j,k]，其中每个 y[i,,] 相等，只是为了看看是否是处理不同大小的矩阵造成了问题。有谁知道为什么会出现这种情况？理论上，通过这个模拟数据，我认为 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
谢谢， 丹
最佳答案
您可以在内循环中添加测试并更改循环的顺序，如下所示:
...
for(j in 1:195){
for(k in 1:6){
if ( !is.na(y(j,k)) ) {
for(i in 1:3000){
...
}
}
}
}
...
为了提高效率，请对内部循环进行矢量化(如上面的评论中所述)。
也可以定义与y具有相同维度的逻辑数组。表示定义位置的子集，例如 subset <- !is.na(y)并使用它来代替。
关于r - 如何处理 R 中 for 循环中的缺失数据 (NA)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/13437257/