<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/monte-carlo-based-prediction-intervals-for-nonlinear-regression/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://rmazing.wordpress.com/2018/05/11/monte-carlo-based-prediction-intervals-for-nonlinear-regression/"> Rmazing</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Calculation of the propagated uncertainty <img alt="sigma_y" title="sigma_y" src="https://s0.wp.com/latex.php?latex=%5Csigma_y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> using <img alt="nabla Sigma nabla^T" title="nabla Sigma nabla^T" src="https://s0.wp.com/latex.php?latex=%5Cnabla+%5CSigma+%5Cnabla%5ET&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (1), where <img alt="nabla" title="nabla" src="https://s0.wp.com/latex.php?latex=%5Cnabla&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is the gradient and <img alt="Sigma" title="Sigma" src="https://s0.wp.com/latex.php?latex=%5CSigma&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> the covariance matrix of the coefficients <img alt="beta_i" title="beta_i" src="https://s0.wp.com/latex.php?latex=%5Cbeta_i&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, is called the “Delta Method” and is widely applied in nonlinear least-squares (NLS) fitting. However, this method is based on first-order Taylor expansion and thus assummes linearity around <img alt="hat{y} = f(x, hat{beta})" title="hat{y} = f(x, hat{beta})" src="https://s0.wp.com/latex.php?latex=%5Chat%7By%7D+%3D+f%28x%2C+%5Chat%7B%5Cbeta%7D%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. The second-order approach can partially correct for this restriction by using a second-order polynomial around <img alt="hat{y}" title="hat{y}" src="https://s0.wp.com/latex.php?latex=%5Chat%7By%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, which is <img alt="nabla Sigma nabla^T + frac{1}{2} rm{tr}(H Sigma H Sigma)" title="nabla Sigma nabla^T + frac{1}{2} rm{tr}(H Sigma H Sigma)" src="https://s0.wp.com/latex.php?latex=%5Cnabla+%5CSigma+%5Cnabla%5ET+%2B+%5Cfrac%7B1%7D%7B2%7D+%5Crm%7Btr%7D%28H+%5CSigma+H+%5CSigma%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (2), where <img alt="rm{tr}(cdot)" title="rm{tr}(cdot)" src="https://s0.wp.com/latex.php?latex=%5Crm%7Btr%7D%28%5Ccdot%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is the matrix trace and <img alt="rm{H}" title="rm{H}" src="https://s0.wp.com/latex.php?latex=%5Crm%7BH%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is the Hessian.<br>
Confidence and prediction intervals for NLS models are calculated using <img alt="t(1 - frac{alpha}{2}, nu) cdot sigma_y" title="t(1 - frac{alpha}{2}, nu) cdot sigma_y" src="https://s0.wp.com/latex.php?latex=t%281+-+%5Cfrac%7B%5Calpha%7D%7B2%7D%2C+%5Cnu%29+%5Ccdot+%5Csigma_y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (3) or <img alt="t(1 - frac{alpha}{2}, nu) cdot sqrt{sigma_y^2 + sigma_r^2}" title="t(1 - frac{alpha}{2}, nu) cdot sqrt{sigma_y^2 + sigma_r^2}" src="https://s0.wp.com/latex.php?latex=t%281+-+%5Cfrac%7B%5Calpha%7D%7B2%7D%2C+%5Cnu%29+%5Ccdot+%5Csqrt%7B%5Csigma_y%5E2+%2B+%5Csigma_r%5E2%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (4), respectively, where the residual variance <img alt="sigma_r^2 = frac{sum_{i=1}^n (y_i - hat{y}_i)^2}{nu}" title="sigma_r^2 = frac{sum_{i=1}^n (y_i - hat{y}_i)^2}{nu}" src="https://s0.wp.com/latex.php?latex=%5Csigma_r%5E2+%3D+%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5En+%28y_i+-+%5Chat%7By%7D_i%29%5E2%7D%7B%5Cnu%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (5).<br>
Now, how can we employ the matrix notation of error propagation for creating Taylor expansion- and Monte Carlo-based prediction intervals?<br>
The inclusion of <img alt="sigma_r^2" title="sigma_r^2" src="https://s0.wp.com/latex.php?latex=%5Csigma_r%5E2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> in the prediction interval can be implemented as an extended gradient and “augmented” covariance matrix. So instead of using <img alt="hat{y} = f(x, hat{beta})" title="hat{y} = f(x, hat{beta})" src="https://s0.wp.com/latex.php?latex=%5Chat%7By%7D+%3D+f%28x%2C+%5Chat%7B%5Cbeta%7D%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (6), we take <img alt="hat{y} = f(x, hat{beta}) + sigma_r^2" title="hat{y} = f(x, hat{beta}) + sigma_r^2" src="https://s0.wp.com/latex.php?latex=%5Chat%7By%7D+%3D+f%28x%2C+%5Chat%7B%5Cbeta%7D%29+%2B+%5Csigma_r%5E2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (7) as the expression and augment the <img alt="i times i" title="i times i" src="https://s0.wp.com/latex.php?latex=i+%5Ctimes+i&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> covariance matrix <img alt="Sigma" title="Sigma" src="https://s0.wp.com/latex.php?latex=%5CSigma&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> to an <img alt="(i+1) times (i+1)" title="(i+1) times (i+1)" src="https://s0.wp.com/latex.php?latex=%28i%2B1%29+%5Ctimes+%28i%2B1%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> covariance matrix, where <img alt="Sigma_{i+1, i+1} = sigma_r^2" title="Sigma_{i+1, i+1} = sigma_r^2" src="https://s0.wp.com/latex.php?latex=%5CSigma_%7Bi%2B1%2C+i%2B1%7D+%3D+%5Csigma_r%5E2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. Partial differentiation and matrix multiplication will then yield, for example with two coefficients <img alt="beta_1" title="beta_1" src="https://s0.wp.com/latex.php?latex=%5Cbeta_1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="beta_2" title="beta_2" src="https://s0.wp.com/latex.php?latex=%5Cbeta_2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and their corresponding covariance matrix <img alt="Sigma" title="Sigma" src="https://s0.wp.com/latex.php?latex=%5CSigma&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">:<br>
<img alt="left[frac{partial f}{partial beta_1}; frac{partial f}{partial beta_2}; 1right] left[ begin{array}{ccc} ;sigma_1^2;;; sigma_1sigma_2;; 0 \ sigma_2sigma_1;; sigma_2^2;;;; 0 \ ;;;0;;;;;; 0;;;;; sigma_r^2 end{array} right] left[ begin{array}{c} frac{partial f}{partial beta_1} \ frac{partial f}{partial beta_2} \ 1 end{array} right]" title="left[frac{partial f}{partial beta_1}; frac{partial f}{partial beta_2}; 1right] left[ begin{array}{ccc} ;sigma_1^2;;; sigma_1sigma_2;; 0 \ sigma_2sigma_1;; sigma_2^2;;;; 0 \ ;;;0;;;;;; 0;;;;; sigma_r^2 end{array} right] left[ begin{array}{c} frac{partial f}{partial beta_1} \ frac{partial f}{partial beta_2} \ 1 end{array} right]" src="https://s0.wp.com/latex.php?latex=%5Cleft%5B%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_1%7D%5C%3B+%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_2%7D%5C%3B+1%5Cright%5D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+%5C%3B%5Csigma_1%5E2%5C%3B%5C%3B%5C%3B+%5Csigma_1%5Csigma_2%5C%3B%5C%3B+0+%5C%5C+%5Csigma_2%5Csigma_1%5C%3B%5C%3B+%5Csigma_2%5E2%5C%3B%5C%3B%5C%3B%5C%3B+0+%5C%5C+%5C%3B%5C%3B%5C%3B0%5C%3B%5C%3B%5C%3B%5C%3B%5C%3B%5C%3B+0%5C%3B%5C%3B%5C%3B%5C%3B%5C%3B+%5Csigma_r%5E2+%5Cend%7Barray%7D+%5Cright%5D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bc%7D+%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_1%7D+%5C%5C+%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_2%7D+%5C%5C+1+%5Cend%7Barray%7D+%5Cright%5D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (8)<br>
<img alt="= left(frac{partial f}{partial beta_1}right)^2sigma_1^2 + 2 frac{partial f}{partial beta_1} frac{partial f}{partial beta_2} sigma_1 sigma_2 + left(frac{partial f}{partial beta_2}right)^2sigma_2^2 + sigma_r^2" title="= left(frac{partial f}{partial beta_1}right)^2sigma_1^2 + 2 frac{partial f}{partial beta_1} frac{partial f}{partial beta_2} sigma_1 sigma_2 + left(frac{partial f}{partial beta_2}right)^2sigma_2^2 + sigma_r^2" src="https://s0.wp.com/latex.php?latex=%3D+%5Cleft%28%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_1%7D%5Cright%29%5E2%5Csigma_1%5E2+%2B+2+%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_1%7D+%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_2%7D+%5Csigma_1+%5Csigma_2+%2B+%5Cleft%28%5Cfrac%7B%5Cpartial+f%7D%7B%5Cpartial+%5Cbeta_2%7D%5Cright%29%5E2%5Csigma_2%5E2+%2B+%5Csigma_r%5E2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> (9)<br>
<img alt="equiv sigma_y^2 + sigma_r^2" title="equiv sigma_y^2 + sigma_r^2" src="https://s0.wp.com/latex.php?latex=%5Cequiv+%5Csigma_y%5E2+%2B+%5Csigma_r%5E2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, which then goes into (4).<br>
The advantage of the augmented covariance matrix is that it can be exploited for creating Monte Carlo-based prediction intervals. This is new from <a rel="nofollow" href="https://cran.r-project.org/web/packages/propagate/index.html" target="_blank">propagate version 1.0-6</a> and is based on the paradigm that we add another dimension by employing the augmented covariance matrix of (8) in the multivariate t-distribution random number generator (in our case <code>rtmvt</code>), with <img alt="mu = 0" title="mu = 0" src="https://s0.wp.com/latex.php?latex=%5Cmu+%3D+0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. All <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> simulations are then evaluated with (7) and the usual <img alt="[1 - frac{alpha}{2}, frac{alpha}{2}]" title="[1 - frac{alpha}{2}, frac{alpha}{2}]" src="https://s0.wp.com/latex.php?latex=%5B1+-+%5Cfrac%7B%5Calpha%7D%7B2%7D%2C+%5Cfrac%7B%5Calpha%7D%7B2%7D%5D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> quantiles calculated for the prediction interval. Using the original covariance matrix with (6) will deliver the MC-based confidence interval.<br>
Application of second-order Taylor expansion or the MC-based approach demonstrates nicely that for the majority of nonlinear models, the confidence/prediction intervals around <img alt="hat{y}" title="hat{y}" src="https://s0.wp.com/latex.php?latex=%5Chat%7By%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> are quite asymmetric, which the classical Delta Method does not capture:<br>
<code><br>
library(propagate)
<p>DNase1 &lt;-subset(DNase, Run == 1)<br>
fm3DNase1 &lt;-nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),<br>
data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))<br>
## first-order prediction interval<br>
set.seed(123)<br>
PROP1 &lt;-predictNLS(fm3DNase1, newdata = data.frame(conc = 2), nsim = 1000000,<br>
second.order = FALSE, interval = &quot;prediction&quot;)<br>
t(PROP1$summary)<br>
</p>
<p><img alt="begin{array}{cc}  Prop.Mean.1 &amp; 0.74804722 \  Prop.sd.1 &amp; 0.02081131 \  Prop.2.5 &amp; 0.70308712 \  Prop.97.5 &amp; 0.79300731 \  end{array}" title="begin{array}{cc}  Prop.Mean.1 &amp; 0.74804722 \  Prop.sd.1 &amp; 0.02081131 \  Prop.2.5 &amp; 0.70308712 \  Prop.97.5 &amp; 0.79300731 \  end{array}" src="https://s0.wp.com/latex.php?latex=%5Cbegin%7Barray%7D%7Bcc%7D++Prop.Mean.1+%26+0.74804722+%5C%5C++Prop.sd.1+%26+0.02081131+%5C%5C++Prop.2.5+%26+0.70308712+%5C%5C++Prop.97.5+%26+0.79300731+%5C%5C++%5Cend%7Barray%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p><code><br>
## second-order prediction interval and MC<br>
set.seed(123)<br>
PROP2 &lt;- predictNLS(fm3DNase1, newdata = data.frame(conc = 2), nsim = 1000000,<br>
                    second.order = TRUE, interval = &quot;prediction&quot;)<br>
t(PROP2$summary)<br>
</code></p>
<p><img alt="begin{array}{cc}  Prop.Mean.1 &amp; 0.74804722 \  Prop.Mean.2 &amp; 0.74815136 \  Prop.sd.1 &amp;  0.02081131 \  Prop.sd.2 &amp;  0.02081520 \  Prop.2.5 &amp;  0.70318286 \  Prop.97.5 &amp; 0.79311987 \  Sim.Mean &amp; 0.74815598 \  Sim.sd &amp; 0.02261884 \  Sim.2.5 &amp; 0.70317629 \  Sim.97.5 &amp; 0.79309874 \  end{array}" title="begin{array}{cc}  Prop.Mean.1 &amp; 0.74804722 \  Prop.Mean.2 &amp; 0.74815136 \  Prop.sd.1 &amp;  0.02081131 \  Prop.sd.2 &amp;  0.02081520 \  Prop.2.5 &amp;  0.70318286 \  Prop.97.5 &amp; 0.79311987 \  Sim.Mean &amp; 0.74815598 \  Sim.sd &amp; 0.02261884 \  Sim.2.5 &amp; 0.70317629 \  Sim.97.5 &amp; 0.79309874 \  end{array}" src="https://s0.wp.com/latex.php?latex=%5Cbegin%7Barray%7D%7Bcc%7D++Prop.Mean.1+%26+0.74804722+%5C%5C++Prop.Mean.2+%26+0.74815136+%5C%5C++Prop.sd.1+%26++0.02081131+%5C%5C++Prop.sd.2+%26++0.02081520+%5C%5C++Prop.2.5+%26++0.70318286+%5C%5C++Prop.97.5+%26+0.79311987+%5C%5C++Sim.Mean+%26+0.74815598+%5C%5C++Sim.sd+%26+0.02261884+%5C%5C++Sim.2.5+%26+0.70317629+%5C%5C++Sim.97.5+%26+0.79309874+%5C%5C++%5Cend%7Barray%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>What we see here is that<br>
i) the first-order prediction interval [0.70308712; 0.79300731] <strong>is symmetric and slightly down-biased</strong> compared to the second-order one [0.70317629; 0.79309874],<br>
and<br>
ii) the second-order prediction interval <strong>tallies nicely up to the 4th decimal</strong> with the new MC-based interval (0.70318286 and 0.70317629; 0.79311987 and 0.79309874).<br>
I believe this clearly demonstrates the usefulness of the MC-based approach for NLS prediction interval estimation…</p>

				
		
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://rmazing.wordpress.com/2018/05/11/monte-carlo-based-prediction-intervals-for-nonlinear-regression/"> Rmazing</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/UU9kE-GjroY" height="1"></a></div></code>