library (binom); library(epiR); library(doBy)

rm(list=ls())

# run dfr & dff functions below !!!

# -------  INPUTs  --------------

AM= read.table("Random & rt.txt", header=T)  
Nt= 5 # transects
n.em= 99 # sampling stations 


# 100 sampling stations (100 x 100m) per scenario are randomly deployed (five randomly deployed transects each). Counts number of healthy & diseased points. Estimates mean & CI for each sample station data set and for simple random sampling.
# 
f1= dff(Nt)
for(i in 1:n.em) {
  f2= dff(Nt)
    f1=rbind(f1,f2)
}
 f1= f1[, c(1,3,2)]
  colnames(f1)= c("size", "diseased", "N"); head(f1)

#  Simple Random Sampling
a= epi.simplesize(N= length(AM$x), Vsq= NA, Py= 0.16, epsilon.r= 0.05, method= "proportion", conf.level= 0.95)
a1=  sample(length(AM$x), a)
a2= AM[a1,]
 srs= data.frame(size= "SRS", sum(a2[,3]), length(a2[,3]))
colnames(srs)= c("size", "diseased", "N"); srs
 f1=rbind(f1,srs)
 
# Calculates IC  
tt= binom.confint(f1$diseased, f1$N, method= "bayes")
 tt$mean= round(tt$mean,2); tt$lower= round(tt$lower,2); tt$upper= round(tt$upper,2)
  CI= cbind(f1, tt[, c(4:6)])
  CI$CIwidth= CI$upper - CI$lower # width of IC

write.table(CI,"R & rt Prev & CI 3.txt", quote=F, sep=" ", row.names=F)

# ------------------

# 1) function dfr.  Generates a 5 transect sampling station within a 100 by 100 sample space. Transect layout is random. Counts number of diseased/healthy points. Transect size is determined by dff function (below).

dfr=function(xt, yt, Nt) { # xt= transect length; yt= transect width; Nt= transect number
  
# --- Determines sampling station (100 * 100m) spatial layout
  xi=as.integer(runif(1,1,1400)); xf= xi + 100 # initial & final x extremes
  yi=as.integer(runif(1,1,100)); yf= yi + 100 # initial & final y extremes
  EM= subset(AM, x>= xi & x<= xf & y>= yi & y<= yf)
# First transect coordinates
   y.ini= sample(100-yt, Nt) + floor(min(EM$y)) 
   x.ini= sample(100-xt, Nt) + floor(min(EM$x)) 
   y1.i= y.ini[1];  y1.f= y1.i + yt 
   x1.i= x.ini[1];  x1.f= x1.i + xt  
# First transecto selection
 t1=subset(EM, x>= x1.i & x<= x1.f & y>= y1.i & y <=y1.f) 
      diseased= sum(t1[,3]==1); N=length(t1$x)
     tran=1 
est=data.frame(tran, xt, yt, N, diseased, x1.i, x1.f, y1.i, y1.f)
# --- Makes the sampling station and count points --- 
 for (i in 2:Nt) {
   tran= tran + 1
   y2= y.ini[i]; y2.f= y2 + yt 
   x2= x.ini[i]; x2.f= x2 + xt
# Selects the other transects
 tn=subset(EM, x>= x2 & x<= x2.f & y>= y2 & y<= y2.f) 
   enf2= sum(tn[,3]==1); nt= length(tn$x)
   est= rbind(est,c(tran, xt, yt, nt, enf2, x2, x2.f, y2, y2.f))
   }
return(est)
colnames(est)= c("tran", "length", "width", "N", "diseased", "Xini", "Xfin", "Yini", "Yfin")
}

# 2)  función dff: Using dfr function generates 5 belt-tran & sum of diseased/healthy points per sampling station. 

dff= function(Nt) { # Nt = número de transectos deseados
BELT.1.10= dfr(10,1,Nt); BELT.1.25=dfr(25,1,Nt); BELT.1.50=dfr(50,1,Nt)  
BELT.2.10= dfr(10,2,Nt); BELT.2.25=dfr(25,2,Nt); BELT.2.50=dfr(50,2,Nt) 
  trans=rbind(BELT.1.10, BELT.1.25, BELT.1.50, BELT.2.10, BELT.2.25, BELT.2.50)
# trans tiene las coordenadas de cada transecto y sirve para checar la selección 
  colnames(trans)= c("tran", "largo", "ancho", "N", "diseased", "Xini", "Xfin", "Yini", "Yfin")
  trans$largo= as.factor(trans$largo)
  trans$ancho= as.factor(trans$ancho)
  trans$size= paste(trans$largo, "x", trans$ancho, sep="")
ESTACION= summaryBy(N + diseased ~ size, data=trans, FUN= sum)
 return(ESTACION)  
}