卡方值v.s. Fisher在2×2表上的卡方检验中的确切概率

如何解决卡方值v.s. Fisher在2×2表上的卡方检验中的确切概率

用R编写的Box1的源代码应符合以下规范；

数据以矩阵形式输出。

选择第2列至第4列中的值，以便当拟合到2×2表格时的边际总数恒定。（边际总数由n_A，n_B，n_T和n_F预先确定。）

从第2-5列的值计算出的卡方值列在第6列。
根据第2-5列的值计算的
Fisher's exact probability列在第6列。

从上述观点来看，有些问题值得关注，但我想首先关注以下“我的问题”。

我的问题

除了上述功能外，我还要添加以下规范；

按第6列的值对所有列进行排序，以使第6列的值随着您的下降而增加。

第8列中的值应该是第7列中从顶部开始加在一起的值。

然后，在x轴上的第6列和y轴上的第7列编写一个图形。

换句话说，我们希望从Box1的代码结果中得到类似于下表1和图1的表格。

现在，每次更改设置时，我都会在Excel中执行此过程，但是可以在R中完成吗？

表1

图1

Box1

#Function to caluculate ln(n!)
ln_fact<-function(n){
  if (n==0){ans=0}else{
    ans=0
    for(i in 1:n) {ans=ans+log(i)}}
  return(ans)
}

#Fuction to caluculate chiq2 value
chiq_2by2<-function(TA,TB,FA,FB){
  nA=TA+FA;nB=TB+FB; ntot=nA+nB
  nF=FA+FB;nT=TA+TB
  ETA=(nT*nA)/ntot;EFA=(nF*nA)/ntot
  ETB=(nT*nB)/ntot;  EFB=(nF*nB)/ntot
  
  ch=((TA-ETA)^2)/(ETA);ch=ch+((TB-ETB)^2)/(ETB)
  ch=ch+((FA-EFA)^2)/(EFA);ch=ch+((FB-EFB)^2)/(EFB)
  return(ch)
}

#main part
##Set marginal total of 2×2.
n_A=14
n_B=6
n_T=13
n_F=n_A+n_B-n_T

##part1 of probability of occurrence
lnop1=ln_fact(n_A)+ ln_fact(n_B)+ln_fact(n_T)+ln_fact(n_F) - ln_fact(n_A+n_B)  


cnt=0;
A_tot=n_A; B_tot=n_B
resul=0
for(i in 0:A_tot){
  for(j in 0:B_tot){
##Calculating the elements of a 2×2 table.
    TA=i;  FA=A_tot-TA
    TB=j;    FB=B_tot-TB

## judging whether or not the elements of a 2×2 are well-defined.
    br1<-(TA+TB==n_T);br2<-(FA+FB==n_F)
    br3<-(TA+FA==n_A);br4<-(TB+FB==n_B)
    br=br1*br2*br3*br4
## To calculate the chi-square value and Fisher's direct probability for the well-defined conditions   
    if (br==1){
      cnt=cnt+1
###ln(probability of occurrence),probability of occurrence is based on the Fisher's direct probability
      lnop=lnop1-(ln_fact(TA)+ ln_fact(TB)+ln_fact(FA)+ln_fact(FB))  
      
      pr=c(cnt,TA,FB,chiq_2by2(TA,FB),exp(lnop),　) #★1
      resul <- rbind(resul,pr)
    }
  }
}

resul

解决方法

这不是一个很聪明的代码，但是我自己弄清楚了。 “此问题的＃答案”下的部分是必不可少的部分；此部分旨在对矩阵的所有行进行排序，以便随着向下移动第六列的值增加。并为第6列和第7列写一个图形其他部分与问题Box1中描述的代码相同。

框

#Function to caluculate ln(n!)
ln_fact<-function(n){
  if (n==0){ans=0}else{
    ans=0
    for(i in 1:n) {ans=ans+log(i)}}
  return(ans)
}


#Fuction to caluculate chiq2 value
chiq_2by2<-function(TA,TB,FA,FB){
  nA=TA+FA;nB=TB+FB; ntot=nA+nB
  nF=FA+FB;nT=TA+TB
  ETA=(nT*nA)/ntot;EFA=(nF*nA)/ntot
  ETB=(nT*nB)/ntot;  EFB=(nF*nB)/ntot
  
  ch=((TA-ETA)^2)/(ETA);ch=ch+((TB-ETB)^2)/(ETB)
  ch=ch+((FA-EFA)^2)/(EFA);ch=ch+((FB-EFB)^2)/(EFB)
  return(ch)
}

#main part
##Set marginal total of 2×2.
n_A=14
n_B=6
n_T=13
n_F=n_A+n_B-n_T

##part1 of probability of occurrence
lnop1=ln_fact(n_A)+ ln_fact(n_B)+ln_fact(n_T)+ln_fact(n_F) - ln_fact(n_A+n_B)  


cnt=0;
A_tot=n_A; B_tot=n_B
resul=0
for(i in 0:A_tot){
  for(j in 0:B_tot){
    ##Calculating the elements of a 2×2 table.
    TA=i;  FA=A_tot-TA
    TB=j;    FB=B_tot-TB
    
    ## judging whether or not the elements of a 2×2 are well-defined.
    br1<-(TA+TB==n_T);br2<-(FA+FB==n_F)
    br3<-(TA+FA==n_A);br4<-(TB+FB==n_B)
    br=br1*br2*br3*br4
    ## To calculate the chi-square value and Fisher's direct probability for the well-defined conditions   
    if (br==1){
      cnt=cnt+1
      ###ln(probability of occurrence),probability of occurrence is based on the Fisher's direct probability
      lnop=lnop1-(ln_fact(TA)+ ln_fact(TB)+ln_fact(FA)+ln_fact(FB))  
      
      pr=c(cnt,TA,FB,chiq_2by2(TA,FB),exp(lnop),0) #★1
      resul <- rbind(resul,pr)
    }
  }
}

#answer of this question
rownames(resul) <- NULL
dat=resul

#↓If you do not want the point (0,0) to appear on the graph,comment-out it.
#dat[1,1]=NA;dat=na.omit(dat)

dat=dat[order(dat[,6]),]

dat[1,8]=dat[1,7]
for(i in 2:length(dat[,8])){dat[i,8]=dat[i-1,8]+dat[i,7]}
dat

plot(dat[,7] ~ dat[,6],col = "red")
curve(dchisq(x,1),col="red",add=T)
#par(new=T) 
#plot(dat[,8] ~ dat[,6])

结果，我们可以获得下面的表和图。

表

图