用于检查基因型 pvalue 的循环卡方检验

Question

以下是我试过的代码

aaxab.res=NULL

chiaaxab <- for (i in 1:length(gt)) {
               observed <- c(aa.count[i], ab.count[i])
               expected <- c(0.5, 0.5)
               chisq.test(observed, p=c(0.5, 0.5))$p.value
               p.value <- ifelse(is.na(T.value ))
               aaxab.res <- c(abxaa.res, p.value)
            }

但我不断收到错误：

Error in chisq.test(observed, p = c(0.5, 0.5)) :
  at least one entry of 'x' must be positive.

有人可以帮忙纠正这个问题，因为我对 R 很陌生

Answer 1

我认为您的意思是比较预期与观察到的情况，在这种情况下您的代码应该是：

               ...
               chisq.test(observed, expected, p=c(0.5, 0.5))$p.value
               ...

当您的观察值都为零时，您将收到该错误。例如：

observed <- c(0, 0)
chisq.test(observed, p=c(0.5, 0.5))

chisq.test 中的错误（观察到，p = c(0.5, 0.5)) ： 至少有一项“x”必须是正数

如果您的数据集为零（有时是正常的），并且您希望循环继续而不抛出错误，您可以添加条件语句来检查零：

chiaaxab <- for (i in 1:length(gt)) {
  observed <- c(aa.count[i], ab.count[i])
  expected <- c(0.5, 0.5)
  if (aa.count[i] == 0 & ab.count[i] == 0) {
    aaxab.res <- c(abxaa.res, NA)
  } else {
    chisq.test(observed, p=c(0.5, 0.5))$p.value
    p.value <- ifelse(is.na(T.value))
    aaxab.res <- c(abxaa.res, p.value) 
    }
}

如果检测到两个零，上面的代码将添加NA。

也不确定你想用这条线完成什么：p.value <- ifelse(is.na(T.value))。语法为：ifelse(test, yes, no)