# Makes ppp's for three distribution patterns and two disease transmission modes

library(spatstat); library(binom)

 rm(list=ls())
  long= 1500; wide= 200 # Sample space dimentions
  dens=.68 # mean density Mahahual 2006
  prev= 16.2  # mean prevalence Mahahual 2006

# INPUTS
  ncol= as.integer(dens * (long * wide))
  nenf= ncol * (prev/100)
  ss=owin(c(0, long), c(0,wide)) 

 # 1. Random 
 set.seed(123)
 Rand= rpoint(n= ncol, win= ss)
 
 # 2. Dispersed  
 set.seed(123)   
 Disp= rSSI(r= 0.5, n= ncol, win= ss, giveup= 100000, x.init= NULL)
    
 # 3. Clustered
 set.seed(123)   
 rad= 2; mu= 15; k1=dens/mu 
 Clust= rMatClust(kappa= k1, r= rad, mu= mu, win= ss)
 
 
 # ----   Ripley´s K normalization (Besaq´s L) function --------     
 Lr=function(PP) { 
   RE= envelope(PP, fun="Lest", nsim=20, nrank=1); L2 = RE 
   L2$obs = L2$obs-L2$r;  L2$theo= L2$theo-L2$r 
   L2$lo= L2$lo-L2$r;  L2$hi= L2$hi-L2$r 
   # Plots Ripley’s L function with 95% CI envelopes
   plot(L2, ylab= "L transf. of Ripley's K", main="", legend=F)
 }
      
# ---  ASSIGN diseased marks to the corresponding ppp:   

PP= Rand   # Choose between Rand, Disp & Clust ppp's
 
# 1)  RANDOM DISEASE TRANSMISSION
   cond=rbinom(PP$n, 1, prev/100) # assign diseased points
   marks(PP)= as.factor(cond)
    df1=as.data.frame(PP); colnames(df1)[3]="cond"
      df1$sim= rep(x, length(df1$x)) # x= ppp´s name
# SCENARIOS DATA: X: type of point distribution; x: type of transmission
    write.table(df1, "X & x DATA.txt", quote=F, sep=" ", row.names=F) 

# 2)  TRANSMISSION  by  DISTANCE ~ PREVALENCE
  nn=sort(nndist(PP, k=1)) # ordering vector by neareast distances
  dd=nn[nn<=nn[nenf]] # as nn is an ordered vector selecting distances <=  nenf. 
# Distance thresholds 
  d1= median(dd); d2=quantile(dd,.80); d3=quantile(dd,.90); d4=quantile(dd,.99) 
  cond=nn
    x1=which(nn<= d1); cond[x1]=1
    x2=which(nn> d1 & nn<= d2); cond[x2]=rbinom(length(x2), 1, 0.9)
    x3=which(nn> d2 & nn<= d3); cond[x3]=rbinom(length(x3), 1, 0.8)
    x4=which(nn> d3 & nn<= d4); cond[x4]=rbinom(length(x4), 1, 0.7) 
    x5=which(nn> d4); cond[x5]=rbinom(length(x5), 1, 0.02) 
   marks(PP)=as.factor(cond)
     df2=as.data.frame(PP); colnames(df2)[3]="cond"
       df2$sim= rep("x", length(df2$x)) # x= ppp´s name

# SCENARIOS DATA: X: type of point distribution; x: type of transmission
   write.table(df2, "X & x DATA.txt", quote=F, sep=" ", row.names=F)

#  --- CONFIRMS DISTRIBUCION type --- 

#  Ripley´s K standarization (Besaq´s L) of the ppp.    
   set.seed(123); Lr(Y)  # Y= ppp (Rand, Disp or Clust) 
   
#  Ripley´s K standarization (Besaq´s L) of random diseased points 
    t1= df1[df1$cond== 1,]; t1= droplevels(t1)
     t2= as.ppp(t1,ss)
  set.seed(321); Lr(t2)
 
#  Ripley´s K standarization (Besaq´s L) of nearest neighbour diseased points   
    t3= df2[df2$cond== 1,]; t3= droplevels(t3)
  t4= as.ppp(t3,ss)
  set.seed(321); Lr(t4)