-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy path3.1_Random_Walk_Fit.R
74 lines (59 loc) · 2.04 KB
/
3.1_Random_Walk_Fit.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
#' Function for the random walk model. Observations are modelled under a Poisson, while the the species is
#' modelled with normal process error
#'
#' @param county.name Name of county of interest
#' @param spp species, one of "albo" or "aegypti"
#' @param data.set data
#' @param n.iter number of iterations, default = 5000
#' @param inits initial conditions, default = NULL
Random_Walk_Fit <- function(county.name, spp, data.set, inits = NULL, ...){
# get county of interest and create a "year-month" column
county.sub <- data.fit %>%
filter(state_county == counties[i])
for(i in 1:nrow(county.sub)){
if(county.sub$month[i]==10){
county.sub$new.month[i] <- 91
} else if(county.sub$month[i]==11){
county.sub$new.month[i] <- 92
} else if(county.sub$month[i]==12){
county.sub$new.month[i] <- 93
} else {
county.sub$new.month[i] <- county.sub$month[i]
}
}
county.sub <- county.sub %>%
unite("year_month", year, month, sep = "-", remove = FALSE) %>%
unite("year_month_seq", year, new.month, sep = "-")
# aggregate counts for each month, as they are separated by trap type
y.albo <- aggregate(county.sub$num_albopictus_collected, by = list(county.sub$year_month_seq), FUN = sum)[,2]
y.aegypti <- aggregate(county.sub$num_aegypti_collected, by = list(county.sub$year_month_seq), FUN = sum)[,2]
# use appropriate data (which species?)
if(spp == "albo"){
y <- y.albo
} else {
y <- y.aegypti
}
# create data list for JAGS
data <- list(y = y,
n.month = length(y))
model <- "
model{
#### Data Model
for(i in 1:n.month){
y[i] ~ dpois(x[i])
}
#### Process Model
for(i in 2:n.month){
x[i] ~ dnorm(x[i-1], tau_proc)
}
#### Priors
x[1] ~ dpois(5)
tau_proc ~ dgamma(0.01,0.01)
}"
j.model <- jags.model(file = textConnection(model),
data = data,
n.chains = 3,
inits = inits,
...)
return(j.model)
}