5_rstan.R

datapath <- '/Users/echellwig/Google Drive/OtherPeople/otterData/'
options(stringsAsFactors = FALSE)
options(mc.cores = parallel::detectCores())
library(rstan)

rstan_options(auto_write = TRUE)

otr <- read.csv(file.path(datapath, 'otterclean2019.csv'))
linear <- readRDS(file.path(datapath, 'models/linear.RDS'))
vl1 <- readRDS(file.path(datapath, 'models/varying1location.RDS'))
vl2 <- readRDS(file.path(datapath, 'models/varying2location.RDS'))
vlP <- readRDS(file.path(datapath, 'models/varyinglocationPOIS.RDS'))



source('functions.R')

#otr$cpop <- scale(otr$pop) #4.62 center; 2.15 scale
#attributes(otr$cpop) <- NULL

#note if a prior is not specified it is uniform.
##########################################################

locfac <- factor(otr$location)

#data for lm stan model
fixedlist <- list(pop=otr$pop,
                  year=otr$year,
                  N=nrow(otr))

multilist <- list(pop=otr$pop,
                  year=otr$year,
                  N=nrow(otr),
                  P=nlevels(locfac),
                  loc=as.integer(locfac))



##################################
#otter population models


#fixed FX model
fixedl <- stan(model_code=linear, data=fixedlist, iter=20000, warmup = 5000,
               chains=1)
saveRDS(fixedl, file.path(datapath, 'models/fixedpost.RDS'))

#Varying FX location model

#identifying any pathologies in the model
multi1 <- stan(model_code=vl1, data=multilist, iter=25000, warmup=5000,
               chains=1, control=list(adapt_delta = 0.99))

saveRDS(multi1, file.path(datapath, 'models/varying1locationpost.RDS'))


#running the model for estimating parameters etc.
multi2long <- stan(model_code=vl2, data=multilist, iter=25000, warmup=5000,
                  chains=1, control=list(adapt_delta = 0.99))
saveRDS(multi2long, file.path(datapath, 'models/varying2locationpost.RDS'))


multiPOIS <- stan(model_code=vlP, data=multilist, iter=25000, warmup=5000,
                   chains=1, control=list(adapt_delta = 0.99))
saveRDS(multiPOIS, file.path(datapath, 'models/varyinglocPOISpost.RDS'))


#############

#use loo to compare models




###########################################################################