fixing some typos in stan code

NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt

@@ -1,157 +1,157 @@
-JAGS model code generated by package mvgam
-
-GAM formula:
-y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs)
-
-Trend model:
-RW
-
-Required data:
-integer n;  number of timepoints per series
-integer n_series;  number of series
-matrix y;  time-ordered observations of dimension n x n_series (missing values allowed)
-matrix ytimes;  time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?)
-matrix X;  mgcv GAM design matrix of dimension (n x n_series) x basis dimension
-matrix S2;  mgcv smooth penalty matrix S2
-matrix S4;  mgcv smooth penalty matrix S4
-vector zero;  prior basis coefficient locations vector of length ncol(X)
-vector p_coefs;  vector (length = 1) of prior Gaussian means for parametric effects
-vector p_taus;  vector (length = 1) of prior Gaussian precisions for parametric effects
-
-
-#### Begin model ####
-model {
-
-## GAM linear predictor
-eta <- X %*% b
-
-## mean expectations
-for (i in 1:n) {
-for (s in 1:n_series) {
-mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s])
-}
-}
-
-## latent factors evolve as time series with penalised precisions;
-## the penalty terms force any un-needed factors to evolve as flat lines
-for (j in 1:n_lv) {
-LV[1, j] ~ dnorm(0, penalty[j])
-}
-
-for (i in 2:n) {
-for (j in 1:n_lv){
-LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j])
-}
-}
-
-## shrinkage penalties for each factor's precision act to squeeze
-## the entire factor toward a flat white noise process if supported by
-## the data. The prior for individual factor penalties allows each factor to possibly
-## have a relatively large penalty, which shrinks the prior for that factor's variance
-## substantially. Penalties increase exponentially with the number of factors following
-## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate
-## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291.
-pi ~ dunif(0, n_lv)
-X2 ~ dnorm(0, 1)T(0, )
-
-# eta1 controls the baseline penalty
-eta1 ~ dunif(-1, 1)
-
-# eta2 controls how quickly the penalties exponentially increase
-eta2 ~ dunif(-1, 1)
-
-for (t in 1:n_lv) {
-X1[t] ~ dnorm(0, 1)T(0, )
-l.dist[t] <- max(t, pi[])
-l.weight[t] <- exp(eta2[] * l.dist[t])
-l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1
-theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[]
-penalty[t] <- max(0.0001, theta.prime[t] * l.var[t])
-}
-
-## latent factor loadings: standard normal with identifiability constraints
-## upper triangle of loading matrix set to zero
-for (j in 1:(n_lv - 1)) {
-for (j2 in (j + 1):n_lv) {
-lv_coefs[j, j2] <- 0
-}
-}
-
-## positive constraints on loading diagonals
-for (j in 1:n_lv) {
-lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1);
-}
-
-## lower diagonal free
-for (j in 2:n_lv) {
-for (j2 in 1:(j - 1)) {
-lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1);
-}
-}
-
-## other elements also free
-for (j in (n_lv + 1):n_series) {
-for (j2 in 1:n_lv) {
-lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1);
-}
-}
-
-## trend evolution depends on latent factors
-for (i in 1:n) {
-for (s in 1:n_series) {
-trend[i, s] <- inprod(lv_coefs[s,], LV[i,])
-}
-}
-
-## likelihood functions
-for (i in 1:n) {
-for (s in 1:n_series) {
-y[i, s] ~ dnegbin(rate[i, s], r[s])
-rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps,
-(r[s] / (r[s] + mu[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_raw[s]
-r_raw[s] ~ dnorm(0, 0.1)T(0, )
-}
-
-## posterior predictions
-for (i in 1:n) {
-for (s in 1:n_series) {
-ypred[i, s] ~ dnegbin(rate[i, s], r[s])
-}
-}
-
-## GAM-specific priors
-## parametric effect priors (regularised for identifiability)
-for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) }
-## prior (non-centred) for s(siteID)...
-for (i in 2:4) {
-b_raw[i] ~ dnorm(0, 1)
-b[i] <- mu_raw1 + b_raw[i] * sigma_raw1
-}
-sigma_raw1 ~ dexp(1)
-mu_raw1 ~ dnorm(0, 1)
-## prior for s(cum_gdd)...
-K2 <- S2[1:4,1:4] * lambda[2]  + S2[1:4,5:8] * lambda[3]
-b[5:8] ~ dmnorm(zero[5:8],K2)
-## prior for s(cum_gdd,siteID)...
-for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) }
-for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) }
-## prior for s(season)...
-K4 <- S4[1:10,1:10] * lambda[6]
-b[24:33] ~ dmnorm(zero[24:33],K4)
-## prior for s(season,series)...
-for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) }
-for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) }
-## smoothing parameter priors...
-for (i in 1:8) {
-lambda[i] ~ dexp(0.05)T(0.0001, )
-rho[i] <- log(lambda[i])
-}
-}
+JAGS model code generated by package mvgam
+
+GAM formula:
+y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs)
+
+Trend model:
+RW
+
+Required data:
+integer n;  number of timepoints per series
+integer n_series;  number of series
+matrix y;  time-ordered observations of dimension n x n_series (missing values allowed)
+matrix ytimes;  time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?)
+matrix X;  mgcv GAM design matrix of dimension (n x n_series) x basis dimension
+matrix S2;  mgcv smooth penalty matrix S2
+matrix S4;  mgcv smooth penalty matrix S4
+vector zero;  prior basis coefficient locations vector of length ncol(X)
+vector p_coefs;  vector (length = 1) of prior Gaussian means for parametric effects
+vector p_taus;  vector (length = 1) of prior Gaussian precisions for parametric effects
+
+
+#### Begin model ####
+model {
+
+## GAM linear predictor
+eta <- X %*% b
+
+## mean expectations
+for (i in 1:n) {
+for (s in 1:n_series) {
+mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s])
+}
+}
+
+## latent factors evolve as time series with penalised precisions;
+## the penalty terms force any un-needed factors to evolve as flat lines
+for (j in 1:n_lv) {
+LV[1, j] ~ dnorm(0, penalty[j])
+}
+
+for (i in 2:n) {
+for (j in 1:n_lv){
+LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j])
+}
+}
+
+## shrinkage penalties for each factor's precision act to squeeze
+## the entire factor toward a flat white noise process if supported by
+## the data. The prior for individual factor penalties allows each factor to possibly
+## have a relatively large penalty, which shrinks the prior for that factor's variance
+## substantially. Penalties increase exponentially with the number of factors following
+## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate
+## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291.
+pi ~ dunif(0, n_lv)
+X2 ~ dnorm(0, 1)T(0, )
+
+# eta1 controls the baseline penalty
+eta1 ~ dunif(-1, 1)
+
+# eta2 controls how quickly the penalties exponentially increase
+eta2 ~ dunif(-1, 1)
+
+for (t in 1:n_lv) {
+X1[t] ~ dnorm(0, 1)T(0, )
+l.dist[t] <- max(t, pi[])
+l.weight[t] <- exp(eta2[] * l.dist[t])
+l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1
+theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[]
+penalty[t] <- max(0.0001, theta.prime[t] * l.var[t])
+}
+
+## latent factor loadings: standard normal with identifiability constraints
+## upper triangle of loading matrix set to zero
+for (j in 1:(n_lv - 1)) {
+for (j2 in (j + 1):n_lv) {
+lv_coefs[j, j2] <- 0
+}
+}
+
+## positive constraints on loading diagonals
+for (j in 1:n_lv) {
+lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1);
+}
+
+## lower diagonal free
+for (j in 2:n_lv) {
+for (j2 in 1:(j - 1)) {
+lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1);
+}
+}
+
+## other elements also free
+for (j in (n_lv + 1):n_series) {
+for (j2 in 1:n_lv) {
+lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1);
+}
+}
+
+## trend evolution depends on latent factors
+for (i in 1:n) {
+for (s in 1:n_series) {
+trend[i, s] <- inprod(lv_coefs[s,], LV[i,])
+}
+}
+
+## likelihood functions
+for (i in 1:n) {
+for (s in 1:n_series) {
+y[i, s] ~ dnegbin(rate[i, s], r[s])
+rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps,
+(r[s] / (r[s] + mu[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_raw[s]
+r_raw[s] ~ dnorm(0, 0.1)T(0, )
+}
+
+## posterior predictions
+for (i in 1:n) {
+for (s in 1:n_series) {
+ypred[i, s] ~ dnegbin(rate[i, s], r[s])
+}
+}
+
+## GAM-specific priors
+## parametric effect priors (regularised for identifiability)
+for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) }
+## prior (non-centred) for s(siteID)...
+for (i in 2:4) {
+b_raw[i] ~ dnorm(0, 1)
+b[i] <- mu_raw1 + b_raw[i] * sigma_raw1
+}
+sigma_raw1 ~ dexp(1)
+mu_raw1 ~ dnorm(0, 1)
+## prior for s(cum_gdd)...
+K2 <- S2[1:4,1:4] * lambda[2]  + S2[1:4,5:8] * lambda[3]
+b[5:8] ~ dmnorm(zero[5:8],K2)
+## prior for s(cum_gdd,siteID)...
+for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) }
+for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) }
+## prior for s(season)...
+K4 <- S4[1:10,1:10] * lambda[6]
+b[24:33] ~ dmnorm(zero[24:33],K4)
+## prior for s(season,series)...
+for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) }
+for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) }
+## smoothing parameter priors...
+for (i in 1:8) {
+lambda[i] ~ dexp(0.05)T(0.0001, )
+rho[i] <- log(lambda[i])
+}
+}