# 30 independent runs for each scenario (counts & proportion models). Results of each validated model run were utilized to estimate the percentage of significan differences between a given scenario survey estimates and Mh real survey estimates. 

rm(list=ls())

library(dispmod); library (vcd); library (statmod); library(MASS); library(binomTools)

tt=read.table("DATA Sims & Mh06 (30 stats, 6 transects each).txt", header= T); head(tt) # File with 25x2m belt transects only, as the purpose is to compare simulated surveys with real Mh surveys (24 transects).Virtual surveys of 24 transects randomly selected from a data set of 180 belt-transects.
 f0= tt[tt$N!= 0,]; summary(f0) # Erase empty belt transects, that do not exist in the Mh data set

# GENERATES DATA SETS for comparisons
# selecting 24 random rows (transects) across all levels of sim; same as in Mh survey.
 
N=24    # Number of rows (ie. transectos) 
temp= lapply(split(f0, f0$sim),  function(z) z[sample(1:nrow(z), N),]) 
t1= do.call('rbind', temp); rownames(t1)= NULL  #  make a data frame

# COUNTS MODEL (NEGBINOMIAL)

distplot(t1$N, type = "nbinomial") # Plots the number of occurrences (counts) against the distribution metameter of the specified distribution. If the distribution fits the data, the plot should show a straight line. 

t1$sim=relevel(t1$sim, ref="Mh") 
m1= glm.nb(N~sim, data=t1); summary(m1)
# validating model by means of randomized quantile residuals plots
qr= qres.nbinom(m1); plot(qr); # qqnorm(qr)


#  BINOMIAL MODEL

t1$sim=relevel(t1$sim, ref="Mh") 
m2= glm(cbind(N,enf) ~ sim, binomial, data=t1)
# validating model by means of randomized quantile residuals plots
qr= qres.binom(m2); plot(qr) # qqnorm(qr)

# model adjusted with no overdispersion
m2d= glm.binomial.disp(m2) 
summary(m2); summary(m2d)