卡方值与Fisher 在 2×2 表上的卡方检验的精确概率

Question

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

• 数据以矩阵形式输出。
• 选择第 2 到第 4 列中的值，以便在拟合到 2×2 表时边际总计保持不变。 （边际总数由 n_A、n_B、n_T 和 n_F 预先确定。）
• 根据第 2-5 列的值计算得出的卡方值列在第 6 列。
• Fisher's exact probability，根据第 2-5 列的值计算得出，列在第 6 列。

从上面的角度来看，有一些关注点，但我想先关注以下我的问题。

我的问题

除了以上功能，我想补充以下规范；

• 将所有列按第 6 列的值排序，使第 6 列的值随着下降而增加。
• 第 8 列中的值应该是第 7 列中的值从顶部加起来的值。
• 然后，编写一个图表，其中 x 轴为第 6 列，y 轴为第 7 列。

换句话说，我们想要从Box1的代码结果中得到如下表1和图1的表格和图形。

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

表 1

图1

盒子1

#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,TB,FA,FB,chiq_2by2(TA,TB,FA,FB),exp(lnop),　) #★1
      resul <- rbind(resul, pr)
    }
  }
}

resul

Answer 1

这不是一个非常聪明的代码，但我自己想出来了。 “#answer of this question”下的部分是必不可少的部分；这部分旨在对矩阵的所有行进行排序，以便随着您的下降而增加第 6 列的值。并为第 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,TB,FA,FB,chiq_2by2(TA,TB,FA,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])

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

表

图表