<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/pk-pd-reserving-models/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/"> R on mages' blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This is a follow-up post on <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank">hierarchical compartmental reserving models</a> using PK/PD models. It will show how differential equations can be used with <a rel="nofollow" href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/claimsprocess" target="_blank">Stan</a>/ <a rel="nofollow" href="https://cran.r-project.org/package=brms" target="_blank"><code>brms</code></a> and how correlation for the same group level terms can be modelled.</p>
<p>PK/ PD is usually short for pharmacokinetic/ pharmacodynamic models, but as <a rel="nofollow" href="https://www.linkedin.com/in/enovik/" target="_blank">Eric Novik</a> of <a rel="nofollow" href="https://www.generable.com/" target="_blank">Generable</a> pointed out to me, it could also be short for <em>Payment Kinetics/ Payment Dynamics Models</em> in the insurance context.</p>
<p>In this post I’d like to discuss an extension to Jake Morris’ <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank"><em>hierarchical compartmental reserving model</em></a>, as described in his original paper (<span>Morris (<a rel="nofollow" href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/#ref-JakeMorris2016" target="_blank">2016</a>)</span>) and my <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank">previous post</a>, which allows for a time varying parameter <span>\(k_{er}(t)\)</span> describing the changing rate of earning and reporting and allows for correlation between <span>\(RLR\)</span>, the reported loss ratio, and <span>\(RRF\)</span>, the reserve robustness factor. A positive correlation would give evidence of a reserving cycle, i.e. in years with higher initial reported loss ratios a higher proportion of reported losses are paid.</p>
<p>I am very grateful for the support Jake Morris and Paul-Christian Bürkner have given me over the last few weeks, answering questions around the model and <a rel="nofollow" href="https://cran.r-project.org/package=brms" target="_blank"><code>brms</code></a> (<span>Bürkner (<a rel="nofollow" href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/#ref-brms" target="_blank">2017</a>)</span>).</p>
<div>
<h1>PK/ PD model</h1>
<p>First of all let’s take a look at my new set of differential equations, describing the dynamics of exposure (<span>\(EX\)</span>), outstanding (<span>\(OS\)</span>) and paid claims (<span>\(PD\)</span>) over time. I replaced the constant parameter <span>\(k_{er}\)</span> with the linear function <span>\(\beta_{er} \cdot t\)</span>.</p>
<p><span>\[<br>
\begin{aligned}<br>
dEX/dt &amp; = -\beta_{er} \cdot t \cdot EX \\<br>
dOS/dt &amp; = \beta_{er} \cdot t \cdot RLR \cdot EX – k_p \cdot OS \\<br>
dPD/dt &amp; = k_{p} \cdot RRF \cdot OS<br>
\end{aligned}<br>
\]</span></p>
<p>The dynamical system is no longer autonomous and initially I can’t be bothered to solve it analytically. Hence, I use an ODE solver instead, but I will get back to integrating the differential equations later.</p>
<div>
<h2>Integrating ODEs with Stan</h2>
<p>Fortunately, an ODE solver is part of the Stan language. The following code demonstrates how I can integrate the differential equations in Stan. You will notice that the model section is empty, as all I want to do for now is run the solver.</p>
<pre><code>functions{
  real[] claimsprocess(real t, real [] y, real [] theta, 
                       real [] x_r, int[] x_i){
    real dydt[3];
    dydt[1] = - theta[1] * t * y[1];
    dydt[2] = theta[1] * t * theta[3] * y[1] - theta[2] * y[2];
    dydt[3] = theta[2] * theta[4] * y[2];
    
    return dydt;
  }
}
data {
  real theta[4];
  int N_t;
  real times[N_t];  
  real C0[3];
}
parameters{
}
transformed parameters{
  real C[N_t, 3];
  C = integrate_ode_rk45(claimsprocess, C0, 0, times, theta,
                         rep_array(0.0, 0), rep_array(1, 1));
}
model{
}
generated quantities{
  real Exposure[N_t];
  real OS[N_t];
  real Paid[N_t];
  Exposure = C[, 1];
  OS = C[, 2];
  Paid = C[, 3];
}</code></pre>
<p>To run the Stan code I provide as an input a list with the parameters, the number of data points I’d like to generate, the time line and initial starting values. Note also that I have to set the <code>algorithm = &quot;Fixed_param&quot;</code> as I run a deterministic model and for that reason one iteration and one chain are sufficient.</p>
<pre><code>library(rstan)
input_data &lt;- list(
  theta=c(Ber=5, kp=0.5, RLR=0.8, RRF=0.9),
  N_t=100, 
  times=seq(from = 0.1, to = 10, length.out = 100),
  C0=c(100, 0, 0))
samples &lt;- sampling(ODEmodel, data=input_data, 
                    algorithm = &quot;Fixed_param&quot;, 
                    iter = 1, chain = 1)</code></pre>
<p>Let’s take a look at the output.</p>
<pre><code>X &lt;- extract(samples)
out &lt;- data.frame(
  Time = input_data$times,
  Exposure = c(X[[&quot;Exposure&quot;]]),
  Outstandings = c(X[[&quot;OS&quot;]]),
  Paid = c(X[[&quot;Paid&quot;]])
)
library(latticeExtra)
xyplot(Exposure + Outstandings + Paid ~ Time,
       data = out, as.table=TRUE, t=&quot;l&quot;,
       ylab=&quot;Amounts ($)&quot;, auto.key=list(space=&quot;top&quot;, columns=3),
       par.settings = theEconomist.theme(with.bg = TRUE, box = &quot;transparent&quot;))</code></pre>
<p><img src="https://i0.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/unnamed-chunk-2-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The exposure declines faster and the outstanding curve has a more pronounced peak as a result of replacing <span>\(k_{er}\)</span> with <span>\(\beta_{er} \cdot t\)</span>.</p>
</div>
</div>
<div>
<h1>Updated model</h1>
<p>My updated model is not much different to the one presented in the earlier <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank">post</a>, apart from the fact that I allow for the correlation between <span>\(RLR\)</span> and <span>\(RRF\)</span> and the mean function <span>\(\tilde{f}\)</span> is the integral of the ODEs above.</p>
<p><span>\[<br>
\begin{aligned}<br>
y(t) &amp; \sim \mathcal{N}(\tilde{f}(t, \Pi, \beta_{er}, k_p, RLR_{[i]}, RRF_{[i]}), \sigma_{y[\delta]}^2) \\<br>
\begin{pmatrix} RLR_{[i]} \\ RRF_{[i]}\end{pmatrix} &amp; \sim<br>
\mathcal{N} \left(<br>
\begin{pmatrix}<br>
\mu_{RLR} \\<br>
\mu_{RRF}<br>
\end{pmatrix},<br>
\begin{pmatrix}<br>
\sigma_{RLR}^2 &amp; \rho \ \sigma_{RLR} \sigma_{RRF}\\<br>
\rho \ \sigma_{RLR} \sigma_{RRF} &amp; \sigma_{RRF}^2<br>
\end{pmatrix}<br>
\right)<br>
\end{aligned}<br>
\]</span></p>
<div>
<h2>Implementation with <code>brms</code></h2>
<p>Let’s load the data back into R’s memory:</p>
<pre><code>library(data.table)
lossData0 &lt;- fread(&quot;https://raw.githubusercontent.com/mages/diesunddas/master/Data/WorkersComp337.csv&quot;)</code></pre>
<p>Jake shows in the appendices of his <a rel="nofollow" href="https://www.casact.org/pubs/forum/16sforum/Morris.pdf" target="_blank">paper</a> how to implement this model in R with the <a rel="nofollow" href="https://cran.r-project.org/package=nlmeODE" target="_blank"><code>nlmeODE</code></a> (<span>Tornoe (<a rel="nofollow" href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/#ref-nlmeODE" target="_blank">2012</a>)</span>) package, together with more flexible models in OpenBUGS (<span>Lunn et al. (<a rel="nofollow" href="https://magesblog.com/post/2018-01-30-pkpd-reserving-models/#ref-lunn2000winbugs" target="_blank">2000</a>)</span>). However, I will continue with <code>brms</code> and Stan.</p>
<p>Using the ODEs with <code>brms</code> requires a little extra coding, as I have to provide the integration function as an additional input, just like I did above, where I defined them as user defined functions for Stan.</p>
<p>To model the group-level terms (<span>\(RLR\)</span>, <span>\(RRF\)</span>) of the same grouping factor (accident year) as correlated I add a unique identifier to the <code>|</code> operator, here I use <code>p</code> (looks similar to <span>\(\rho\)</span>), i.e. <code>~ 1 + (1 | p | accident_year)</code>.</p>
<p>Below is my code for the new updated model. I am using again Gamma distributions as my priors for the parameters and <a rel="nofollow" href="http://stla.github.io/stlapblog/posts/StanLKJprior.html" target="_blank">LKJ(2)</a> as a prior for the correlation coefficient.</p>
<pre><code>myFuns &lt;- &quot;
real[] ode_claimsprocess(real t, real [] y,  real [] theta, 
                         real [] x_r, int[] x_i){
  real dydt[3];
  
  dydt[1] = - theta[1] * t * y[1]; // Exposure
  dydt[2] = theta[1] * t * theta[2] * y[1] - theta[3] * y[2]; // OS
  dydt[3] = theta[3] * theta[4] * y[2]; // Paid
  
  return dydt;
}
real claimsprocess(real t, real premium, real Ber, real kp, 
                   real RLR, real RRF, real delta){
  real y0[3];
  real y[1, 3];
  real theta[4];
  theta[1] = Ber;
  theta[2] = RLR;
  theta[3] = kp;
  theta[4] = RRF;
  y0[1] = premium;
  y0[2] = 0;
  y0[3] = 0;
  y = integrate_ode_rk45(ode_claimsprocess, 
                        y0, 0, rep_array(t, 1), theta,
                        rep_array(0.0, 0), rep_array(1, 1),
                        0.001, 0.001, 100); // tolerances, steps
      return (y[1, 2] * (1 - delta) + y[1, 3] * delta);
}
&quot;
library(brms) 
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
b4 &lt;- brm(
  bf(loss_train ~ claimsprocess(dev, premium, Ber, kp, RLR, RRF, delta),
     RLR ~ 1 + (1 | p | accident_year), # 'p' allow for correlation with RRF 
     RRF ~ 1 + (1 | p | accident_year), # 'p' allow for correlation with RLR
     Ber ~ 1,
     kp ~ 1,
     sigma ~ 0 + deltaf, # different sigma for OS and paid
     nl = TRUE), 
  stan_funs = myFuns, # defintion of 'claimsprocess' function
  data = lossData0[cal &lt;= max(accident_year) &amp; dev > 0], 
  family = brmsfamily(&quot;gaussian&quot;, link_sigma = &quot;log&quot;),
  prior = c(prior(gamma(4, 5), nlpar = &quot;RLR&quot;, lb=0),
            prior(gamma(4, 5), nlpar = &quot;RRF&quot;, lb=0),
            prior(gamma(12, 3), nlpar = &quot;Ber&quot;, lb=0),
            prior(gamma(3, 4), nlpar = &quot;kp&quot;, lb=0),
            set_prior(&quot;lkj(2)&quot;, class = &quot;cor&quot;)),
  control = list(adapt_delta = 0.999, max_treedepth=15),
  seed = 1234, iter = 500)</code></pre>
<p>While the original model with the analytical solution ran in about 3 minutes, this code took about 9 hours per chain for 500 samples. Any suggestions how to speed this up will be much appreciated. The frequentist model using <code>nlmeODE</code> runs in seconds.</p>
<pre><code>b4</code></pre>
<pre><code>## Warning: There were 235 divergent transitions after warmup. Increasing adapt_delta above 0.999 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup</code></pre>
<pre><code>##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: loss_train ~ claimsprocess(dev, premium, Ber, kp, RLR, RRF, delta) 
##          RLR ~ 1 + (1 | p | accident_year)
##          RRF ~ 1 + (1 | p | accident_year)
##          Ber ~ 1
##          kp ~ 1
##          sigma ~ 0 + deltaf
##    Data: lossData0[cal &lt;= max(accident_year) &amp; dev > 0] (Number of observations: 110) 
## Samples: 4 chains, each with iter = 500; warmup = 250; thin = 1; 
##          total post-warmup samples = 1000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Group-Level Effects: 
## ~accident_year (Number of levels: 10) 
##                                  Estimate Est.Error l-95% CI u-95% CI
## sd(RLR_Intercept)                    0.18      0.06     0.11     0.32
## sd(RRF_Intercept)                    0.15      0.05     0.09     0.27
## cor(RLR_Intercept,RRF_Intercept)     0.54      0.24    -0.04     0.88
##                                  Eff.Sample Rhat
## sd(RLR_Intercept)                       197 1.01
## sd(RRF_Intercept)                       397 1.00
## cor(RLR_Intercept,RRF_Intercept)        326 1.01
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## RLR_Intercept        0.86      0.06     0.74     0.99        161 1.01
## RRF_Intercept        0.83      0.05     0.73     0.93        179 1.02
## Ber_Intercept        5.66      0.32     4.97     6.21        359 1.01
## kp_Intercept         0.40      0.01     0.38     0.41        446 1.00
## sigma_deltafos      -3.61      0.11    -3.81    -3.40        593 1.00
## sigma_deltafpaid    -4.99      0.11    -5.20    -4.79        587 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).</code></pre>
<p>The output above looks very similar to the output of <code>nlmeODE</code>. That’s good, also that <code>Rhat</code> is close to 1 for all parameters. However, there is a warning message of 235 divergent transitions after warm-up, which I will ignore for the time being.</p>
<p>The correlation coefficient <span>\(\rho\)</span> between <span>\(RLR_{[i]}\)</span> and <span>\(RRF_{[i]}\)</span> is estimated as 0.54, but with a wide 95% credible interval from -0.04 to 0.88. Therefore there is moderate evidence of a correlation.</p>
<p>As Jake puts it, a positive correlation between the reported loss ratio and reserve robustness factor parameters by accident year is indicative of a case reserving cycle effect, i.e. more conservative case reserves (low <span>\(RRF_{[i]}\)</span>) in a hard market (low <span>\(RLR_{[i]}\)</span>) to create cushions for the future.</p>
<pre><code>stanplot(b4, 
         pars= c(&quot;b_RLR_Intercept&quot;, &quot;b_RRF_Intercept&quot;, &quot;b_Ber_Intercept&quot;,
                 &quot;b_kp_Intercept&quot;, &quot;b_sigma_deltafos&quot;, &quot;b_sigma_deltafpaid&quot;,
                 &quot;sd_accident_year__RLR_Intercept&quot;, 
                 &quot;sd_accident_year__RRF_Intercept&quot;, 
                 &quot;cor_accident_year__RLR_Intercept__RRF_Intercept&quot;), 
         type=&quot;dens_overlay&quot;)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/densityPlot1-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The density plots look OK, all chains seem to have behaved similarly. The last chart in the bottom right shows the distribution of the correlation parameter <span>\(\rho\)</span>. Moderate evidence of a correlation might be an understatement. Note also that <span>\(\sigma_{\delta[OS]}\)</span> and <span>\(\sigma_{\delta[PD]}\)</span> are on a log-scale.</p>
</div>
<div>
<h2>Predict future claims developments</h2>
<p>To validate my model against the test data set I can use the <code>predict</code> function in <code>brms</code>. Before I can apply the method I have to expose the Stan functions I wrote above, namely <code>claimsprocess</code> to R.</p>
<pre><code>expose_functions(b4, vectorize = TRUE) # requires brms >= 2.1.0</code></pre>
<p>With the <code>claimsprocess</code> function now in the R’s memory I can predict and plot the model in the same way as I did in the <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank">previous post</a>.</p>
<pre><code>predClaimsPred &lt;- predict(b4, newdata = lossData0, method=&quot;predict&quot;)
plotDevBananas(`2.5%ile`/premium*100 + `97.5%ile`/premium*100 + 
                 Estimate/premium*100 + loss_train/premium*100 + 
                 loss_test/premium*100 ~ 
                 dev | factor(accident_year), ylim=c(0, 100),
               data=cbind(lossData0, predClaimsPred)[delta==0], 
               main=&quot;Outstanding claims developments&quot;, 
               ylab=&quot;Outstanding loss ratio (%)&quot;)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/predData-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<pre><code>plotDevBananas(`2.5%ile`/premium*100 + `97.5%ile`/premium*100 + 
                 Estimate/premium*100 + loss_train/premium*100 + 
                 loss_test/premium*100 ~  
                 dev | factor(accident_year), ylim=c(0, 120),
               data=cbind(lossData0, predClaimsPred)[delta==1], 
               main=&quot;Paid claims developments&quot;, 
               ylab=&quot;Paid loss ratio (%)&quot;)</code></pre>
<p><img src="https://i0.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/unnamed-chunk-4-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The graphs look very similar to the output from the <a rel="nofollow" href="https://magesblog.com/post/hierarchical-compartmental-reserving-models/" target="_blank">previous model</a>. Yet, the prediction for 1997 improved significantly, perhaps because of the correlation between RLR and RRF no longer being assumed to be zero. But this model has tested my patience, as I used a remote 2-core machine on <a rel="nofollow" href="https://m.do.co/c/3a02238dca6d" target="_blank">Digital Ocean</a> and waited 18 hours for the results.</p>
</div>
</div>
<div>
<h1>Analytical solution</h1>
<p>The long run time did in the end motivate me to look for an analytical solution of the ODEs, with initial values <span>\(\mbox{EX}(0) = \Pi\)</span> (premium), <span>\(\mbox{OS}(0) = 0\)</span>, <span>\(\mbox{PD}(0) = 0\)</span>.</p>
<p>According to <a rel="nofollow" href="https://www.wolframalpha.com/input/?i=solve+%7Bf%27(x)+%3D+-x++b++f(x),+g%27(x)+%3D+b++x+++L++f(x)+-+k++g(x),+h%27(x)+%3D+k++R++g(x),+f(0)%3DP,+g(0)%3D0,+h(0)%3D0%7D" target="_blank">Wolfram Alpha the analytical solution</a> is:</p>
<p><span>\[<br>
\begin{aligned}<br>
\mbox{EX}(t) &amp; = \Pi \exp \left(-\frac{\beta_{er} t^2}{2} \right) \\<br>
\mbox{OS}(t) &amp; = – \frac{\Pi \cdot RLR}{2 \sqrt{\beta_{er}}}<br>
\exp\left(-\frac{\beta_{er} t^2}{2} – k_p t \right)<br>
\left[ \sqrt{2 \pi} k_p \mbox{erf}\left(\frac{k_p}{\sqrt{2 \beta_{er}}}\right)<br>
\left(-\exp\left(\frac{k_p^2}{2 \beta_{er}} + \frac{\beta_{er} t^2}{2} \right) \right) \right. – \\<br>
&amp; \left. \qquad \sqrt{2 \pi} k_p \exp\left(\frac{k_p^2}{2 \beta_{er}} + \frac{\beta_{er} t^2}{2}\right)<br>
\mbox{erf}\left(\frac{\beta_{er} t – k_p}{\sqrt{2 \beta_{er}}}\right) + 2 \sqrt{\beta_{er}}<br>
\exp\left(k_p t\right) – 2 \sqrt{\beta_{er}}<br>
\exp\left(\frac{\beta_{er} t^2}{2}\right)<br>
\right] \\<br>
\mbox{PD}(t) &amp; = \frac{\Pi \cdot RLR \cdot RRF}{2 \sqrt{\beta_{er}}}<br>
\exp\left(-k_p t\right) \left[-\sqrt{2 \pi} k_p<br>
\exp\left(\frac{k_p^2}{2 \beta_{er}}\right)<br>
\mbox{erf}\left(\frac{\beta_{er} t – k_p}{\sqrt{2 \beta_{er}}}\right) \right. + \\<br>
&amp; \left. \qquad \sqrt{2 \pi} k_p<br>
\left(-\exp\left(\frac{k^2}{2 \beta_{er}}\right)\right)<br>
\mbox{erf}\left(\frac{k}{\sqrt{2 \beta_{er}}} \right) +<br>
2 \sqrt{\beta_{er}} \exp\left(k_p t\right) –<br>
2 \sqrt{\beta_{er}}<br>
\right]<br>
\end{aligned}<br>
\]</span></p>
<p>Now I can write down <span>\(\tilde{f}(t)\)</span> in a closed-form: <span>\[<br>
\tilde{f}(t, \delta,\dots) = (1 – \delta) \cdot \mbox{OS}(t, \dots) + \delta \cdot \mbox{PD}(t, \dots)<br>
\]</span></p>
<p>The <a rel="nofollow" href="https://en.wikipedia.org/wiki/Error_function" target="_blank">error function</a> is part of the Stan language (not directly in R, although it is just <code>erf &lt;- function(x) 2 * pnorm(x * sqrt(2)) - 1</code>). Thus, I create a Stan function for the analytical claims process again and feed this into <code>brm</code>.</p>
<pre><code>library(brms)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

myFun &lt;- &quot;
real anaclaimsprocess(real t, real premium, real Ber, real kp, 
                      real RLR, real RRF, real delta){

 real os;
 real paid;

 os = -(RLR * exp(-(Ber * t^2)/2.0 - kp * t) * 
          (sqrt(2.0 * pi()) * kp * erf(kp/(sqrt(2.0) * sqrt(Ber))) * 
             (-exp(kp^2/(2.0 * Ber) + (Ber * t^2)/2.0)) - 
             sqrt(2.0 * pi()) * kp * exp(kp^2/(2.0 * Ber) + (Ber *t^2)/2.0) * 
             erf((Ber*t - kp)/(sqrt(2.0) * sqrt(Ber))) + 
             2.0 * sqrt(Ber) * exp(kp * t) - 2.0 * sqrt(Ber) * 
             exp((Ber*t^2)/2.0)))/(2.0 * sqrt(Ber));
 paid = (RLR * RRF * exp(-kp * t) * 
           (-sqrt(2.0 * pi()) * kp * exp(kp^2/(2.0 * Ber)) * 
              erf((Ber*t - kp)/(sqrt(2.0) * sqrt(Ber))) + 
              sqrt(2.0 * pi()) * kp * (-exp(kp^2/(2.0 * Ber))) * 
              erf(kp/(sqrt(2.0) * sqrt(Ber))) + 2.0 * sqrt(Ber) * exp(kp * t) - 
              2.0 * sqrt(Ber)))/(2.0 * sqrt(Ber));

 return (premium * (os * (1 - delta) + paid * delta)); 
}
&quot;

b5 &lt;- brm(bf(loss_train ~ anaclaimsprocess(dev, premium, Ber, kp, 
                                           RLR, RRF, delta),
             RLR ~ 1 + (1 | p | accident_year),
             RRF ~ 1 + (1 | p | accident_year),
             Ber ~ 1,
             kp ~ 1,
             sigma ~ 0 + deltaf,
             nl = TRUE),
          stan_funs = myFun, # defintion of 'anaclaimsprocess' function
          data = lossData0[cal &lt;= max(accident_year)], 
          family = brmsfamily(&quot;gaussian&quot;, link_sigma = &quot;log&quot;),
          prior = c(prior(gamma(4, 5), nlpar = &quot;RLR&quot;, lb=0),
                    prior(gamma(4, 5), nlpar = &quot;RRF&quot;, lb=0),
                    prior(gamma(12, 3), nlpar = &quot;Ber&quot;, lb=0),
                    prior(gamma(3, 4), nlpar = &quot;kp&quot;, lb=0),
                    set_prior(&quot;lkj(2)&quot;, class = &quot;cor&quot;)),
          control = list(adapt_delta = 0.999, max_treedepth=15),
          seed = 1234, iter = 1000)</code></pre>
<p>Wonderful, the code ran in about 5 minutes, with twice the number of samples and without any warnings messages. Let’s take a look at the output:</p>
<pre><code>b5</code></pre>
<pre><code>##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: loss_train ~ anaclaimsprocess(dev, premium, Ber, kp, RLR, RRF, delta) 
##          RLR ~ 1 + (1 | p | accident_year)
##          RRF ~ 1 + (1 | p | accident_year)
##          Ber ~ 1
##          kp ~ 1
##          sigma ~ 0 + deltaf
##    Data: lossData0[cal &lt;= max(accident_year)] (Number of observations: 110) 
## Samples: 4 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 2000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Group-Level Effects: 
## ~accident_year (Number of levels: 10) 
##                                  Estimate Est.Error l-95% CI u-95% CI
## sd(RLR_Intercept)                    0.18      0.06     0.11     0.32
## sd(RRF_Intercept)                    0.15      0.04     0.09     0.26
## cor(RLR_Intercept,RRF_Intercept)     0.52      0.26    -0.08     0.89
##                                  Eff.Sample Rhat
## sd(RLR_Intercept)                       852 1.01
## sd(RRF_Intercept)                       860 1.00
## cor(RLR_Intercept,RRF_Intercept)        885 1.00
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## RLR_Intercept        0.86      0.06     0.73     0.97        480 1.00
## RRF_Intercept        0.83      0.05     0.72     0.93        723 1.00
## Ber_Intercept        5.69      0.35     5.03     6.42       1968 1.00
## kp_Intercept         0.40      0.01     0.38     0.41       1701 1.00
## sigma_deltafos      -3.61      0.11    -3.80    -3.40       1337 1.00
## sigma_deltafpaid    -4.99      0.11    -5.19    -4.75       1547 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).</code></pre>
<p>Perfect, the estimates are very much the same as from the ODE model and the plots haven’t changed much either.</p>
<pre><code>stanplot(b5, 
         pars= c(&quot;b_RLR_Intercept&quot;, &quot;b_RRF_Intercept&quot;, &quot;b_Ber_Intercept&quot;,
                 &quot;b_kp_Intercept&quot;, &quot;b_sigma_deltafos&quot;, &quot;b_sigma_deltafpaid&quot;,
                 &quot;sd_accident_year__RLR_Intercept&quot;, 
                 &quot;sd_accident_year__RRF_Intercept&quot;, 
                 &quot;cor_accident_year__RLR_Intercept__RRF_Intercept&quot;), 
         type=&quot;dens_overlay&quot;)</code></pre>
<p><img src="https://i0.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/densityPlot2-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The density plots look smoother and more consistent across chains with twice the number of samples compared to the previous model.</p>
<pre><code>expose_functions(b5, vectorize = TRUE) # requires brms >= 2.1.0</code></pre>
<pre><code>predClaimsPred2 &lt;- predict(b5, newdata = lossData0, method=&quot;predict&quot;)
plotDevBananas(`2.5%ile`/premium*100 + `97.5%ile`/premium*100 + 
                 Estimate/premium*100 + loss_train/premium*100 + 
                 loss_test/premium*100 ~ 
                 dev | factor(accident_year), ylim=c(0, 100),
               data=cbind(lossData0, predClaimsPred2)[delta==0], 
               main=&quot;Outstanding claims developments&quot;, 
               ylab=&quot;Outstanding loss ratio (%)&quot;)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/predData2-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<pre><code>plotDevBananas(`2.5%ile`/premium*100 + `97.5%ile`/premium*100 + 
                 Estimate/premium*100 + loss_train/premium*100 + 
                 loss_test/premium*100 ~  
                 dev | factor(accident_year), ylim=c(0, 120),
               data=cbind(lossData0, predClaimsPred2)[delta==1], 
               main=&quot;Paid claims developments&quot;, 
               ylab=&quot;Paid loss ratio (%)&quot;)</code></pre>
<p><img src="https://i1.wp.com/magesblog.com/post/2018-01-30-pkpd-reserving-models_files/figure-html/unnamed-chunk-6-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>I am much happier now. The model runs in an acceptable time, allowing me to play around with my assumptions further. I have yet to understand why the integration routine in Stan took so long.</p>
</div>
<div>
<h1>Session Info</h1>
<pre><code>sessionInfo()</code></pre>
<pre><code>## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] methods   stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
## [1] brms_2.1.0          Rcpp_0.12.15        data.table_1.10.4-3
## [4] latticeExtra_0.6-28 RColorBrewer_1.1-2  lattice_0.20-35    
## [7] rstan_2.17.3        StanHeaders_2.17.2  ggplot2_2.2.1      
## 
## loaded via a namespace (and not attached):
##  [1] mvtnorm_1.0-7        gtools_3.5.0         zoo_1.8-1           
##  [4] assertthat_0.2.0     rprojroot_1.3-2      digest_0.6.15       
##  [7] mime_0.5             R6_2.2.2             plyr_1.8.4          
## [10] backports_1.1.2      stats4_3.4.3         evaluate_0.10.1     
## [13] coda_0.19-1          colourpicker_1.0     blogdown_0.5        
## [16] pillar_1.1.0         rlang_0.1.6          lazyeval_0.2.1      
## [19] curl_3.1             miniUI_0.1.1         DT_0.3.3            
## [22] Matrix_1.2-12        rmarkdown_1.8        labeling_0.3        
## [25] shinythemes_1.1.1    shinyjs_1.0          stringr_1.2.0       
## [28] htmlwidgets_1.0      loo_1.1.0            igraph_1.1.2        
## [31] munsell_0.4.3        shiny_1.0.5          compiler_3.4.3      
## [34] httpuv_1.3.5         xfun_0.1             pkgconfig_2.0.1     
## [37] base64enc_0.1-3      rstantools_1.4.0     htmltools_0.3.6     
## [40] tibble_1.4.2         gridExtra_2.3        bookdown_0.6        
## [43] codetools_0.2-15     threejs_0.3.1        matrixStats_0.53.0  
## [46] dplyr_0.7.4          grid_3.4.3           nlme_3.1-131        
## [49] xtable_1.8-2         gtable_0.2.0         magrittr_1.5        
## [52] scales_0.5.0         stringi_1.1.6        reshape2_1.4.3      
## [55] bindrcpp_0.2         dygraphs_1.1.1.4     xts_0.10-1          
## [58] tools_3.4.3          glue_1.2.0           markdown_0.8        
## [61] shinystan_2.4.0      crosstalk_1.0.0      rsconnect_0.8.5     
## [64] abind_1.4-5          parallel_3.4.3       yaml_2.1.16         
## [67] inline_0.3.14        colorspace_1.3-2     bridgesampling_0.4-0
## [70] bayesplot_1.4.0      knitr_1.18           bindr_0.1           
## [73] Brobdingnag_1.2-4</code></pre>
</div>
<div>
<h1>References</h1>
<div>
<div>
<p>Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” <em>Journal of Statistical Software</em> 80 (1): 1–28. doi:<a rel="nofollow" href="https://doi.org/10.18637/jss.v080.i01" target="_blank">10.18637/jss.v080.i01</a>.</p>
</div>
<div>
<p>Lunn, David J, Andrew Thomas, Nicky Best, and David Spiegelhalter. 2000. “WinBUGS-a Bayesian Modelling Framework: Concepts, Structure, and Extensibility.” <em>Statistics and Computing</em> 10 (4). <a rel="nofollow" href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.471.604&rep=rep1&type=pdf" target="_blank">http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.471.604&amp;rep=rep1&amp;type=pdf</a>; Springer: 325–37.</p>
</div>
<div>
<p>Morris, Jake. 2016. “Hierarchical Compartmental Models for Loss Reserving.” In. Casualty Actuarial Society Summer E-Forum; <a rel="nofollow" href="https://www.casact.org/pubs/forum/16sforum/Morris.pdf" target="_blank">https://www.casact.org/pubs/forum/16sforum/Morris.pdf</a>.</p>
</div>
<div>
<p>Tornoe, Christoffer W. 2012. <em>NlmeODE: Non-Linear Mixed-Effects Modelling in Nlme Using Differential Equations</em>. <a rel="nofollow" href="https://cran.r-project.org/package=nlmeODE" target="_blank">https://CRAN.R-project.org/package=nlmeODE</a>.</p>
</div>
</div>
</div>
	<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://magesblog.com/post/2018-01-30-pkpd-reserving-models/"> R on mages' blog</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/dWlvd1RhyI0" height="1"></a></div>