Skip to content

Commit

Permalink
Merge pull request #3 from farrmt/working
Browse files Browse the repository at this point in the history 
Working
  • Loading branch information
@farrmt
farrmt authored Jun 9, 2020
2 parents 6f94b22 + df07c5f commit 2230432
Show file tree
Hide file tree
Showing 7 changed files with 745 additions and 320 deletions.
139 changes: 117 additions & 22 deletions DataAnalysis/HMSNO.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,10 @@ for(i in 1:5){
}

CovData <- bind_rows(CovData, read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 1) %>%
select(`deployment ID`, days, elevation, edge, Year)) %>%
bind_rows(., read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>%
select(`deployment ID`, days, elevation, edge, Year)) %>%
distinct()
select(`deployment ID`, days, elevation, edge, Year)) %>%
bind_rows(., read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>%
select(`deployment ID`, days, elevation, edge, Year)) %>%
distinct()

Data <- full_join(Data, CovData, by = c("deployment ID", "Year")) %>% drop_na()

Expand All @@ -67,8 +67,11 @@ y <- Data2$value

#Indices
nyrs <- Data %>% group_by(`deployment ID`) %>% summarize(min(Year), max(Year))
nstart <- as.numeric(nyrs$`min(Year)`) - 2009 + 1
nend <- 9 - (2017 - as.numeric(nyrs$`max(Year)`))
nstart <- as.numeric(nyrs$`min(Year)`) - as.numeric(min(nyrs$`min(Year)`)) + 1
nend <- 1 + as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) -
(as.numeric(max(nyrs$`max(Year)`)) - as.numeric(nyrs$`max(Year)`))

nyrs <- as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) + 1

nsites <- as.numeric(Data %>% summarize(n_distinct(`deployment ID`)))
nspec <- as.numeric(Data %>% summarize (n_distinct(species)))
Expand All @@ -95,42 +98,134 @@ days <- t(as.matrix(days[,-1]))
density <- Data %>% distinct(`deployment ID`, density)
density <- as.numeric(density$density)

elev[is.na(elev)] <- 0
edge[is.na(edge)] <- 0
days[is.na(days)] <- 0

#--------------#
#-Compile Data-#
#--------------#

#Data
jags.data <- list(y = y, nobs = nobs, nstart = nstart, nend = nend, nsites = nsites, nspec = nspec,
jags.data <- list(y = y, nobs = nobs, nstart = nstart, nend = nend, nsites = nsites, nspec = nspec,
yr = yr, site = site, spec = spec, elev = elev, edge = edge, density = density, days = days)

#Initial values
Nst <- array(NA, dim = c(9, nsites, nspec))
Sst <- array(NA, dim = c(9-1, nsites, nspec))
Gst <- array(NA, dim = c(9-1, nsites, nspec))
Nst <- array(NA, dim = c(nyrs, nsites, nspec))
Sst <- array(NA, dim = c(nyrs-1, nsites, nspec))
Gst <- array(NA, dim = c(nyrs-1, nsites, nspec))
for(j in 1:nsites){
Nst[nstart[j],j,] <- 10
Sst[nstart[j]:(nend[j]-1),j,] <- 5
Gst[nstart[j]:(nend[j]-1),j,] <- 2
Nst[nstart[j],j,] <- 10
Sst[nstart[j]:(nend[j]-1),j,] <- 5
Gst[nstart[j]:(nend[j]-1),j,] <- 2
}

alpha0.fun <- function(){
alpha0 <- NULL
alpha0[1] <- runif(1, -0.2, 0.2)
alpha0[2] <- runif(1, -4, -2.75)
alpha0[3] <- runif(1, -2, -1.5)
alpha0[4] <- runif(1, -5, -3)
alpha0[5] <- runif(1, -1, -0.6)
alpha0[6] <- runif(1, -1.5, -1)
alpha0[7] <- runif(1, -2, -1.5)
alpha0[8] <- runif(1, -4, -2)
alpha0[9] <- runif(1, -2, -1.5)
return(alpha0)
}

beta0.fun <- function(){
beta0 <- NULL
beta0[1] <- runif(1, 0.4, 0.8)
beta0[2] <- runif(1, 1.5, 3.5)
beta0[3] <- runif(1, 0, 1)
beta0[4] <- runif(1, -7, -4)
beta0[5] <- runif(1, 0.5, 1)
beta0[6] <- runif(1, -0.5, 0)
beta0[7] <- runif(1, 1, 3)
beta0[8] <- runif(1, 0, 3)
beta0[9] <- runif(1, 0.5, 1)
return(beta0)
}

inits <- function()list(N=Nst, S=Sst, G=Gst)
beta4.fun <- function(){
beta4 <- NULL
beta4[1] <- runif(1, -0.5, 0.5)
beta4[2] <- runif(1, 0, 0.5)
beta4[3] <- runif(1, -0.5, 0)
beta4[4] <- runif(1, -2, 2)
beta4[5] <- runif(1, 0, 0.6)
beta4[6] <- runif(1, 1, 1.5)
beta4[7] <- runif(1, 0, 2)
beta4[8] <- runif(1, 1, 2)
beta4[9] <- runif(1, 0.4, 0.6)
return(beta4)
}

omega0.fun <- function(){
omega0 <- NULL
omega0[1] <- runif(1, 1.4, 2.2)
omega0[2] <- runif(1, -5, -3)
omega0[3] <- runif(1, 0.5, 1)
omega0[4] <- runif(1, -2, -0.4)
omega0[5] <- runif(1, 0.6, 1)
omega0[6] <- runif(1, 2, 4)
omega0[7] <- runif(1, 0, 1)
omega0[8] <- runif(1, 0, 3)
omega0[9] <- runif(1, 2, 3)
return(omega0)
}

# gamma0.fun <- function(){
# gamma0 <- NULL
# gamma0[1] <- runif(1, 0, 1)
# gamma0[2] <- runif(1, -0.2, 0.2)
# gamma0[3] <- runif(1, -1.5, -0.5)
# gamma0[4] <- runif(1, -2.5, -1.5)
# gamma0[5] <- runif(1, -0.6, -0.2)
# gamma0[6] <- runif(1, -3, -2)
# gamma0[7] <- runif(1, -1.5, -0.5)
# gamma0[8] <- runif(1, -3, -1)
# gamma0[9] <- runif(1, -3, -2)
# return(gamma0)
# }

inits <- function()list(N=Nst, S=Sst, G=Gst,
mu.a0 = runif(1, 0.1, 0.25), tau.a0 = runif(1, 0, 5), alpha0 = alpha0.fun(), alpha1 = runif(1, -0.25, 0.25),
mu.b0 = runif(1, 0, 3), tau.b0 = runif(1, 0, 1), beta0 = beta0.fun(),
mu.b1 = runif(1, -1, 1), tau.b1 = runif(1, 0, 10), beta1 = runif(nspec, -1, 1),
mu.b2 = runif(1, -1, 1), tau.b2 = runif(1, 0, 10), beta2 = runif(nspec, -0.5, 0.5),
mu.b3 = runif(1, -1, 1), tau.b3 = runif(1, 0, 10), beta3 = runif(nspec, -1, 1),
mu.b4 = runif(1, 0, 2), tau.b4 = runif(1, 0, 2), beta4 = beta4.fun(),
mu.o0 = runif(1, 0.5, 0.9), tau.o0 = runif(1, 0, 1), omega0 = omega0.fun(),
mu.o1 = runif(1, 0, 4), tau.o1 = runif(1, 0, 1), omega1 = runif(nspec, -5, 5),
mu.o2 = runif(1, -1, 1), tau.o2 = runif(1, 0, 10), omega2 = runif(nspec, -1, 1),
mu.o3 = runif(1, -1, 1), tau.o3 = runif(1, 0, 10), omega3 = runif(nspec, -1, 1),
gamma0 = runif(1, -4, 0)
)

#Parameters to save
params <- c("mu.a0", "tau.a0",
"mu.b0", "sig.b0", "mu.b1", "tau.b1", "mu.b2", "tau.b2", "mu.b3", "tau.b3", "mu.b4", "tau.b4",
"mu.o0", "sig.o0", "mu.o1", "sig.o1", "mu.o2", "sig.o2", "mu.o3", "sig.o3",
"mu.g0", "sig.g0",
"mu.b0", "tau.b0", "mu.b1", "tau.b1", "mu.b2", "tau.b2", "mu.b3", "tau.b3", "mu.b4", "tau.b4",
"mu.o0", "tau.o0", "mu.o1", "tau.o1", "mu.o2", "tau.o2", "mu.o3", "tau.o3",
"alpha0", "alpha1", "beta0", "beta1", "beta2", "beta3", "beta4",
"omega0", "omega1", "omega2", "omega3", "gamma0")


#MCMC settings
ni <- 115000
# ni <- 115000
# na <- 1000
# nb <- 100000
# nc <- 3
# nt <- 5

na <- 1000
nb <- 100000
nc <- 3
nt <- 5
ni <- 100000
nb <- 20000
nc <- 1
nt <- 40

#Run JAGS model
out <- jagsUI(jags.data, inits, params, model.file = "../MSdyn_cov.txt",
n.adapt = na, n.chains = nc, n.iter = ni, n.burnin = nb, n.thin = nt,
parallel = TRUE)
parallel = FALSE)
Loading

0 comments on commit 2230432

Please sign in to comment.