gistfile1.txt

#! Simulações de mortalidade de aves
#! Maria Dias e Paulo Eduardo Cardoso
setwd("C:/...")
visf<-read.table("visf.csv", sep=";", header=T) #tabela inicial, que tem os dias do ano em que foram feitas as prospecções, segundo cada um dos esquemas a analisar (MAp,MAs,MAm,Maq1 e Maq2); esta tabela tem também os valores de taxa de permanencia (r) associada a cada visita
load("txperm.fun")


######## 1. SIMULAÇÃO DE CENÁRIOS DE MORTALIDADE

### 1.1. criação de tabelas com diferentes cenários de mortalidade real; nrep: número de simulações para cada valor de mortalidade real

p=0.21 #probabilidade de detecção
nrep=999 #número de repetições
cens<-c(5,50,100,200) #define os cenários de mortalidade real

Mapf<-data.frame()



for (g in 1:length(cens)) #--------faz variar a mortalidade real entre os valores definidos em "cens"
{

y<-cens[g]


mapc1<-numeric()
mapc2<-numeric()



for (x in 1:nrep) # -----------faz repetir a simulação, para o valor de mortalidade real definido acima (em y), "nrep" vezes
{

mort2=data.frame(pp=sort(trunc(runif(y,1,336))), mort=rep(1,y)) #cria uma tabela com os uma simulação dos dias em que morreram animais, assumindo uma distribuição uniforme das mortes ao longo do ano, e assumindo que um ano tem 336 dias (i.e., 12 meses*4 semanas*7 dias)



###1.2 cálculo, para cada um dos esquemas de prospecção, da estimativa de mortalidade para o cenário criado acima




#-------------------------Esquema MAp------------
tab=visf[visf$map>0,] #faz uma sub-selecção da tabela visf, apenas com os dias nos quais foram feitos prospecções segundo o esquema em análise
tab$pe<-rep(9999,nrow(tab)) #cria uma coluna que irá ser a probabilidade de encontro (permanecia*detecção*não ter detectado em visitas anteriores)
tab$pe2<-rep(9999,nrow(tab)) 
tab$pe[1]<-mean(txperm(0:tab$mapd[1]))*p #probabilidade de encontrar um cadáver na 1ª visita
tab$pe[2]<-mean(txperm(0:(tab$mapd[2]-tab$mapd[1])))*p+mean(txperm((tab$mapd[2]-tab$mapd[1]):(tab$mapd[2])))*p*(1-p)


for (j in 3:nrow(tab)) #
{
a<-mean(txperm((tab$mapd[j]-tab$mapd[1]):(tab$mapd[j])))*p*(1-p)^(j-1) #prob de encontro de um cadáver que tenha morrido na primeira semana
for (i in 0:(j-2)) #define as semanas anteriores às ultimas visitas
{
a<-c(a,mean(txperm((tab$mapd[j]-tab$mapd[j-i]):(tab$mapd[j]-tab$mapd[j-i-1])))*p*((1-p)^(i))) #calcula a probabilidade de encontro para as semanas anteriores à visita
a
}
tab$pe[j]<-sum(a)
}


tab$pe2[1]<-mean(txperm(0:tab$mapd[1]))*p

for (j in 2:nrow(tab)) tab$pe2[j]<-mean(txperm(0:(tab$mapd[j]-tab$mapd[j-1])))*p




denc<-numeric()
for (z in 1:nrow(mort2)) #faz variar os dias em que ocorreram mortes (i.e, vai percorrer todas as linhas da tabela mort2, analisando as mortes uma a uma)
{
prob<-numeric()
pos=tab$mapd[tab$mapd>=mort2$pp[z]] #dias possíveis p encontrar o cadáver (i.e., dias das visitas que foram feitas após a morte em análise)
if(length(pos)>0) for(w in 1:length(pos)) prob[w]<-rbinom(mort2$mort[z],1,txperm(pos[w]-mort2$pp[z])*p*(1-p)^(w-1)) #determina se o cadáver foi encontrado ou não, em cada dia em que tal é possível
else (prob=0) #no caso de não haver nenhuma visita após a morte
if (sum(prob)==0) denc[z]=0 # no caso de o cadáver nunca vir a ser encontrado
else(denc[z]<-pos[prob==1][1]) #dia em que o cadáver foi encontrado
}
mort2$dencMAp<-denc #cria uma nova coluna na tabela mort2, com os dias em que foram encontrados os cadáveres; no caso do cadáver não ser encontrado, toma o valor 0
if(sum(mort2$dencMAp>0)) #no caso de pelo menos um cadáver ter sido encontrado...
{
mort3<-aggregate(mort2$mort[mort2$dencMAp>0], list(dencMAp=mort2$dencMAp[mort2$dencMAp>0]), sum) #cria uma tabela com o número total de cadáveres encontrados (coluna "cadáveres") por dia de visita (coluna "dencMAp")
mort3$dencMAp<-as.numeric(as.character(mort3$dencMAp))
names(mort3)[2]<-"cadaveres"
mort3$pe<-tab$pe[which(tab$mapd%in%mort3$dencMAp)] #no caso do equema MAp, os valores assumidos para a taxa de permanencia variam consoante o dia da prospoecção; nesta operação são adicionados os valores de r da tabela tab (q é uma sub-selecção da visf) à tabela mort3
mort3$pe2<-tab$pe2[which(tab$mapd%in%mort3$dencMAp)] #no caso do equema MAp, os valores assumidos para a taxa de permanencia variam consoante o dia da prospoecção; nesta operação são adicionados os valores de r da tabela tab (q é uma sub-selecção da visf) à tabela mort3


#estimativa da mortalidade segundo o esquema em análise
MAp=sum(mort3$cadaveres/(mort3$pe)) 
MAp2=sum(mort3$cadaveres/(mort3$pe2))


}
else #no caso de nenhum cadáver ter sido encontrado
{
mort3=0
MAp=0
MAp2=0
}

mapc1<-c(mapc1,MAp) #guarda o valor de mortalidade estimada, para o cenário de mortalidade real e p a simulação em causa (das "nrep" simulações)
mapc2<-c(mapc2,MAp2)
}

###1.3.acumulação dos valores das várias simulações das estimativas de mortalidade

Mapf<-rbind(Mapf,data.frame(mortreal=rep(y,nrep),Map=mapc1, Mapnc=mapc2))

}

Mapf$erro<-(Mapf$Map-Mapf$mortreal)*100/Mapf$mortreal
Mapf$erronc<-(Mapf$Mapnc-Mapf$mortreal)*100/Mapf$mortreal



#1.4. boxplots para cada um dos cenários de mortalidade acima definidos
par(las=1, mfrow=c(1,2))
boxplot(erro~mortreal, data=Mapf, main="erro estimativa (com correcção)", xlab="mortalidade real", ylab="erro na mortalidade estimada (%)")
abline(h=c(-50,0,50), col="red", lty=c(3,2,3))
boxplot(erronc~mortreal, data=Mapf, main="erro estimativa (sem correcção)", xlab="mortalidade real", ylab="erro na mortalidade estimada (%)")
abline(h=c(-50,0,50), col="red", lty=c(3,2,3))