-
Notifications
You must be signed in to change notification settings - Fork 0
/
9_parameter_post_analysis.R
133 lines (90 loc) · 4.17 KB
/
9_parameter_post_analysis.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
125
126
127
128
129
130
131
132
133
datapath <- '/Users/echellwig/Drive/OtherPeople/otterData'
options(stringsAsFactors = FALSE)
library(rstan)
library(ggplot2)
library(rethinking)
library(reshape2)
library(loo)
source('functions.R')
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'))
DPbr <- readRDS(file.path(datapath, 'models/DeclinePBetaRegPost.RDS'))
Amod <- readRDS(file.path(datapath, 'models/AlphaFixedModel.RDS'))
Aint_mod <- readRDS(file.path(datapath, 'models/AlphaInterceptModel.RDS'))
Bmod <- readRDS(file.path(datapath, 'models/BetaFixedModel.RDS'))
Bint_mod <- readRDS(file.path(datapath, 'models/BetaInterceptModel.RDS'))
DPmod <- readRDS(file.path(datapath, 'models/DeclinePFixedModel.RDS'))
DPbrmod <- readRDS(file.path(datapath, 'models/DeclinePBetaRegModel.RDS'))
DPbIntmod <- readRDS(file.path(datapath,
'models/DeclinePBetaInterceptModel.RDS'))
Afit <- modfit(Amod)
Bfit <- modfit(Bmod)
DPfit <- modfit(DPbrmod, betareg = TRUE)
####
respvars <- c('alpha','beta','declineP')
char <- c('PopSize', 'PopTrend','DeclineProb')
###########################################################################
###########################################################################
##alpha
Apost <- sapply(1:2, function(column) Afp$beta[,column])
Adf <- data.frame(char='IGS', par=c('intercept','Latitude'),
mu=apply(Apost, 2, mean),
lower=apply(Apost, 2, HPDI, prob=0.975)[1,],
upper=apply(Apost, 2, HPDI, prob=0.975)[2,])
Ares <- processFixedMod(Amod, av$latitude, av$alpha, rname='alpha')
###beta
Bpost <- sapply(1:2, function(column) Bfp$beta[,column])
Bdf <- data.frame(char='PGR', par=c('intercept','Latitude'),
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
DPpost <- sapply(1:2, function(column) DPbr$beta[,column])
DPdf <- data.frame(char='Decline', par=c('intercept','Latitude'),
mu=apply(DPpost, 2, mean),
lower=apply(DPpost, 2, HPDI, prob=0.95)[1,],
upper=apply(DPpost, 2, HPDI, prob=0.95)[2,])
DPres <- processFixedMod(DPbrmod, av$latitude, av$declineP, rname='declineP')
#######################################################
avars <- c('ID','loc','region','latitude','Latitude','Longitude')
modres <- cbind(Ares, Bres[,2:4], DPres[,2:4])
ResidDF <- merge(av[,avars], modres, by.x='latitude', by.y='lat')
write.csv(ResidDF, file.path(datapath, 'Residuals.csv'), row.names = FALSE)
#######################################################
#37.7-38.3 centered: 37.98501, scaled:0.10711
brvars <- c('declineP','alpha','beta','Latitude','latitude')
avr <- av[,brvars]
avr$type <- 'obs'
LatX <- seq(37.8, 38.22, by=0.01)
latX <- (LatX - 37.98501)/0.10711
avfit <- data.frame(Latitude=LatX,
latitude=latX,
type='fit',
alpha=Afit(latX),
beta=Bfit(latX),
declineP=DPfit(latX))
AV <- rbind(avr, avfit)
avm <- melt(AV, id.vars=c('Latitude','latitude','type'),
variable.name = 'response', value.name = 'value')
write.csv(avm, file.path(datapath, 'LatitudeModelResults.csv'),
row.names = FALSE)
#######################################################
#crossvalidation
estimates <- rbind(Adf, Bdf, DPdf)
write.csv(estimates, file.path(datapath, 'ParameterEstimates.csv'),
row.names = FALSE)
Aloglike <- extract_log_lik(Amod, parameter='loglik')
Aintlike <- extract_log_lik(Aint_mod, parameter='loglik')
Bloglike <- extract_log_lik(Bmod, parameter='loglik')
Bintlike <- extract_log_lik(Bint_mod, parameter='loglik')
DPloglike <- extract_log_lik(DPbrmod, parameter='log_lik')
DPintlike <- extract_log_lik(DPbIntmod, parameter='log_lik')
Aslope <- loo(Aloglike)
Aint <- loo(Aintlike)
Bslope <- waic(Bloglike)
Bint <- waic(Bintlike)
DPslope <- loo(DPloglike)
DPint <- loo(DPintlike)