####################################
# Chi-square test
####################################

Chi.t=function(a,b){
	n1=sum(a)
	n2=sum(b)
	N=length(a)

  ###Expected number###
      E1=n1*n2/N
      E2=(N-n1)*(N-n2)/N
      E3=n1*(N-n2)/N
      E4=(N-n1)*n2/N
      EE=c(E1,E2,E3,E4)
 ####Observed number########
	O1=sum(a*b)
	O2=sum((1-a)*(1-b))
	O3=sum(a*(1-b))
	O4=sum((1-a)*b) 
        OO=c(O1,O2,O3,O4)
 ##### Chi-sq stat.#########
    temp.j=sum(EE<10)   
    
    chi.s1 = sum( (OO-EE)^2/EE ) #without Yate's correction

    temp=abs(OO-EE)-0.5
    temp1=temp[temp>0]	
    chi.s2=sum(temp1^2/EE[temp>0]) #with Yate's correction

    chi.s = chi.s1*(temp.j==0)+chi.s2*(temp.j>0)
    p=pchisq(chi.s,1,lower.tail=F)
    return(p)
}


##################################
# Bin. and Poi. test
##################################    
BP.t=function(a,b,posi.test=TRUE){
	n1=sum(a)
	n2=sum(b)
	N=length(a)

  ###Expected number###
      E1=n1*n2/N
      E2=(N-n1)*(N-n2)/N
      E3=n1*(N-n2)/N
      E4=(N-n1)*n2/N
 
 ####Observed number########
	O1=sum(a*b)
	O2=sum((1-a)*(1-b))
	O3=sum(a*(1-b))
	O4=sum((1-a)*b) 

 ####Binomial################
  p1.p=pbinom(O1-1, N, E1/N, lower.tail = FALSE) #P1 for testing positive association
  p1.n=pbinom(O1, N, E1/N, lower.tail = TRUE)    #P1 for tesing negative association
  p1=p1.p*posi.test+p1.n*(1-posi.test)           #P1

  p22.p=pbinom(O2-1, N, E2/N, lower.tail = FALSE) 
  p22.n=pbinom(O2, N, E2/N, lower.tail = TRUE)    
  p2.p=p1.p*p22.p                           #P2 for testing positive association
  p2.n=p1.n*p22.n                           #P2 for tesing negative association
  p2=p2.p*posi.test+p2.n*(1-posi.test)      #P2 

  p31.p=pbinom(O3 , N, E3/N, lower.tail = TRUE)
  p32.p=pbinom(O4 , N, E4/N, lower.tail = TRUE)
  p3.p=p31.p*p32.p             #P3 for testing positive association
  p31.n=pbinom(O3-1, N, E3/N, lower.tail = FALSE) 
  p32.n=pbinom(O4-1, N, E4/N, lower.tail = FALSE) 
  p3.n=p31.n*p32.n             #P3 for tesing negative association
  p3=p3.p*posi.test+p3.n*(1-posi.test)      #P3 

 ####Poisson################
  p4.p=ppois(O1-1, E1, lower.tail = FALSE) #P4 for testing positive association
  p4.n=ppois(O1, E1, lower.tail = TRUE)    #P4 for tesing negative association
  p4=p4.p*posi.test+p4.n*(1-posi.test)           #P4

  p5.p=ppois(O1+O2-1,E1+E2, lower.tail = FALSE) #P5 for testing positive association
  p5.n=ppois(O1+O2, E1+E2, lower.tail = TRUE)    #P5 for tesing negative association
  p5=p5.p*posi.test+p5.n*(1-posi.test)           #P5

  p6.p=ppois(O3+O4, E3+E4, lower.tail = TRUE) #P6 for testing positive association
  p6.n=ppois(O3+O4-1, E3+E4 , lower.tail =FALSE )    #P6 for tesing negative association
  p6=p6.p*posi.test+p6.n*(1-posi.test)           #P6
  
  p.all=c(p1,p2,p3,p4,p5,p6)
  return(round(p.all,digits=5))
}



############# using data(lansing) for demonstration############
library(letsR)
library(spatstat)
library(cooccur)

data(lansing)
data.t=split(lansing)
data.names=names(data.t)
i.end=length(data.names)
w<-1   #width of the study area
l<-1  #length of the study area
g<-0.05     #width of each cell
r<-g^2   #area of each cell

XX=NULL
N.all=NULL
for (i in 1:i.end){
	test.d=lansing[lansing$marks==data.names[i],]
	species=test.d$marks
	xy=cbind(test.d$x,test.d$y)
        PAM<- lets.presab.points(xy, species, xmn =0, xmx =l,ymn =0, ymx =w,remove.cells=FALSE,resol=g)
  	T=PAM$Presence_and_Absence_Matrix[,3]
  	XX=rbind(XX,T) #PAM 
	N.all=rbind(N.all,test.d$n)
}

N <- dim(XX)[2] #number of the cells




#########Table : summary of lansing data #########
data.frame(names=data.names,N.all,presence.cells=apply(XX,1,sum),presence.rate=apply(XX,1,sum)/N)


#########Table: p-values for the Chi-sq test and testing positive association #########
oupt=NULL
i.end=dim(XX)[1]-1
j.end=dim(XX)[1]
for (i in 1:i.end){
 for (j in (i+1):j.end){
      T1=XX[i,]
	T2=XX[j,]
	p.chi=Chi.t(T1,T2)
        p.bp=BP.t(T1,T2)
	temp=c(i,j,p.chi,p.bp)
	oupt=rbind(oupt,temp)
	}
}
cooccur.lansing <- cooccur(mat = XX, type = "spp_site", thresh = TRUE, spp_names = TRUE)
p.veech=prob.table(cooccur.lansing)[,9]
oupt=cbind(oupt,p.veech)

oupt
#########Table: p-values for the Chi-sq test and testing  negative association #########
oupt=NULL
i.end=dim(XX)[1]-1
j.end=dim(XX)[1]
for (i in 1:i.end){
 for (j in (i+1):j.end){
      	T1=XX[i,]
	T2=XX[j,]
	p.chi=Chi.t(T1,T2)
      	p.bp=BP.t(T1,T2,FALSE)
      	temp=c(i,j,p.chi,p.bp)
	oupt=rbind(oupt,temp)
	}
}

cooccur.lansing <- cooccur(mat = XX, type = "spp_site", thresh = TRUE, spp_names = TRUE)
p.veech=prob.table(cooccur.lansing)[,8]
oupt=cbind(oupt,p.veech)
oupt