Empirical-Bayes-Resampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9d975f1f",
   "metadata": {},
   "source": [
    "# Empirical Bayes uncertainty estimation for (U–Th)/He data\n",
    "\n",
    "This Jupyter notebook demonstrates the use of a form of empirical Bayes resampling (also known as hierarchical Bayes) for more robust (U-Th)/He data uncertainty estimation for use in thermal-history inversions. See [Malinverno and Briggs (2004)](https://doi.org/10.1190/1.1778243) for description of the hierarchical Bayesian method. As described and implemented here, it is particularly suited for zircon datasets (high n) that span a broad range of effective uranium (eU). The aim is to expand uncertainty accounting where the hyperparameters (i.e., observed He dates) will have a prior distribution that expresses their initial uncertainty and a posterior distribution that is determined by the data directly. The individual date errors are treated as hyperparameters drawn from a probability distribution and the variance is used to infer the 'empirical' date uncertainty. This weighted 1σ uncertainty is inferred from the scatter of the data as determined by the standard deviation of the data weighted by a Gaussian kernel in eU space (σeU = 100 ppm). The 100 ppm eU kernel is taken to represent the range over which zircon grains with similar eU should have similar ages. This value was chosen as a good balance because smaller values tend to converge on the internal uncertainty and larger values begin to reduce the influence of individual dates during thermal modelling. The empirical uncertainty is estimated by summing the internal and external uncertainties in quadrature.\n",
    "\n",
    "Helium data are typically overdispersed with respect to analytical uncertainties (approx. 2-5% for apatite and zircon). The age reproducibility of helium age standards such as the Durango (DUR) apatite and Fish Canyon Tuff (FCT) zircon suggest uncertainties are greater, on the order of 6-7% for apatite and 8-10% for zircon. For more information regarding date reproducibility and ICPMS measurements, see Reiners and Nicolescu (2006) (https://www.geo.arizona.edu/~reiners/arhdl/arhdlrep1.pdf), [Guenthner et al. (2016)](https://doi.org/10.1002/2016GC006311) for U–Th/He, and [Gleadow et al. (2015)](https://doi.org/10.1016/j.epsl.2015.05.003) for FCT specifically. \n",
    "\n",
    "These are still likely minimum estimates, since those grains are laboratory standards with well-behaved diffusion behavior demonstrated experimentally. Natural, 'wild' samples typically have total uncertainties that exceed 10% due various reasons such as (unaccounted for) isotopic zoning, imperfect grain measurement/grain morphology characterization for alpha-loss correction, or unidentified mineral or fluid inclusions (to name a few). There is also the possibility for trapping of pre- and post-'closure' 4He in various grain sinks and imperfections that further complicate the age reproducibility of slowly cooled apatite; see [Zeitler et al. (2017)](https://doi.org/10.1016/j.gca.2017.03.041), [McDannell et al. (2018)](https://doi.org/10.1016/j.gca.2017.11.031), and [Guo et al. (2021)](https://doi.org/10.1016/j.gca.2021.07.015) for recent experiments and discussion of He age dispersion and diffusion systematics in apatite.\n",
    "\n",
    "Note: some publications bin single-grain dates by eU and then average the dates in each bin to create 'synthetic' He dates, but this approach lacks rigor and quickly becomes an arbitrary exercise depending on data quality and quantity. This is even more problematic for apatite dates that normally span a narrow eU range. This is usually carried out for time-temperature (t-T) modelling, however, it removes temporal information and artificially increases date errors, and as a result, makes it easier to meet p-value statistical thresholds of acceptance for pure Monte Carlo modelling. Thus its only purpose is to 'hack' the t-T search and allow more paths to be found more easily. Do not do this. Since most age dispersion is explainable to first-order by the effects of radiation damage on diffusivity, each grain is a separate thermochronometer (that is the advantage of using radiation-damage models in the first place!). Averaging dates results in a loss of valuable t-T information and decreased t-T resolution. Importantly, the empirical Bayes approach retains the observed date calculated from the measured isotopic data but the uncertainties are estimated from the data directly. Error resampling is helpful during t-T modelling when errors are underestimated and prevents the t–T search algorithm from becoming trapped in local minima by reducing overprecision (in the case of learning algorithms such as Markov chain Monte Carlo). A notable point is that the empirical uncertainties only act to assist in improving the search, as the observed age is still the target datum, therefore, it is different from a simple Monte Carlo approach where increasing the uncertainties (or reducing the number of precise data by averaging) actually results in more paths being found and accepted during the random sifting of t-T space."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b25cd7ad",
   "metadata": {},
   "source": [
    "## Load required Julia packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f09b2f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## --- Load (and install if neccesary) the StatGeochem package\n",
    "\n",
    "using Plots, Statistics, StatGeochem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a25d7011",
   "metadata": {},
   "source": [
    "## Enter sample information\n",
    "\n",
    "First we read in our sample information for completeness and record-keeping purposes for output. This is of course optional, and we could instead just read in ages, actinide concentrations, and uncertainties for simplicity.\n",
    "\n",
    "sample name; equivalent spherical radius (ESR) or grain halfwidth; U-Th-(Sm) are used for calculating eU;\n",
    "and the observed uncorrected He age and its corresponding 1σ uncertainty. A row can be added for Sm if required and the eU equation will have to be updated accordingly. See [Cooperdock et al. (2019)](https://gchron.copernicus.org/articles/1/17/2019/), for the most recent and complete eU equation(s) and derivation: eU = U + 0.238 x Th + 0.0012 x Sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0788464",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ds = Dict{String, Union{Vector{Float64}, Vector{String}}}()\n",
    "ds[\"Name\"] = [\"04RF1zA\",\"04SC1zB\",\"04SC1zA\",\"04GF1zB\",\"04GF1zA\",\"04RF1zB\",\"04MT1zA\",\"04R1zB\",\"04MT1zB\",\"04R1zA\",\"04MT1zF\",\"04EQ1zA\",\"04EQ1zF\",\"04SH1zB\",\"04EQ1zB\",\"04MT1zD\",\"04MT1zC\",\"04MT1zE\",\"04SG1zB\",\"04EQ1zC\",\"04EQ1zE\",\"04EQ1zD\",\"04SH1zA\",]\n",
    "ds[\"ESR\"] = ESR = [64,71,47,63,56,72,67,51,68,42,54,41,43,56,42,80,45,49,55,36,35,31,50.,]\n",
    "ds[\"U [ppm]\"] = U = [84,415,422,254,247,155,362,580,346,656,567,781,647,883,748,739,980,974,1738,894,1107,866,1011.,]\n",
    "ds[\"Th [ppm]\"] = Th = [67.3,111,110,104,88.4,177,88.8,174,78,167,75.48,203,284,184,218,140.59,251.7,197.01,1171,246,351,166,228,]\n",
    "ds[\"Age [Ma]\"] = Age = [770,659,649,638.4,619.8,557,545,500,493,357,329,253,240.8,224.6,224.5,217,193,190,72,57,42,29.5,11.4,]\n",
    "ds[\"Age_sigma [Ma]\"] = Age_sigma = [15.0,51.0,13.0,12.6,18.4,10.0,12.0,10.0,43.0,7.0,7.0,5.2,4.6,4.6,6.0,5.0,4.0,4.0,2.0,1.1,1.1,0.8,0.2,]\n",
    "ds[\"eU [ppm]\"] = eU = U+.238*Th"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b26684b",
   "metadata": {},
   "source": [
    "#### optional: set minimum 1σ error to 5% for zircon date, if not already so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7873dbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#t = Age_sigma .< (0.05.*Age)\n",
    "#Age_sigma[t] = (0.05.*Age)[t]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0837f02",
   "metadata": {},
   "source": [
    "## Calculate the empirical standard deviation on the single-grain dates and plot the results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffced5f0",
   "metadata": {},
   "source": [
    "Standard deviation of the Gaussian kernel in eU space representing the range over which grains with similar eU should have similar ages. Plotting x,y with both internal and empirical date uncertainties. We then export the results and our sample information to a .csv file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68f525e",
   "metadata": {},
   "outputs": [],
   "source": [
    "σeU = 100\n",
    "\n",
    "# Calculate errors\n",
    "ds[\"Age_sigma_empirical [Ma]\"] = Age_sigma_empirical = zeros(length(Age))\n",
    "for i = 1:length(Age)\n",
    "    W = normpdf.(eU[i], σeU, eU)\n",
    "    σ_external = nanstd(Age, W)\n",
    "    σ_internal = Age_sigma[i]\n",
    "    Age_sigma_empirical[i] = sqrt(σ_external^2 + σ_internal^2)\n",
    "end\n",
    "\n",
    "h = plot(xlabel=\"eU (ppm)\", ylabel=\"Date (Ma)\", framestyle=:box)\n",
    "plot!(eU, Age, yerror=Age_sigma_empirical, seriestype=:scatter, c=:black, msc=:black, label=\"empirical\")\n",
    "plot!(eU, Age, yerror=Age_sigma, seriestype=:scatter, c=:red, msc=:red, label=\"internal\")\n",
    "xq = range(minimum(eU), maximum(eU), length=300)\n",
    "savefig(h,\"Age_sigma_empirical.pdf\")\n",
    "display(h)\n",
    "\n",
    "println(\"Age_sigma_empirical:\")\n",
    "println(round.(Age_sigma_empirical, sigdigits=4))\n",
    "exportdataset(ds, \"sample-emp.csv\", ',')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}