extracted posterior estimates for parameter parameters and calculated…

… residuals. Then saved those values for use in the autocorrelation measurements

9_parameter_post_analysis.R

@@ -3,94 +3,55 @@ options(stringsAsFactors = FALSE)
 library(rstan)
 library(ggplot2)
 library(rethinking)
-otr <- read.csv(file.path(datapath, 'otterclean.csv'))
-vl1 <- readRDS(file.path(datapath, 'models/varying1locationpost.RDS'))
-vl2 <- readRDS(file.path(datapath, 'models/varying2locationpost.RDS'))
+av <- read.csv(file.path(datapath, 'allvars.csv'))
+Afp <- readRDS(file.path(datapath, 'models/AlphaFixedPost.RDS'))
+Bfp <- readRDS(file.path(datapath, 'models/BetaFixedPost.RDS'))
+DPfp <- readRDS(file.path(datapath, 'models/DeclinePFixedPost.RDS'))
 
+Amod <- readRDS(file.path(datapath, 'models/AlphaFixedModel.RDS'))
+Bmod <- readRDS(file.path(datapath, 'models/BetaFixedModel.RDS'))
+DPmod <- readRDS(file.path(datapath, 'models/DeclinePFixedModel.RDS'))
 
-source('functions.R')
 
-locs <- levels(factor(otr$location))
+source('functions.R')
 
 ###########################################################################
-otr2 <- otr[,c('location','year','pop')]
-names(otr2) <- c('Site','YearID','O_Otters')
-
 ###########################################################################
 
-##param analysis
-
-b0fixed <- extractpar(vl2, 'beta', location=NA, rows=1)
-b1fixed <- extractpar(vl2, 'beta', location=NA, rows=2)
-
-b0samples <- data.frame(sapply(1:14, function(loc) {
-    extractpar(vl2, 'u', location=loc, rows=1) + b0fixed
-}))
-names(b0samples) <- locs
-
-
-
-b1samples <- data.frame(sapply(1:14, function(loc) {
-    extractpar(vl2, 'u', location=loc, rows=2) + b1fixed
-}))
-names(b1samples) <- locs
-
-b1quants <- t(round(apply(b1samples, 2, quantile, probs=c(0.025, .10, .25,
-                                                          .5, .75,.9, .975)),3))
-
-
-locpredict <- sapply(1:14, function(i) {
-    locpost(vl2, i, 'u', response='function')
-})
-
-
-
-yrs <- 0:10
-df <- expand.grid(1:14, yrs)
-attributes(df)$out.attrs <- NULL
-df$pop <- predictpop(df, locpredict, median)
-allpredictions <- predictpop(df, locpredict, same)
-hpdints <- t(apply(allpredictions, 2, HPDI, prob=0.95))
-attributes(hpdints)$dimnames <- NULL
-
-dfci <- cbind(df, hpdints)
+##alpha
 
-dfci$location <- rep(locs, length(yrs))
-names(dfci) <- c('ID','YearID','P_Otters','Lower95','Upper95','Site')
-dfci$Year <- dfci$YearID + 2012
-dfm <- merge(dfci, otr2, all.x = TRUE)
+Apost <- sapply(1:2, function(column) Afp$beta[,column])
+ACIs <- data.frame(mu=apply(Apost, 2, mean),
+                 lower=apply(Apost, 2, HPDI, prob=0.95)[1,],
+                 upper=apply(Apost, 2, HPDI,  prob=0.95)[2,])
+ACIs
 
-write.csv(dfm, file.path(datapath, 'vis/popplot.csv'), row.names = FALSE)
+Ares <- processFixedMod(Amod, av$latitude, av$alpha, rname='alpha')
 
 
-predict10 <- sapply(locpredict, function(lfun) {
-    do.call(lfun, list(10, median))
-})
+###beta
 
-extinctp10 <- sapply(locpredict, function(lfun) {
-    mean(do.call(lfun, list(10, same))<0)
-})
+Bpost <- sapply(1:2, function(column) Bfp$beta[,column])
+Bdf <- data.frame(mu=apply(Bpost, 2, mean),
+                  lower=apply(Bpost, 2, HPDI, prob=0.95)[1,],
+                  upper=apply(Bpost, 2, HPDI, prob=0.95)[2,])
 
+Bres <- processFixedMod(Bmod, av$latitude, av$beta, rname='beta')
+###DeclineP
 
-pdat <- data.frame(loc=locs,
-                   declineP=sapply(1:14, function(i) {
-                       mean(b1samples[,i]<0)
-                   }),
-                   extinct10=extinctp10,
-                   year5=round(sapply(locpredict, function(lfun) {
-                       do.call(lfun, list(5, median))
-                   }),1),
-                   year10=round(sapply(locpredict, function(lfun) {
-                       do.call(lfun, list(10, median))
-                   }),1))
+DPpost <- sapply(1:2, function(column) DPfp$beta[,column])
+DPdf <- data.frame(mu=apply(DPpost, 2, mean),
+                  lower=apply(DPpost, 2, HPDI, prob=0.95)[1,],
+                  upper=apply(DPpost, 2, HPDI, prob=0.95)[2,])
 
-write.csv(pdat, file.path(datapath,'vis/likelihoods.csv'), row.names = FALSE)
+DPres <- processFixedMod(DPmod, av$latitude, av$declineP, rname='declineP')
 
+#######################################################
+avars <- c('ID','loc','region','latitude','Latitude','Longitude')
 
-paramdf <- data.frame(alpha=apply(b0samples, 2, mean),
-                      beta=apply(b1samples, 2, mean))
-responsedf <- cbind(pdat, paramdf)
+modres <- cbind(Ares, Bres[,2:4], DPres[,2:4])
+ResidDF <- merge(av[,avars], modres, by.x='latitude', by.y='lat')
 
-write.csv(responsedf, file.path(datapath,'responsevars.csv'), row.names = FALSE)
+write.csv(ResidDF, file.path(datapath, 'Residuals.csv'), row.names = FALSE)