Chi squared value v.s. Fisher’s exact probability on the Chi squared test on 2×2 table

It’s not a very smart code, but I figured it out on my own. The part under the “#answer of this question” is the essential part; this part is intended to sort all rows of the matrix, so as to the 6th column’s value increase as you go down. And write a graph for column 6 vs. column 7 Other parts are as same as the code described in the Box1 of the question.

Box

#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])

As the result, we can get the following Table and graph.

Table enter image description here

Graph enter image description here

CLICK HERE to find out more related problems solutions.

Leave a Comment

Your email address will not be published.

Scroll to Top