first addition of code

.Rbuildignore

@@ -1,2 +1,3 @@
 ^mvgam\.Rproj$
 ^\.Rproj\.user$
+^LICENSE\.md$

DESCRIPTION

@@ -1,15 +1,21 @@
 Package: mvgam
-Title: What the Package Does (One Line, Title Case)
+Title: Multivariate Bayesian Generalised Additive Models for discrete time series
 Version: 0.0.0.9000
 Authors@R: 
-    person(given = "First",
-           family = "Last",
+    person(given = "Nicholas",
+           family = "Clark",
            role = c("aut", "cre"),
-           email = "[email protected]",
-           comment = c(ORCID = "YOUR-ORCID-ID"))
-Description: What the package does (one paragraph).
-License: `use_mit_license()`, `use_gpl3_license()` or friends to
-    pick a license
+           email = "[email protected]")
+Description: Multivariate Bayesian Generalised Additive Models for discrete time series
+Maintainer: Nicholas J Clark <[email protected]>
+License: MIT + file LICENSE
+Depends: R (>= 4.0.0), mgcv, rjags
+Imports: parallel,
+    coda,
+    MCMCpack,
+    runjags,
+    MCMCvis,
+    dplyr
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2021
+COPYRIGHT HOLDER: Nicholas Clark

LICENSE.md

@@ -0,0 +1,21 @@
+# MIT License
+
+Copyright (c) 2021 Nicholas Clark
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

NAMESPACE

@@ -1,2 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(hpd)
+export(mvjagam)
+export(plot_mvgam_fc)
+export(plot_mvgam_trend)

R/all_neon_tick_data.R

@@ -0,0 +1,15 @@
+#' NEON Amblyomma and Ixodes tick abundance survey data
+#'
+#' A dataset containing timeseries of Amblyomma americanum and Ixodes scapularis nymph abundances at NEON sites
+#'
+#' @format A tibble/dataframe containing covariate information alongside the main fields of:
+#' \describe{
+#' \item{Year}{Year of sampling}
+#' \item{epiWeek}{Epidemiological week of sampling}
+#'   \item{plot_ID}{NEON plot ID for survey location}
+#'   \item{siteID}{NEON site ID for survey location}
+#'   \item{amblyomma_americanum}{Counts of A. americanum nymphs}
+#'   \item{ixodes_scapularis}{Counts of I. scapularis nymphs}
+#' }
+#' @source \url{https://www.neonscience.org/data}
+"all_neon_tick_data"

R/hpd.R

@@ -0,0 +1,58 @@
+# Calculate the highest posterior density interval
+#'
+#'This function uses estimated densities to calculate HPD intevals. Code originally supplied by Martyn Plummer
+#'
+#'@param x \code{vector} of values representing the distribution to be summarised
+#'@param coverage \code{numeric} value specifying the width of the HPD interval. Default is 0.95
+#'@return A \code{list} containing the reconciled forecast distributions for each series in \code{y}. Each element in
+#'the \code{vector} with three values: lower estimate, median estimate and upper estimate of the HPD interval
+#'
+#'@export
+#'
+hpd <- function(x, coverage = 0.95)
+{
+  x <- as.matrix(x)
+  out <- matrix(NA, nrow = ncol(x), ncol = 3)
+  rownames(out) <- dimnames(x)[[2]]
+  colnames(out) <- c("mode", "lower", "upper")
+
+  f <- function(p) {
+    if (p == density.range[2]) {
+      set.coverage <- 0
+    }
+    else {
+      p.upper <- min(y.density$y[y.density$y > p])
+      p.lower <- max(y.density$y[y.density$y <= p])
+      cov.upper <- sum(y.counts[y.density$y >= p.upper]) / sum(y.counts)
+      cov.lower <- sum(y.counts[y.density$y >= p.lower]) / sum(y.counts)
+      c <- (p.upper - p)/(p.upper - p.lower)
+      set.coverage <- c * cov.upper + (1 - c) * cov.lower
+    }
+    return(set.coverage - coverage)
+  }
+
+  for (i in 1:ncol(x)) {
+    y <- unclass(x[,i])
+    y.density <- stats::density(y, n=1024)
+    m <- length(y.density$x)
+
+    ## Find the midpoint
+    out[i,1] <- y.density$x[which.max(y.density$y)]
+    dx <- diff(range(y.density$x)) / m
+    breaks <- c(y.density$x[1] - dx / 2, y.density$x + dx /2)
+    y.counts <- hist(y, breaks = breaks, plot = FALSE)$counts
+    density.range <- range(y.density$y)
+    uniroot.out <- stats::uniroot(f, density.range)
+
+    ## Assuming that we have a single interval, find the limits
+    out[i,2:3] <- range(y.density$x[y.density$y > uniroot.out$root])
+
+    ## Check!
+    if (sum(abs(diff(y.density$y > uniroot.out$root))) != 2) {
+      warning("HPD set is not a closed interval")
+    }
+  }
+  out <- c(out[2], out[1], out[3])
+  names(out) <- c('lower', 'median', 'upper')
+  return(out)
+}