# published on: # http://www.r-statistics.com/2010/02/barnards-exact-test-a-powerful-alternative-for-fishers-exact-test-implemented-in-r/ Barnardextest<-function(Ta,Tb =NULL,Tc =NULL,Td =NULL, to.print = F, to.plot = T) { # The first argument (Ta) can be either a table or a matrix of 2X2. # Or instead, the values of the table can be entered one by one to the function #Barnard's test calculates the probabilities for contingency tables. It has been shown that for 2x2 tables, Barnard's test #has a higher power than Fisher's Exact test. Barnard's test is a non-parametric test that relies upon a computer to generate #the distribution of the Wald statistic. Using a computer program, one could find the nuisance parameter that maximizes the #probability of the observations displayed from a table. #Despite giving lower p-values for 2x2 tables, Barnard's test hasn't been used as often as Fisher's test because of its #computational difficulty. This code gives the Wald statistic, the nuisance parameter, and the p-value for any 2x2 table. #The table can be written as: # Var.1 # --------------- # a b r1=a+b # Var.2 # c d r2=c+d # --------------- # c1=a+c c2=b+d n=c1+c2 #One example would be # Physics # Pass Fail # --------------- # Crane 8 14 # Collage # Egret 1 3 # --------------- # #After implementing the function, simply call it by the command: #Barnardextest(8,14,1,3) #This will display the results: #"The contingency table is:" # [,1] [,2] #[1,] 8 14 #[2,] 1 3 #"Wald Statistic:" #0.43944 #"Nuisance parameter:" #0.9001 #"The 1-tailed p-value:" #0.4159073 #On a side note, I believe there are more efficient codes than this one. For example, I've seen codes in Matlab that run #faster and display nicer-looking graphs. However, this code will still provide accurate results and a plot that gives the #p-value based on the nuisance parameter. I did not come up with the idea of this code, I simply translated Matlab code #into R, occasionally using different methods to get the same result. The code was translated from: # #Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez, L. Rodriguez-Cardozo. Probability Test. A MATLAB file. URL #http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6198 # #My goal was to make this test accessible to everyone. Although there are many ways to run this test through Matlab, I hadn't #seen any code to implement this test in R. I hope it is useful for you, and if you have any questions or ways to improve #this code, please contact me at calhoun.peter@gmail.com. # Tal edit: choosing if to work with a 2X2 table or with 4 numbers: if(is.null(Tb) | is.null(Tc) | is.null(Td) ) { # If one of them is null, then Ta should have an entire table, and we can take it's values if(is.table(Ta) | is.matrix(Ta)) { if(sum(dim(Ta) == c(2,2))==2) { Tb <- Ta[1,2] Tc <- Ta[2,1] Td <- Ta[2,2] Ta <- Ta[1,1] } else {stop("The table is not 2X2, please check it again...") } } else {stop("We are missing value in the table") } } c1<-Ta+Tc c2<-Tb+Td n<-c1+c2 pao<-Ta/c1 pbo<-Tb/c2 pxo<-(Ta+Tb)/n TXO<-abs(pao-pbo)/sqrt(pxo*(1-pxo)*(1/c1+1/c2)) n1<-prod(1:c1) n2<-prod(1:c2) P<-{} for( p in (0:99+.01)/100) { TX<-{} S<-{} for( i in c(0:c1)) { for( j in c(0:c2)) { if(prod(1:i)==0){fac1<-1} else {fac1<-prod(1:i)} if(prod(1:j)==0){fac2<-1} else {fac2<-prod(1:j)} if(prod(1:(c1-i))==0){fac3<-1} else {fac3<-prod(1:(c1-i))} if(prod(1:(c2-j))==0){fac4<-1} else {fac4<-prod(1:(c2-j))} small.s<-(n1*n2*(p^(i+j))*(1-p)^(n-(i+j))/(fac1*fac2*fac3*fac4)) S<-cbind(S,small.s) pa<- i/c1 pb<-j/c2 px <- (i+j)/n if(is.nan((pa-pb)/sqrt(px*(1-px)*((1/c1)+(1/c2))))) { tx<-0 } else { tx <- (pa-pb)/sqrt(px*(1-px)*((1/c1)+(1/c2))) } TX<-cbind(TX,tx) } } P<-cbind(P,sum(S[which(TX>=TXO)])) } np<-which(P==max(P)) p <- (0:99+.01)/100 nuisance<-p[np] pv<-P[np] if(to.print) { print("The contingency table is:") print(matrix(c(Ta,Tc,Tb,Td),2,2)) print("Wald Statistic:") print(TXO) print("Nuisance parameter:") print(nuisance) print("The 1-tailed p-value:") print(pv) } if(to.plot) { plot(p,P,type="l",main="Barnard's exact P-value", xlab="Nuisance parameter", ylab="P-value") points(nuisance,pv,col=2) } return( list( contingency.table = as.table(matrix(c(Ta,Tc,Tb,Td),2,2)), Wald.Statistic = TXO, Nuisance.parameter = nuisance, p.value.one.tailed = pv ) ) } Barnardextest(matrix(c(8,1,14,3),2,2)) fisher.test(matrix(c(8,1,14,3),2,2)) Convictions <- matrix(c(2, 10, 15, 3), nrow = 2, dimnames = list(c("Dizygotic", "Monozygotic"), c("Convicted", "Not convicted"))) Convictions fisher.test(Convictions, alternative = "less") Barnardextest(Convictions) library(highlight) highlight( "C:\\R\\aaaa good code\\Barnardtest.R", renderer = renderer_html(), output = "C:\\R\\aaaa good code\\Barnardtest.R.html" )