-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathHierarchicalBayesianQuantileRegressionModels.R
158 lines (147 loc) · 5.5 KB
/
HierarchicalBayesianQuantileRegressionModels.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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
# ----------------------------------------------------------------------------------------------------------------------------------- #
# ----------------------------------------------------------------------------------------------------------------------------------- #
# ----------------------------------------------------------------------------------------------------------------------------------- #
# Citation: Precipitation Scaling with Temperature in the Northeast US: Variations by Weather Regime, Season, and Precipitation Intensity
# Geophysical Research Letters, 2021
# All rights reserved: Nasser Najibi, Sudarshana Mukhopadhyay, Scott Steinschneider
# Department of Biological and Environmental Engineering, Cornell University, Ithaca 14853, NY, United States, [email protected]
# -----------------------------------------------------------------------------------------------------------------------------------
# Version 1.0 August 2021
# Hierarchical Bayesian quantile Regression Model #
# for Td
bayes_cc_qreg_wr_multilevel <- function() {
for(i in 1:nsites)
{
for(t in 1:nevents[i])
{
# Level 1 #
precip[i,t] ~ dnorm(me[i,t],pe[i,t])
me[i,t] <- (1-2*p)/(p*(1-p))*w[i,t] + mu[i,t]
mu[i,t] <- alpha.wrs[i,t] + beta.wrs[i,t]*td[i,t]
pe[i,t] <- (tau.q[i]*p*(1-p))/(2*w[i,t])
w[i,t] ~ dexp(tau.q[i])
# Level 2 #
beta.wrs[i,t] <- b1[i]*wr1[i,t]+b2[i]*wr2[i,t]+b3[i]*wr3[i,t]+b4[i]*wr4[i,t]
alpha.wrs[i,t] <- a1[i]*wr1[i,t]+a2[i]*wr2[i,t]+a3[i]*wr3[i,t]+a4[i]*wr4[i,t]
}
# Prior #
lsigma[i] ~ dunif(-10,10)
sigma[i] <- exp(lsigma[i]/2)
tau.q[i] <- pow(sigma[i],-2)
a1[i] ~ dnorm(0,0.01)
a2[i] ~ dnorm(0,0.01)
a3[i] ~ dnorm(0,0.01)
a4[i] ~ dnorm(0,0.01)
b1[i] ~ dnorm(mu.b1,tau.b1)
b2[i] ~ dnorm(mu.b2,tau.b2)
b3[i] ~ dnorm(mu.b3,tau.b3)
b4[i] ~ dnorm(mu.b4,tau.b4)
}
# Prior #
mu.b1 ~ dnorm(0,1)
mu.b2 ~ dnorm(0,1)
mu.b3 ~ dnorm(0,1)
mu.b4 ~ dnorm(0,1)
tau.b1 ~ dgamma(1,0.1)
tau.b2 ~ dgamma(1,0.1)
tau.b3 ~ dgamma(1,0.1)
tau.b4 ~ dgamma(1,0.1)
sigma.b1 <- 1/sqrt(tau.b1)
sigma.b2 <- 1/sqrt(tau.b2)
sigma.b3 <- 1/sqrt(tau.b3)
sigma.b4 <- 1/sqrt(tau.b4)
}
# -----------------------------------------------------------------------------------------------------------------------------------
# Hierarchical Bayesian Quantile Regression Model #
# for Ta
bayes_cc_qreg_wr_multilevel <- function() {
for(i in 1:nsites)
{
for(t in 1:nevents[i])
{
# Level 1 #
precip[i,t] ~ dnorm(me[i,t],pe[i,t])
me[i,t] <- (1-2*p)/(p*(1-p))*w[i,t] + mu[i,t]
mu[i,t] <- alpha.wrs[i,t] + beta.wrs[i,t]*ta[i,t]
pe[i,t] <- (tau.q[i]*p*(1-p))/(2*w[i,t])
w[i,t] ~ dexp(tau.q[i])
# Level 2 #
beta.wrs[i,t] <- b1[i]*wr1[i,t]+b2[i]*wr2[i,t]+b3[i]*wr3[i,t]+b4[i]*wr4[i,t]
alpha.wrs[i,t] <- a1[i]*wr1[i,t]+a2[i]*wr2[i,t]+a3[i]*wr3[i,t]+a4[i]*wr4[i,t]
}
# Prior #
lsigma[i] ~ dunif(-10,10)
sigma[i] <- exp(lsigma[i]/2)
tau.q[i] <- pow(sigma[i],-2)
a1[i] ~ dnorm(0,0.01)
a2[i] ~ dnorm(0,0.01)
a3[i] ~ dnorm(0,0.01)
a4[i] ~ dnorm(0,0.01)
b1[i] ~ dnorm(mu.b1,tau.b1)
b2[i] ~ dnorm(mu.b2,tau.b2)
b3[i] ~ dnorm(mu.b3,tau.b3)
b4[i] ~ dnorm(mu.b4,tau.b4)
}
# Prior #
mu.b1 ~ dnorm(0,1)
mu.b2 ~ dnorm(0,1)
mu.b3 ~ dnorm(0,1)
mu.b4 ~ dnorm(0,1)
tau.b1 ~ dgamma(1,1)
tau.b2 ~ dgamma(1,1)
tau.b3 ~ dgamma(1,1)
tau.b4 ~ dgamma(1,1)
sigma.b1 <- 1/sqrt(tau.b1)
sigma.b2 <- 1/sqrt(tau.b2)
sigma.b3 <- 1/sqrt(tau.b3)
sigma.b4 <- 1/sqrt(tau.b4)
}
# -----------------------------------------------------------------------------------------------------------------------------------
# Bayesian Quantile Regression Model #
# for Td
bayes_cc_qreg_wr_singlelevel <- function() {
for(i in 1:nsites)
{
for(t in 1:nevents[i])
{
# single Level #
precip[i,t] ~ dnorm(me[i,t],pe[i,t])
me[i,t] <- (1-2*p)/(p*(1-p))*w[i,t] + mu[i,t]
mu[i,t] <- alpha.wrs[i] + beta.wrs[i]*td[i,t]
pe[i,t] <- (tau.q[i]*p*(1-p))/(2*w[i,t])
w[i,t] ~ dexp(tau.q[i])
}
# Priors #
alpha.wrs[i] ~ dnorm(0,0.01)
beta.wrs[i] ~ dnorm(0,0.01)
lsigma[i] ~ dunif(-10,10)
sigma[i] <- exp(lsigma[i]/2)
tau.q[i] <- pow(sigma[i],-2)
}
}
# -----------------------------------------------------------------------------------------------------------------------------------
# Bayesian Quantile Regression Model #
# for Ta
bayes_cc_qreg_wr_singlelevel <- function() {
for(i in 1:nsites)
{
for(t in 1:nevents[i])
{
# single Level #
precip[i,t] ~ dnorm(me[i,t],pe[i,t])
me[i,t] <- (1-2*p)/(p*(1-p))*w[i,t] + mu[i,t]
mu[i,t] <- alpha.wrs[i] + beta.wrs[i]*ta[i,t]
pe[i,t] <- (tau.q[i]*p*(1-p))/(2*w[i,t])
w[i,t] ~ dexp(tau.q[i])
}
# Priors #
alpha.wrs[i] ~ dnorm(0,0.01)
beta.wrs[i] ~ dnorm(0,0.01)
lsigma[i] ~ dunif(-10,10)
sigma[i] <- exp(lsigma[i]/2)
tau.q[i] <- pow(sigma[i],-2)
}
}
# ----------------------------------------------------------------------------------------------------------------------------------- #
# ----------------------------------------------------------------------------------------------------------------------------------- #
# ----------------------------------------------------------------------------------------------------------------------------------- #