Fisher-b in R

From InterSciWiki
Jump to: navigation, search


Pat: the 3-way probability test would be to look at the row containing the tab1 cell 1,1 value. In this case (for the original fisherb.R below) that's row 6 in either the dougtest$nicemat[6,]=.00032 (rounded) or the dougtest$possmat[6,] (scientific notation) attributes: the value reported in the original article, (p=.00032761 with rounding error of R versus Fortran).

In Fisher-b.pdf or [Fisher-b.rtf] Its the vocabulary I probably dont get:

#No library needed
tab1<-matrix(c(9,8,9,5),nrow=2,ncol=2)
rownames(tab1)<-c("bwabsent","bwpresent")
colnames(tab1)<-c("patriabsent","patripresent")
tab2<-matrix(c(76,11,22,46),nrow=2,ncol=2)
rownames(tab2)<-rownames(tab1)
colnames(tab2)<-colnames(tab1)
#Assume that you have done the following commands:
tab1  #works
tab2  #works created tab1 and tab2
source('fisherb.r')
source("http://intersci.ss.uci.edu/wiki/R/fisherb.R")
dougtest<-fisherb(tab1,tab2)
ls()  #This will give you a listing of the objects in memory.  You should have at least tab1, tab2, fisherb, and dougtest. Now issue the command:
dougtest
dougtest$possmat[6,]  #p=.0003168398
dougtest$possmat
     a11 a12 a21 a22         pval
[1,]   4  14  13   0 3.132032e-15
[2,]   5  13  12   1 3.737766e-12
[3,]   6  12  11   2 1.278710e-09
[4,]   7  11  10   3 1.728531e-07
[5,]   8  10   9   4 1.059343e-05
[6,]   9   9   8   5 3.168398e-04
[7,]  10   8   7   6 4.815964e-03
[8,]  11   7   6   7 3.788783e-02
[9,]  12   6   5   8 1.542035e-01
[10,]  13   5   4   9 3.182411e-01
[11,]  14   4   3  10 3.182411e-01
[12,]  15   3   2  11 1.415431e-01
[13,]  16   2   1  12 2.372610e-02
[14,]  17   1   0  13 1.013685e-03
dougtest$nicemat[6,] #p=.00032
dougtest$nicemat
     a11 a12 a21 a22    pval
[1,]   4  14  13   0 0.00000
[2,]   5  13  12   1 0.00000
[3,]   6  12  11   2 0.00000
[4,]   7  11  10   3 0.00000
[5,]   8  10   9   4 0.00001
[6,]   9   9   8   5 0.00032
[7,]  10   8   7   6 0.00482
[8,]  11   7   6   7 0.03789
[9,]  12   6   5   8 0.15420
[10,]  13   5   4   9 0.31824
[11,]  14   4   3  10 0.31824
[12,]  15   3   2  11 0.14154
[13,]  16   2   1  12 0.02373
[14,]  17   1   0  13 0.00101
dougtest$outmat
         P value   Gamma
Table One  0.7168 -0.2308
Table Two  0.0000  0.8705


Doug 13:18, 14 June 2012 (PDT) I still dont get it, sorry to be so dense

do I need to install a library?


source("http://intersci.ss.uci.edu/wiki/R/fisherb.R")
foo<-xtabs(tab1,tab2) 
 source("http://intersci.ss.uci.edu/wiki/R/fisherb.R")
 fisherb<-function(tab1,tab2) {
 	###################################################################
 	# Function to do exact significant test for three-way interaction
 	# interaction effects.
 	#	Version 1.0 J. Patrick Gray June 11, 2012
 	# Article:	Douglas R. White, Robert Pesner, and Karl P. Reitz
 	#						"An Exact Significance Test for Three-Way Interaction
 	#						Effects" Behavior Science Research 18(2):103-122. 1983.
 	# Input:	tab1 = 2x2 matrix
 	#					tab2 = 2x2 matrix
 	# Suggest making the table with the fewest observations tab1
 	# Output:	possmat = a matrix with five columns running from:
 	#         underline a111 (the smallest possible cell 111 for tab1
 	#														given marginal constraints)
 	#         overline a111 (the largest possible cell 111 for tab1
 	#														given marginal constraints)
 	#					The first four columns are the cells of table 1
 	#					The last column is the point probability of the table
 	#################################################################  
 	
 	###############################################
 	# Calculate marginals needed in the denominator
 	A.11<-tab1[1,1]+tab2[1,1]
 	A.12<-tab1[1,2]+tab2[1,2]
 	A.21<-tab1[2,1]+tab2[2,1]
 	A.22<-tab1[2,2]+tab2[2,2]
 	#################################################
 	
 	##########################################
 	# Find mina111 & maxa111
 	maxa111<-upa(tab1,tab2)
 	maxa<-maxa111[1,1]
 	mina111<-doa(tab1,tab2)
 	mina<-mina111[1,1]
 	##########################################
 	
 	#############################################################
 	# Now set up matrix of possibilities from maxa111 to mina111
 	# use seq based on range of mina to maxa
 	gap<-(maxa-mina)
 	a11<-mina111[1,1]:(mina111[1,1]+gap)
 	a12<-mina111[1,2]:(mina111[1,2]-gap)
 	a21<-mina111[2,1]:(mina111[2,1]-gap)
 	a22<-mina111[2,2]:(mina111[2,2]+gap)
 	pval<-rep(0,gap+1)
 	possmat<-as.matrix(cbind(a11,a12,a21,a22,pval))
 	###############################################################
 
 	##############################################################
 	# Calculate denominator in equations 3 & 4
 	denom<-0
 	# target is the last row in possmat
 	target<-nrow(possmat)
 	for (ii in 1:target) {
 		n1<-lchoose(A.11,possmat[ii,1])
 		n2<-lchoose(A.12,possmat[ii,2])
 		n3<-lchoose(A.21,possmat[ii,3])
 		n4<-lchoose(A.22,possmat[ii,4])
 		# Remember these are logs so addition, not muliplication
 		tmp<-n1+n2+n3+n4
 		# Use exp() to convert from logs
 		denom<-denom+exp(tmp)
 	}
 	################################################################
 	
 	
 	#################################################################
 	# Do probability for each table in possmat: The numerator in 
 	#    equation 3
 	for (ii in 1:length(pval)) {
 		tmp<-0
 		n1<-lchoose(A.11,possmat[ii,1])
 		n2<-lchoose(A.12,possmat[ii,2])
 		n3<-lchoose(A.21,possmat[ii,3])
 		n4<-lchoose(A.22,possmat[ii,4])
 		tmp<-n1+n2+n3+n4
 		possmat[ii,5]<-exp(tmp)/denom
 	}
 	#################################################################
 	
 	##############################################
 	# Make the probabilities look nice
 	nicemat<-possmat
 	nicemat[,5]<-round(nicemat[,5],digit=5)
 	##############################################
 	
 	##########################################################################################
 	# Create display for individual table independence and gamma
 	outmat<-matrix(rep(0,4),nrow=2)
 	rownames(outmat)<-c("Table One","Table Two")
 	colnames(outmat)<-c("P value","Gamma")
 	tmp<-fisher.test(tab1)
 	outmat[1,1]<-round(tmp$p.value,4)
 	gamma<-(tab1[1,1]*tab1[2,2]-tab1[1,2]*tab1[2,1])/(tab1[1,1]*tab1[2,2]+tab1[1,2]*tab1[2,1])
 	outmat[1,2]<-round(gamma,4)
 	tmp<-fisher.test(tab2)
 	outmat[2,1]<-round(tmp$p.value,4)
 	gamma<-(tab2[1,1]*tab2[2,2]-tab2[1,2]*tab2[2,1])/(tab2[1,1]*tab2[2,2]+tab2[1,2]*tab2[2,1])
 	outmat[2,2]<-round(gamma,4)
 	###########################################################################################
 	
 	return=list(possmat=possmat,nicemat=nicemat,outmat=outmat)
 }
 
 upa<-function(x1,x2) {
 	# Function to find the largest 111 cell given constraints 
 	ori<-x1[1,1]
 	out<-x1
 	max11<-x1[1,1]+x2[1,1]
 	for (ii in ori:max11) {
 		x1[1,1]<-x1[1,1]+1
 		x2[1,1]<-x2[1,1]-1
 		x1[1,2]<-x1[1,2]-1
 		x2[1,2]<-x2[1,2]+1
 		x1[2,1]<-x1[2,1]-1
 		x2[2,1]<-x2[2,1]+1
 		x1[2,2]<-x1[2,2]+1
 		x2[2,2]<-x2[2,2]-1
 		tmp1<-all(x1>=0)
 		tmp2<-all(x2>=0)
 		#break condition if any cell below zero
 		if (tmp1==FALSE | tmp2==FALSE) {
 			break
 		}else {
 			out<-x1
 		}
 	}
 	return(out)
 }
 
 doa<-function(x1,x2) {
 		# Function to find the smallest 111 cell given constraints 
 	ori<-x1[1,1]
 	out<-x1
 	for (ii in ori:0) {
 		x1[1,1]<-x1[1,1]-1
 		x2[1,1]<-x2[1,1]+1
 		x1[1,2]<-x1[1,2]+1
 		x2[1,2]<-x2[1,2]-1
 		x1[2,1]<-x1[2,1]+1
 		x2[2,1]<-x2[2,1]-1
 		x1[2,2]<-x1[2,2]-1
 		x2[2,2]<-x2[2,2]+1
 		tmp1<-all(x1>=0)
 		tmp2<-all(x2>=0)
 		#break condition if any cell below zero
 		if (tmp1==FALSE | tmp2==FALSE) {
 			break
 		}else {
 			out<-x1
 		}
 	}
 	return(out)
 }