Work on #11

- Add simple code to see how process rates change as a function of BD

docs/introduction.html

@@ -261,7 +261,7 @@ <h3 data-number="2.2.1" class="anchored" data-anchor-id="what-do-the-various-fun
  }
  return(results)
 }
-&lt;bytecode: 0x7f986a1ac640&gt;
+&lt;bytecode: 0x7fee7dc14240&gt;
 &lt;environment: namespace:PRE&gt;</code></pre>
 </div>
 </div>

docs/search.json

@@ -18,7 +18,7 @@
  "href": "introduction.html#installation-of-the-r-package-pre",
  "title": "2  Introduction",
  "section": "2.2 Installation of the R package PRE",
- "text": "2.2 Installation of the R package PRE\nAll functions that are used repeatedly throughout this projec are implemented in an R package called PRE. The source code of this package is stored in a publically accessible way on GitHub. You can install the package directly from GitHub using the install_github function from the remotes package.\n\ninstall.packages(\"remotes\")\nremotes::install_github(\"https://github.com/Damian-Oswald/PRE\")\n\nAlternatively, you could download the GitHub repository and manually install the PRE package locally. Once the package is installed, you can attach it to the search path in the usual way.\n\nlibrary(PRE)\n\n\n2.2.1 What do the various functions do?\nThe R code shown on this project help page is extremely simplified. If you want want to understand what a specific function of the PRE package does, simple visit the corresponding R help page. For example, you could read on the function crossValidate by using one of the the following commands.\n\nhelp(\"crossValidate\", package = \"PRE\")\n?crossValidate\n\nIf you want to jump really deep into a function, you could print it out in your console.\n\nprint(crossValidate)\n\nfunction (FUN, x, y, k = nrow(x), r = 1, hyperparameter = NULL) \n{\n results <- data.frame()\n for (R in 1:r) {\n I <- matrix(c(sample(1:nrow(x)), rep(NA, k - nrow(x)%%k)), \n ncol = k, byrow = TRUE)\n for (K in 1:k) {\n i <- na.omit(I[, K])\n result <- cbind(cost = FUN(x_train = x[-i, ], y_train = y[-i, \n ], x_test = x[i, ], y_test = y[i, ], hyperparameter = hyperparameter), \n r = R, k = K)\n results <- rbind(results, result)\n }\n }\n return(results)\n}\n<bytecode: 0x7f986a1ac640>\n<environment: namespace:PRE>\n\n\nNote this function source code output is lacking any comments. Those were automatically removed when compiling the R package. To see the commented source code, you’ll have to jump into the package files, or visit the GitHub page. Simply use the following URL:\nhttps://github.com/Damian-Oswald/PRE/blob/main/R/crossValidate.R\nNotice how the base name (the last part of the URL) matches the function name. You can replace that by any function name of your liking."
+ "text": "2.2 Installation of the R package PRE\nAll functions that are used repeatedly throughout this projec are implemented in an R package called PRE. The source code of this package is stored in a publically accessible way on GitHub. You can install the package directly from GitHub using the install_github function from the remotes package.\n\ninstall.packages(\"remotes\")\nremotes::install_github(\"https://github.com/Damian-Oswald/PRE\")\n\nAlternatively, you could download the GitHub repository and manually install the PRE package locally. Once the package is installed, you can attach it to the search path in the usual way.\n\nlibrary(PRE)\n\n\n2.2.1 What do the various functions do?\nThe R code shown on this project help page is extremely simplified. If you want want to understand what a specific function of the PRE package does, simple visit the corresponding R help page. For example, you could read on the function crossValidate by using one of the the following commands.\n\nhelp(\"crossValidate\", package = \"PRE\")\n?crossValidate\n\nIf you want to jump really deep into a function, you could print it out in your console.\n\nprint(crossValidate)\n\nfunction (FUN, x, y, k = nrow(x), r = 1, hyperparameter = NULL) \n{\n results <- data.frame()\n for (R in 1:r) {\n I <- matrix(c(sample(1:nrow(x)), rep(NA, k - nrow(x)%%k)), \n ncol = k, byrow = TRUE)\n for (K in 1:k) {\n i <- na.omit(I[, K])\n result <- cbind(cost = FUN(x_train = x[-i, ], y_train = y[-i, \n ], x_test = x[i, ], y_test = y[i, ], hyperparameter = hyperparameter), \n r = R, k = K)\n results <- rbind(results, result)\n }\n }\n return(results)\n}\n<bytecode: 0x7fee7dc14240>\n<environment: namespace:PRE>\n\n\nNote this function source code output is lacking any comments. Those were automatically removed when compiling the R package. To see the commented source code, you’ll have to jump into the package files, or visit the GitHub page. Simply use the following URL:\nhttps://github.com/Damian-Oswald/PRE/blob/main/R/crossValidate.R\nNotice how the base name (the last part of the URL) matches the function name. You can replace that by any function name of your liking."
  },
  {
  "objectID": "data-preparation.html#description-of-the-measurements",
@@ -109,7 +109,14 @@
  "href": "sensitivity-analysis.html#how-do-the-process-rates-correlate",
  "title": "6  Sensitivity Analysis",
  "section": "6.3 How do the process rates correlate?",
- "text": "6.3 How do the process rates correlate?\n\n\n\n\nFigure 6.2: Combinations of scatterplots of the estimated process rates for nitrification, denitrification and reduction. Each point represents one mean process rate for a given depth, column and set of randomly drawn parameters according to table 6.1. As to be seen, clusters are formed based on combinations of depth and column. The clusters themselves can each be well described as multivariate Gaussians (colored area). The scatterplots are drawn with constrained x- and y-axis limits (A) and unconstrained ones (B).\n\n\n\nTwo things are striking when looking at the clusters in figure 6.2: Both the within-group variance as well as the total variance (i.e. the variance of the cluster location) seems to be lower for the deeper depths. Additionally, one of the cluster centers constitutes a massive outlier – it’s column 12 at 7.5 cm depth (Figure 6.2 B).\n\n\n\n\nBates, D., M. Mächler, B. Bolker, and S. Walker. 2014. Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:1406.5823.\n\n\nDecock, C., and J. Six. 2013. On the potential of \\(\\delta\\) 18O and \\(\\delta\\) 15N to assess N2O reduction to N2 in soil. European journal of soil science 64:610–620.\n\n\nDenk, T. R., J. Mohn, C. Decock, D. Lewicka-Szczebak, E. Harris, K. Butterbach-Bahl, R. Kiese, and B. Wolf. 2017. The nitrogen cycle: A review of isotope effects and isotope modeling approaches. Soil Biology and Biochemistry 105:121–137.\n\n\nLewicka-Szczebak, D., J. Augustin, A. Giesemann, and R. Well. 2017. Quantifying N₂O reduction to N₂ based on N₂O isotopocules – validation with independent methods (helium incubation and ¹⁵N gas flux method). Biogeosciences 14:711–732.\n\n\nSobol’, I. M. 1990. On sensitivity estimation for nonlinear mathematical models. Matematicheskoe modelirovanie 2:112–118.\n\n\nWell, R., and H. Flessa. 2008. Isotope fractionation factors of N₂O diffusion. Rapid Communications in Mass Spectrometry: An International Journal Devoted to the Rapid Dissemination of Up-to-the-Minute Research in Mass Spectrometry 22:2621–2628.\n\n\nWickham, H. 2011. The split-apply-combine strategy for data analysis. Journal of statistical software 40:1–29."
+ "text": "6.3 How do the process rates correlate?\n\n\n\n\nFigure 6.2: Combinations of scatterplots of the estimated process rates for nitrification, denitrification and reduction. Each point represents one mean process rate for a given depth, column and set of randomly drawn parameters according to table 6.1. As to be seen, clusters are formed based on combinations of depth and column. The clusters themselves can each be well described as multivariate Gaussians (colored area). The scatterplots are drawn with constrained x- and y-axis limits (A) and unconstrained ones (B).\n\n\n\nTwo things are striking when looking at the clusters in figure 6.2: Both the within-group variance as well as the total variance (i.e. the variance of the cluster location) seems to be lower for the deeper depths. Additionally, one of the cluster centers constitutes a massive outlier – it’s column 12 at 7.5 cm depth (Figure 6.2 B)."
+ },
+ {
+ "objectID": "sensitivity-analysis.html#what-effect-does-the-bulk-density-have",
+ "href": "sensitivity-analysis.html#what-effect-does-the-bulk-density-have",
+ "title": "6  Sensitivity Analysis",
+ "section": "6.4 What effect does the bulk density have?",
+ "text": "6.4 What effect does the bulk density have?\nFirst, we’ll determine a range of bulk density to be searched.\nBD <- getParameters()$BD + seq(-0.15, 0.15, 0.05)\nNext, we have to construct a function to return the process estimates for a give bulk density x.\nf <- function(x) {\n parameters <- getParameters(BD = x)\n data <- subset(measurements, column==1) |>\n getN2ON(parameters = parameters) |>\n getMissing() |>\n calculateFluxes(parameters = parameters)\n return(longPRE(data, 1, 7.5, n = 15)$processes)\n}\nFinally, we’ll apply said function to all bulk densities in the vector BD.\nresult <- t(sapply(BD, f))\nThis process will return us the data to be analyzed.\n\n\n\nFigure 6.3: Illustration to the sensitivity of estimated process rates to a change in prior bulk density.\n\n\n\n\n\n\n\n\n\n\n\n\nDependent variable:\n\n\n\n\n\n\n\n\n\n\n\n\nNitrification\n\n\nDenitrification\n\n\nReduction\n\n\n\n\n\n\n(1)\n\n\n(2)\n\n\n(3)\n\n\n\n\n\n\n\n\nBD\n\n\n-20.907***\n\n\n-4.672***\n\n\n31.267***\n\n\n\n\n\n\n(0.993)\n\n\n(0.619)\n\n\n(0.926)\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nConstant\n\n\n47.404***\n\n\n26.665***\n\n\n-27.614***\n\n\n\n\n\n\n(1.678)\n\n\n(1.045)\n\n\n(1.565)\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nObservations\n\n\n7\n\n\n7\n\n\n7\n\n\n\n\nR2\n\n\n0.989\n\n\n0.919\n\n\n0.996\n\n\n\n\nAdjusted R2\n\n\n0.987\n\n\n0.903\n\n\n0.995\n\n\n\n\nResidual Std. Error (df = 5)\n\n\n0.263\n\n\n0.164\n\n\n0.245\n\n\n\n\nF Statistic (df = 1; 5)\n\n\n442.914***\n\n\n57.032***\n\n\n1,139.130***\n\n\n\n\n\n\n\n\nNote:\n\n\np<0.1; p<0.05; p<0.01\n\n\n\n\n\n\n\n\nBates, D., M. Mächler, B. Bolker, and S. Walker. 2014. Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:1406.5823.\n\n\nDecock, C., and J. Six. 2013. On the potential of \\(\\delta\\) 18O and \\(\\delta\\) 15N to assess N2O reduction to N2 in soil. European journal of soil science 64:610–620.\n\n\nDenk, T. R., J. Mohn, C. Decock, D. Lewicka-Szczebak, E. Harris, K. Butterbach-Bahl, R. Kiese, and B. Wolf. 2017. The nitrogen cycle: A review of isotope effects and isotope modeling approaches. Soil Biology and Biochemistry 105:121–137.\n\n\nLewicka-Szczebak, D., J. Augustin, A. Giesemann, and R. Well. 2017. Quantifying N₂O reduction to N₂ based on N₂O isotopocules – validation with independent methods (helium incubation and ¹⁵N gas flux method). Biogeosciences 14:711–732.\n\n\nSobol’, I. M. 1990. On sensitivity estimation for nonlinear mathematical models. Matematicheskoe modelirovanie 2:112–118.\n\n\nWell, R., and H. Flessa. 2008. Isotope fractionation factors of N₂O diffusion. Rapid Communications in Mass Spectrometry: An International Journal Devoted to the Rapid Dissemination of Up-to-the-Minute Research in Mass Spectrometry 22:2621–2628.\n\n\nWickham, H. 2011. The split-apply-combine strategy for data analysis. Journal of statistical software 40:1–29."
  },
  {
  "objectID": "references.html",

docs/sensitivity-analysis.html

@@ -209,6 +209,7 @@ <h2 id="toc-title">Table of contents</h2>
  <li><a href="#parameter-priors" id="toc-parameter-priors" class="nav-link active" data-scroll-target="#parameter-priors"><span class="header-section-number">6.1</span> Parameter priors</a></li>
  <li><a href="#using-a-mixed-effects-model-as-an-emulator" id="toc-using-a-mixed-effects-model-as-an-emulator" class="nav-link" data-scroll-target="#using-a-mixed-effects-model-as-an-emulator"><span class="header-section-number">6.2</span> Using a mixed effects model as an emulator</a></li>
  <li><a href="#how-do-the-process-rates-correlate" id="toc-how-do-the-process-rates-correlate" class="nav-link" data-scroll-target="#how-do-the-process-rates-correlate"><span class="header-section-number">6.3</span> How do the process rates correlate?</a></li>
+ <li><a href="#what-effect-does-the-bulk-density-have" id="toc-what-effect-does-the-bulk-density-have" class="nav-link" data-scroll-target="#what-effect-does-the-bulk-density-have"><span class="header-section-number">6.4</span> What effect does the bulk density have?</a></li>
  </ul>
 </nav>
  </div>
@@ -681,6 +682,240 @@ <h2 data-number="6.3" class="anchored" data-anchor-id="how-do-the-process-rates-
 </div>
 </div>
 <p>Two things are striking when looking at the clusters in <a href="#fig-pairs-processes">figure&nbsp;<span>6.2</span></a>: Both the within-group variance as well as the total variance (i.e.&nbsp;the variance of the cluster location) seems to be lower for the deeper depths. Additionally, one of the cluster centers constitutes a massive outlier – it’s column 12 at 7.5 cm depth (<a href="#fig-pairs-processes">Figure&nbsp;<span>6.2</span></a> <strong>B</strong>).</p>
+</section>
+<section id="what-effect-does-the-bulk-density-have" class="level2" data-number="6.4">
+<h2 data-number="6.4" class="anchored" data-anchor-id="what-effect-does-the-bulk-density-have"><span class="header-section-number">6.4</span> What effect does the bulk density have?</h2>
+<p>First, we’ll determine a range of bulk density to be searched.</p>
+<div class="sourceCode" id="cb4"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a>BD <span class="ot">&lt;-</span> <span class="fu">getParameters</span>()<span class="sc">$</span>BD <span class="sc">+</span> <span class="fu">seq</span>(<span class="sc">-</span><span class="fl">0.15</span>, <span class="fl">0.15</span>, <span class="fl">0.05</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<p>Next, we have to construct a function to return the process estimates for a give bulk density <code>x</code>.</p>
+<div class="sourceCode" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>f <span class="ot">&lt;-</span> <span class="cf">function</span>(x) {</span>
+<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a> parameters <span class="ot">&lt;-</span> <span class="fu">getParameters</span>(<span class="at">BD =</span> x)</span>
+<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a> data <span class="ot">&lt;-</span> <span class="fu">subset</span>(measurements, column<span class="sc">==</span><span class="dv">1</span>) <span class="sc">|&gt;</span></span>
+<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a> <span class="fu">getN2ON</span>(<span class="at">parameters =</span> parameters) <span class="sc">|&gt;</span></span>
+<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a> <span class="fu">getMissing</span>() <span class="sc">|&gt;</span></span>
+<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a> <span class="fu">calculateFluxes</span>(<span class="at">parameters =</span> parameters)</span>
+<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span>(<span class="fu">longPRE</span>(data, <span class="dv">1</span>, <span class="fl">7.5</span>, <span class="at">n =</span> <span class="dv">15</span>)<span class="sc">$</span>processes)</span>
+<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<p>Finally, we’ll apply said function to all bulk densities in the vector <code>BD</code>.</p>
+<div class="sourceCode" id="cb6"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a>result <span class="ot">&lt;-</span> <span class="fu">t</span>(<span class="fu">sapply</span>(BD, f))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
+<p>This process will return us the data to be analyzed.</p>
+<div id="fig-BD" class="quarto-figure quarto-figure-center anchored">
+<figure class="figure">
+<p><img src="scripts/sensitivity-analysis/output/bulk-density-sensitiviy.svg" class="img-fluid figure-img" style="width:100.0%"></p>
+<figcaption class="figure-caption">Figure&nbsp;6.3: Illustration to the sensitivity of estimated process rates to a change in prior bulk density.</figcaption>
+</figure>
+</div>
+
+<table style="text-align:center">
+<tbody><tr>
+<td colspan="4" style="border-bottom: 1px solid black">
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td colspan="3">
+<em>Dependent variable:</em>
+</td>
+</tr>
+<tr>
+<td>
+</td>
+<td colspan="3" style="border-bottom: 1px solid black">
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+Nitrification
+</td>
+<td>
+Denitrification
+</td>
+<td>
+Reduction
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+(1)
+</td>
+<td>
+(2)
+</td>
+<td>
+(3)
+</td>
+</tr>
+<tr>
+<td colspan="4" style="border-bottom: 1px solid black">
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+BD
+</td>
+<td>
+-20.907<sup>***</sup>
+</td>
+<td>
+-4.672<sup>***</sup>
+</td>
+<td>
+31.267<sup>***</sup>
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+(0.993)
+</td>
+<td>
+(0.619)
+</td>
+<td>
+(0.926)
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+</td>
+<td>
+</td>
+<td>
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+Constant
+</td>
+<td>
+47.404<sup>***</sup>
+</td>
+<td>
+26.665<sup>***</sup>
+</td>
+<td>
+-27.614<sup>***</sup>
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+(1.678)
+</td>
+<td>
+(1.045)
+</td>
+<td>
+(1.565)
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+</td>
+<td>
+</td>
+<td>
+</td>
+<td>
+</td>
+</tr>
+<tr>
+<td colspan="4" style="border-bottom: 1px solid black">
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+Observations
+</td>
+<td>
+7
+</td>
+<td>
+7
+</td>
+<td>
+7
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+R<sup>2</sup>
+</td>
+<td>
+0.989
+</td>
+<td>
+0.919
+</td>
+<td>
+0.996
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+Adjusted R<sup>2</sup>
+</td>
+<td>
+0.987
+</td>
+<td>
+0.903
+</td>
+<td>
+0.995
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+Residual Std. Error (df = 5)
+</td>
+<td>
+0.263
+</td>
+<td>
+0.164
+</td>
+<td>
+0.245
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+F Statistic (df = 1; 5)
+</td>
+<td>
+442.914<sup>***</sup>
+</td>
+<td>
+57.032<sup>***</sup>
+</td>
+<td>
+1,139.130<sup>***</sup>
+</td>
+</tr>
+<tr>
+<td colspan="4" style="border-bottom: 1px solid black">
+</td>
+</tr>
+<tr>
+<td style="text-align:left">
+<em>Note:</em>
+</td>
+<td colspan="3" style="text-align:right">
+<sup><em></em></sup><em>p&lt;0.1; <sup><strong></strong></sup><strong>p&lt;0.05; <sup></sup></strong></em>p&lt;0.01
+</td>
+</tr>
+
+</tbody></table>
 
 
 <div id="refs" class="references csl-bib-body hanging-indent" data-line-spacing="2" role="list" style="display: none">