update portal example

NEON_manuscript/portal_example.R

@@ -57,7 +57,7 @@ test <- mvjagam(formula =  y ~ s(season, bs = "cc", k = 12) +
                 data_test = data_test,
                 family = 'poisson',
                 chains = 4,
-                burnin = 5000,
+                burnin = 2000,
                 trend_model = 'AR1')
 
 plot_mvgam_fc(test, series = 1, data_test = data_test, ylim = c(0, 100))
@@ -69,6 +69,58 @@ plot_mvgam_smooth(test, series = 1, smooth = 3)
 plot_mvgam_trace(test, 'trend')
 summary_mvgam(test)
 
+# For visualising how covariate functions change with different lags, use
+# the predict_mvgam function
+# Set up prediction data by zeroing out all covariates apart from the covariate of
+# interest
+newdata <- data_test
+newdata$year <- rep(0, length(newdata$year))
+newdata$season <- rep(0, length(newdata$season))
+newdata$precip <- matrix(0, ncol = ncol(newdata$precip),
+                         nrow = nrow(newdata$precip))
+
+# Set up plot colours and initiate plot window
+cols <- viridis::inferno(5)
+plot(1, type = "n",
+     xlab = 'Mintemp',
+     ylab = 'Predicted response function',
+     xlim = c(min(data_train$mintemp), max(data_train$mintemp)),
+     ylim = c(-0.6, 0.6))
+
+for(i in 1:5){
+  # Set up prediction matrix for mintemp with lag i as the prediction sequence
+  newdata$mintemp <- matrix(0, ncol = ncol(newdata$precip),
+                            nrow = nrow(newdata$precip))
+  newdata$mintemp[,i] <- seq(min(data_train$mintemp),
+                             max(data_train$mintemp),
+                             length.out = length(newdata$year))
+
+  # Predict on the link scale (intercept still included in prediction by default)
+  preds <- predict_mvgam(test, series = 1, newdata = newdata, type = 'link')
+
+  # Calculate prediction quantiles
+  probs = c(0.05, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.95)
+  cred <- sapply(1:NCOL(preds),
+                 function(n) quantile(preds[,n],
+                                      probs = probs))
+
+  # Plot expected function posterior median in varying colours per lag; subtract out the
+  # intercept so that the function is roughly centred on zero
+  pred_upper <- cred[4,] - test$mgcv_model$coefficients[1]
+  pred_lower <- cred[6,] - test$mgcv_model$coefficients[1]
+  pred_vals <- seq(min(data_train$mintemp),
+                   max(data_train$mintemp),
+                   length.out = length(newdata$year))
+  polygon(c(pred_vals, rev(pred_vals)), c(pred_upper, rev(pred_lower)),
+          col = scales::alpha(cols[i], 0.6), border = scales::alpha(cols[i], 0.7))
+  lines(pred_vals, cred[5,] - test$mgcv_model$coefficients[1],
+        col = scales::alpha(cols[i], 0.8), lwd = 2.5)
+}
+abline(h = 0, lty = 'dashed')
+legend('topleft',legend=paste0('lag', seq(0, 4)),
+       bg = 'white', bty = 'n',
+       col=cols,lty=1,lwd=6)
+
 
 # Initiate particle filter
 pfilter_mvgam_init(test, data_assim = data_test)

R/mvjagam.R

@@ -62,7 +62,7 @@
 #'@param ar_prior \code{character} specifying (in JAGS syntax) the prior distribution for the AR terms
 #'in the latent trends
 #'@param r_prior \code{character} specifying (in JAGS syntax) the prior distribution for the Negative Binomial
-#'overdispersion parameters. Note that this prior acts on the INVERSE of \code{r}, which is convenient
+#'overdispersion parameters. Note that this prior acts on the SQUARE ROOT of \code{r}, which is convenient
 #'for inducing a complexity-penalising prior model whereby the observation process reduces to a Poisson
 #'as the sampled parameter approaches \code{0}. Ignored if \code{family = 'poisson'}
 #'@param tau_prior \code{character} specifying (in JAGS syntax) the prior distributions for the independent gaussian
@@ -337,8 +337,8 @@ for (i in 1:n) {
 ## 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] ~ dexp(0.5)
+ r[s] <- pow(r_raw[s], 2)
+ r_raw[s] ~ dexp(0.05)
 }
 
 
@@ -508,7 +508,7 @@ for (i in 1:n) {
          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] <- theta.prime[t] * l.var[t]
+         penalty[t] <- max(min_eps, theta.prime[t] * l.var[t])
         }
 
         ## Latent factor loadings: standard normal with identifiability constraints
@@ -558,8 +558,8 @@ for (i in 1:n) {
         ## 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] ~ dexp(0.5)
+         r[s] <- pow(r_raw[s], 2)
+         r_raw[s] ~ dexp(0.05)
         }
 
         ## Posterior predictions
@@ -770,6 +770,10 @@ for (i in 1:n) {
         ss_jagam$jags.data$n_lv <- min(2, floor(ss_jagam$jags.data$n_series / 2))
       } else {
         ss_jagam$jags.data$n_lv <- n_lv
+        ss_jagam$jags.ini$LV <- matrix(rep(0.1, n_lv), nrow = NROW(ytimes),
+                                       ncol = n_lv)
+        ss_jagam$jags.ini$X1 <- rep(1, n_lv)
+        ss_jagam$jags.ini$X2 <- 1
       }
       if(ss_jagam$jags.data$n_lv > ss_jagam$jags.data$n_series){
         stop('Number of latent variables cannot be greater than number of series')

R/plot_mvgam_uncertainty.R

@@ -59,7 +59,7 @@ plot_mvgam_uncertainty = function(object, series, data_test, legend_position = '
   }
 
   n_samples <- NROW(trend)
-  size <- MCMCvis::MCMCsummary(object$jags_output, 'r')[,series]$mean
+  size <- MCMCvis::MCMCsummary(object$jags_output, 'r')[series,]$mean
 
   # Full uncertainty interval for the mean
   if(class(data_test) == 'list'){

R/predict_mvgam.R

@@ -49,26 +49,26 @@ predict_mvgam = function(object, series = 1, newdata, type = 'link'){
 
     if(object$use_lv){
       lv_preds <- do.call(cbind, lapply(seq_len(object$n_lv), function(lv){
-        rnorm(NROW(newdata), 0, sqrt(1 / taus[x,lv]))
+        rnorm(length(newdata$series), 0, sqrt(1 / taus[x,lv]))
       }))
 
       if(type == 'link'){
         as.vector(((Xp[which(as.numeric(newdata$series) == series),] %*% betas[x,])) +
                            ( lv_preds %*% lv_coefs[[series]][x,]))
       } else {
-        rnbinom(n = NROW(newdata), size = sizes[x, series],
+        rnbinom(n = length(newdata$series), size = sizes[x, series],
                 mu = exp(((Xp[which(as.numeric(newdata$series) == series),] %*% betas[x,])) +
                            ( lv_preds %*% lv_coefs[[series]][x,])))
       }
 
     } else {
       if(type == 'link'){
         as.vector(((Xp[which(as.numeric(newdata$series) == series),] %*% betas[x,])) +
-                           (rnorm(NROW(newdata), 0, sqrt(1 / taus[x,series]))))
+                           (rnorm(length(newdata$series), 0, sqrt(1 / taus[x,series]))))
       } else {
-        rnbinom(n = NROW(newdata), size = sizes[x, series],
+        rnbinom(n = length(newdata$series), size = sizes[x, series],
                 mu = exp(((Xp[which(as.numeric(newdata$series) == series),] %*% betas[x,])) +
-                           (rnorm(NROW(newdata), 0, sqrt(1 / taus[x,series])))))
+                           (rnorm(length(newdata$series), 0, sqrt(1 / taus[x,series])))))
       }
 
     }

mvgam.Rproj

@@ -2,7 +2,7 @@ Version: 1.0
 
 RestoreWorkspace: No
 SaveWorkspace: No
-AlwaysSaveHistory: Default
+AlwaysSaveHistory: No
 
 EnableCodeIndexing: Yes
 UseSpacesForTab: Yes