-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstartup.R
124 lines (119 loc) · 4.71 KB
/
startup.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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
message("Cargando los packages...")
library(TMB)
library(ggplot2)
library(reshape2)
library(plyr)
message("Cargando los datos...")
esfuerzo <- read.csv('datos/Effort.txt')
df <- read.csv("datos/CH_Matrix.csv")
names(df)[1:2] <- c("numero", "evento")
df$evento <- as.factor(df$evento)
individuo <- read.csv("datos/Green_tags_events.csv")
individuo <- individuo[!is.na(individuo$Size),]
100*mean(individuo$Size>115)
individuo$length <- individuo$Size #(individuo$Size-1.469)/1.1578
sum(duplicated(unique(individuo$Tag_num)))
sum(duplicated(unique(df$numero)))
which.missing <- which(! df$numero %in% individuo$Tag_num)
## which.missing <- which(! individuo$Tag_num %in% df$numero )
df[which.missing,]
df <- df[-which.missing,]
nrow(df); nrow(individuo)
df <- merge(x=individuo, y=df, by.x='Tag_num', by.y='numero')
df <- droplevels(subset(df, Sex %in% c('M', 'F')))
## df <- subset(df, length > 100)
CH <- unname(as.matrix(df[, -(1:10)]))
last <- as.numeric(apply(CH, 1, function(xx)
tail(which(!is.na(xx) & xx>0), n=1)))
first <- as.numeric(apply(CH, 1, function(xx)
head(which(!is.na(xx) & xx>0), n=1)))
counts <- unname(t(apply(CH,1, function(xx){
xx[is.na(xx)]=0
cumsum(xx)
})))
data <- list(I=nrow(CH), K=ncol(CH), CH=CH, last=last, counts=counts,
effort=esfuerzo$Trap_haul/1000, lengths=df$length,
first=first, sexo=as.numeric(df$Sex)-1,
evento=as.numeric(df$Event)-1,
lengths_pred=seq(50, 110, len=30),
esfuerzo_pred=seq(.05, 4, len=20))
message("Cargando las funciones...")
##### Una funcion temporaria para hacer una simulacion
simulator <- function(make.plots){
## La entrada de los datos
I <- 7000 # numero de los individuos
K <- 29 # numero de los periodos
## El primero event de las capturas es 1 (2nd octubre hasta 18 octubre
## 2008) con 3000 individuos
## El segundo es en periodo 10 con 2000 individuos
## El ultimo es en 14 con 2000 individuos
##
## los parametros de usar en la simulacion, por ahora son iguales que el
## adjuste de los datos reales.
M <- .09/24
r <- .01
k <- .0001
## Selectividad
a <- .4
b <- 90
## componente de observacion
counts <- matrix(0, I,K)
alive <- p <- CH <- matrix(NA, I,K)
first <- last <- rep(NA, len=I)
p <- matrix(NA, nrow=I, ncol=K) # capture probability
phi <- matrix(NA, I, K); # survival probability
last <- rep(NA, len=I)
effort <- runif(K, min=500, max=2000)
effort[16:24] <- 0
## las longitudes de los individuos de los datos reales
lengths <- sample(na.omit(df$length), size=I, replace=TRUE)
events <- rep(NA, len=I)
for(i in 1:I){
if(i <= 3000) t0 <- 1
if(i > 3000 & i <= 5000) t0 <- 10
if(i > 5000) t0 <- 13
## inicializacion
## estan vivo en periodo uno
first[i] <- t0
counts[i,t0] <- CH[i,t0] <- alive[i,t0] <- 1
phi[i,t0] <- exp(-M-counts[i,t0]*r)
p[i,t0] <- 1 # captura
for(t in (t0+1):K) {
## calcula sobrevivencia
phi[i,t] <- exp(-M-counts[i,t-1]*r)
## calcula probabilidad de recaptura para todos los periodos except el
## primero porque conocido (captura)
p[i,t] <- 1*(1-exp(-k*effort[t]))/(1+exp(-a*(lengths[i]-b)))
## actualiza esta probabilidad por la longitud
## vida ahora solo si vida y sobrevivido en el periodo anterior
alive[i,t] <- alive[i,t-1]*rbinom(n=1, size=1, prob=phi[i,t-1])
## recapturado solo si vida
CH[i,t] <- alive[i,t]*rbinom(n=1, size=1, prob=p[i,t])
counts[i,t] <- counts[i,t-1]+ ifelse(!is.na(CH[i,t]), CH[i,t],0)
}
## el ultimo periodo en que un individuo fue visto
last[i] <- tail(which(CH[i,]==1), n=1)
}
CH[, 16:24] <- NA
simdata <- list(I=I, K=K, CH=CH, last=last, counts=counts, effort=effort,
lengths=lengths, first=first, lengths_pred=seq(60,115, by=1))
simpars <- list(logM=log(M), logr=log(r), logk=log(k), a=a, b=b)
if(make.plots){
par(mfcol=c(2,2))
x <- seq(70,120, len=1000)
plot(x, (1-exp(-k*max(effort)))/(1+exp(-a*(x-b))), ylab='Prob. captura',
type='l', xlab='Longitud')
legend('topleft', legend=c('Max. Esfuerzo', 'Min. Esfuerzo'), lty=1, col=c(1,2))
lines(x, (1-exp(-k*min(effort[effort>0])))/(1+exp(-a*(x-b))), col='red')
hist(lengths, xlim=c(70,120), xlab='Longitud')
plot(0,0, xlim=c(1,K), ylim=c(0,.25), xlab='Periodo',
ylab='Pr. captura de 5 individuos',
type='n')
trash <- sapply(1:5, function(i) lines(2:K, p[i,-1], col=1))
plot(0,0, xlim=c(1,K), ylim=c(.8,1), xlab='Periodo',
ylab='Pr. sobrevivencia de 5 individuos',
type='n')
trash <- sapply(1:5, function(i) lines(2:K, phi[i,-1], col=1))
}
return(list(simdata=simdata, simpars=simpars, p=p, phi=phi))
}