求助：如何一次性计算出二项分布检验的p值

数据框testdata中有若干个表格的数据，其中的8个表格数据如下：

a，b，c，d，a1，b1，c1，d1

1 22，1，5，3，13.2266667，1.03333333，8.99，7.75

2 8，0，0，1，3.6906475，0.29136691，2.6546763，2.3633094

3 0，1，0，0，0.3925926，0.03703704，0.3037037，0.2666667

4 0，1，0，0，0.3925926，0.03703704，0.3037037，0.2666667

5 0，1，0，0，0.3925926，0.03703704，0.3037037，0.2666667

6 16，1，31，19，24.327381，1.99404762，22.5327381，18.1458333

7 2，0，7，2，4.2428571，0.35357143，3.4964286，2.9071429

8 0，0，0，5，1.9343066，0.16423358，1.4963504，1.4051095

说明：上面8个表格中各有8个数值，每个表格中的前四个值为观测频次（a,b,c,d），后面四个值为前面四个观测频次对应的期望频次(a1,b1,c1,d1)。例如，第一个表格中的四个观测频次为22，1，5，3，它们各自对应的期望频次为13.2266667，1.03333333，8.99，7.75。

---

我现在需要对每个表格中的每个观测频次进行二项分布检验（binomial test），计算出p值。

计算使用的公式为: binom.test(观测频次, 观测频次总数, 期望频次/（即除以）观测频次总数, alternative="??")$p.value，其中alternative选项部分由观测频次与其对应的期望频次大小关系来决定，当观测频次大于期望频次时，选用alternative="greater"；当观测频次小于期望频次时，选用alternative="less".

例如第一个表格中四个观测频次的总数为31（即22+1+5+3），对第一个值22与其对应的期望频次13.2266667进行二项分布检验，计算过程为binom.test(22,31,13.2266667/31,alternative="greater")$p.value

最后的p值为0.001356598；对第二个值1与其对应的期望频次1.03333333进行二项分布检验，计算过程为binom.test(1,31,1.03333333/31,alternative="less")$p.value，最后的p值为 0.723323。

---

问题：1. 要对每个表格中的的每个观测频次与其期望频次进行二项分布检验，如何一次行计算出其所有p值？

      2. 再将计算出的所有p值进行对数转换（以10为底）？

1. s=testdata
   
2. bt=g=matrix(rep(0,nrow(s)*4),ncol=4)
   
3. for(i in 1:nrow(s)){
   
4. g[i,]=unlist(lapply(1:4,function(x) ifelse(s[i,x]<s[i,x+4],"less","greater")))
   
5. bt[i,]=unlist(lapply(1:4,function(x)binom.test(s[i,x],sum(s[i,1:4]),s[i,x+4]/sum(s[i,1:4]),alternative=g[i,x])$p.value))
   
6. }
   
7. g
   
8. bt
   
9. lbt=log10(bt)

复制代码

自己顶一下！

write a function

你已经描述的这么清楚了， 也应该知道这么写了吧

用EXCEL可以很快做到（我不会用R）：
在空白EXCEL sheet的第A2～H9（这个8格×8格）的范围内输入你的数据(8×8)，例如A2输22，B2输入1，..., A3输入8,...，然后在J列输入求和公式，J2格输入=SUM(A2:D2),COPY该公式到J3～J9.然后在K2输入公式=IF(A2>E2,1-BINOMDIST(A2,$J2,E2/$J2,TRUE),BINOMDIST(A2,$J2,E2/$J2,TRUE))，COPY该公式到K2~N9的范围内。那么K2到N9范围内显示的就是所有你要的p值。
如果你要取LOG，可以在K11输入公式=IFERROR(LOG10(K2),"N/A")，然后COPY公式到K11～N18的范围内，注：IFERROR是为了不出错，因为有p值为0，取LOG无意义。就全部做到你要的结果了。
不过，我计算出来的第一个p值是0.0003746，不是你举的例子0.001356598(虽然第2个跟你的一样0.723323)，不知为何，你可以检查一下，如我有错，请指正。

求证1加1 发表于 2014-2-24 22:13
write a function

就烦劳版主大人帮忙写个函数吧。

jmpamao 发表于 2014-2-24 22:26
你已经描述的这么清楚了， 也应该知道这么写了吧

可恨自己根本不会写代码啊

TimeT 发表于 2014-2-24 23:03
用EXCEL可以很快做到（我不会用R）：
在空白EXCEL sheet的第A2～H9（这个8格×8格）的范围内输入你的数据( ...

多谢多谢啊！另外，第一个p值就是我算出的那个值。
> binom.test(22,31,13.2266667/31,alternative="greater")$p.value
[1] 0.001356598
我还不知道问题出在什么地方。

灰常感谢9楼的yywan0913，这个代码相当管用，一次性解决了我上面的两个问题！