-
Notifications
You must be signed in to change notification settings - Fork 0
/
PositiveSim.R
81 lines (64 loc) · 2.18 KB
/
PositiveSim.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
library(spatstat)
library(letsR)
source("D://Co-Occ.R")
##########Simulation function for generating two species with positive association##########
SimPos = function(c1,c2,M){
pv=NULL
ncell=NULL
npoint=NULL
j1 <- function(x,y){ c1*x }
j2 <- function(x,y){ c2*x }
########## Simulation of M times##########
for (i in 1:M){
set.seed(i)
jh1 <- rpoispp(j1)
jh2 <- rpoispp(j2)
j1xy <- cbind(jh1$x, jh1$y)
j2xy <- cbind(jh2$x, jh2$y)
a=rep('species1',jh1$n)
b=rep('species2',jh2$n)
##########PAM##########
PAM1 <- lets.presab.points(j1xy,a,xmn=0,xmx=1,ymn=0,ymx=1,remove.cells=FALSE,resol=0.05)
PAM2 <- lets.presab.points(j2xy,b,xmn=0,xmx=1,ymn=0,ymx=1,remove.cells=FALSE,resol=0.05)
T1=PAM1$Presence_and_Absence_Matrix[,3]
T2=PAM2$Presence_and_Absence_Matrix[,3]
########P-values for testing positive association under Binomial and Possison test###############
p=BP.t(T1,T2)
p.vee=Veech.t(T1,T2)
pv=rbind(pv,c(p,p.vee)) #p-value
n1n2=c(sum(T1),sum(T2))
ncell=rbind(ncell,n1n2)
N1N2=c(jh1$n,jh2$n)
npoint=rbind(npoint,N1N2)
}
output=cbind(npoint,ncell,pv)
return(output)
}
##################################################################################################
#######output for table 2###########################################################################################
c.intensity=c((1:10)*100,(2:5)*1000)
M=500
k=length(c.intensity)
outp=NULL
for (i in 1:k){
c12=c.intensity[i]
temp=SimPos(c12,c12,M)
temp11=temp[,5:11]<0.05
temp1=c(round(apply(temp[,1:4],2,mean)),apply(temp11,2,sum)/500)
outp=rbind(outp,temp1)
}
outp
##################################################################################################
#######output for table 3###########################################################################################
c.intensity=c((1:10)*100,(2:5)*1000)
M=500
k=length(c.intensity)
outp=NULL
for (i in 1:k){
c12=c.intensity[i]
temp=SimPos(400,c12,M)
temp11=temp[,5:11]<0.05
temp1=c(round(apply(temp[,1:4],2,mean)),apply(temp11,2,sum)/500)
outp=rbind(outp,temp1)
}
outp