Skip to content

Commit

Permalink
harmonise parameter names; begin adding support for other exponential…
Browse files Browse the repository at this point in the history 
… families
  • Loading branch information
@nicholasjclark
nicholasjclark committed Mar 28, 2023
1 parent 77fc921 commit 7b8472c
Show file tree
Hide file tree
Showing 51 changed files with 1,204 additions and 1,104 deletions.
40 changes: 20 additions & 20 deletions R/add_base_dgam_lines.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,23 +223,23 @@ add_base_dgam_lines = function(use_lv, stan = FALSE, offset = FALSE){
}
for (j in 1:n_lv) {
LV_raw[2, j] ~ dnorm(phi[j] + ar1[j]*LV_raw[1, j], penalty[j])
LV_raw[2, j] ~ dnorm(drift[j] + ar1[j]*LV_raw[1, j], penalty[j])
}
for (j in 1:n_lv) {
LV_raw[3, j] ~ dnorm(phi[j] + ar1[j]*LV_raw[2, j] + ar2[j]*LV_raw[1, j], penalty[j])
LV_raw[3, j] ~ dnorm(drift[j] + ar1[j]*LV_raw[2, j] + ar2[j]*LV_raw[1, j], penalty[j])
}
for (i in 4:n) {
for (j in 1:n_lv) {
LV_raw[i, j] ~ dnorm(phi[j] + ar1[j]*LV_raw[i - 1, j] +
LV_raw[i, j] ~ dnorm(drift[j] + ar1[j]*LV_raw[i - 1, j] +
ar2[j]*LV_raw[i - 2, j] + ar3[j]*LV_raw[i - 3, j], penalty[j])
}
}
## AR components
for (s in 1:n_lv) {
phi[s] ~ dnorm(0, 10)
drift[s] ~ dnorm(0, 10)
ar1[s] ~ dnorm(0, 10)
ar2[s] ~ dnorm(0, 10)
ar3[s] ~ dnorm(0, 10)
Expand Down Expand Up @@ -318,24 +318,24 @@ add_base_dgam_lines = function(use_lv, stan = FALSE, offset = FALSE){
## likelihood functions
for (i in 1:n) {
for (s in 1:n_series) {
y[i, s] ~ dnegbin(rate[i, s], r[s])T(, upper_bound[s]);
rate[i, s] <- ifelse((r[s] / (r[s] + mus[i, s])) < min_eps, min_eps,
(r[s] / (r[s] + mus[i, s])))
y[i, s] ~ dnegbin(rate[i, s], phi[s])T(, upper_bound[s]);
rate[i, s] <- ifelse((phi[s] / (phi[s] + mus[i, s])) < min_eps, min_eps,
(phi[s] / (phi[s] + mus[i, s])))
}
}
## complexity penalising prior for the overdispersion parameter;
## where the likelihood reduces to a 'base' model (Poisson) unless
## the data support overdispersion
for (s in 1:n_series) {
r[s] <- 1 / r_inv[s]
r_inv[s] ~ dexp(5)
phi[s] <- 1 / phi_inv[s]
phi_inv[s] ~ dexp(5)
}
## posterior predictions
for (i in 1:n) {
for (s in 1:n_series) {
ypred[i, s] ~ dnegbin(rate[i, s], r)T(, upper_bound[s])
ypred[i, s] ~ dnegbin(rate[i, s], phi[s])T(, upper_bound[s])
}
}
Expand All @@ -361,22 +361,22 @@ add_base_dgam_lines = function(use_lv, stan = FALSE, offset = FALSE){
}
for (s in 1:n_series) {
trend[2, s] ~ dnorm(phi[s] + ar1[s]*trend[1, s], tau[s])
trend[2, s] ~ dnorm(drift[s] + ar1[s]*trend[1, s], tau[s])
}
for (s in 1:n_series) {
trend[3, s] ~ dnorm(phi[s] + ar1[s]*trend[2, s] + ar2[s]*trend[1, s], tau[s])
trend[3, s] ~ dnorm(drift[s] + ar1[s]*trend[2, s] + ar2[s]*trend[1, s], tau[s])
}
for (i in 4:n) {
for (s in 1:n_series){
trend[i, s] ~ dnorm(phi[s] + ar1[s]*trend[i - 1, s] + ar2[s]*trend[i - 2, s] + ar3[s]*trend[i - 3, s], tau[s])
trend[i, s] ~ dnorm(drift[s] + ar1[s]*trend[i - 1, s] + ar2[s]*trend[i - 2, s] + ar3[s]*trend[i - 3, s], tau[s])
}
}
## AR components
for (s in 1:n_series){
phi[s] ~ dnorm(0, 10)
drift[s] ~ dnorm(0, 10)
ar1[s] ~ dnorm(0, 10)
ar2[s] ~ dnorm(0, 10)
ar3[s] ~ dnorm(0, 10)
Expand All @@ -387,24 +387,24 @@ add_base_dgam_lines = function(use_lv, stan = FALSE, offset = FALSE){
## likelihood functions
for (i in 1:n) {
for (s in 1:n_series) {
y[i, s] ~ dnegbin(rate[i, s], r[s])T(, upper_bound[s]);
rate[i, s] <- ifelse((r[s] / (r[s] + mus[i, s])) < min_eps, min_eps,
(r[s] / (r[s] + mus[i, s])))
y[i, s] ~ dnegbin(rate[i, s], phi[s])T(, upper_bound[s]);
rate[i, s] <- ifelse((phi[s] / (phi[s] + mus[i, s])) < min_eps, min_eps,
(phi[s] / (phi[s] + mus[i, s])))
}
}
## complexity penalising prior for the overdispersion parameter;
## where the likelihood reduces to a 'base' model (Poisson) unless
## the data support overdispersion
for (s in 1:n_series) {
r[s] <- 1 / r_inv[s]
r_inv[s] ~ dexp(5)
phi[s] <- 1 / phi_inv[s]
phi_inv[s] ~ dexp(5)
}
## posterior predictions
for (i in 1:n) {
for (s in 1:n_series) {
ypred[i, s] ~ dnegbin(rate[i, s], r[s])T(, upper_bound[s])
ypred[i, s] ~ dnegbin(rate[i, s], phi[s])T(, upper_bound[s])
}
}
Expand Down
2 changes: 1 addition & 1 deletion R/add_poisson_lines.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#' @return A modified `JAGS` model file
add_poisson_lines = function(model_file, upper_bounds){

odis_begin <- grep('r\\[s\\] <- ', model_file) - 4
odis_begin <- grep('phi\\[s\\] <- ', model_file) - 4
odis_end <- odis_begin + 7
model_file <- model_file[-c(odis_begin:odis_end)]

Expand Down
16 changes: 8 additions & 8 deletions R/add_stan_data.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,16 +35,16 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE,

if(family == 'nb'){
stan_file[grep('// raw basis', stan_file) + 2] <-
'\n// negative binomial overdispersion\nvector<lower=0>[n_series] r_inv;\n'
'\n// negative binomial overdispersion\nvector<lower=0>[n_series] phi_inv;\n'

stan_file[grep('// priors for smoothing', stan_file) + 2] <-
paste0('\n// priors for overdispersion parameters\n',
'r_inv ~ student_t(3, 0, 0.1);\n')
'phi_inv ~ student_t(3, 0, 0.1);\n')

to_negbin <- gsub('poisson_log', 'neg_binomial_2',
stan_file[grep('y[i, s] ~ poisson', stan_file, fixed = T)])
stan_file[grep('y[i, s] ~ poisson', stan_file, fixed = T)] <-
gsub(');', ', inv(r_inv[s]));', to_negbin)
gsub(');', ', inv(phi_inv[s]));', to_negbin)

add_exp_open <- gsub('\\(eta', '(exp(eta',
stan_file[grep('y[i, s] ~ neg_binomial', stan_file, fixed = T)])
Expand All @@ -55,16 +55,16 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE,

stan_file[grep('matrix[n, n_series] ypred;', stan_file, fixed = T)] <-
paste0('matrix[n, n_series] ypred;\n',
'matrix[n, n_series] r_vec;\n',
'vector[n_series] r;\n',
'r = inv(r_inv);\n',
'matrix[n, n_series] phi_vec;\n',
'vector[n_series] phi;\n',
'phi = inv(phi_inv);\n',
'for (s in 1:n_series) {\n',
'r_vec[1:n,s] = rep_vector(r[s], n);\n}\n')
'phi_vec[1:n,s] = rep_vector(phi[s], n);\n}\n')

to_negbin <- gsub('poisson_log_rng', 'neg_binomial_2_rng',
stan_file[grep('ypred[i, s] = poisson_log_rng', stan_file, fixed = T)])
stan_file[grep('ypred[i, s] = poisson_log_rng', stan_file, fixed = T)] <-
gsub(');', ', r_vec[i, s]);', to_negbin)
gsub(');', ', phi_vec[i, s]);', to_negbin)

add_exp_open <- gsub('\\(eta', '(exp(eta',
stan_file[grep('ypred[i, s] = neg_binomial', stan_file, fixed = T)])
Expand Down
Loading

0 comments on commit 7b8472c

Please sign in to comment.