Merge branch 'master' of https://github.com/SPATIAL-Lab/isoWater

.github/workflows/r.yml

@@ -45,10 +45,7 @@ jobs:
         run: |
           brew install jags
       - name: Install dependencies
-        run: |
-          install.packages(c("remotes", "rcmdcheck"))
-          remotes::install_deps(dependencies = TRUE)
-        shell: Rscript {0}
+        uses: r-lib/actions/setup-r-dependencies@v2
       - name: Check
         run: rcmdcheck::rcmdcheck(args = "--no-manual", error_on = "error")
         shell: Rscript {0}

.github/workflows/test-coverage.yaml

@@ -25,27 +25,9 @@ jobs:
           sudo apt-get install libcurl4-openssl-dev
           sudo apt-get install jags
           
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Restore R package cache
-        uses: actions/cache@v2
-        with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
       - name: Install dependencies
-        run: |
-          install.packages(c("remotes"))
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("covr")
-        shell: Rscript {0}
-
+        uses: r-lib/actions/setup-r-dependencies@v2
+
       - name: Test coverage
         run: covr::codecov()
         shell: Rscript {0}

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: isoWater
 Type: Package
 Title: Discovery, Retrieval, and Analysis of Water Isotope Data 
-Version: 1.1.2
+Version: 1.1.2.9000
 Authors@R: person("Gabe", "Bowen", email = "gabe.bowen@utah.edu",
                   role = c("aut", "cre"))
 Maintainer: Gabe Bowen <gabe.bowen@utah.edu>

NEWS.md

@@ -1,5 +1,8 @@
 # isoWater news
 
+## isoWater 1.1.2.9000
+* Added gamma distribution as an option for the prior on the evaporation effect in both mixSource() and mwlSource().
+
 ## isoWater 1.1.2
 * Bug fixes
 * mixSource function now allows optional prior on O-isotope evaporation effect

R/watercomp.R

@@ -4,7 +4,8 @@
 
 #takes values of observed water (type 'iso'), MWL (see below), hypothesized EL slope value
 #and number of parameter draws
-mwlSource = function(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1){
+mwlSource = function(obs, MWL = NULL, slope, stype = 1, edist = "unif", 
+                     eprior = NULL, ngens=1e4, ncores = 1){
 
   if(is.null(MWL)){
     data("GMWL", envir = environment())
@@ -59,10 +60,13 @@ mwlSource = function(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1){
   #Find se range for d18O
   o_keep = o_eval[h_el >= h_MWL.min & h_el <= h_MWL.max]
   o_var = diff(range(o_keep))^2
+
+  #check evap arguments
+  eprior = check_evap(edist, eprior)
 
   d = list(o_cent = o_cent, o_var = o_var,
            obs = obs, obs.vcov = obs.vcov, MWL = MWL,
-           slope = slope, ndat = ndat)
+           slope = slope, eprior = eprior, ndat = ndat)
 
   p = c("source_d2H", "source_d18O", "S", "E")
 
@@ -73,15 +77,17 @@ mwlSource = function(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1){
   if(ncores > 1){
     tmf = tempfile(fileext = ".txt")
     tmfs = file(tmf, "w")
-    cat(mwlModel, file = tmfs)
+    switch(match(edist, c("unif", "gamma")), cat(mwlModel, file = tmfs),
+           cat(mwlModel.gam, file = tmfs))
     close(tmfs)
     post = jags.parallel(data = d, parameters.to.save = p, 
                          model.file = tmf, 
                          n.iter = n.iter, n.chains = ncores, 
                          n.burnin = n.burnin, n.thin = n.thin)    
   } else{
+    model = switch(match(edist, c("unif", "gamma")), mwlModel, mwlModel.gam)
     post = jags(data = d, parameters.to.save = p, 
-                model.file = textConnection(mwlModel), n.iter = n.iter, 
+                model.file = textConnection(model), n.iter = n.iter, 
                 n.burnin = n.burnin, n.thin = n.thin)    
   }
 
@@ -102,8 +108,9 @@ mwlSource = function(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1){
 
 #takes values of observed and hypothesized endmember source waters (each type 'iso'),hypothesized EL slope,
 #prior (as relative contribution of each source to mixture), and number of parameter draws
-mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)), 
-                   shp=1, eprior = c(0, 15), ngens=1e5, ncores = 1){
+mixSource = function(obs, sources, slope, prior = rep(1,nrow(sources)), 
+                   shp = 1, edist = "unif", eprior = NULL, ngens = 1e5, 
+                   ncores = 1){
 
   if(!inherits(obs, "iso")){
     warning("Expecting iso object for obs, this argument may be
@@ -145,24 +152,8 @@ mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)),
   #dirchlet priors
   alphas = prior/min(prior) * shp
 
-  #evap priors
-  if(inherits(eprior, "numeric")){
-    if(length(eprior) != 2){
-      stop("eprior must be length 2")
-    }
-    if(any(eprior < 0)){
-      stop("eprior values must be equal to or greater than zero")
-    }
-    if(any(eprior > 15)){
-      message("eprior values greater than 15 are very unlikely in most systems")
-    }
-    eprior = sort(eprior)
-    if(eprior[2] <= eprior[1]){
-      eprior[2] = eprior[1] + 1e-3
-    }
-  } else{
-    stop("eprior must be numeric")
-  }
+  #check evap arguments
+  eprior = check_evap(edist, eprior)
 
   #data
   d = list(nsource = nsource, ndat = ndat, 
@@ -181,15 +172,17 @@ mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)),
   if(ncores > 1){
     tmf = tempfile(fileext = ".txt")
     tmfs = file(tmf, "w")
-    cat(mixModel, file = tmfs)
+    switch(match(edist, c("unif", "gamma")), cat(mixModel, file = tmfs),
+           cat(mixModel.gam, file = tmfs))
     close(tmfs)
     post = jags.parallel(data = d, parameters.to.save = p, 
                          model.file = tmf, 
                          n.iter = n.iter, n.chains = ncores, 
                          n.burnin = n.burnin, n.thin = n.thin)
   } else{
+    model = switch(match(edist, c("unif", "gamma")), mixModel, mixModel.gam)
     post = jags(data = d, parameters.to.save = p, 
-                model.file = textConnection(mixModel), 
+                model.file = textConnection(model), 
                 n.iter = n.iter, n.burnin = n.burnin, n.thin = n.thin)
   }
 
@@ -210,3 +203,36 @@ mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)),
 
   return(wcout)
 }
+
+check_evap = function(edist, eprior){
+  #evap distribution
+  if(!(edist %in% c("unif", "gamma"))){
+    stop("edist must be unif or gamma")
+  }
+
+  #evap prior
+  if(is.null(eprior)){
+    eprior = switch(edist, unif = c(0, 15), gamma = c(1, 1))
+  }
+  if(inherits(eprior, "numeric")){
+    if(length(eprior) != 2){
+      stop("eprior must be length 2")
+    }
+    if(edist == "unif"){
+      if(any(eprior < 0)){
+        stop("eprior values must be equal to or greater than zero")
+      }
+      if(any(eprior > 15)){
+        message("eprior values greater than 15 are very unlikely in most systems")
+      }
+      eprior = sort(eprior)
+      if(eprior[2] <= eprior[1]){
+        eprior[2] = eprior[1] + 1e-3
+      }
+    }
+  } else{
+    stop("eprior must be numeric")
+  }
+
+  return(eprior)
+}

data-raw/prepData.R

@@ -12,7 +12,30 @@ mwlModel = "model{
   o_pred = source_d18O + E
   
   #evap prior
-  E ~ dunif(0, 15)
+  E ~ dunif(eprior[1], eprior[2])
+  
+  #EL slope prior
+  S ~ dnorm(slope[1], 1 / slope[2] ^ 2)
+  
+  #MWL source prior
+  source_d2H ~ dnorm(source_d18O * MWL[1] + MWL[2], 1 / sy ^ 2) 
+  sy = MWL[5] * sqrt(MWL[7] + (source_d18O - MWL[3])^2 / MWL[4])
+  source_d18O ~ dnorm(o_cent, 1 / o_var)
+}"
+
+mwlModel.gam = "model{
+  #Data model
+  for(i in 1:ndat){
+    obs[i, 1:2] ~ dmnorm.vcov(c(h_pred, o_pred), 
+                              obs.vcov[(1 + (i-1) * 2):(2 + (i-1) * 2),])
+  }
+  
+  #Process model
+  h_pred = source_d2H + E * S
+  o_pred = source_d18O + E
+  
+  #evap prior
+  E ~ dgamma(eprior[1], eprior[2])
   
   #EL slope prior
   S ~ dnorm(slope[1], 1 / slope[2] ^ 2)
@@ -52,7 +75,37 @@ mixModel = "model{
   }
 }"
 
-use_data(mwlModel, mixModel, internal = TRUE, overwrite = TRUE)
+mixModel.gam = "model{
+  #Data model
+  for(i in 1:ndat){
+    obs[i, 1:2] ~ dmnorm.vcov(c(h_pred, o_pred), 
+                              obs.vcov[(1 + (i-1) * 2):(2 + (i-1) * 2),])
+  }
+  
+  #Process model
+  h_pred = mixture_d2H + E * S
+  o_pred = mixture_d18O + E
+  
+  #evap prior
+  E ~ dgamma(eprior[1], eprior[2])
+  
+  #EL slope prior
+  S ~ dnorm(slope[1], 1 / slope[2] ^ 2)
+  
+  #mixture
+  mixture_d18O = sum(fracs * h_s[, 2])
+  mixture_d2H = sum(fracs * h_s[, 1])
+  fracs ~ ddirch(alphas)
+  
+  #sources prior
+  for(i in 1:nsource){
+    h_s[i, 1:2] ~ dmnorm.vcov(sources[i, 1:2], 
+                              sources.vcov[(1 + (i-1) * 2):(2 + (i-1) * 2),])
+  }
+}"
+
+use_data(mwlModel, mwlModel.gam, mixModel, mixModel.gam, 
+         internal = TRUE, overwrite = TRUE)
 
 lsl = wiDB_sites(countries = "US", types = "Lake_or_pond")
 psl = wiDB_sites(countries = "US", types = "Precipitation")

man/mixSource.Rd

@@ -11,8 +11,8 @@ Given isotopic compositions of two or more potential sources, generate a posteri
 }
 
 \usage{
-mixSource(obs, sources, slope, prior = rep(1,nrow(sources)), shp = 1, 
-            eprior = c(0, 15), ngens = 1e5, ncores = 1)
+mixSource(obs, sources, slope, prior=rep(1,nrow(sources)), 
+  shp = 1, edist = "unif", eprior = NULL, ngens = 1e5, ncores = 1)
 }
 
 \arguments{
@@ -31,8 +31,11 @@ mixSource(obs, sources, slope, prior = rep(1,nrow(sources)), shp = 1,
   \item{shp}{
   numeric. Shape parameter constant used in specifying prior estimates of source contributions (see Details). 
   }
+  \item{edist}{
+  character. One of \code{"unif"} or \code{"gamma"}, specifying whether the evaporation prior is modeled using a uniform or gamma distribution. 
+  }
   \item{eprior}{
-  numeric. Vector of length 2 giving prior estimates of maximum and minimum oxygen isotope evaporation effect.
+  numeric. Vector of length 2 giving prior parameter estimates for the oxygen isotope evaporation effect. For \code{edist = "unif"} these are maximum and minimum values and the default values are \code{c(0, 15)}. For \code{edist = "gamma"} these are shape and rate parameters, and the defaults are \code{c(1, 1)}.
   }
   \item{ngens}{
   integer. Number of posterior samples to obtain (per chain).

man/mwlSource.Rd

@@ -11,7 +11,8 @@ Given parameters describing a meteoric water line in H-O isotope space, generate
 }
 
 \usage{
-mwlSource(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1)
+mwlSource(obs, MWL = NULL, slope, stype = 1, edist = "unif", 
+  eprior = NULL, ngens=1e4, ncores = 1)
 }
 
 \arguments{
@@ -27,6 +28,12 @@ mwlSource(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1)
   \item{stype}{
   integer. Line statistic used to constrain the source prior: 1 = confidence interval, 2 = prediction interval (see Details).
   }
+  \item{edist}{
+  character. One of \code{"unif"} or \code{"gamma"}, specifying whether the evaporation prior is modeled using a uniform or gamma distribution. 
+  }
+  \item{eprior}{
+  numeric. Vector of length 2 giving prior parameter estimates for the oxygen isotope evaporation effect. For \code{edist = "unif"} these are maximum and minimum values and the default values are \code{c(0, 15)}. For \code{edist = "gamma"} these are shape and rate parameters, and the defaults are \code{c(1, 1)}.
+  }
   \item{ngens}{
   integer. Number of posterior samples to obtain (per chain).
   }