Talk:Fisher-b in R

From InterSciWiki
Jump to: navigation, search


This is a silly example, but if you read the sccs files in and were planning to run fisherb without having to input the tables by hand, this is how you would do it. The bipaint.tab file should go in your working directory.

Pat - thanks for helping me along. Now we can edit down to the nub so that the merest simpleton like myself can get the results. (Pat): the 3-way probability test would be to look at the row containing the tab1 cell 11 value (????). That is 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:

"When using the call sign the results to a new object ...."

i.e., foo <- what?
foo<-xtabs(tab1,tab2) #?
(sorry: just never had a course in R)

Doug: Let me walk you through the doc file with some vocabulary.

  1. With what you did below you have two objects in memory (tab1 and tab2). I assume you realize they could be named anything (e.g., foo1 and foo2). YES-GOT THAT
  2. Distinguish between the function (named fisherb) and the script file that contains it (named fisherb.r). The fisherb.r script file should be in your working directory (usually Rdata). If it is in another directory change to that one. YES
  3. Bring the fisherb function into memory by the command: source('fisherb.r'). If you get an error is it almost always because the file isn't in your working directory. YES-I use source("http://intersci.ss.uci.edu/wiki/R/fisherb.R")
  4. Now the fisherb function is in memory (and continues to be so until you close the session). To call it you need to assign the results an object name. In the doc file I called the object test (you can name it anything). That object contains three attributes: test$outmat, test$possmat, test$nicemat. If you want to see them all at once just type test. Or you can view them one at a time using the objectname$attribute format as I did in the doc (test$outmat, etc.) WELL- test at command line doesnt show anything- am doing it wrong
  5. Each attribute of the object can now be used for further calculations. For example, to find the probability of cell A111 being 6 or less you could just do: sum(test$possmat[1:3,5]) or for the probability of cell A111 being 12 or higher: sum(test$possmat[9:14,5]).

Does this help? If not try to run an example and send me screen shots of any errors.

Pat Gray

Doug: I think I may have confused you by thinking you knew more R than you did. I should not have mentioned xtabs in the context of you having tables in memory. Let me break the issue into two task--one for this e-mail and in another.

Assume that you have done to following: 1. created tab1 and tab2 2. source('fisberb.r')

issue the command:

dougtest<-fisherb(tab1,tab2)

When you do this you should just get the R prompt >

issue the command:

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?
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)
tab1  #works
tab2  #works
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)
 }