File: glmer.html

package info (click to toggle)
r-cran-rstanarm 2.21.1-1
links: PTS, VCS
area: main
in suites: bullseye
size: 7,964 kB
sloc: cpp: 47; sh: 18; makefile: 2
file content (544 lines) | stat: -rw-r--r-- 142,154 bytes
<!DOCTYPE html>

<html>

<head>

<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />

<meta name="viewport" content="width=device-width, initial-scale=1" />

<meta name="author" content="Jonah Gabry and Ben Goodrich" />

<meta name="date" content="2020-07-19" />

<title>Estimating Generalized (Non-)Linear Models with Group-Specific Terms with rstanarm</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css" data-origin="pandoc">
a.sourceLine { display: inline-block; line-height: 1.25; }
a.sourceLine { pointer-events: none; color: inherit; text-decoration: inherit; }
a.sourceLine:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode { white-space: pre; position: relative; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
code.sourceCode { white-space: pre-wrap; }
a.sourceLine { text-indent: -1em; padding-left: 1em; }
}
pre.numberSource a.sourceLine
  { position: relative; left: -4em; }
pre.numberSource a.sourceLine::before
  { content: attr(title);
    position: relative; left: -1em; text-align: right; vertical-align: baseline;
    border: none; pointer-events: all; display: inline-block;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {  }
@media screen {
a.sourceLine::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */

</style>
<script>
// apply pandoc div.sourceCode style to pre.sourceCode instead
(function() {
  var sheets = document.styleSheets;
  for (var i = 0; i < sheets.length; i++) {
    if (sheets[i].ownerNode.dataset["origin"] !== "pandoc") continue;
    try { var rules = sheets[i].cssRules; } catch (e) { continue; }
    for (var j = 0; j < rules.length; j++) {
      var rule = rules[j];
      // check if there is a div.sourceCode rule
      if (rule.type !== rule.STYLE_RULE || rule.selectorText !== "div.sourceCode") continue;
      var style = rule.style.cssText;
      // check if color or background-color is set
      if (rule.style.color === '' && rule.style.backgroundColor === '') continue;
      // replace div.sourceCode by a pre.sourceCode rule
      sheets[i].deleteRule(j);
      sheets[i].insertRule('pre.sourceCode{' + style + '}', j);
    }
  }
})();
</script>



<style type="text/css">body {
background-color: #fff;
margin: 1em auto;
max-width: 700px;
overflow: visible;
padding-left: 2em;
padding-right: 2em;
font-family: "Open Sans", "Helvetica Neue", Helvetica, Arial, sans-serif;
font-size: 14px;
line-height: 1.35;
}
#TOC {
clear: both;
margin: 0 0 10px 10px;
padding: 4px;
width: 400px;
border: 1px solid #CCCCCC;
border-radius: 5px;
background-color: #f6f6f6;
font-size: 13px;
line-height: 1.3;
}
#TOC .toctitle {
font-weight: bold;
font-size: 15px;
margin-left: 5px;
}
#TOC ul {
padding-left: 40px;
margin-left: -1.5em;
margin-top: 5px;
margin-bottom: 5px;
}
#TOC ul ul {
margin-left: -2em;
}
#TOC li {
line-height: 16px;
}
table {
margin: 1em auto;
border-width: 1px;
border-color: #DDDDDD;
border-style: outset;
border-collapse: collapse;
}
table th {
border-width: 2px;
padding: 5px;
border-style: inset;
}
table td {
border-width: 1px;
border-style: inset;
line-height: 18px;
padding: 5px 5px;
}
table, table th, table td {
border-left-style: none;
border-right-style: none;
}
table thead, table tr.even {
background-color: #f7f7f7;
}
p {
margin: 0.5em 0;
}
blockquote {
background-color: #f6f6f6;
padding: 0.25em 0.75em;
}
hr {
border-style: solid;
border: none;
border-top: 1px solid #777;
margin: 28px 0;
}
dl {
margin-left: 0;
}
dl dd {
margin-bottom: 13px;
margin-left: 13px;
}
dl dt {
font-weight: bold;
}
ul {
margin-top: 0;
}
ul li {
list-style: circle outside;
}
ul ul {
margin-bottom: 0;
}
pre, code {
background-color: #f7f7f7;
border-radius: 3px;
color: #333;
white-space: pre-wrap; 
}
pre {
border-radius: 3px;
margin: 5px 0px 10px 0px;
padding: 10px;
}
pre:not([class]) {
background-color: #f7f7f7;
}
code {
font-family: Consolas, Monaco, 'Courier New', monospace;
font-size: 85%;
}
p > code, li > code {
padding: 2px 0px;
}
div.figure {
text-align: center;
}
img {
background-color: #FFFFFF;
padding: 2px;
border: 1px solid #DDDDDD;
border-radius: 3px;
border: 1px solid #CCCCCC;
margin: 0 5px;
}
h1 {
margin-top: 0;
font-size: 35px;
line-height: 40px;
}
h2 {
border-bottom: 4px solid #f7f7f7;
padding-top: 10px;
padding-bottom: 2px;
font-size: 145%;
}
h3 {
border-bottom: 2px solid #f7f7f7;
padding-top: 10px;
font-size: 120%;
}
h4 {
border-bottom: 1px solid #f7f7f7;
margin-left: 8px;
font-size: 105%;
}
h5, h6 {
border-bottom: 1px solid #ccc;
font-size: 105%;
}
a {
color: #0033dd;
text-decoration: none;
}
a:hover {
color: #6666ff; }
a:visited {
color: #800080; }
a:visited:hover {
color: #BB00BB; }
a[href^="http:"] {
text-decoration: underline; }
a[href^="https:"] {
text-decoration: underline; }

code > span.kw { color: #555; font-weight: bold; } 
code > span.dt { color: #902000; } 
code > span.dv { color: #40a070; } 
code > span.bn { color: #d14; } 
code > span.fl { color: #d14; } 
code > span.ch { color: #d14; } 
code > span.st { color: #d14; } 
code > span.co { color: #888888; font-style: italic; } 
code > span.ot { color: #007020; } 
code > span.al { color: #ff0000; font-weight: bold; } 
code > span.fu { color: #900; font-weight: bold; } 
code > span.er { color: #a61717; background-color: #e3d2d2; } 
</style>




</head>

<body>




<h1 class="title toc-ignore">Estimating Generalized (Non-)Linear Models with Group-Specific Terms with rstanarm</h1>
<h4 class="author">Jonah Gabry and Ben Goodrich</h4>
<h4 class="date">2020-07-19</h4>


<div id="TOC">
<ul>
<li><a href="#introduction">Introduction</a></li>
<li><a href="#glms-with-group-specific-terms">GLMs with group-specific terms</a></li>
<li><a href="#priors-on-covariance-matrices">Priors on covariance matrices</a><ul>
<li><a href="#overview">Overview</a></li>
<li><a href="#details">Details</a></li>
</ul></li>
<li><a href="#comparison-with-lme4">Comparison with <strong>lme4</strong></a><ul>
<li><a href="#advantage-better-uncertainty-estimates">Advantage: better uncertainty estimates</a></li>
<li><a href="#advantage-incorporate-prior-information">Advantage: incorporate prior information</a></li>
<li><a href="#disadvantage-speed">Disadvantage: speed</a></li>
</ul></li>
<li><a href="#relationship-to-glmer">Relationship to <code>glmer</code></a></li>
<li><a href="#relationship-to-gamm4">Relationship to <code>gamm4</code></a></li>
<li><a href="#relationship-to-nlmer">Relationship to <code>nlmer</code></a></li>
<li><a href="#conclusion">Conclusion</a></li>
</ul>
</div>

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{stan_glmer: GLMs with Group-Specific Terms}
-->
<div class="sourceCode" id="cb1"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb1-1" title="1"><span class="kw">library</span>(ggplot2)</a>
<a class="sourceLine" id="cb1-2" title="2"><span class="kw">library</span>(bayesplot)</a>
<a class="sourceLine" id="cb1-3" title="3"><span class="kw">theme_set</span>(bayesplot<span class="op">::</span><span class="kw">theme_default</span>())</a></code></pre></div>
<div id="introduction" class="section level1">
<h1>Introduction</h1>
<p>This vignette explains how to use the <code>stan_lmer</code>, <code>stan_glmer</code>, <code>stan_nlmer</code>, and <code>stan_gamm4</code> functions in the <strong>rstanarm</strong> package to estimate linear and generalized (non-)linear models with parameters that may vary across groups. Before continuing, we recommend reading the vignettes (navigate up one level) for the various ways to use the <code>stan_glm</code> function. The <em>Hierarchical Partial Pooling</em> vignette also has examples of both <code>stan_glm</code> and <code>stan_glmer</code>.</p>
</div>
<div id="glms-with-group-specific-terms" class="section level1">
<h1>GLMs with group-specific terms</h1>
<p>Models with this structure are refered to by many names: multilevel models, (generalized) linear mixed (effects) models (GLMM), hierarchical (generalized) linear models, etc. In the simplest case, the model for an outcome can be written as <span class="math display">\[\mathbf{y} = \alpha + \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \mathbf{b} + \boldsymbol{\epsilon},\]</span> where <span class="math inline">\(\mathbf{X}\)</span> is a matrix predictors that is analogous to that in Generalized Linear Models and <span class="math inline">\(\mathbf{Z}\)</span> is a matrix that encodes deviations in the predictors across specified groups.</p>
<p>The terminology for the unknowns in the model is diverse. To frequentists, the error term consists of <span class="math inline">\(\mathbf{Z}\mathbf{b} + \boldsymbol{\epsilon}\)</span> and the observations within each group are <em>not</em> independent conditional on <span class="math inline">\(\mathbf{X}\)</span> alone. Since, <span class="math inline">\(\mathbf{b}\)</span> is considered part of the random error term, frequentists allow themselves to make distributional assumptions about <span class="math inline">\(\mathbf{b}\)</span>, invariably that it is distributed multivariate normal with mean vector zero and structured covariance matrix <span class="math inline">\(\boldsymbol{\Sigma}\)</span>. If <span class="math inline">\(\epsilon_i\)</span> is also distributed (univariate) normal with mean zero and standard deviation <span class="math inline">\(\sigma\)</span>, then <span class="math inline">\(\mathbf{b}\)</span> can be integrated out, which implies <span class="math display">\[\mathbf{y} \thicksim \mathcal{N}\left(\alpha + \mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}+\mathbf{Z}^\top \boldsymbol{\Sigma} \mathbf{Z} \right),\]</span> and it is possible to maximize this likelihood function by choosing proposals for the parameters <span class="math inline">\(\alpha\)</span>, <span class="math inline">\(\boldsymbol{\beta}\)</span>, and (the free elements of) <span class="math inline">\(\boldsymbol{\Sigma}\)</span>.</p>
<p>Consequently, frequentists refer to <span class="math inline">\(\mathbf{b}\)</span> as the <em>random effects</em> because they capture the random deviation in the effects of predictors from one group to the next. In contradistinction, <span class="math inline">\(\alpha\)</span> and <span class="math inline">\(\boldsymbol{\beta}\)</span> are referred to as <em>fixed effects</em> because they are the same for all groups. Moreover, <span class="math inline">\(\alpha\)</span> and <span class="math inline">\(\boldsymbol{\beta}\)</span> persist in the model in hypothetical replications of the analysis that draw the members of the groups afresh every time, whereas <span class="math inline">\(\mathbf{b}\)</span> would differ from one replication to the next. Consequently, <span class="math inline">\(\mathbf{b}\)</span> is not a “parameter” to be estimated because parameters are unknown constants that are fixed in repeated sampling.</p>
<p>Bayesians condition on the data in-hand without reference to repeated sampling and describe their <em>beliefs</em> about the unknowns with prior distributions before observing the data. Thus, the likelihood in a simple hierarchical model in <strong>rstarnarm</strong> is <span class="math display">\[\mathbf{y} \thicksim \mathcal{N}\left(\alpha + \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{b}, \sigma^2 \mathbf{I}\right)\]</span> and the observations are independent conditional on <span class="math inline">\(\mathbf{X}\)</span> and <span class="math inline">\(\mathbf{Z}\)</span>. In this formulation, there are</p>
<ul>
<li>intercept(s) and coefficients that are <em>common across groups</em></li>
<li>deviations in the intercept(s) and / or coefficients that <em>vary across groups</em></li>
</ul>
<p>Bayesians are compelled to state their prior beliefs about all unknowns and the usual assumption (which is maintained in <strong>rstanarm</strong>) is that <span class="math inline">\(\mathbf{b} \thicksim \mathcal{N}\left(\mathbf{0},\boldsymbol{\Sigma}\right),\)</span> but it is then necessary to state prior beliefs about <span class="math inline">\(\boldsymbol{\Sigma}\)</span>, in addition to <span class="math inline">\(\alpha\)</span>, <span class="math inline">\(\boldsymbol{\beta}\)</span>, and <span class="math inline">\(\sigma\)</span>.</p>
<p>One of the many challenges of fitting models to data comprising multiple groupings is confronting the tradeoff between validity and precision. An analysis that disregards between-group heterogeneity can yield parameter estimates that are wrong if there is between-group heterogeneity but would be relatively precise if there actually were no between-group heterogeneity. Group-by-group analyses, on the other hand, are valid but produces estimates that are relatively imprecise. While complete pooling or no pooling of data across groups is sometimes called for, models that ignore the grouping structures in the data tend to underfit or overfit (Gelman et al.,2013). Hierarchical modeling provides a compromise by allowing parameters to vary by group at lower levels of the hierarchy while estimating common parameters at higher levels. Inference for each group-level parameter is informed not only by the group-specific information contained in the data but also by the data for other groups as well. This is commonly referred to as <em>borrowing strength</em> or <em>shrinkage</em>.</p>
<p>In <strong>rstanarm</strong>, these models can be estimated using the <code>stan_lmer</code> and <code>stan_glmer</code> functions, which are similar in syntax to the <code>lmer</code> and <code>glmer</code> functions in the <strong>lme4</strong> package. However, rather than performing (restricted) maximum likelihood (RE)ML estimation, Bayesian estimation is performed via MCMC. The Bayesian model adds independent prior distributions on the regression coefficients (in the same way as <code>stan_glm</code>) as well as priors on the terms of a decomposition of the covariance matrices of the group-specific parameters. These priors are discussed in greater detail below.</p>
</div>
<div id="priors-on-covariance-matrices" class="section level1">
<h1>Priors on covariance matrices</h1>
<p>In this section we discuss a flexible family of prior distributions for the unknown covariance matrices of the group-specific coefficients.</p>
<div id="overview" class="section level3">
<h3>Overview</h3>
<p>For each group, we assume the vector of varying slopes and intercepts is a zero-mean random vector following a multivariate Gaussian distribution with an unknown covariance matrix to be estimated. Unfortunately, expressing prior information about a covariance matrix is not intuitive and can also be computationally challenging. When the covariance matrix is not <span class="math inline">\(1\times 1\)</span>, it is often both much more intuitive and efficient to work instead with the <strong>correlation</strong> matrix and variances. When the covariance matrix is <span class="math inline">\(1\times 1\)</span>, we still denote it as <span class="math inline">\(\boldsymbol{\Sigma}\)</span> but most of the details in this section do not apply.</p>
<p>The variances are in turn decomposed into the product of a simplex vector (probability vector) and the trace of the implied covariance matrix, which is defined as the sum of its diagonal elements. Finally, this trace is set equal to the product of the order of the matrix and the square of a scale parameter. This implied prior on a covariance matrix is represented by the <code>decov</code> (short for decomposition of covariance) function in <strong>rstanarm</strong>.</p>
</div>
<div id="details" class="section level3">
<h3>Details</h3>
<p>Using the decomposition described above, the prior used for a correlation matrix <span class="math inline">\(\Omega\)</span> is called the LKJ distribution and has a probability density function proportional to the determinant of the correlation matrix raised to a power of <span class="math inline">\(\zeta\)</span> minus one:</p>
<p><span class="math display">\[ f(\Omega | \zeta) \propto \text{det}(\Omega)^{\zeta - 1}, \quad \zeta &gt; 0. \]</span></p>
<p>The shape of this prior depends on the value of the regularization parameter, <span class="math inline">\(\zeta\)</span> in the following ways:</p>
<ul>
<li>If <span class="math inline">\(\zeta = 1\)</span> (the default), then the LKJ prior is jointly uniform over all correlation matrices of the same dimension as <span class="math inline">\(\Omega\)</span>.</li>
<li>If <span class="math inline">\(\zeta &gt; 1\)</span>, then the mode of the distribution is the identity matrix. The larger the value of <span class="math inline">\(\zeta\)</span> the more sharply peaked the density is at the identity matrix.</li>
<li>If <span class="math inline">\(0 &lt; \zeta &lt; 1\)</span>, then the density has a trough at the identity matrix.</li>
</ul>
<p>The <span class="math inline">\(J \times J\)</span> covariance matrix <span class="math inline">\(\Sigma\)</span> of a random vector <span class="math inline">\(\boldsymbol{\theta} = (\theta_1, \dots, \theta_J)\)</span> has diagonal entries <span class="math inline">\({\Sigma}_{jj} = \sigma^2_j = \text{var}(\theta_j)\)</span>. Therefore, the trace of the covariance matrix is equal to the sum of the variances. We set the trace equal to the product of the order of the covariance matrix and the square of a positive scale parameter <span class="math inline">\(\tau\)</span>:</p>
<p><span class="math display">\[\text{tr}(\Sigma) = \sum_{j=1}^{J} \Sigma_{jj} = J\tau^2.\]</span></p>
<p>The vector of variances is set equal to the product of a simplex vector <span class="math inline">\(\boldsymbol{\pi}\)</span> — which is non-negative and sums to 1 — and the scalar trace: <span class="math inline">\(J \tau^2 \boldsymbol{\pi}\)</span>. Each element <span class="math inline">\(\pi_j\)</span> of <span class="math inline">\(\boldsymbol{\pi}\)</span> then represents the proportion of the trace (total variance) attributable to the corresponding variable <span class="math inline">\(\theta_j\)</span>.</p>
<p>For the simplex vector <span class="math inline">\(\boldsymbol{\pi}\)</span> we use a symmetric Dirichlet prior, which has a single <em>concentration</em> parameter <span class="math inline">\(\gamma &gt; 0\)</span>:</p>
<ul>
<li>If <span class="math inline">\(\gamma = 1\)</span> (the default), then the prior is jointly uniform over the space of simplex vectors with <span class="math inline">\(J\)</span> elements.</li>
<li>If <span class="math inline">\(\gamma &gt; 1\)</span>, then the prior mode corresponds to all variables having the same (proportion of total) variance, which can be used to ensure that the posterior variances are not zero. As the concentration parameter approaches infinity, this mode becomes more pronounced.</li>
<li>If <span class="math inline">\(0 &lt; \gamma &lt; 1\)</span>, then the variances are more polarized.</li>
</ul>
<p>If all the elements of <span class="math inline">\(\boldsymbol{\theta}\)</span> were multiplied by the same number <span class="math inline">\(k\)</span>, the trace of their covariance matrix would increase by a factor of <span class="math inline">\(k^2\)</span>. For this reason, it is sensible to use a scale-invariant prior for <span class="math inline">\(\tau\)</span>. We choose a Gamma distribution, with shape and scale parameters both set to <span class="math inline">\(1\)</span> by default, implying a unit-exponential distribution. Users can set the shape hyperparameter to some value greater than one to ensure that the posterior trace is not zero. In the case where <span class="math inline">\(\boldsymbol{\Sigma}\)</span> is <span class="math inline">\(1\times 1\)</span>, this shape parameter is the cross-group standard deviation in the parameters and its square is the variance.</p>
</div>
</div>
<div id="comparison-with-lme4" class="section level1">
<h1>Comparison with <strong>lme4</strong></h1>
<p>There are several advantages to estimating these models using <strong>rstanarm</strong> rather than the <strong>lme4</strong> package. There are also a few drawbacks. In this section we briefly discuss what we find to be the two most important advantages as well as an important disadvantage.</p>
<div id="advantage-better-uncertainty-estimates" class="section level3">
<h3>Advantage: better uncertainty estimates</h3>
<p>While <strong>lme4</strong> uses (restricted) maximum likelihood (RE)ML estimation, <strong>rstanarm</strong> enables full Bayesian inference via MCMC to be performed. It is well known that (RE)ML tends to underestimate uncertainties because it relies on point estimates of hyperparameters. Full Bayes, on the other hand, propagates the uncertainty in the hyperparameters throughout all levels of the model and provides more appropriate estimates of uncertainty for models that consist of a mix of common and group-specific parameters.</p>
</div>
<div id="advantage-incorporate-prior-information" class="section level3">
<h3>Advantage: incorporate prior information</h3>
<p>The <code>stan_glmer</code> and <code>stan_lmer</code> functions allow the user to specify prior distributions over the regression coefficients as well as any unknown covariance matrices. There are various reasons to specify priors, from helping to stabilize computation to incorporating important information into an analysis that does not enter through the data.</p>
</div>
<div id="disadvantage-speed" class="section level3">
<h3>Disadvantage: speed</h3>
<p>The benefits of full Bayesian inference (via MCMC) come with a cost. Fitting models with (RE)ML will tend to be much faster than fitting a similar model using MCMC. Speed comparable to <strong>lme4</strong> can be obtained with <strong>rstanarm</strong> using approximate Bayesian inference via the mean-field and full-rank variational algorithms (see <code>help(&quot;rstanarm-package&quot;, &quot;rstanarm&quot;)</code> for details). These algorithms can be useful to narrow the set of candidate models in large problems, but MCMC should always be used for final statistical inference.</p>
</div>
</div>
<div id="relationship-to-glmer" class="section level1">
<h1>Relationship to <code>glmer</code></h1>
<p>In the <strong>lme4</strong> package, there is a fundamental distinction between the way that Linear Mixed Models and Generalized Linear Mixed Models are estimated. In Linear Mixed Models, <span class="math inline">\(\mathbf{b}\)</span> can be integrated out analytically, leaving a likelihood function that can be maximized over proposals for the parameters. To estimate a Linear Mixed Model, one can call the <code>lmer</code> function.</p>
<p>Generalized Linear Mixed Models are appropriate when the conditional mean of the outcome is determined by an inverse link function, <span class="math inline">\(\boldsymbol{\mu} = g\left(\alpha + \mathbf{X} \boldsymbol{\beta} + \mathbf{Z}\mathbf{b}\right)\)</span>. If <span class="math inline">\(g\left(\cdot\right)\)</span> is not the identity function, then it is not possible to integrate out <span class="math inline">\(\mathbf{b}\)</span> analytically and numerical integration must be used. To estimate a Generalized Linear Mixed Model, one can call the <code>glmer</code> function and specify the <code>family</code> argument.</p>
<p>In the <strong>rstanarm</strong> package, there is no such fundamental distinction; in fact <code>stan_lmer</code> simply calls <code>stan_glmer</code> with <code>family = gaussian(link = &quot;identity&quot;)</code>. Bayesians do not (have to) integrate <span class="math inline">\(\mathbf{b}\)</span> out of the likelihood and if <span class="math inline">\(\mathbf{b}\)</span> is not of interest, then the margins of its posterior distribution can simply be ignored.</p>
</div>
<div id="relationship-to-gamm4" class="section level1">
<h1>Relationship to <code>gamm4</code></h1>
<p>The <strong>rstanarm</strong> package includes a <code>stan_gamm4</code> function that is similar to the <code>gamm4</code> function in the <strong>gamm4</strong> package, which is in turn similar to the <code>gamm</code> function in the <strong>mgcv</strong> package. The substring <code>gamm</code> stands for Generalized Additive Mixed Models, which differ from Generalized Additive Models (GAMs) due to the presence of group-specific terms that can be specified with the syntax of <strong>lme4</strong>. Both GAMs and GAMMs include nonlinear functions of (non-categorical) predictors called “smooths”. In the example below, so-called “thin-plate splines” are used to model counts of roaches where we might fear that the number of roaches in the current period is an exponentially increasing function of the number of roaches in the previous period. Unlike <code>stan_glmer</code>, in <code>stan_gamm4</code> it is necessary to specify group-specific terms as a one-sided formula that is passed to the <code>random</code> argument as in the <code>lme</code> function in the <strong>nlme</strong> package.</p>
<div class="sourceCode" id="cb2"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb2-1" title="1"><span class="kw">library</span>(rstanarm)</a>
<a class="sourceLine" id="cb2-2" title="2"><span class="kw">data</span>(roaches)</a>
<a class="sourceLine" id="cb2-3" title="3">roaches<span class="op">$</span>roach1 &lt;-<span class="st"> </span>roaches<span class="op">$</span>roach1 <span class="op">/</span><span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb2-4" title="4">roaches<span class="op">$</span>log_exposure2 &lt;-<span class="st"> </span><span class="kw">log</span>(roaches<span class="op">$</span>exposure2)</a>
<a class="sourceLine" id="cb2-5" title="5">post &lt;-<span class="st"> </span><span class="kw">stan_gamm4</span>(</a>
<a class="sourceLine" id="cb2-6" title="6">  y <span class="op">~</span><span class="st"> </span><span class="kw">s</span>(roach1) <span class="op">+</span><span class="st"> </span>treatment <span class="op">+</span><span class="st"> </span>log_exposure2,</a>
<a class="sourceLine" id="cb2-7" title="7">  <span class="dt">random =</span> <span class="op">~</span>(<span class="dv">1</span> <span class="op">|</span><span class="st"> </span>senior),</a>
<a class="sourceLine" id="cb2-8" title="8">  <span class="dt">data =</span> roaches, </a>
<a class="sourceLine" id="cb2-9" title="9">  <span class="dt">family =</span> neg_binomial_<span class="dv">2</span>, </a>
<a class="sourceLine" id="cb2-10" title="10">  <span class="dt">QR =</span> <span class="ot">TRUE</span>,</a>
<a class="sourceLine" id="cb2-11" title="11">  <span class="dt">cores =</span> <span class="dv">2</span>,</a>
<a class="sourceLine" id="cb2-12" title="12">  <span class="dt">chains =</span> <span class="dv">2</span>, </a>
<a class="sourceLine" id="cb2-13" title="13">  <span class="dt">adapt_delta =</span> <span class="fl">0.99</span>,</a>
<a class="sourceLine" id="cb2-14" title="14">  <span class="dt">seed =</span> <span class="dv">12345</span></a>
<a class="sourceLine" id="cb2-15" title="15">)</a></code></pre></div>
<div class="sourceCode" id="cb3"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb3-1" title="1"><span class="kw">plot_nonlinear</span>(post)</a></code></pre></div>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAu4AAAHPCAIAAACDZ8iTAAAACXBIWXMAABcRAAAXEQHKJvM/AAAgAElEQVR4nOzdd3gc5bUw8Hf6bG/qlovkKhv3isENGwPGGEyvTugJubnhpt0QciE3CSQkge9CAoQQSkzoxWATMMUG94abbFmyLEtW3ZW0vc7stO+PldayJGxJ3tXsrM7v4eGRRmvtsbW7Ovu+5z0HUxQFAQAAAABoE652AAAAAAAAAwepDAAAAAA0DFIZAAAAAGgYpDIAAAAA0DBIZQAAAACgYZDKAAAAAEDDIJUBAAAAgIZBKgMAAAAADYNUBgAAAAAaBqkMAAAAADQMUhkAAAAAaNgQSmUikcgvfvGLV199Ve1AAAAAAJAy2NAZJ+l2u3Nzc1euXLlhwwa1YwEAAABAagyhVRkAAAAAZB9IZQAAAACgYZDKAAAAAEDDIJUBAAAAgIZBKgMAAAAADYNUBgAAAAAaBqkMAAAAADQMUhkAAAAAaBikMgAAAADQMEhlAAAAAKBhkMoAAAAAQMMglQEAAACAhkEqAwAAAAANg1QGAAAAABpGqh0AAAAA0D+KorhDfGsg2hqIxeJi4iKOYQSOnXE7DFFE93fsJI5jZ96KwLHkJQzDcs1saZ45TZGDdIBUBgAAgDaIstIeiLkC0bZATJDkbl+VFUWWlG4XBbH7zc6pPRijCbzYYRx4oGBwQSoDAAAgo/GC1BqIuQJRd4iT5e7JSjpUNPtyzTqGIgbhvsD5g1QGAABAJorwossfbQ1EfWG+1/wFxzE9TTIUkdgckhUFoTO2jhSkKN3+pKLI3S8hBXXcTlFQhBcRQoIoH230zizNTdXfBaQVpDIAAAAyCC9ITd5Ikycc4oReb0ARuI4mdAxJkwTW6y3Og4K4KC8ihJz+qNMfLbTqU30PIPUglQEAAKA+RUHtwViDJ9waiPVYS0EYQhRF6ClCx5A9K3lTyG5gOEFKbGMdbfTmGFmKhKO+mQ5SGQAAAGriBKnZG6lvD0U7zyKdhiGGJPQ0qWfI7qeT0gPHMZuB8YQ4hBAvSBVN3mmjcgbhfsH5gFQGAACACmRFcfmjDe6wJ8R1W4TBEGJpQs9QeprAsMHIYLoyMGSUJ2JxCSHU5I0MsxtyzbpBjgH0C6QyAAAABlWYExo94UZPJC5K3b5E4JiBIY1sOspg+sFuZJ2+aKJA+HC9Z/HEIjKdu1rgPEEqAwAAYDBIstIWiNW7Q+4Q1/1rGGJJwshSeppEauYwHQgcsxoYb5hDCHGCVNXiv2C4Xe2gwLeCVAYAAEDaNXrCFU0+sUdfO4rADSxlGKxSmL4zsmSUJzhBQgjVt4eKbHq7kVU7KNA7WDEDAACQRoqCjjX5Dtd7uuYxGIYMDJlv0Rfa9GYdlWl5TILdyCYKdRSEDtd7pEHpzgcGAFZlAAAApIsoyQdOudsCseQVmsSNLGVgyMGv5+0vksAsetofiSOEIrxY7fSXDbOpHRToBaQyAAAA0iIaF/fVtCU73eEY5jCxOlpL0wDMLB3jJV6UEEK1bcFCm8Gqp9UOCnQHG0wAAABSzxfhd1S5knkMgWN5Fp228hiEEMKQ3cR0bDMp6HC9u+fcA6A6SGUAAACkWLM3sutEK9951pohiQKrntZm21yKwC26jpWYUEyocQXVjQf0pMkHFgAAgIxV7fQfPOVOjrA2MGSeRZeZhb19ZNbRyTysxhUIxXofDgXUAqkMAACA1BBl5ZuT7dXOQPKKRU87TGzGF/ieC4YcJjYxdltWlMMNHthlyiiQygAAAEgBTpB2VbtcgWjiUxxDOWbWki1FshSBm3UdB2X8Eb6uDbaZMgikMgAAAM6XL8Jvq3IGovHEpwSO5Vn0ejqrDslaDExyKPdxpz/C9xh+CVQCqQwAAIDz0uyL7DrRygvZUOR7FhhCdhOT2CuTZKW83qNyQKBTtj3UAAAADKZqp/9g3ekiX732i3zPgiEJI0slPvaEuXp3SN14QAKkMgAAAAai1yLfnCwo8j0rq4FJTu2ubPbH4rDNpD5IZQAAAPRbL0W+puwp8j0LDEPJuZKiJB9p9KobD0CQygAAAOiv3ot8mawq8j0Lljq9zdQWiDV7I+rGAyCVAQAA0A8tvsjuIVDke3Y2A52sB6po8ib/NYAqhtaDDwAAwPmodvoP1LmloVHkexYYhiW3meKiXNEE20xqglQGAADAuYmy8k3tkCvyPQsdTSQb57T4oi5/VN14hjJIZQAAAJxDR5Fv529rDEO5WdTJd8DsRgbvXJE60ugVJFndeIYsSGUAAACcTc8i33yLXpddnXwHBscxm4FJfMwL0rEmn7rxDFmQygAAAPhWUOR7dgaG1NFE4uNGT7g9GFM3nqEJHo4AAAB616PIlxqaRb5nZzeyeGfFUHmDV4RtpkEHqQwAAIDupN6LfJkhW+R7FgSOWTvLhmJx8bjTr248QxCkMgAAAM7ACdJOKPLtD6OOYqmObaZTbSFvmFc3nqEGUhkAAACnQZHvwNg7l6wUhMobPMn5mmAQQCoDAACgAxT5DhiJ4xZ9x2mmMCdUuwJnvz1IIXiAAgAAQApCVS1nFPkahmon3wEzs1Qy7TvZGkiubIF0g1QGAACGOklW9te217jOKPJ1DOFOvgOEIYeRxRCGEFIUVN7gURTYZhoMkMoAAMCQFouLO45DkW9qUCRu1ncMzQ5E4ydbg+rGM0RAKgMAAEOXP8JvP+4Kxk4X+eZBke/5MespqnObqdoZCHGCuvEMBZDKAADAENXii+7qUuRLEXi+RcdAke/5wRDmMDAIQwghWVEO13tglynd4CELAABD0cnW4IG69i6dfMkCq54k4JdCCtAUYWI7tpn8Ef5Ue0jdeLIePGoBAGBokWXl4Cl3ZfPp2YdmHQVFvqll1TMk0fEPWtXii/KiuvFkN0hlAABgCEl08m32RhKfYhhymFirgYE0JrUwDNmNbOJjSVYON3jUjSe7QSoDAABDRTAW33Hc5e/sd4LjWL5Fb2CgyDctWIowdm4zeUJcgyesbjxZDFIZAAAYEpz+6I7jrli8Y6eDJvACqw46+aaVzUAny4+ONfm4zgprkFrwIAYAgOx3sjV4oLZLkS9N5lv1JA6/AtILwzCbsaNDjyjJ5bDNlB7wOAYAgGyWLPJNngg26yiHGYp8B4mOIvWdW3htgViLL6JuPFkJUhkAAMhavCDtPNF6usgXQZGvCmwGJjnK6mijlxdhmynFIJUBAIDsFIzFtx93+SN84lMcx/ItOijyHXwEjtkMHUOz46J8rMl39tuD/oJUBgAAslD3Il+SKLTqaIpQN6ohS8+QyXEQzd6IKxBVN54sA6kMAABkm25FvjqayLewBBT5qspuZPDOAqWjDV5BktWNJ5vAIxsAALKHrCiHzizyNemoHLMOgypftRH46dNMnCB17bYMzpPmUxlRFB988MEvvvhC7UAAAEBlcVHefaK16cwiXxsU+WYMA0OxndtMDe5we5BTN56soflU5vXXX6+trVU7CgAAUFkwFt9W5fSGTxf55kGRb+ZxGOnkNtORBo8ow9TsFNB2KlNRUVFTU6N2FAAAoLK2QGxndWuyyJci8QKrjoEi38xD4LjV0LHNFI2L1S1+dePJDhpOZaLR6HvvvbdmzRq1AwEAADXVtYX2nWwTO8tIdRRZYNFBJ9+MZWAphur46dS1BZMLaWDANPxYf/nll2+77TaGYdQOBAAA1CEryqF6T0WT94wiXwsLRb6ZDEPIYezotqwgVN7gkWGb6fxoNZXZuXNnfn7+mDFj1A4EAADUERfl3SfamjrnLWMIsxsZKPLVBJLALbqObaYwJ5xoDagbj9ZpsiLM5/Nt3rz5l7/85Tlv5nK5kp/6/X6WZdMcGgAADIZQTNh3si3aWRxD4JjDxLJQHKMdZh0djYtxUUYI1bgCBRa9RU+rHZRWaTKVefHFF++99178XDvBX3755RNPPNH1SklJSTrjAgCAwdAWjB2ocyeLYygSzzWzUByjMRiyGxmXn0NIURR0uN6zYEIB7AwOjPZSmU8++WTGjBn5+fnnvOWKFSvmz5+f/NTn882cOXP06NHpjA4AANKrri10rEtxDEuTuSYGfgVqEU0SFj0ViMYRQsFYvLYtNDrfrHZQmoQpisaqjX7xi18cO3as1y/de++9V1111bf9QbfbnZubu3Llyg0bNqQtOgAASBdZUY40eBs7i2MQQiYdZdMzCNIYzVIQcvmiiSEGOIYtLCs0spTaQWmP9lZl7rrrrnD49DPZ6/U+88wzq1atmjFjRnFxsYqBAQBA+sRFeX9duyfU2R8WQzY9Y9LBrz1twxByGBlXMIYUJCvK4XrP/PEFkJr2l/ZSmXHjxnX9NFHYO3LkyBkzZqgUEQAApFfPIt8cEwsd8LIDTREmlgrFBISQL8LXt4dG5ZrUDkpjoEwMAAAyWnswtqPalcxjSALPt0An36xi1TMk0bEWU9nsi/KiuvFoDqQyAACQuRrc4b1dOvmyNFlg1ZEEvHRnFQxDDmNHrxBJVsobPOrGozna22DqpqCgYP369WpHAQAAKaYoqKLJe6o9lLxiZEm7gYUi36zEUISRocK8gBByh7gmT7jYYVQ7KM2A1B4AADKOIMq7a1pP5zEYshkZuxHymGxmM9IE3vEDPtrk4wRJ3Xg0BFIZAADILBFe2FHtSh5WwnEsz6wzwRndbIdhWHKbSZTkI7DN1GeQygAAQAZpD8a2V7nCnJD4lCSwfIsOJhIMESxN6JmOwo/WQMzpj6obj1ZAKgMAAJkiUeQrdBb56miiwKqnoMh3KLEbGKKzd/ORBm9chG2mc4NnCAAAqE9R0NFGb3mDJ9mA3ciSuWYWh4kEQwyOY1Yjk/g4LkrHmnzqxqMJkMoAAIDKBFHec2aRrz1R5AtVvkOSgSGT20xN3khrIKZuPJkPUhkAAFBThBd3VLvcySJfDMsz62AQzxBnMzD46W0mT7KxEOgVpDIAAKCa9iC3vcqZLPKlCCzfCkW+ABE4ZjV0bDNxglTV4lc3ngwHqQwAAKijW5EvQ+H5FijyBR2MLKmjO5La+vaQN8yd/fZDGTxnAABgsCU6+ZY3eJTOKl8jS+ZbdDgOxTHgNJuRSTwiFIQO13skWTnXnxiiIJUBAIBBJUryvtq2urbTRb5WAxT5gl6QOG7W04mPI7xY7YRtpt5BKgMAAIMnwovbj7vaOs+k4BiWa9KZdVDkC3pn0tHJKei1bUF/NK5uPJkJUhkAABgk3jC/4/iZnXytumQ9BAA9YQjZjUziMJOioMP1blmBbabuIJUBAIDB0OAO7zrRmmzeylBEART5gj6gCNyi69hmCsWEGldQ3XgyEDyLAAAgvXhBOtLQs8iXhSJf0EdmPUV3Zr01rkCoc2EPJJBqBwAAAFnLHeLq3SGXP5ZMYhCGrHoGimNAP2EOM+v0xRBSZEU5XO+5aFwBzLRIglQGAABSzxvmyxs84TPfPWMYlmtiWBpeeEG/UQRu1pHBmIAQ8kf4uvZgaZ5Z7aAyBTyjAAAglSRZqWrxn2oLdi3OxBBiaMJmYKA4BgyYRc9E41JiiMHxFn+BRZ8c1TTEwZMKAABSxhfht1U567rkMQSOmXVUoV2fZ9ZBHgPOB4Yhh4lJbCtJsnK43qNyQBkDEjoAAEgBSVaOO/11bcHTVTEYsuoZk46E3ncgVRiSMLJUouzXE+bq3aGROSa1g1IfpDIAAHC+fBH+cP0ZlTEMSdhNsJ0EUs9ioKNxSZJlhFBlsz/foof5o/A0AwCAgeMEqbzBs7NL4zsMQ1YDk2+B7SSQFjiGOUwdQ7NFSS5vgG0mWJUBAIABkWXlZFuwxhXoOuSPpgiHERZjQHqxFGFgqQgnIITaArFmb2SY3aB2UGqCVAYAAPrN6Y9WNvmicTF5BceQWc+YWQoKY1IiLkqKrDBwcP1b2PQ0FxcTaXRFkzfHzDLk0N1mgkcJAAD0Q5gTKpp87cFY14t6hrQaaBKHxZjzEuOEw9UN+ytPfXOsvqKmWZQkhJBRz9IkoWdpnY6hKcKoY1iGoknSbGApitTRVGGu9bKLLjDpWbXDH1Q4jtmNbOJxGBflikbfjJIctYNSDaQyAADQJ6KsVDX76t2hruP8aIqw62l6yNddDhjHC4eONxyoqt9XcSqZvnQVjnIIIW8wcpZv8tS/Prt8/uQbLp1dVlqYxlgzjI4m9DSZWBps8UWKbPoCq17toNQBqQwAAJybO8SV13u67igROGY1MAboUdZ/HC8crm48cOzU3mOnjp1sFsQz0hccx8eNzJ81cZTNYohx8XhcjMT4GB+PC1IoyvGCxMfjkSgfF6VoLB7heEmSOV748KsDH351YNKYYTcsm33Z/AuGyM6UzchwfkmWFYTQ0Uavw8QOzTqtIfHDBgCAAZNkpdrpr2093fUOQ5hRR1r1NAZTcPrjQGX9nvKT+yrrK2qaeqQv2LiRBbMmjpo1sWR62Yi+7xYpirLnSO27X3yz9cBxSZIraporapr/378+W7Vo2nXLZo8otKfh75FBCByzGRhPiEMIcYJ0rMk3daRD7aBUgCldl0qzmtvtzs3NXbly5YYNG9SOBQCgDe3BWHmDN9ZlMYalCLuRJQlIYvpKlpWvv6l66cOtlbXOrtdxHBs7In/2pJKZZaOml40wG3Tncy9t3uAHmw58sOkbtz+cuIJh2NzJpTdcOnvRzHF4VpcxtQe4mNDxEJ07Jj/XPLTKhhCkMgAA0CtRko81+xrc4eQVDEMWPWPSUZDF9JEsK5v2Hnvhva9rm9qTF4vzbHMnl86ZXDp7UonVlOLaDlmWtx+seXPj7r1H65K/3XKsxqsWTrth+eyCHEtq7y5DSLLi9EVlRUEI6Why0cQiEh9aD1JIZQAAoLsWX+RYk48TTm+CsBRhNzFwRqmP4qK0/quD/9ywo7nNl7hiNelvvnzOjZfOtpoHowNKXYv7vc/3fbz1cCjKJa6QBLFkTtmNl86aUTYy+3YGwzHBG+ETH5fkmSYVZ/nOWjeQygAAwGm+CH+syefr/K2AEMIwzGZgjCxUFvZJjBPWfbV/7Yadbd5g4ordbLhh+ezbVswzDvp5aV4Qt+4//sYnuw9XNyYvjixyrFo8/dpLZlqM57WllVEUhNoCMV6QEEIYQvPHF9gMjNpBDR5IZQAAACGEOEGqavY1ec849MvSpN0IDWP6JBLj3/li32sbdvpD0cSVolzrrSvmXbd0lurniSprnR9s+ubf28s5vmO+BEORy+ZNuv3KC8ePKlA3tlQRZdnpiyZ+pRtZauGEQnzIbDNBKgMAGOpESa5pDda2BeUuIwhIHLcYaDhr3Re+YPSdz/e+8cnu5G7O8HzbmlUXX7N4OpFJZ4PDUe6zXRVvfrq7a+1OWWnhtUtnXXnxFJahVIwtJYIxwd+5oDimwDKhyKpuPIMGUhkAwNAlycqp9lBNa0AQ5eRFDEMmljLrmSHznnbgnJ7AvzbsXPfVgeRqx7gR+bevnL/i4skZe2hIUZS9R+s+2LR/875KSer4uZv07PILJ92yYl7psFx1wzsvCnIFonFRRghhGLp4fKFFT6sd02CAVAYAMBQlkpiTrYF4lyQGYcjIUBY9TUAWcy6nnJ5XP9r+ybbyZH/e6RNG3H3NgvnTxqobWN+53IF1m/av++pA8vw2jmPzpoz5j5sumVCi1a7Bgii7/DEFKQghi56+eHxh1pU49wJSGQDA0CLJSn176GRrkD+zSxtLk1Y9RQ/hmXx9VN3gemnd9k17KpL7cfOnjb37mgXTJ4xQN7CBEUTpq72V7365b/+x+sQVk559988P5NnN6gY2YIFoPBCNJz6eUGQdU5CdR9C7gm1gAMBQIStKkydS7fR3PWWNEGIowqKnWZijdC6HqxtfWrd1x6GaxHtgHMeWzCq7e/UC7a5hIIQoklg+/4Ll8y+obWp/+/O9736+LxTl/veF9X/9xW0aPbNt1lPRuJjYM612BvKtehOr+TKgs4NUBgCQ/b41iSEJs57W0ZDEnMOhqvpX1u/YdqA68SmOY0vnTLz/+sWlxVquLDlTaXHuQ3ddyceF9V8f2nW45oNNB65bNlPtoAYCQ5jDwLiCMaQgWVEO13suGlegzaysryCVAQBkM0VRmn3Raqc/yotdr0MS0xeKomw7cOKldVuP1DQlrlAksfzCC+69dlG2jjf66ZrL9x2pc3oCT/1r4+wLSkYUaPKvSVOEiaVCMQEh5I/wp9pDJXkmtYNKI0hlAADZ6duSGJokzDpKD6esz0qWle0HT/ztva+q6joGJ+lY6prFM76z6iLtFpH0hVHPPvr9q7//2GsxTnj0+Q9fevROjXZnseoZLi4KkoIQqmrx5Vt0WfyYz9q/GABgyFIQcvmiVS2+yJlJDEXiFh2dxS/oKSGI0mc7j768busppydxxaBjbrx09h1XzU/5yKTMNOeC0hsunfXO5/sOH294c+Pu21ZcqHZEA4FhyGZk2wIxhJAkK4cbPBeOzVc7qHSBpzQAIHt0JjH+SGebkwSKwC16SGLOIS5K732x77WPd7Z6OmYOOCyG266cf8Py2QZ2SLQnSfrRrct3Ha5pbPX99a1NF00fN6rQoXZEA8FShJGlwpyAEPKEuAZPeITDqHZQaQFPbABANkgkMced/sQLdxJF4GY9baBJpMldgkEiy8qnO8qfe/srp9ufuFKYY12zav41i2eoPnNAFTqW+vX3V9/7m1f4uPjIc+te/d+7Mrbj39nZDDQnSKIkI4SONfnyzLqsPKk3FB+jAIBskhikd7zFH4zFu16HJKaPdhw68cybm07UuxKfjixy3LlqwYoFk0kiC3/n9d30CSNuu2Lea//edfRE06vrd9x1zQK1IxoIDMNsRro9wCGEREk+0uCZPTpP7aBSD1IZAICGtQZi1U5/siFYAonjZj1tZCCJOYdjJ5uffuPLfRV1iU9tZv3dqxfetHxORg1OUtEPblm2s/zkyca2F977esH0sWNHanLwpI4i9QyZKH5vDcRafJEim0HtoFIMUhkAgCa5Q1xlsw+SmIFxuv3Pvf3VJ9vLE83u9Cx90/I5d127cKjVxJwdTRK/fWD1ml/9QxCl/3l23WuP30dpsxm0zcDwgiTJCkLoaKPXYWIZbf5Fvg2kMgAAjek1iSFw3AJJTB/4g5HXPt71+ie74qKEEKJIYtWiad+78RKHJdveqafEhJLCO1dd9OK6rdUNrf/4YOv3b1yidkQDQeCYzcC4QxxCKC7Kx5p800flqB1UKkEqAwDQBllRnL5obVuw95UYlkCQxZxVjBPe+nzPy+u2RWI8QgjDsGVzJ/7w1mXFeTa1Q8to912/ePvhE5W1zpc/2rZwxrhJY4apHdFA6BlSH+/YZmr2RopshnyLTu2gUgZSGQBApuMEqb491OAJ82eOHSBwzKKjDToSgyTmrERJ+ujrQy+8+1VyBPTcyaUP3rZ8/ChNFn8MMoLAf/vAtbc99AIviI88t+6NP3xPo6e67AaGi0uyoiCEjjR4HBOLyGwpitLkzwMAMEQEovG6tmCzL5oo6UgiMMykp00sld2TZc6foihf7jn27FubGlzexJWJo4f96NZlsyeVqBuYtpQW595z7aJn395U1+J+/p3ND96+XO2IBgLHMZuR9oR4hBAnSMeafVNGaLJfTk+QygAAMo4kK83eyKn2ULfz1QghisCNLGWEJKYP9hypffqNL5KTBwpzrN+/ccmVC6ZodOCzuu68+uIdh2sOVdX/65NdC2aMnzlxpNoRDYSBoaK8GItLCKEGd3iYzeAwsWoHlQJYt/c6Wcztdufm5q5cuXLDhg1qxwIA6J0kKw3ucE1roNteEoYQQxEmHQyA7JOTjW3PvPllcpC11Wy4Y+WFt624kM6ucyuD7JTTc8t/P8/HxWF5trf++H2NnvaSZNnpiyW2mfQ0uXBiEanNIVNdwaoMACAjSLJS7w6dbA12S2JwHDOylIklCW22Wx1kTk/g5Q+2rvvqgCwrCCEdS928fO5dqxcYdIzaoWneqELHD25a+tRrnzW3+f7yxhe/uOtKtSMaCALHrQbaG+YRQtG4WN3in1is+bpvSGUAACqTZaW+t5UYmsSNLGVgYC+pT/yh6GsbdiZPWZMEcfXiafffsCTHmp1jd1Rx6xXztuw/vv/YqXe/+GbhjHHzp41VO6KBMLBUhBcTT7e69mCRTW81aDvThQ0mAIBqZEVp8kSqnX7uzCSGIQmzHvaS+ipxyvqVD7eHoxzqPGX9g5uXjiiwqx1aFmpp99/0s+ciXDzPbn7nT983GzR5pFmQZJc/mvj9b2SphRMKcS1vM8GqDABABYqiNPui1U5/otFFEk3gZhhh3We9nrL+z1suLSstVDewLFaUa33w9sse+8eGNm/wqdc+//X3rlY7ooGgCNyio/3ROEIozAknWgPjC61qBzVw8HoBABhU35bEUCRu0UES0w9b91c/9a+NDc6OU9alxXn3X7/o0nmT1I1qKLh26Yyvv6nacejE+q8PLpwx7pI5ZWpHNBBmHR2Ni3FRRgjVuAIFFr1Fr8lCZgSpDABg0MiK0uSN1LgC3ZMYArcYaL02246p4nB149NvfHmoqj7xaaHDcte1C1cvmY5DZfSgwDDskftX3fDT54KR2GP/2DBtwgi7WYNjHzDkMLFOXwwhRVFQeYPn4vGFGq1Lg9cOAEDayYrS4A6fbA3G4t1XYqx6WgdJTJ/VNrW/8N7XX+yuSHxqNenvuGr+rVfMYyj4NxxUuTbTT79z+SPPrfMFo3989dM//Of1akc0EBSBW/RUYhJIIBqvbQuOzjerHdRAwKMfAJBG31bYSxG4WU8baJj+2Fcud+CldVvXfXVQlmWEEMtQt1w2985rLjbqs6HFmRatXDh1yzfHN+099vnOo0tnT7j0wgvUjmggzHo6youCJCOEjrf48y06I0upHVS/QSoDAEgLUVbq20O1rUFePDOJIXGLjtYxJOQwfRQIx9au3/HGp7t5QUQIEQR+zeLpcMo6E/zynlJBVuIAACAASURBVJWHjtd7ApHfv/zJ9LJRWvyJYAjZTUybP6YgJCtKeYPnwnEFmntuQioDAEgxUZIbPZGefWIoEjfrYCWmHzheePOzPa9+uD0U5RJXFswY9+PbLxtZlCWjc7TOZtY/fO9VP/7zW/5Q9Hcvrv+/n92qdkQDwZCEUUeFYgJCyBvm69tDo3JNagfVP5DKAABSJswJTd5IfXsosV6dRFOEBWYO9IcsK59sL//Lm1+2+0KJK1PHDf/RrcumTdDk6J8stnjWhMsvumDjjqNb91dv2HLoqkXT1I5oIKx6JsZLoiwjhKpa/PkWnbYq2NIbqyiKR48edTqdgUDAZDIVFBRMmTKForS3DwcAOAtOkJy+qNMfSXRD74qhCDMkMf2050jtk699VtPQmvi0tDj3/usXwynrjPXQ3SsPVjW0eoJ/+ufGWReUFDosakfUbxiGHCamNRBDCImSfLjeM29svtpB9UO6Upna2tqnn376s88+C4VCXa+zLLts2bIHHnhg6tSpabprAMDgECXZFYg5fZG2YKxn23BIYgbgSE3T069/caCy45R1QY7l7tVwyjrTmfTso/df/YPf/ysc5f73+Y+ef/gOLc4eZyjCyFBhXkAIuUNckydc7NBM6U/qUxlJkp566qm//OUvBEFMmTIlPz/farUajcZIJOL3+71e75YtWz7++OM1a9Y88sgjOp0mWz4DMJTJitIWiDV7I63BWGJmYVcYhnQUadJTDAxh7o/a5vYX3j19ytpi1K1ZdRGcstaKeVNGr1o07aOvD+49WvveF9/csHy22hENhNVIxwRRkhWEUEWzL9esYyhtPItT/CQJh8N33323Tqd7++23Z8yYwTC9TKgSRXH//v2vv/769ddf//rrr1utGm6WDMDQoSDkC/NOf6TZG0l0CO0KQ4imCD1NGhhS08NcBl+rJ/iPD7Z0O2X93WsuNsEpa0356Xcu33u0zun2/98bn8+dPHpEofYGYOEY5jCybcEYQkgQ5fIGz+zReWoH1SepHCcZi8Xuv//+H/zgB3Pnzu3L7Xft2vXcc8+98sorJDkYbztgnCQAAxPiBKcv0uSJRM9scJdAEbieIQ0sScImSD8FI7F/frTjjY27+biIEMJxbOmciQ/efmlhDrzB06R9FXXf+91aRVGmjhv+0q/v1Oi2oDvEJftxzyzNLbTq1Y2nL1KZQ/z5z3/+8Y9/PG1aX+u3L7zwQoPB8Nxzz/3nf/5nCsMAAKRELC66/LEmbzjRDLQbEsf1DGFgKYrQ5Ou1ugRRWr/l0HNvb/IFo4krcyeX/vj25WNHFqgbGDgfsyeV3LR89luf7T1c3fjGJ7tvXzlf7YgGwm5g+LgkKQpC6EiD12Fk6IzfLE5lKrNmzZqRI/t3UHDKlCkWi/aKvQHIYoIotwZjTZ6wO8T1/CqOY3qaMDCUVjbRM40sK5v2HnvmjS+b23yJK5PHFP/otktnlMEp62zwn7deuuvIyfoWz1/f2jR/2pjSYm1s0HSF45jVyHhCHEIoLkrHmnzTRuWoHdQ5pHKDKcPBBhMAZyHJSlsg1uQNtwW5ni8LGIZ0NGlgKB2NI+hwN1B7jtT+v9c+q+48ZV1SlPO9G5fAKessU36i8a5HX5FleUJJ4drf3UMSmkz6u24zzR6dl2/J6DM6UBsPwJCmKMgT5po8YZc/KvY4joQwxJKEgaV0NAm1vOdj/7H6597ZdLCqIfFpvsP8vRuWXLVwqkbLKcBZTBk7/I4r5/1zw86qOuerH+2459qFakc0EDYDw8UlWVEQQkcbvQ5jIZnBW8nqpDLr169ftWqVKncNAEjwRfhmb6TF18txJIQQQ+EGhtLRJAEpzHmIxPh/byt/94t9JxvbElfMBt13r774lsvnMprqpgr65Xs3XrLtYE1tU9uLH2xZOGvcuBHaK4EicMxqYLxhDiEUi4tVLf4LhmfumSx1nkuHDh2CVAYAVSgIOX3Rk62BXot5KRI30KSBJQlYLTg/Tk9g7Uc7Nmw9FOU6/p0Zirz58jnfvXqBxZjRa/Xg/DEU+dsHVq/5nxcFUXrkr+tee/w+KuMrZ3sysmSUJxIz7evbQ0U2vd2YoQ0CUpnK/PCHP/T5fOe8Gc/z+/bte+SRR1J41wCAs4vGxVZ/rC0YC0T5nsswHceRGIoiIYM5X81tvlc/2r5+yyGhcyR4cZ7t2mWzrl483WbWwLlWkBJlpYV3X7Pw7+9/Xd3Q+vf3v/7BTUvVjmgg7EbW6Y8oClIQOlzvXVhWmJnLtKlMZcaOHfvEE0+k8BsCAM6HgpA/wrcGYq2BaGLsbTc4jiWa2tEUkYmvT1rT2Op7ed3Wf28rF6WOJObiaWNvvnzOvCljoG3gEHTPtQu27K86fsr16vodi2dOmDRmmNoR9RtJYBY97Y/EEUIRXqh2+suG2dQOqhepPMEUDoeXLFny2GOP5ebmftttFEVxu92/+93vtm7dmqr77SM4wQSGCElW2oOx1kCsLRDjOxcGuiFx3KynDAylwVkxmajB5f3HB1s+3XFEkmSEEIZhi2aOv/faRWWlhWqHBtR0srHttodeiItSSVHOG3/4niYLpBTU2vlKgmHoovGFVj2tdkzdpfKf1Wg03nvvvYsXL6bpc/w9KysrU3i/AACEECdIrYFoayDmDnE9RyMhhBCGGJLQUQRN4ixNwJnqlHD7wy9+sGXdpgOJlRgcx5bMKrv3uoXjoNkdQGj08Lz7b1zylze+rGtxP/v2ph/fcZnaEfUfhuwmxuWLKggpCiqv9yyYUJBp8zJTnCHeeuutbre7qKjo7DdbvXp1au8XgCErxAltgVhrIOoN873eAMMQSxF6htLRBJ5hL0CaFo5yaz/e+fonu2KcgBDCcWzZvEn3rl44erj2uqKB9PnOyvlb9h0vP9H4xqe7F8+aoMVeiBSBW/S0PxpHCAVj8ROu4LjCzOptCy3yANAeWVHcIa41EGv1RxPnC3oiCUxHUzqaYCloapdiibEDz7+92RuMJK7MnVz64G3Lx4+ClRjQi6Y2380/fz7KxYfl2d764/cNbMZt0PSB4vLF4pKMEMIxbEFZoYml1A7pNA3u2wEwVJ19NjVCCEOIInEdTepokoazSGkgy/KGrYdfeO9rlzuQuDJl7PAf3rJs5kTtvdUGg6Y4z/aDm5f+6dVPm9t8T7/++S/vXql2RAOAOcys0xdDSJEV5XC956JxBZmzyKtOKvPkk0/+5Cc/UeWuAdAcTpA8Yc4T4loDMb63NRgMQyxF6mhCT5NwUiZNREnauOPoy+u2nnJ6ElcSYweWzZ2YaXUDIAPdfNmcbQeqd5effP/L/Ytmjr9o2li1I+o3isDNOjIYExBC/ghf1x4szTOrHVSH9G4wlZeXb9682el0ynLHO0hFUVpbW3fv3n3ixIn03W+vYIMJaAgnSJ4Q5wlznhAf4Xs5R40QInBMR5M6hmBJEn6Zpk9clDZsOfTSuq3JlZh8h/meaxdds3g6kcGt3EGmaWn33/Tz5yMxPtdmevfPD5gN2uuUqCjI6Y+KkowQInBsUVmRnsmIvZ00BrF27dqHHnoofd8fgGwiSnIgFg9E48Fo3BfhI52D3HoicEzPkAaapGE2dZpFuPg7n+99/eNdyZqYwhzrd1ZddPWS6QyVEa/gQEOKcq3/dcdlv/v7+nZf6M//3PibB7R3/AXDkMPItgWiCkKSrBxu8Fw4Nl/toBBKayrz7LPPFhQUPPDAA6NGjSI6R4PKstzQ0LB27dr03S8AmhAXpUA0HojGA7F4MBo/S+6SQBE4SxMsReooAqp4080bjLzz2b63Nu4JRmKJKyOLHN9ddfGVC6ZodNAxyASrl0z/am/ljkMnPt56+JI5ZYtnTVA7on5jKNzIUiFOQAh5Qly9OzQyx6R2UOlMZTiOe/zxx6+88sqeXzpLDz0AsgwvSGFO8IQ5T5iPixJNEhhCYU74tpNHXVEEzlIEQxEMRWRmv/Dsc8rpef3jnR9vPcwLHcnl2JEFd1198aXzJsIUa3CeMAz7n/uuuvFnzwcjscde3DBt3HCr2aB2UP1mMdDRuCTJMkKostmfb9Gzaq8QpzGVmTt3rij2/kbzsss02CYIgL4RRDkQi/sivC/C+yM9Bx71XviSgGMYRWI02Zm+QAnMINp/rP61f+/YfvBEssHgBWOL7756wcKZ46CwF6RKnt388ztX/Oqv73sCkd+//MkTD96gdkT9hmOYw0i3BTmEkCjJ5Q2eOaNV7qWUxlTm4Ycf/v3vf3/llVeSZPd7eeWVV+6999703TUAg0NRlAgvhjmh8/9CmBN6PSb9bQgcowicJhP/EbD2MvgkSd6059jaj3ceq21JXMFxbMGMcXdcOV+L3cxA5ltx8eTNe49t3lv5xe6KS3aWXTb/ArUj6jeWJg0sFeEEhFBbINbsjQyzq7m8lMpU5sknn4zH412v8Dx/5513Tpw4MXlFUZSGhoYtW7ZAKgM0R5KVMCck/wtxQoQX+3gGEMcwliJIAlMUpCCFwHGKwGkKJ2HPQj2RGL9u0/43N+51uv2JKwxNrlww9baV80cVOtSNDWS3X96z8mBVvS8YfeKVT2ZOHJVjNaodUb/Z9DQXFyVZQQhVNHlzzCxDqrbNlMpU5tChQ5s3b+55vdeLAGQ4UZIT+Ury/zFe7HvrAhzHGBJnSIImcZKElCWzON3+tz7d88HmA5FYx7QHu9lww/LZNy6fYzPr1Y0tKykK4gWRF2VJVhDqeBp1HRSmKKffFMidHyXWLCmSoAicIrNqwdJuNjx8z1U/feptfyj627+vf/rnt6odUb/hOGY3Mu1BDiEUF+WKRt+Mkhy1gkllKrNixYpQKHTTTTcxDPNtt1EUpaWl5Z133knh/QJw/uKiHIrFw7wQinUsuvSlLDcp8ZpLEhhJEBSBJz5OX7RgwI6fcv3r37s+23k0Mf0RITSiwH7jZXOuvWQmy2RQI/YsoCiIE6S4KHGCxAtyMoPplxiSOsrLMEThOEXiiQ1ZkiBIjZeSXTKnbMXFkz/ZfmTbgeqPvj549eLpakfUbzqa1NNkNC4ihFp8kdH5ZotKQ7NTmcpceumlubm5y5YtO+ctS0pKUni/AAxMXJTDXDwxzCgQjZ/7DyRgiMIxsiNfwenE+0Vtv6hmP1lWth888ebG3XuO1CYvTh0/4tYr5i6dUwZHk1JFURAvSnxc4kQpLsjKgNKXb/nWSJBkQTpdiIZhKLFmQxO4kaW0+BT877uuPFDV4HIH/vzPjXMmlxY6MmtGY1/YjAznl2RZ5WGOgzpO0ufz2Wy2Qbu7bqDb7xCXSFxCXMe6S4gTeh0C0A2GMJLAKAInSZwmMJLEKQLWW7QkLkqf7zz66kfba5vbE1dIglgye8IdV154wdhidWPLDon0hRMkXjhH+sJQhJE5nXOQXSa1kzh+xnWEEELRuBiKCRFeOOevKQLHbAYmQzrP9svu8pM/+P2/FEWZPankb79ao8U3RXFRDnHCmHzzcIdqFT/p7Svz4IMPtrW1vf/++4kfT0tLy69//euHH344L0/lg1tgKIiLcuI4tC/Ch2ICL/YhcUm8z0v8R+IkgVM4Dv3oNMobjLz7+b63P9vrD0UTVww6ZtWiaXesnF+Qo723vxlFUVBi54gTxLhwtrfkDEnYTYzDyDpM7MBmKcuKEooJIU4IxeKJD2Lx7m0+JFlxhzgmRtiMjLYGqc6bMnr1kukfbD6wr6LunS/23bR8jtoR9RtN4g4jo2Ieg9Kayjz99NMbNmwoKiqSJClxHnvSpEn333//dddd99FHH9nt9vP55jU1NWvXrq2pqTGZTPPnz7/llltoWotj00GKyYrSFoi5/NGz9/7vCscQQ5F6mmBoAopzs0N9i+edz/d+sHk/3/k7b1ie7dplM69fNsukZ9WNTbu6pC9SXJDOkr7QJOEwMg7TwNOXrnAMs+hpi55GqOO4ryDJoZgQ4uKhmOCL8MndYV6UWv1RA0tZ9bSG5qr+ZM0V+yrqGlt9T7/+xYWTx4woPK9fjkNTGjeYli1bdvPNN990000m0xldja+44orp06c//vjjA/7OBw8efP755+fNm8dx3K5duwKBwMyZMx999NGz/ynYYMoysqL4wnwgFucFSVFQlBcRhjxhTjhrW5fEuktnjQtBETj0ocsmh6rqX1m/Y/vBE8lXtrLSwlsun3fFRZNh9ONAKCguyZwgcXGRF6Wz/LogCdyqp3PNuhwTa9bTg/mcag3EKhq90S5LNRiGzDrarKO18tw+VFV/z29elWVlytjhL//vnVos3poyQs3+BWlclRk1atQ999zT8/qwYcO+/PLLAacyiqK8//77Tz31lNFoRAjdfvvtP/7xj/fv319ZWVlWVnZeEQMtiIuyOxRz+WNtwZgonbsZHYFjNEkwJCQu2UwQpc92Hl378c6ahtbEFRzHLpo29q6rL546foS6sWmPFtKXrvItulxzUX17+LjTn3hNUBQUiMajvGDRa6OAZtqEkTctn/Pmxj3lJxr/9cnuNSvnqx2RxqTxZ9zrjo8oigcPHuR5fsDf9uTJk8uXL0/kMQghs9l8zTXX/P3vf6+trYVUJltJsuINc+4Q5w5xwWj8nAuJJIHpaIohMZqCPaMsF45y67ccWrthZ5s3mLiiZ+krLpoMbe76S1ZQjBeicZETzpa+0CRuN7IOI+swMSadaulLNziGleSZCm36yiZfs69jjLkgKe4Qx/KkzUBTGb8m98Nblu0sr6lv8Tz31qaLpo4ZPRwqSvshjalMWVnZm2++ecsttySvyLL86KOPulyuJUuWDPjbFhUVlZaWdr1SWFiIEDIYtDeUC5yFoiB/lE+kL74wL3/LiyuGYRSBIQwjMAzDEIFheoZk1J5tBgZBY6vvrU93f/j1gRjXMdYqx2q8btmsmy+fazHq1I1NQxSEeEGKcEI0Ln5bBkORuMPAOkyZlb70xFLE9JKckbmmiiZvsoCGi4tOQTSxlEVP4xm8Jssy1G8eWH3nIy/HRelXz37w2mP3wgz2vktjKnPfffetWLFi06ZNixcvttvt9fX177//fmVlJUmSP/3pTwf8bfX67r0429vbcRyfOnXq+cULMkKYE9whzh3k3GHuLPtHFIHrGFJPETRkLUPPoar6Nzfu3bS3UpY7HiHjRhbcfuWFl190Abz6950gyVFejHCiKPfyRCNxzGpgVN88GgC7kVkwobDJG6ls9nU0XFBQKCZEeNGio00slbFnEiePKf7OyvmvrN9+/JTr5Q+33XfdYrUj0ow0pjIMw7z55ps/+clP/vu//zt5cfjw4Y899ti0adNSeEd79uxZsWJFz441TU1NVVVVyU9DoZDZbE7h/YIUkhXFHeSOO/1naVVH4BhLESxNsBRBwLbR0OMPRTfuOPLhVwer612JK4mCmDtWzp81cZSqoWmJLCsRXozExXhvfZUserrIZtBc+tJTsd1QYNGdcAXq2kKJNV1ZVnwRPsIJNiOTsQu399+4ZOvB6pONbf/4YNuimRPGjypQOyJtGIwWebW1tQcOHBAEYcyYMVOnTk3tqemqqqq//e1vf/zjH3t+23ffffeJJ57oeiUajY4ePRpOMGUORUHtwViTN9wWiIm9NafAMURTBEsROoqkNNUuAqSKJMk7DtWs33Jw24FqobM5EEOTKxdOve1KKIjpMwVx4rduJDEUUWTTF9uNajWeT58IL1a1+Jy+aNeLOpqwGZnMLKSrqnOu+dU/REkaMyL/X4/fR6s3o7Ff1D3BNKjdfhP+/ve/T5gwYeHChef/rTiO+/nPf/6zn/1s+PDh57wxHMbOKFFebPSEGz3h3kYdYQyFsRTJUgRD4Shjl4NBmp1sbFv/9cFPtpd7ApHkxeI821VLpl+/dBbMfewjQZIjnBDhO4YYd4VjWK6ZLbYbC6w6LfaZ7Tt3iKto8oZiQvIKhjCjjszMApoX3vv6hfe+RgjdefXFP7zl3LOAMkHWHsZGCImiWFFR4ff7kwmTKIpHjx795z//uWPHjvP//mvXrl2zZk1f8hiQIRSE2gKx2ragJ8T1/CqJ40aWNOpo7XS3Aqnn9oe/2F3x762Hj9W2JC/qWXrZvIlXLZw+o2xEdv/STZWOjSReiPfWacmip4vtxmF2g7Z64w5YjoldOKGw0ROpavEl/kEUpIRiQpQXrQbGkGEHtu9evWDrgeOVtc5/bti5ePaEyWNgwsY5pPHn19jYuHr1aqfT2fNLo0ePPv/v/+mnn44dO3bWrFnn/63AIFAUpdkXrXEFwpzQ7UskgRsYUkeTQ+RVFfQqHOU2763auPPI3qN1yXpeDMNmlI1YtWj60rkT9Wy27X2kQ2IedYQXonGx5zQkliKG2Q0jcoyGoTcGHMOwETnGQqu+2uU/1R5KvL+WZMUT4kIx3GZkmIzZyiEJ4rcPrL7tob/zgvjos+ve+MP3YGz72aUxlfnlL385evToq666yufzORwdS08ulysUCv32t789z2++bds2QRCuuOKK5BWXy+V0OqdP196c9KwXiMZd/miTN9JtcgqGkI4hjQzJ0CS8yx6yeEHccfDEp9uPbD9YzQunHyHFebYrF05duXDqsDzVZtBqS1yUI7wQ7W0jicCxPItuZI4pxzTU5zZQJD6p2D7cYTza6POGO9aG46Lc6o8ZGNJqYIjMWBMuLc773o2XPP3656ecnr++temn37lc7YgyWhpTGZvN9swzzyCE/H5/OBwuLu5YIvuv//ovhmHO5zsfOnTotddeW7Ro0euvv564Eo1Gy8vL//SnP51nzCCFYnGxti3o8sd6zn4jCdzEUnqGzJBXDTD4ZFkpr274YvexT3ce9QdPl8JYTfqlc8quXDBl6njYSOoTUZKjcSnC9TIWCUPIZmSK7cYiu4GE51oXZh09f1x+t4kHEV6MxkWzjs6Q01t3XHnhtgPHD1TWv/XZniWzJ8yEY3rfLo2pTHIlxmq1Hj58OJnKzJo16/HHH09kOQNQV1f3+OOPcxz39ttvd71+xRVXsOxQf8OROVp8kfIGb8/GMBSBm3WUgSWhmHdoUhTl6MnmjduPfLG7wu0PJ68b9ewlcyZcPn/ynAtKtDiAZvApihKJixFO4IVeSmFMLFXsMA6zG9hMPXWcCfItuhxT0cnW4MnWQGIpq2PiQVzMt+hULwfGceyR+66+5aHnY5zw67999NYfv2+APdZvkcZUprm5effu3QUFBaNGjaJp+t13373hhhs4jtu/f/+WLVsG/G1LSkreeeedFMYJUkuU5IomX6Mn3PUijmEsTehpUsfAXtJQJMtyeXXjF7uPbd5X2eoJJq/TJDF3yuhL501aOmei7rxHKA8JCuIEMbF+0PP4KU3iRTZDscNozboz1WlC4Ni4Qstwh6Gy2d+SnHggyt4Qn2NW/73xiEL7D2+59I+vfNLc5vu/1z57+N6r1I4oQ6UxlVm1atV1112HENq8efOFF154xx13PPHEE8FgMBKJjBo1Kn33C1TkCXOHT3m6jqjVM6SRoRgKpjgORXxc3HOk9ovdFVv3Hw9FT59Zw3Fs8tjhl86buOKiyVYzjBzpk7Ofqc6z6IrthrwMWEvQIh1NzijJGZVrOtroDcbiCKFoXIzwQiYUR9+0fPa2A9W7Dte8v2n/otkTLp42Vu2IMlF6+8q8/fbb27dv//3vf280Gl0u14033njy5EmTyfTcc89dcskl6bvfXkFfmbSSZeW401/bFkw+oDAMcxi1MZYWpJY3GNnyTdXX3xzfe6S2ayUvRRKzJpUsmTVh8ewJOVajihFqiKwoUV6K8ALfW3PexEZSscOQOadvNC0uSlsqnYl/agxDBVZ9JsyhdLr9N/7s+UiMz7Ob3/nT982GTBwxNoRa5CmK0tzcnJeXl9qGv30EqUz6+CP8oXpP11PWDEU4TBnaTBOknChJJxvbK2pbKk82V9S2VNe3yl3G+hhYev70cZfMnnDRtDFGvfqL9pqgIIWLyxFeiPGS0uNQNUsRhTb9cIfRrIONpBRrD8b21LQlPqZJPN+qwzKgsO/Drw785oX1CKEVF0/53X9cq3Y4vcjmFnkIoXA4vG3btqNHj1qt1gsuuGDevHlwKiGbKIpS2xaqajndBRHDkEXPmHQU/JizmCzLdc3uipMtlbUtFXUtJ065uq6+JORYjYtmjl88e8LsC0q10nw9E3RsJHGi1ON9Jo5j+RZdsd2YZ9bB62ia5Jp1JXmmurYQQiguyoGokAmFR9csmbF5b+X2gyc+2V6+dG7ZktllakeUWdKbyuzatetHP/pRc3Nz8srs2bOffPLJlLTIA6oLROOH6z2JreUEmiIcRiYTlmRBaimKUu/yHkvkLrUtx0+1xHq0OkQIGVh6QmnhlLHDF84cP3lMMQ4HgPtMkpUoL0Y4Id7bQHiLnh6ZYyqy6Ul4cqVf2TCbJ8QnXtlC0ThLEZlwEOx/7lt148+eC4Rjj7348bTxI2FwR1dp3GBqaWlZsmRJJBKZN2/euHHj8vPzvV7vvn373G73pk2bTCZTmu7328AGUwopinLCFTzhCpxejEGYRU+ZdXQGrMWC1Gj3hSprnZV1LcdqW46caPKHoj1vQxLEiEL79PEjpo4fMbG0cFRRLqQv/aIoKBYXI7wQE6SezXl1NFlk04/MMUHN2SALccL2KmeiwprAsQKrPhOaYH2648jDf3kfIbR41oSnfnqz2uGcIWs3mP70pz+VlJQ8++yz3dZg3nnnnaeeeurRRx9N312DtApxwqFT7kD09GIMReIOIwtjB7Sua+5ytKbJF/zW3KWspKistHBiSeGkMcUUbB4NSFyUw1w8yktyj/eTJIEXWHTFDiM051WLiaXKhtmONnoRQpKseMN8bgaczb7ioslf76v6YnfF199UfbrjyBUXTVY7okyRxlSmurr67bfftlgs3a7feOONDz30UPruF6TVydbg8RZ/lxdfzKKnLHoKUZ1VagAAIABJREFUWt5pUTASq6x1HjrecKy2pbK2pWvPuiSCwEcWOpK5S9noYQwFKwQDJ8pylJfCXFzsrTmvw8QWO4yFmbEGMMSNyjW1B2OtgRhCKBYXw5xgzIDWRw/dteJgVb3bH/7DS/+eMWFkvsOsdkQZIY0vScXFxT3zmAS3252++wVpwgvSoXpPezCWvAKLMRoSinLtvpA3EKlv8ZxocNU2tZ9qcfeau+A4XjIsZ9LoorLSokklRWNHFUDucv5kJVEKI/Jib2eqddRwu3GY3cBkQE0GSJo60pE8m+2L8AxFqF4IaDUbHrp75U+efCsU5R57ccMzv7hN3XgyRBpfob5tjEBVVZXH40nf/YJ08IS4A6fcpztbYMjMUhYDA+8cM4QoSd5AxO0Pe/xhtz/U7gv5AtE2X9AbiLT7Qh5/uOcJoyQMw0YW2CcmcpfSovGjiqDxbqqIsszF5Vhc4OJyzzPVFIkXWvXFdqPdeF5j6UCa0CQxbWTOnppWhJCiIE+Iy4Sz2UtmT7hywdR/bzu8/dCJD786cM2SGerGkwnSmMoUFBQ8/vjjP/vZzyiq42XR7XZ/8MEHTz/99J133pm++wUpV9cWPNbsS+4pETieY2Lg7WOvREmKcnGEkCwr4RgviXKUiwcjHUtZgiQFw7FAOMbFezn+E43x8pmbDixLkb1VooiiFI3xnkDEE4i0+0K+QNgbjPa9hF/P0qOG5Y4Znjt6eH5ZSWFZSaFBB79KU0lBKMoJIV6M99bXDsewfIuu2GHMM7PQnCLD5ZrZ0jxzbVsQJc5mR+JWg/pPlp/fecW+iro2b/DJtZ/NnVxamGNVOyKVpfEEUywWW7p0aSQSKSsrQwg1Nzc3NDSIojh27NjPPvvsPIdjDwCcYBoASVbKGzzN3tODi/U0aTcxQ7M5ujcYafMEWz1Bpyfg9oVaPcFWT8DtC3mDUYRQOMZ37QunOotR57AaHRZjrs1ktxrzrEa71ZRnNxXmWIpyrfAbNE0kWQlzQjgm9OwKgxCyGZhiu6HIblB9nwL0nawo26tcibPZGEK5Fl0mnM3edbjmP/7wuqIosyeV/O1Xa1R/RmftCSadTvfhhx8+9NBDGzduTF5cvXr1b37zm8HPY8AAROPi/tr20yeVMGTTMyZdlm89RGJ8qyfQ4g60eoJtnqDT7W/1BNt8QVd74Cx7NIOPIgm7xZBjMzkshlyrKcdmspsNuXaT3WLItZsdFiN0pRtkvCCFOCHKS+jMjSQMQ3YDm2/VFVj0cKZai3AMm1GSs63KmThs5glxhVa96k0HLpw65tpLZry/af++irq3P9t78+Vz1Y1HXel9XuXl5b300kvt7e2VlZU4jk+cONFut6f1HkGqtAZiB0+5xc5uXTiO5ZrYbNpUkmW50eU70dha3+Jxefwud9DpDrR6ApEY35c/ThJErtWYl2PJs5vNeiYxFYWmSZrueE6ZDToMw0x6hiZJo57Vsx0NQw06xmxkE6+DFEH2vSolLkqCIMJOUAbiBCkQifes5zXpqFG5pkKrAUrjtc7IUhOLbUcaOs9mR/hMOCf/X3dctvtIbXOb75k3v5w/deyIwqH76zWNqcw333xz9OjR73znO7m5ubm5uem7I5ByLn90f117coGcJvEcE6v1NqOeQKSmofVEvaumqT1xhIePn2OVBcfxHKuxMMeS57AU2E35OZbCHEu+3ZznsDgshkFe0aVJAhZaMg0vSP5ovNugRwxD+RZ9Sa7JkQG/7UCqjMwxtQc4VyCKEIryYphS/2y2nqV//b1r7v/dqxwvPPr8upd+fSc+VMfepTGVue+++9xu96WXXjps2LD03QtIubZg7MApdzKPMbKUzcCovQ/bb3FRqm1sq65vrWlsra5vPdHg6rXhW4LVbChwmAsclo6sxWEucJjzHZZcm4nQeAIH0iQuSsGYEOXPyIZJAh/uMJTkmfU0bCRloSkjHf5KnkuezSYJSu31tpkTR95y+bzXP9l1uLpx7ce7vrvqInXjUUsan2+TJk1avHhxr3nMkSNHJk+GNoWZqC0Q+6a2PdkBz2bQRnFMXJSC4Vh9i+d4vbOy1llZ56x3uqXeZtkghGxm/dgRBWOG540ZkTduRP6oYbnJ3R8AzkkQ5UAs3i2JoUl8dL5lVK4JWttlMZrEp43K2XOiVUFIUZA7xBXY9Kr/vP/j5qU7Dlafcnr+9s7mBdPHjh6ep3ZEKkhjKvOHP/xh/fr1iqL0XId/6623IJXJQJ4w1zWPsWZqHuP2h4+fch4/5TrR0OpsD7S0+3pt9ZZAk0Tp8LxxI/NHD88fOzxvzIj8HKtxMKMFWSMuysFoPBYXu5b1UgRemm8uyTVpfQcW9EWOiS3NN59sDSKEBEn2R+I2g8pvhBia/M0Dq7/76MtxUfqf59at/d09JDHkdqLTmMq8+OKLbrd76dKls2fP7prNBIPBjRs3PvbYY+m7azAAIU745uQZeYw5Y/IYp9tfVeesqnNV1rVU1TnPkrgghGxm/cTSYeNG5o8bWTB2RP7IQgdsEoHzxMXFYEzgzqyJIXGsJM9cmm+Gk9VDyvgiqzvEJY52hmJxHYWzau8nXjC2+LurLnr5w21Vdc6X1m27//rF6sYz+NL4A6itrd20aRNC6Pjx4+m7F5ASnCDtrWkTOndk1M1jFEVpbPUdr2upqnNVnnJW1Tl7ncmMECp0WIbl24ryrIU5VpvZUJxnLS3OK8jpfVwGAP2lIMTFpUAs3q3THY5hxQ7D+EJrNp3pA32EY9iMktxtlS2irCCEPCG+wEaovrF4//WLtx08caLe9dK6bYtmjp9QUqhuPIMsjanM8uXLCYK4/vrru00waG9vf+aZZ9J3v6C/JFn55mRbrPM4j5GlBj+PCUe5b46dOlTVUFHbfPxUazjK9bwNjuOjCh0TSgomjCocX1I4bmSBxagb5DjBEKEgJcKJwVj3oY8kgY/MMZbmmSGJGcoMDDmx2F7e4EEISYriDXO5ZpVfiyiS+O0D19zx8IuCKP3q2XVv/OH+IXXgMY2pzLJly4YPH75o0aKeX4pGv/UsCRhkCkIH69z+zj54Opq0D2Jb7lNOz85DJ7YfPLH/2CmhR1sOkiBKh+VMKC2aMKqgrKRw3MhCmA0EBkEsLvki3P9v777DoyrTv4Hfp0yfZDLphSQkkEJPABUQpSi4oIK6KoLY2EWx89rW1bWtuu7aFrGtPxQRVxZUFCwgRVFUUFARkV4jLaRPnzNzyvvHmRmGgBBgkjOTfD/XXnudOZlyJ8TMd57z3M/TLMQYeK4ww1qUictJQERUkG6tc/kONHqJyBeQ3L6gVesr8qWF2X++7NxX312xa1/Nf9794o4JI7Stpy217h5M2dnZx/zSpEmTWu914aRs2d+orpRARHodl55kbO290mRZ/n7Drq9+3Lrq5x37ahqjv8RzXGlhVreinPKi3G5F2V0KsrAnM7SlQFBqPGqdGIuB75Jl65Rq0XyBV4grPfPTGtzh3myvYNBrv2/2pLHnfPXjtk0798/+ZPWQ/uV9SvO1rafNxPJ94o033rjhhpNbokdRlNdff33y5MkxLANabn+DR52KT0Qcy2YkGVt1/Zide2sWf7vh06/XH6p3Rp/PSU8Z3LdkWP/yirICowHjLqABUZabPM1brJNN+pJsWzw03EIcatabXe/0Z2n9q8Jx7N9vueTq+18TguIjr3449583d5C/qLGMMv369XvkkUcef/zxlj/koYceuuCCC2JYA7Sc2x/c8Fu9eswyTGaysZVmrjU5PYu+3fDRinXbfjsUOannucpuhYMqSs7u07W4ExaDBs1IsuL0Bd3+QPT+j0YdV56bkpdmRYiB44juzQ7ER292cV7GLVed9++3l/x2sGH6/5bfd/0obetpG7GMMhUVFRs3brz66qufeOKJoqKi4995y5YtTz/99OjRo88555wY1gAtJCvKT7vr1Bn4RJSaZIj5spWyrHz3y46FX/781Q9bAuF5MAzD9C7pNGZoxchBvSxYmA40Jcmywxf0+ILRk2J4lumSZSvOSta8JwUSQhz2Zl89asCKtVt+3lI1b8maYWeUn9HjBG/H7UCMf+JXX321wWAYNWpUv379hg4d2q9fv9TU1JSUFLPZ7Ha7nU5ndXX12rVrV61atXHjxmnTpg0fPjy2BUALbdrXqO5ZT0TJJn1s11lvcnk/+ern95at3Xvo8FSYgpzUC8+tGHV2r06Z9hi+FsApEGXZ5Qu6jwwxDFFeqqVbnh3dSdByod7sLQfVzXfjoTebZZnHbr7kqr+84vMHH3114bvP3Nzut6GNfXi8/PLLKyoqHnzwwUcfffSYd2AY5qKLLvrss89yc3Nj/urQEgcaPXtqXeqxnmdtlphdTN2wY9+7S9YuW/1rZBjGaNCdf1b3S4b1rSwvaOP9FwGOFpRkpzfgCYgUlWIYohy7pSTbFp/LW0Ocsxj47p3sv1SFerPr3f5MrXuz87Psd04Y+c+Znx6sa3pu9pKHbxqjbT2trVXGwbp27Tpv3ryNGzcuWrRozZo11dXVbrc7OTk5Ozv77LPPHjVqVElJSWu8LrSERxB/+a1BPWYYJi3JxJx2z5IQFBd9vf7dpWu37qmOnOyck3blBWdePKSi3X8ggIQQECWXL3h0iMm2m8tyUjTf5RgSWkGatc7pP9DoISJ/QHL5g0la/0ZdMaL/irWbv9+wa8GKn4af1W1wRXt+22UURTnxvdqFurq6jIyMiy666OOPP9a6Fs0oRKu3VTe4BfVmepLRbDitOKsoyuJvf335f8sP1jvUMzzHDT+r2xXn9+/brRDDMBAPAqLsPGoDSJZhcu3mkhybpWO0eEBrC0ryys0H1bVGGaKsFJPmi9QdrHeMu/dVt9efYU9695lbWnVN0d4Faa335CeERTs6lt01zkiOSTLpTjPHbNix7/nZS9Zv26veTE+xXnxuxRUjz8DWARAnhKDk9AV8gWNsO1CSbTNpPT0T2hMdx1Z2Tl+9vVpRSFEnzaSYtf00l5Nmu/vaCx77z8LaRtfTby568vY/allNa2rd/5JFUfT5fBaLhWXZ6urq6dOnb9u27ZJLLpk4cWKrvi4ck9sf3HqgST3mOcZmPvXrPnsO1k9/Z9mXP2xRb6anWG++ctiYIZXYuBHihBCUHN5Asw0gOZYpSLd2ybIZMbEXWkGq1dAly7aj2kHqvtlewd6Gi6cf09ihlSvWbl7547bF324Yfla3887srm09raQVo8zkyZPXrVt3/fXX33jjjQzD/PnPf163bl1mZub06dMlSbruuuta76XhaIpCP1fVS7JCRAxRmtV4alPsXV7/K3O/mP/5j6IkEZHJqLv2orOvuWiQGZ3VEB98x9oAkufY/DRL1ywbupOgVZXlpNS7/I0egYhcvqBRx5v0Gv/K/W3ymCu3vdLk8v7j9U8qywtTky3a1tMaWvEztMViWbFixW233abX61966aV169YNGTLkhx9+WLVq1cqVK1vvdeGYtlc3NXlCl5asJt2p/UH/9uftV9zzyryla0RJYln2kmF9F/z7jpsuH4ocA5pTiHwBqbrJW+v0RecYPc+W5tjO65nXo1Mqcgy0Noahys7pfHh8ut7tl2SNJ6Smp1jvv2E0ETU6vU/O+ETbYlpJK47KFBQUJCUlEVFNTc0rr7yi1+ufffZZjuOIqKysrPVeF47W6BG2V4em5eo4NuXkLy15/IHnZy/58Isf1Ztn9Sq+a+LIksJj77EF0JYUIp8gOryBoCRHn1c3gCzOTOZx3RPakNnA9+hkX19VT0SyrNS7hcxko7YljRzU8/O1W5at/nXF2s2Lvvll9ODe2tYTc60YZaqqqtSDxx57zOv13nXXXZGFZNRAA21DlJWf99SFO9WYtJPfaOnHTXseeXXBgdomIrKYDHdfe8HYoZXoTgLNKaR4BcnhFZrtYm3S88WZyQXpVqzYC5rIT7PWHu7NFl2+oOZLFv31hlHrNu+pa3I//ebi/t07Z6Yma1tPbLVilMnMzJwyZYrT6fzqq6/OPPPMqVOnquc//PDDzZs3t97rQjO/7m3whNtQUyx6/clsUBAQpRfnLP/fZ9/JskJEA3p3eeimMTlpaFACjSlEHr/o8ArNRu/Ner4oM7kw3YpdrEFbvQpSGz2C2pvd5BEMelav6Wf4lGTLg3+++P89+z+nx/f3//v4xb9MaE8fR1tx3PW+++4rLCysr68fO3bsa6+9xnHchg0bbrzxxmnTpvE8eiDbSHWTd1+9Wz026Ljkk1m1qcnlvfmJt95ZtFqWFZNRd/+kC1/+60TkGNCWrJDLFzzQ4Gk4chaC2cD3Lkgb1iOvKDMJOQY0p/Zmq2lBIap3Cpov4jakf9nFQyqIaNXP2z/44ieNq4kpLJHXnvmD0lebDwRFmYhYhsm2m3i2peF1x2+H7nz6fwfrmoioV9dOj992WUF2aivWCnAiikJuf9DhC8hHjsQkm/Ql2bZsuxn5BeLN1gNNkXmKSUad3apxb7bb67/yvler6xwWo37u0zfnxW5HPG2XyMNsuPbs5z11ao4holSroeU55uuftt3wyEw1x4w6u9eMR29AjgEtKeT2Bw82ehs9QnSOSbEYzuiSeW63nBzkGIhLpTkpKeGlZVz+YLPVGtue1Wx8+MYxDMN4/IFH/7NQ1rq7KlZwoafdOtDoqXP51WOL8SQW9v3vJ6umzVkuyzLLMrePH3HdxYNarcY2o0gyyYoiy4qsqAeypJAiK7KiyDIppMiKwjIMMQzLEEPEMBS+yTBEDMOwDBFD6k1FUUSZJFkWJVmU5aNHNhliOJbhOIZjWVZ9lPrMDKOQoijEsYyOZU9786sOwReQmrxCJJSrUq2Gkmxbhtab9gEcH8NQ387pKyP7Zrv9OSlmbWejD+jd5fIR/d9buvbHTXvmLvl+wqgBGhYTK4gy7ZMoK5v2NarHHMukWlq07ouiKM/MWjx3yRoiMhv1T972xyH9471tXlZIFCVRURhiWJbUrTEVRRFlOSjKoqQEJTkoKURx9+GDYUjPc3qeVf8fyeZogig1eQQheESISTHry/Ps6Ukad7cCtJDZwPfMT/15Tx3FTW/2nRNGrP55x76axpfmfj6ooqRzjpbXhmICUaZ92lHtiCzZbrfoWzJTXZLkx15b+MnK9USUnW6bdt/40oK4WDZGVhRJVmRZkWRFVhRRVgdXFCKSZcUvSvGXUlpEUUgISkJQIgrSUcmG79htxMfcANJi4Mty7bl2s1ZVAZyaTqmWGocv0pvt9AWSTVouK2o26h+75dLJf3/TLwQfefnDN/8+iW3x9IP4hCjTDnmE4K5DTvXYoOPMLdj4NyhKf53+/hdrNhNRcV7GKw9e06arDigkykpQkgKirCiKrJAoy6HscoyrN6eLZRgdz+o5VsezOo7VcSzLMjzLSrKiKEow/NJqZhIP3zyiFI5lzHreqOdNek531Apskqz4g5I/KPkDYlCSW7LcZ7NkwzKMgecMOlav4/QdJtcEJdkriF5BbL7YnY4rzbYVpFvbU/sodCjRvdkOT8Ck54/+u9GWKssLJowe+N9PVm3Yse+tj1fdMHawhsWcPkSZdmjTvsbQuy5DqS3YzCwgSvc+P+/rn7YRUXlRzst/vcae3IoffBWFfEFRvfojybIoy6JEMb8AxDKM2cBbjTqrUWcx8GY9rwYXPc+dci5QiERJlmWFYZiTWp6HiIKSrCiKKCkMQxzLuv1Bh1do8gYc3oDHHzz6m5cVxRcUfUEiImJIz7F6njPqOL2Obfn07UQhyrJXkDz+YLMEQ0Q8x3bJSi7OTO4oaQ7aqWb7Ztc5/Zrvm33rlcNXrdu+a3/ta++tOKeypGtBlpbVnB5Emfam1uk75PCpx0lGne5E77hCULz3uXnf/LydiPqUFbz4lwlWc6tcx1WI/AHJFxC9gnjKQy16njPwrF7H6TiWYYhjGJvZkJ5stBh4dVRDIYUhhmHIqONj/meCIdJxLJ3SMlfqJzB9+D+4VKshNdyWKcqK0xtweAMOr+DwBtxCsPmPR6GAKAdE2e0PEhHHMgae0+tYg47Tc1zijlPIsuIRRK8gCuIx2jp4lslPTyrJtp1sagSIT+pc9W0HQ/tmN3qEVE17sw16/rGbL7n+kTcCovTQKx++/eRkPmEX4keUaVdkRdkYnu3LsswJL8cKAfHu5+et+nk7EVWWF7x4/8SYbwypJhivEPQFpBMmGHXsRJ0sYuA5g44z6Dg9zxp1nJ7nDPxx37YZanmXVlzhWSacbJKISJQVh1dodAuNHqHB07xzh4gkWfEGRG9AvcXoefVSVGiGTVtXf/JkWfEGJK8QPOY8J5ZlMpNNeXZLps2EkRhoZ0qyU2qdoX2z3f6gUc+Z9Vr+1erRNe+GsYNf/2Dl1j3VM+avvPnKYRoWczoS8k8//J5dh5zqB3ciSjHrj/9O4PEH7vzXOz9triKift0Lp9830XQyawGfUCAouQTRK4jHXIbRqOMsRp1Zz5v0vFnPqReD9HyifiaIIZ5l0qzGNGtobMztDzZ4Qskm8o8bRQmISkCUXeEBG3XisDp2xcbNiI2kKIGg5A/KQlAMSPLRCYZhmIxkY67dkm0zJUQgAzgFDEOVRekrN4d6sxvcgiHl1C95x8Tky4as/GHrtt8OzVz49ZD+Zd2LczUs5pQhyrQf3oAYWVZSz7PW4872FSXp3ufnqTlmQO8uz911VQxzjD8gOX2BSAtVNJOez7KZ8lIt9hZM4gEiUqf7FKRZiSggyk2e0GiNwyOIR80mlmTFF5Aiy3DpOFav4wwcq9dxeq6tm71FWRaCsjqX+ehJMCqGKDXJmGe3ZKeYcSEJOgKznu+Vn7ou3Jvd4PanJ5s0zDI6nvv7bZde88CMoCg9/MqH7zx1k0GXeMEg8SqG37NpX6PaKcMQpVoMx3/f+sfrn373y04iGlxZ8uzdV8VqOMQXkBxeIXDUNRGzns+xm3PslhSzli2IiU7Ps5k2U6bNRESKQi5/oMkTaPQITR7Bfay5w0FJDkqyh4hCzd6sgef0Os5wMnOf1TX91P8RkXqVUFGH2pTQTTn6S7IiyrI/KB2/b8tuMeTaLTl2s1GHoTjoWPJSLTVO3/4GDxH5ApLLF0zWdN/s0oLsG/849OV5n+/aV/vqvC+mThypYTGnBlGmnWhwC9VNXvXYYtTpj/v28PoHKxes+ImIunfJe3rqlTHJMd6A6PQGmoUYHcfmp1lz7eYUjMHEGsNQskmfbNIXpFuJSJTkJm8o1jR5AkdPpFUUEoKyEJTJF7oUpeM5liEllEoioSQcTZRwXokds55PtRpSrcb0ZKO2UwQAtNUrP7XRLXjDvdlGHaftqOT1Y87+8sctG3fs/++i1UP7l1WUF2pYzCnAX5N2YsuB0GxfhmGOP/KxdNWvr763gohyM1JeuG+CsQWrzhyfIEqNnkDgyMtJRh1XlJlcmG7FvIe2wXNsepIxsgauVxAbPUKTV2jyBBzewNETriVZkQLiUU8TYwyR1aRLsxrVBIMBGAAVz7GVRemrtqm92Uq9S+PebI5j/37LpRPu/48QEB/9z8L//fPm2E6dbG2IMu1Bncvf4BbUY5tJx/7+tYMNO/Y98uoCRVGSzMbp91+dZrOczutKstLkETxHLslq1HFds20F6db4mXPaAZkNvNnA56VaiEhWFIc30OQJNHmFRo/QbAndmFBbzflwh7yO5+wWfarVmGo1aLsOGEDcslsMJdkp2w42Uag3259q1XJDg6Lc9FvHnff820t+q26YPmfZXyaN1rCYk4Uo0x5sPdCkHrAsY/39BuwDtU3/7+n/CUGR57hn7rqyOC/jlF9RIXJ5A05fIHo6hEnPd8lKRoiJNyzD2C0GuyXU7C2IUpMn4PQGPEJQUhSWYTiWYYjhOYaIdBzLMAzHMizDsAxxLMswxEeSClHoS6z6KMKoG8ApK8m21bl86gdRt1806UWTphdeJ4wa8OUPW37aXPXusrVDzyg/q1exhsWcFESZhFfr9KmrFBBRskn/eyMyQlC869m5DU4PEf31T6PP7Hnqv6PegNjkEUTpcIrhWaZrtq04M/k4A0IQJww8l2UzZdmwpzSAxhiGKjof0ZudrWlvNssyj958yVX3ver1B/7+n4Xznrm5lVZMjTl8okp4W8JDMhzLJP3+1c2n31y8raqaiK65cOClw/ud2mvJslLr9Nc5/dE5Ji/VMrRHXtdsG3IMAMBJMev5XgWp6rEkK/Uuv7bb43bKtE+9eiQRHax3PDt7iaa1nAREmcRW3eR1hJd9TTbpf+/CzuJvN3z4xY9E1K974R0TRpzaa/kD4oEmry9qrqjNrD+7LLuyczqmcwIAnJo8u6VTamjaoj8ouXxHr4TZpv54fr+zK0qI6KMv161Yu1nbYloIUSaBKURbDx4ekrH+zpDMzr01j//fR0SUZrP84/bLuVOZ3KA0eoQap18OT40x6Lg+hWmDy3Ow0h0AwGnqmZ8a2XfF4Wm+qkUbYxjm4ZvGJFtMRPTkjE/UaQlxDlEmgR1o8ETyu8187CEZrz/wl2nv+YUgyzKP33pZhj3pZF9FlOVqhz/6g0K2zTykW05+mhXXkwAATh/PsZWd09W/4QopdS5/jNd0OkkZ9qS7r72AiBqcnmdmLdaylJZBlElUikLbwkMyPMdYjceewf2P1z/etb+WiG65cviA3l1O9lW8gniw0RdZM4Zjmd4Faf27ZGCzJACAGLJbDKU5KeqxKMmN4fU1tHLxkIrzzuxOREtW/bps9a/aFnNCiDKJan+jJ7Kgi81soGPtUzD3s+8XfbOBiAZXlFw/ZvBJPb9C1OgW6lz+yIqvySb9OeU56tqyAAAQW12zbKnW0CV7txD0tv4ilsd3/58uTEm2ENFTby6ud8T1ZSZEmYSkKLQjvHOkjmMtx1qKYPOug9P+u5SIctJTHr/tspNqL5Jk5VCT1xUWL+RnAAAgAElEQVS1D3NhetLgsuzfm44DAACniWGosnN6ZFXJBrdw/I3MWluazfLXSaOJqMnpeXLGxxpWckKIMgnpQJPHHc4ZyWb90SMyXn/ggZfeD4gSx7H/uOOPNutJLCLiC4rVTd7IvDOeY/sWZfQqSEWvNQBAqzJF9WbLslKndW/2iAE9LhjUk4i+/GHLx1/9rGktx4Mok3gUhbYfjAzJMMccknlyxsdVB+qJ6Pbx5/cpzW/pMxM1eYRahz/yUcBm1p9TnpNrN8eicAAAOIFcuyU/LXQdXwhKrvByG1r5yw2j01OsRPTsW58dqndqW8zvQZRJPPsa3FFDMoajh2S+/mnb4m83ENHgypJrLhzYwqdVFKp3+Z1RnUp5dsug0myLAUtCAwC0nR6d7JE/vE3eoLa92SlJ5oduHENELq//sdcWKtr2Vv0ORJkEc8JZMh5/4KmZnxKR1Wz82+SLmZZthyTLyiGHN7LRIMcyfQrTKovSNVxCGwCgYwr3Zqt/fpU6l0/TOTN0Tt/SsUMriei7X3bOX/6jlqX8DkSZBFPt8EYal445S+bl/y2vrnMQ0dSJIzNTk1vynKIkH3Icnhxj0vNnl2VHRjgBAKCNpVgMZTk29ViUlMhGe1q5+9oLctJsRDTtnaX7ahq1LeZoiDIJZueh0KVK/lizZH7dvu/dZWuJqF/3wkuHVbbkCUVZrnH4g9LhyTGDy7KTf397bQAAaANdsm1p1tBujh5/MDJqrgmr2fjwlLEMw3j9gUdfXSBrO0x0FESZRFLv9jeFs3mSqfmQjChJT7z+iSwrBh3/tz+PacmlpaAk1zj8ohwaj8lINg0qzTZgQyUAAK0xRBWd06J7syN/qzVxVq/iK0ecQUQ/ba6as/g7DSs5GqJMItlZHRqSYVnGami+xMuM979S976+fuzgwty0Ez6bEJQOOXzq5vJElGUzndElA5NjAADihEnP9y4M/TGXFaXeJZCmoyF3TBhRkJ1KRC/P+3z3gTotSzkSokzCcPoCNU6fepxs1DUbc9m4Y//Mj74hotKCrD9des4Jn80fEGucvsggYXaKuV9xBtuyOcIAANA2clLM0b3ZDp+Wvdkmo+6xmy9hWUYIiA+/8qEkaTlKFA1RJmHsCA/JMAxZjlx1NyhKj762UJJkPc89cdtlPHeCK0S+gFTrPLxdWUG6tV8RcgwAQDzqkZ8a6c12eIORTfE00aesYOLogUS0ccf+Nxd+o2El0RBlEoNXEA82hbbAsBp1zS4DvfXRtzv31hDR5D8O7VqQdfyn8gWkOufhFSRLc2y9C9IQYwAA4hPPMpVFGYd7s91+WdPFXW4eN7y4UwYRzfjgq22/VWtYSQSiTGLYWeNUf3UZYpq1F+2vaXxjwUoi6pKfed2YQcd/HiEo1bl8Svhya7c8e2QvVgAAiE8pZn15bmTfbKXRo+VlJoOO//stl3IcGxSlh1/6MChqOUqkQpRJAAFR3lfvVo/NRr7ZkMzzby8RAiLDMH+bfPHxLy0FRLkm6rpSj072LlktWngGAAC0VZyVnJYUL73Z3YtzJ409h4i2/XZoxgdfaViJClEmAVTVuiKbIiUZj1hL5tuft69Yu4WILh7S5/h7LQUlucbpi6w5XZabUpSJHAMAkBgYosrO6To+et9sLWfdTv7juWWds4nozYXfbNyxX8NKCFEm/smKsqfOpR4b9byePzzuIgTEp99cRERWs/H28SOO8ySiJNc4DvcrFWcml2TbWq1kAACIPaOO61NwuDe7TtPebJ7jHr/lUj3PSZL8yKsLfELwxI9pNYgy8W5/g0cIz1dPOrJxaeaCr/ceaiSiW68anmaz/N4ziJJyyOGLjOsUpFm7d7K3Wr0AANBaslPMBemHe7OdmvZmdy3IuumKYUS0a3/t399YpGEliDLxbndNaEhGx7Mm/eEhmaoD9W999A0RdS/OveL8/r/3cElWaqJyTK7d3KvgxKvnAQBAfOrRKdUa/ljb5AsIms66ve7iQb26diKij7/eoGEZiDJxrd7lj4TuZkMy/5q1KCBKLMs88KcLWfbY/45qjomsdZ2dYq7snIG+awCAxMWxTGXn9NBKYArVOwVFu95slmX/NfWKK0eeMf3uK7SqgRBl4tzu2tCQDMsylqidCj5fs+m7X3YS0aXD+3bvknfMx6o5JhhejTHbZu5XlI4cAwCQ6GxmfVmkN1uW691a7pudnW67f9KFw/uXaVgDokz88gXEQw6vepwUtVOBEBD//fZSIkpJttx21fnHfGyzHJOZbOpbnN6SDSYBACD+FWclp4d7s72C6NG0N1tziDLxq6rOHRo1ZMgadXXp7U9WHahtIqJbrxxus5qOfqCiKDXOwzkmPcmI/ZUAANoThqiic7qej5d9s7WFKBOnZFn5rS68LJ7u8LJ4dU3uWeq2kYXZlw6vPMYjFap1CUEx9DudlmQ8o0sm9rsGAGhnjDqud7iNQ1GUeqfG+2ZrCFEmTu1v9ATC89KtpsNDMi/P/dzrDxDR3ddecMzZvo1ewR8IjTSmWg1nIscAALRT2SnmwvQk9VgQNd43W0OIMnHqcA82xxp1oR7sLbsPfrzyZyIa0r/sjB5FRz/KIwRdvtA6RWYD378YOQYAoD3r3skemYHg8AUETffN1gqiTDxqcAuHe7CjhmSenf2ZLCs6nps6YeTRjwqIckN4HjvPsWd2yYxcRgUAgHapeW+2S8vebK3grS4eVR2rB3vZdxt/2lxFROP/cFZhbvNl7iRZqXH6wrtnU9/O6dYj16EBAIB2yRa9b7Z8+DNtx4EoE3cConSwKdSDbTXwatQWguL0OcuJyJ5s/tNl5zZ7iKJQrfPwFkslOSmZtmN0NgEAQLtUnJWckRzeN7vj9WYjysSdqjq3HB4etBj16sG8JWv21zQS0ZTLhyWZjdH3V4jqXP5AuGUp124pzcFWkQAAHUtFYXpkv+EGtyBKHegyE6JMfFEU2hvuwTbpOR3HEJHL639zwddEVNwp44/n92v2EKc34Au3LKWY9X0KscUSAECHY9Bxkb//iqLUu/wdpzcbUSa+HHJ6veFcYg0Pyby54GuH20dEt191frMGbF9AcnpDE4SNOg5LyAAAdFhZNlNhxuHe7KYO05uNKBNfqmpDQzIcy5p0HBHVNDjnLllDRL26djq3X2n0nYOSUu/yq7GbZZkzumQadBwBAEBH1T3PHtl72OkLCMEOsQQwokwc8QhindOnHicZeWKIiF559wu/ECSiu64ZGb2JkqJQncsXmVXTKz/VZta3dcUAABBPOJapLDrcm13n8ssdoDcbUSaO7Kl1hfdcIotRR0Tbqqo/WbmeiIadUd6nrCD6zg2ew7sTFGYk5adZ27haAACIQ8kmfbe8UG+2JMsN7vZ/mYnXuoBT5HA4Zs2axTCMIAgcx02ZMsVsNmtd1GmRZGVvfXjTJUNo06Vp7yyTZYXnuDvGj4i+s0cIevyhVX3tFkOPTvY2rhYAAOJWUWZyjdNf6/QRkVcIenSspV2vNJaQozIej+eee+5JS0u744477r333vT09EcffdTn82ld12nZ3+ARw3tZJ5n0RPTdLzu/+2UnEV0+on/0mniBoNTgCq2ApOfZvpGxRAAAACIiqihMM4R7sxs9gWC77s1OyCgzd+7cpqamcePGqTfHjRu3c+fOhQsXalvVaaqqC63wq+dZPc8qivLSvC+IyGIyTL5sSORusqLUuoTIdaiKwnSTPlGH1gAAoJVE92bLilLv9lP7bc5OvCgTDAaXLVtWXl6u04WGywwGQ2lp6dKlS2U5UadqN3kDjnBPtbrhwJc/bN20cz8RTbxwoD358LWzepcghb/N0lys6gsAAMeWaTN1DvdmB4JSkzeobT2tJ/GizLZt27xeb3FxcfTJ4uLiurq6HTt2aFXVafotvOkSwzAWg06W5ZfnfUFEKUnmiRcOjNzN7Q9GVsPLTDaVZGNVXwAA+F3d8uyRPYmdvoC/ne6bnXhRZs+ePURktR7RsKPerK2t1aSk0yRK8v5Gj3psMfAMQ5+s/GXXvhoium7MYIvJoH4pKMmNntAUGYOOq+iMVX0BAOB4OJbp2zmDZUO92Q1uoV32ZifeNAun00lEzfqVLBYLHRVltm7d+t1330VuejyetLR4fPvf3+iRwjtBJhl1QVF6/YOviCjDnjRu5BnqeYWozumP/Ab2KUiL7LUBAADwe5JMum659o37GohIlOR6lxDZeLLdSLwoYzKZiKjZ+v3q2nHB4BEXAnfu3PnBBx9EbkqSZLfHY9Py3vrQkIyB53Q8+86i1ftqGonopiuGGQ2hgUGHJxAM9zcVZyVjigwAALRQUWZSrctX4/ARkS8geoSgxdCuerMTL8pkZ2cTkcfjiT6p3kxOTo4+OXr06NGjR0du1tXVZWRklJeXt0mZLeXwBprCl40sRt7jD7y54BsiKsxNGzukQj3vD4jO8FYaNrO+PDdFk1IBACBB9SlIW7nloBCUiKjRLRh4jucSb4bJ70m870SNMl6vN/qkerOoqEibmk5DpAebZchi4Ocs/q7B6SGiW64YxnEsEUmyUu8OZR2OZSo7YxUZAAA4OUf2ZlOdq131ZidelMnLyzObzVVVVdEnq6qqzGZzly5dtKrq1IiSfKAhNLxkNupcXv9/P15FRKWF2ecP6KGeb3AJkZk0PTqlWtv1io0AANBKMpNNRZnh3mxRbk+92YkXZXQ63dChQzdu3BhZRUaSpM2bN48cOZLjEmwm7L4GjxiOKVYjP/vjVS6vn4huu+o8dfaPyxf0BUPd1zkp5oJ0bLQEAACnqFuuPdkU2njY5Q0I7aU3O/GiDBGNHz/ebDYvWLBAvTl//nybzRZZ/DeBRK4uGXSc1yvM++x7Iupdkj+4soSIAqLU5A1dWjLp+d6F8dh+BQAAiYJlmcqidHWPP4WozuWX5fZwmSnxpv0Skc1m+/e//z1r1qzXXntNFMVgMPjMM8+o/dgJpN7td/lC43tJRt3MD7/y+ANENOWKoUSkKFTvEtTua4ahys5punY0RQsAADSRZNR1y7P/ureBiCRZaXAL6Ynfm52QUYaIbDbbnXfeqXUVp6UqvMIvyzJCIBAZkhnQuwsRNXqESPd11yxbqjXhf9UAACAedM5IqnX6Djl8ROQNiG5/MNFnYeKDvjaEoFTdFNrK22rg3/l0dfSQjD8guv2hARu7xVCag+5rAACImT6FaQZdZN/sw5+cExSijDb21rvV1aMZomBQnPvZ90TUp6xgQO8usnK4+5rn2MrO6Wi+BgCAGNLzXEV4/qWiUL3LryRybzaijAYUharq3OqxQc//b9Eqrz9ARDdfMZSIGt2Hu6+759nNhkS9CAgAAHErI9lUlBlaVzYgyg5PQNt6TgeijAaqHd7IBtdSMPju0rVE1Kes4Myexb6A5BFCX8pINqH7GgAAWkm3vJTDvdm+YOLum40oo4HIhF+eY+YtXq0Oydw6brikKA1uv/olHcf2Qfc1AAC0GpY5oje73uWXErM3G1Gmrbn8wTpXKK8EhMC7y9YSUb/uhf27d26I+jXqkZ9q1CXYin8AAJBYkoy67nmhjZYlWWkM7wmYWBBl2tqemtCQDMPQe5995xeCRHTzFcPdftEXCA3u5aSYO6Um2DI5AACQiAozkrJsJvXYKxzun00giDJtKijJ+xpCE369Xt8Hn/9IRGf2LO5dlt/kCQ3VGHVcrwJcWgIAgDbSpzA9oXuzEWXa1N56d+QS0ruLVwtBkYhuvWp4gzsQuUDZpzBNz+PfBQAA2oieZysK09V1P8K92YkEb5ltR1FoT3jCb6PD/fFX64no3H6lXfKz/OGGpsL0pIxkk2YlAgBAh5SRbCzKOtyb3ZRQk2YQZdrOgSaPN9xo/cHS70VJYhjm5iuGNYS7+Q06rjwPC/sCAIAGynNTbOZQb7bbF4xM34x/iDJtZ0e1Qz1odLo/+/YXIhravywzIzWyMWmPTqnYMxIAADTBMkxl56jebHfC7JuNN842Ut3kjeyD/cGS74OixDDM9WPPcYXnimckm3LtZu0KBACAjs5q1PXolKoey7JSH17qLM4hyrSR7eEhmZoGx+JvfiGic/uVZWSkqrtesAzTM9+uYXkAAABEVJBuzU4Jfa72BSS3LwF6sxFl2kKNw+fwhibEvLd4lShJLMtcP2ZwILxKdNfsZIshsfdYBwCA9qF3QVpkjdZGbwL0ZiPKtIUdh0JDMvsP1X/+/SYiGjGgR1p6aBjGYuC7Ztk0Kw4AACCKnmcrOkf1ZjvjvTcbUabV1bn8De5QV9v7S76XZYXj2PGjB0WmU/XMT2VZRrsCAQAAjpCeZCyO9GZLclN875uNKNO6FIU27m1Qj/fXNHz1w2YiumBQzxR76FckO8WMhWQAACDelEX1Zrt8gcj6Z3EIUaZ17WtwR3qUFi5bI8sKy7KXjRwQme0b2ccLAAAgfkT3ZhNRvVuQlDi90IQo04pEWdlyoEk9rmtwfr5mExGdd1a39PCQTGGG1WzgNasPAADg90X3Zkuy0uCK095sRJlWtKPaIYR7lD5c9r0kyQzDXDbiLPUMz7El2ZjtCwAA8asg3ZoT1Zvtist9sxFlWotHCO465FSPHS7P0tW/EtHQ/mU5maGE2yUrWc9zmtUHAADQAr0L00z60AWEJo8QEONuQwNEmdaycV+jrF5WZOjDZaHlfS85LzQkY9RxxZnJWtYHAADQAjqOrShMO9yb7Rbibc4MokyrqHZ4axw+9djn8X26cj0RnV1RUtgpQz1ZmpPCoQEbAAASQVqSsUt4RkRQlJu88bVvNqJM7AVEecNvoQZslmUWfL42KEpEdPkFoSEZq1GXn2bVrD4AAICTVJZjSzncmx30BeOoNxtRJva27G+MzPZVRHHBip+I6KxeXTp3ylJPlmTbGIzIAABA4mAYprIonQ9fT2hwCVLc7JuNKBNjjR5hb71bPTbpuQ+XrxECIhFd8YcB6kmjjsvBDtgAAJBoLAZdj/zDvdn17ni5zIQoE0uSrKyvqldjKsMQR/J7S9cSUUV5YdfCHPU+RZnJLMZkAAAgAeWnWXPtFvXYHxDd8dGbjSgTS1sONEb+XW0m/XtL13r8ASK6alRoSIbn2IJ0zJIBAIBE1asgNdKb3egWApL2vdmIMjFT6/TtrnGpxwae4xia+9n3RNS9OK+8S4F6viDNquPwMwcAgESl49jKzunq1QWFqN6pfW823lZjIyBKP1fVq8cMQ6lJhvmf/+Bw+4joqlFnESnq+aLMJC2rBAAAOG2pVkPXrHBvtiQ3eTSeNIMoExvrq+ojXUt2i1GW5Xc+XU1ERXkZfcqL1fM5KZbIoBwAAEDiKs2xpVgM6rHLHzwUXkpNE4gyMfBbnTvyr2jS81Yjv3DFutpGFxFNuHCgEp7jiyEZAABoHxiG6ds5nQ9Pmfilqj4oyVoVgygTA9urHeoBxzJpVoMoSW8t/IaI8rPsA/qUql+ymfX2cIAFAABIdGYD37OTXT0WREnDy0yIMjFg0nNERAylWQ0sy3z01fqD9Q4imnjR2VJ4MlTnDAzJAABAu9IpzVqYnkRESSad3WrUqgxM3YiB/sWZ+xs8dS6/QcdKkvzWwq+JKCc95dz+3QVJJiIdz0Ya8QEAANqNXgWpZbk2HcdpuGIaRmViQM+zRZlJBh1LRIu/3bD3UCMRXXvxoIAcunCYn2rF5pEAANAu6XktcwwhysSWLMuvf/gVEWWmJp83sGek1R7L4gEAALQSRJlY+nzN5t8ONhDR9WPODoancqclGa1GnZZlAQAAtF+IMjGjKMrMBV8TUUqyZfQ5FUExlGUK0jAkAwAA0FoQZWLmm3Xbt+6pJqLxo84Swyd5js1OwT7YAAAArQVRJjYURXnt/S+JyGo2jht5pscfCjO5djMm/AIAALQeRJnY+PibDZt2HSCiq0cPYDlOCc/4VRvuAQAAoJUgysSAoiiPz/yMiGxW08QLB3qF0JCMzay3mfWalgYAANDOIcrEwJLvNq/bupeIJoweYDDo/OF9JfMx4RcAAKCVIcrEwPod+4koyWy86g9nRWbJMAyTY8eEXwAAgNaFjQti4JY/nitJcn6nLKvZeKDBo57MspkMPKdtYQAAAO0eokwMJJkND1x/wS+/1fuDkiSHJvxihV8AAIA2gAtMseQJT/jV82xGkknbYgAAADoCRJmYkRXyhaNMXqpF2721AAAAOghEmZjxBoJyeDmZTqm4ugQAANAWEGVixu0LqgdYTgYAAKDNIMrEhssfDIT3j8QKvwAAAG0GUSY2Ij3YHMvkYjkZAACAtoIoExv7G0NRJtNm4jn8VAEAANoI3nRjoNEjRPZdyrNbtC0GAACgQ0GUiYED4SEZHcdm2rCcDAAAQNtBlImB6iafepCdYmaxngwAAEAbQpSJAT0f+jFiK2wAAIA2hj2YYqB/ccbeeneySZ9qNWhdCwAAQMeCKBMDJj1fmpOidRUAAAAdES4wAQAAQAJDlAEAAIAEhigDAAAACQxRBgAAABIYogwAAAAkMEQZAAAASGCIMgAAAJDAEGUAAAAggSHKAAAAQAJDlAEAAIAEhigDAAAACQxRBgAAABIYogwAAAAkMEZRFK1raCN1dXUZGRk2m61r165a1wIAAAAnMHHixKlTp574fkqHIYrixx9/3Eo/7u7du9vt9lZ6cjgpPM/37NnTYrFoXQgQEZlMpp49e+r1eq0LASIim83Ws2dPrauAkPT09LKyMq2riGv33HNPS97fea3rbDscx1100UVK64xCnXPOOU899dSYMWNa48nhpDQ2No4YMeLLL7/s37+/1rUAbdq06dprr929e3dubq7WtQAtX778/vvvl2WZYRitawGaPXv2e++9t2XLFq0LSXiYKwMAAAAJDFEGAAAAEhiiDAAAACQwRBkAAABIYNyjjz6qdQ3tQVFRUe/evZOSkrQuBIjjuJKSku7du5tMJq1rATIYDN26dSstLdXpdFrXAmS1Wnv16tWlSxetCwEiopSUlF69ehUUFGhdSMLrQOvKAAAAQPuDC0wAAACQwBBlAAAAIIEhygAAAEAC60Cr/bYSh8Mxa9YshmEEQeA4bsqUKWazWeuiOro1a9bMmzfvueee07qQDm3Hjh2zZ8/esWNHUlLSoEGDxo8fj+0LNLRr167Zs2dv3brVZDINGjTommuuMRgMWhcFJIriPffcc+GFF44YMULrWhIYRmVOi8fjueeee9LS0u6444577703PT390Ucf9fl8WtfVcW3ZsuXll19+8sknq6qqtK6lQ1u3bt3TTz/duXPnwYMHe73e+fPnP/XUU1oX1XFt3br1jTfeuOCCC/7yl7+Ul5d/9NFHM2fO1LooICJ65513du3apXUVCQ9R5rTMnTu3qalp3Lhx6s1x48bt3Llz4cKF2lbVkRUVFd16663oNdWWoijz589//vnnJ02adMstt7z88suZmZk//vjj5s2btS6tg1q+fPlDDz00cODAioqKe++9Nz8/f+3atVoXBbRx48YdO3ZoXUV7gChz6oLB4LJly8rLyyMLZhgMhtLS0qVLl8qyrG1tHZY6Zo4VZbS1c+fOkSNHWq1W9WZycvIll1xCRPj0qQlZlidMmGA0GtWbDMN06tSpU6dO2lYFXq/3/fffv/baa7UupD1AlDl127Zt83q9xcXF0SeLi4vr6uoQtKEjy83NHTx4cPSZnJwcIrJYLBpV1KGxLGu32yM3BUHYt2/f+PHjNSwJiGjmzJlXX301ZizFBKLMqduzZw8RRT56qtSbtbW1mpQEEA/MZjPLHvG3pba2lmXZPn36aFUSqHw+33PPPXfRRRd169ZN61o6tFWrVmVlZXXt2lXrQtoJdDCdOqfTSUTN+pXUz52IMgDRvv/++9GjR0ePDUDbe//99xcuXOhwOL7//vvq6uobbrhB64o6qMbGxi+++OKBBx7QupD2A6Myp06dkNHs0yfDMEQUDAa1qQkg/mzZsqWxsfH666/XupCO7vLLL3/ppZduv/12nU734Ycfrl+/XuuKOqgZM2ZMnjy52XsHnA78KE9ddnY2EXk8nuiT6s3k5GRtagKIM36//5VXXrn77ruxqEw8sNlsI0aMmDJlChEhymhi0aJFffv2zcrK0rqQdgUXmE6dGmW8Xm/0SfVmUVGRNjUBxJnZs2dfe+21+fn5WhcChw0fPvyNN95o9jEM2sbKlSs3bdo0ffr06JMvvvjiiy++OHny5IsvvlirwhIaosypy8vLM5vNzZZiq6qqMpvNWNcEgIgWL15cUlLSv39/rQuBI7Asm5SU1Kz7EtrGpEmT3G535GZDQ8P06dPHjBnTt29fdMifMkSZU6fT6YYOHfrll1/Ksqxe9ZQkafPmzX/4wx84jtO6OgCNff3118FgcNSoUZEz1dXVBw8erKys1LAqICKHw+Hz+c444wytC+mISktLo29WV1cTUWFhYd++fTWqqD1AlDkt48ePX7NmzYIFCy677DIimj9/vs1miyz+C1rxeDySJGldRYf2888/v/3220OGDHnnnXfUM16v95dffnnmmWe0LawDEgThySef7NWr16hRo6xWqyAIb7zxxp133pmamqp1aQCxgShzWmw227///e9Zs2a99tproigGg8FnnnkG64BpaMuWLd99993u3bsVRfm///u/Xr16DRw4UOuiOpzdu3f/4x//8Pv98+bNiz4/atSoyJqz0GZ4nk9JSVm8ePFnn302ePDgpKSkP/3pTzabTeu6AGKGURRF6xoAAAAAThGasQEAACCBIcoAAABAAkOUAQAAgASGKAMAAAAJDFEGAAAAEhiiDAAAACQwRBkAAABIYIgyAAAAkMAQZQAAACCBIcoAAABAAkOUAQAAgASGKAMAiWf37t1Tp0697777tC4EALSHKAMAiaS2tvaKK64oKSl54YUXXC5X9Jd27do1bdq0uro6ratEhyQAAAeHSURBVGoDAE3wWhcAAHAS0tPTZ8yY0djY+Pnnnzf70k033bR8+fI9e/ZMmzZNk9oAQBMYlQGARMIwTEpKyllnnXX0ly688MLc3Nxhw4ad1BMqivLggw/GqDoA0ACiDAAkHpY9xt+uqVOn7t+/f+zYsSf1VA888MAHH3wQo7oAQAOIMgDQcf3rX//65z//qXUVAHBaEGUA4NTV1tbOmTPnuuuuKyoqIqK5c+eOHDnSbrcPGjTov//9b+Ruhw4deuedd66//vrCwkIimjVrVnZ2dp8+fRoaGohIluVp06aNHDmyoKCga9eu11133b59+5q90IoVK6688so//OEPF1988V133dVswq/f7//iiy/+9re/DRo06Pnnn2/22LVr115++eWjRo2qrKwcOHDgzJkz1fOzZ8+eMWMGEe3evbuysrKysnLRokWRRy1btmzcuHHjxo3r27fvOeec8/DDD3s8nhZ+RwDQphQAgFO1cOHCyy67jIgMBsMTTzzRr1+/22+//cwzz1T/vPzzn/9U7/bRRx9NmjSJiHieX7t27Zw5c9QZLcuWLXO5XGefffarr74qiqKiKEuWLNHpdCkpKevXr4+8yvPPP5+bm/vjjz9Gnk2n0xHRlClT1DPr16+fM2dOZmZm9Iuq7rvvvr59++7cuVO9OXDgQCK6//771ZvLly8novLy8mbf1y233JKfn79hwwb15qeffmo0GgsLCyNVHec7itGPFgBaClEGAE7LunXriIhhmA8++CBy8tZbbyUik8l04MAB9czmzZvVu3366aeKojQ1Nc2ZMycQCNx9992XXHJJ9BNOnDiRiIYNG6beXL16NcMwM2fOjL7PhAkToqOMSp0lEx1lpk+frtPptm/fHjkzf/58IurZs6d685hR5q233iKit99+O/rkq6++qj5QEITjf0ct/skBQGwgygDAadmxYwcRGY3G6JMejyclJYWInn76afWMes2I4zhJkiJ3czgcPM/n5ub2j5KXl6cO86ihYeDAgRzHRQKESu05ahZl1AwUiTKBQCA1NfW8885rVvDSpUt37dqlHh8zyuTk5BBRTU1N9Emv15uamkpEs2bNOs53BABtD3NlACD2zGazOkaydevW6PMMw0Q3H23atEkUxdtvv31tlH379imK4vf79Xq90+lcvXp1fn6+Xq9v9jwnrGHVqlUNDQ0lJSXNzo8YMUKd2XNMBw4cOHjwIMMwaWlp0edNJpPaAf7DDz8c5zsCgLaH/wIBoFV07tyZiGRZPs59du3aRUR79+79vTts3LjxlAtQR01Odh6uWpKiKIcOHWr2JXWG74EDB065JABoDYgyANCKunXrdpyv5ubmEpF6ledo1dXVfr+fiPbt2yeK4sm+dFZWFhGtW7dOUZSWPyoyirN79+5mXzIYDETUq1evk60EAFoVogwAtIrt27czDHP++ecf5z5lZWVEtG3btk8++aTZl15//fVffvmlZ8+eRCSK4tF3OCH1sdu3b58zZ06zLz355JPV1dWRm8FgMHKclZWljr4sXLiw2aN+++03Iho8ePDJVgIArQpRBgBiQJKk6IGThoaGBQsWXHPNNZWVleoZdWik2QBJTk7OpZdeSkSTJk1auXJl5KlefvnlhQsXjhw5MiMj46KLLiKihx56qL6+PvJAda5xs9VlJEmiqEta2dnZaqPTjTfe+Oabb6rlCYLw17/+NRAIZGdnE1FycjIRVVVVCYJAROp9nnvuOSKaNWuWelLldDo///zzsWPHRsLZMb8jANCAhlOOAaAdUFMFEU2ZMmXv3r2KojQ0NAwbNmzw4MG1tbWRu23YsEG926FDh6Ifvm3bNrVjiIi6d+8+ZMiQ7Ozsvn37NjU1qXfYvXu32tPUtWvXp556as6cOddcc43NZiOirKysF154Yf/+/eo9hw8fTkR33XVX5MkPHDgQuWBkNBp79Ohhs9muvPJKdQ0bRVEcDofZbCaiu+++e8aMGY899ph6/r777mMYZuzYsWrnVCAQGD9+/IABAw4ePHjC7wgA2hiiDACclkgz9vz583v06NG7d+8hQ4Y88cQT0SusPPvss6Wlpeob/4ABA1544QWPxxP56qFDhyZNmlRWVmaxWCoqKp599lm/3x/9EnV1dTfddFN5eTnLsjzP33jjjQ899JD6KmvWrJEkaeXKlXfeeafaSZSSkvLYY49VVVWpj3W5XFOnTu3Vq5fVaq2srHzjjTdkWY5+8pkzZ2ZmZpaWlj744IPRbdVffvnlZZdd1qdPn0svvfSqq6567rnnIgHohN8RALQlRsHoKACchp07d3bt2tVoNPp8vtZ+rUAgIMuy0WiUJInjuNZ+OQBICLzWBQAAtFRkdRnkGACIwLRfADgtCma/AoCmEGUA4LSoa9AJguBwOLSuBQA6IkQZADh1jz/+uNrwTEQXXnjhk08+Gd0yDQDQBjDtFwAAABIYRmUAAAAggSHKAAAAQAJDlAEAAIAEhigDAAAACQxRBgAAABIYogwAAAAkMEQZAAAASGCIMgAAAJDAEGUAAAAggSHKAAAAQAJDlAEAAIAEhigDAAAACez/A4W0rPaF9K+kAAAAAElFTkSuQmCC" width="60%" style="display: block; margin: auto;" /></p>
<p>Here we see that the relationship between past and present roaches is estimated to be nonlinear. For a small number of past roaches, the function is steep and then it appears to flatten out, although we become highly uncertain about the function in the rare cases where the number of past roaches is large.</p>
</div>
<div id="relationship-to-nlmer" class="section level1">
<h1>Relationship to <code>nlmer</code></h1>
<p>The <code>stan_gamm4</code> function allows designated predictors to have a nonlinear effect on what would otherwise be called the “linear” predictor in Generalized Linear Models. The <code>stan_nlmer</code> function is similar to the <code>nlmer</code> function in the <strong>lme4</strong> package, and essentially allows a wider range of nonlinear functions that relate the linear predictor to the conditional expectation of a Gaussian outcome.</p>
<p>To estimate an example model with the <code>nlmer</code> function in the <strong>lme4</strong> package, we start by rescaling the outcome and main predictor(s) by a constant</p>
<div class="sourceCode" id="cb4"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb4-1" title="1"><span class="kw">data</span>(<span class="st">&quot;Orange&quot;</span>, <span class="dt">package =</span> <span class="st">&quot;datasets&quot;</span>)</a>
<a class="sourceLine" id="cb4-2" title="2">Orange<span class="op">$</span>age &lt;-<span class="st"> </span>Orange<span class="op">$</span>age <span class="op">/</span><span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb4-3" title="3">Orange<span class="op">$</span>circumference &lt;-<span class="st"> </span>Orange<span class="op">$</span>circumference <span class="op">/</span><span class="st"> </span><span class="dv">100</span></a></code></pre></div>
<p>Although doing so has no substantive effect on the inferences obtained, it is numerically much easier for Stan and for <strong>lme4</strong> to work with variables whose units are such that the estimated parameters tend to be single-digit numbers that are not too close to zero. The <code>nlmer</code> function requires that the user pass starting values to the ironically-named self-starting non-linear function:</p>
<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" title="1">startvec &lt;-<span class="st"> </span><span class="kw">c</span>(<span class="dt">Asym =</span> <span class="dv">2</span>, <span class="dt">xmid =</span> <span class="fl">7.25</span>, <span class="dt">scal =</span> <span class="fl">3.5</span>)</a>
<a class="sourceLine" id="cb5-2" title="2"><span class="kw">library</span>(lme4)</a>
<a class="sourceLine" id="cb5-3" title="3">nm1 &lt;-<span class="st"> </span><span class="kw">nlmer</span>(circumference <span class="op">~</span><span class="st"> </span><span class="kw">SSlogis</span>(age, Asym, xmid, scal) <span class="op">~</span><span class="st"> </span>Asym<span class="op">|</span>Tree,</a>
<a class="sourceLine" id="cb5-4" title="4">             <span class="dt">data =</span> Orange, <span class="dt">start =</span> startvec)</a>
<a class="sourceLine" id="cb5-5" title="5"><span class="kw">summary</span>(nm1)</a></code></pre></div>
<pre><code>Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX</code></pre>
<pre><code>Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX</code></pre>
<pre><code>Nonlinear mixed model fit by maximum likelihood  ['nlmerMod']
Formula: circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree
   Data: Orange

     AIC      BIC   logLik deviance df.resid 
   -49.2    -41.4     29.6    -59.2       30 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9170 -0.5421  0.1754  0.7116  1.6820 

Random effects:
 Groups   Name Variance Std.Dev.
 Tree     Asym 0.100149 0.31646 
 Residual      0.006151 0.07843 
Number of obs: 35, groups:  Tree, 5

Fixed effects:
     Estimate Std. Error t value
Asym   1.9205     0.1558   12.32
xmid   7.2791     0.3444   21.14
scal   3.4807     0.2631   13.23

Correlation of Fixed Effects:
     Asym  xmid 
xmid 0.384      
scal 0.362 0.762</code></pre>
<p>Note the warning messages indicating difficulty estimating the variance-covariance matrix. Although <strong>lme4</strong> has a fallback mechanism, the need to utilize it suggests that the sample is too small to sustain the asymptotic assumptions underlying the maximum likelihood estimator.</p>
<p>In the above example, we use the <code>SSlogis</code> function, which is a lot like the logistic CDF, but with an additional <code>Asym</code> argument that need not be one and indicates what value the function approaches for large values of the first argument. In this case, we can interpret the asymptote as the maximum possible circumference for an orange. However, this asymptote is allowed to vary from tree to tree using the <code>Asym | Tree</code> syntax, which reflects an assumption that the asymptote for a randomly-selected tree deviates from the asymptote for the population of orange trees in a Gaussian fashion with mean zero and an unknown standard deviation.</p>
<p>The <code>nlmer</code> function supports user-defined non-linear functions, whereas the <code>stan_nlmer</code> function only supports the pre-defined non-linear functions starting with <code>SS</code> in the <strong>stats</strong> package, which are</p>
<pre><code> [1] &quot;SSasymp&quot;     &quot;SSasympOff&quot;  &quot;SSasympOrig&quot; &quot;SSbiexp&quot;     &quot;SSfol&quot;      
 [6] &quot;SSfpl&quot;       &quot;SSgompertz&quot;  &quot;SSlogis&quot;     &quot;SSmicmen&quot;    &quot;SSweibull&quot;  </code></pre>
<p>To fit essentially the same model using Stan’s implementation of MCMC, we add a <code>stan_</code> prefix</p>
<div class="sourceCode" id="cb10"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb10-1" title="1">post1 &lt;-<span class="st"> </span><span class="kw">stan_nlmer</span>(circumference <span class="op">~</span><span class="st"> </span><span class="kw">SSlogis</span>(age, Asym, xmid, scal) <span class="op">~</span><span class="st"> </span>Asym<span class="op">|</span>Tree,</a>
<a class="sourceLine" id="cb10-2" title="2">                    <span class="dt">data =</span> Orange, <span class="dt">cores =</span> <span class="dv">2</span>, <span class="dt">seed =</span> <span class="dv">12345</span>, <span class="dt">init_r =</span> <span class="fl">0.5</span>)</a></code></pre></div>
<div class="sourceCode" id="cb11"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb11-1" title="1">post1</a></code></pre></div>
<pre><code>stan_nlmer
 family:       gaussian [inv_SSlogis]
 formula:      circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree
 observations: 35
------
     Median MAD_SD
Asym 1.9    0.1   
xmid 7.2    0.4   
scal 3.4    0.3   

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.1    0.0   

Error terms:
 Groups   Name Std.Dev.
 Tree     Asym 0.310   
 Residual      0.089   
Num. levels: Tree 5 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg</code></pre>
<p>In <code>stan_nlmer</code>, it is not necessary to supply starting values; however, in this case it was necessary to specify the <code>init_r</code> argument so that the randomly-chosen starting values were not more than <span class="math inline">\(0.5\)</span> away from zero (in the unconstrained parameter space). The default value of <span class="math inline">\(2.0\)</span> produced suboptimal results.</p>
<p>As can be seen, the posterior medians and estimated standard deviations in the MCMC case are quite similar to the maximum likelihood estimates and estimated standard errors. However, <code>stan_nlmer</code> produces uncertainty estimates for the tree-specific deviations in the asymptote, which are considerable.</p>
<div class="sourceCode" id="cb13"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb13-1" title="1"><span class="kw">plot</span>(post1, <span class="dt">regex_pars =</span> <span class="st">&quot;^[b]&quot;</span>)</a></code></pre></div>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAu4AAAHPCAIAAACDZ8iTAAAACXBIWXMAABcRAAAXEQHKJvM/AAAgAElEQVR4nOzdeVxTV94/8BM2kRBERFYVqqwCrigSLFA31ILa1qK/cQG1UkbptO44VupGXarj1nmkaqMiLY8zI24l6oiAlFgVFYGyiMNOANlJIJBAkt8ft08eHpQYMBKuft5/+Lrcc+/JlyDhk3PuuWHI5XICAAAAQE9ami4AAAAAoPcQZQAAAIDGEGUAAACAxhBlAAAAgMYQZQAAAIDGEGUAAACAxhBlAAAAgMYQZQAAAIDGEGUAAACAxhBlAAAAgMYQZQAAAIDGEGXeTiUlJeHh4fHx8ZouBOCtIpfLBQKBQCCQyWSargUA/oAo83YqLy/fv39/UlKSpgsBeKvI5fLm5ubm5mZEGYD+A1EGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGhMR9MFALzrhC2tN5If3bzzqLC0qq5ROIhlMNxy6DSvsfNmegw1GaTp6gAA+juMygBoTGubZN+Jf9p4Bi9au++XpIdiqdza2pyhrcN7nLt6y7HhHkHrd52qaxBqukwAgH4NozIAmlFeWbtg9e707EK/aZOXLZoz2tGWwWAoWkvKqi5cuv336Ph/XeddOvn1RDc7DZYKANCfYVQGQAP4VXWeH20oLHv+w+HN+75Z4+L0XuccQwixGW6x+S9LLnB2aWlr+wRueZCRr6lSAQD6OZVGZa5fv3727NnW1lZCyMGDBx0cHF48Zv/+/TweT09Pz9HR0cjIaPPmzV1emgkhR48enTBhwvvvv//6datFamrqgQMHCCEmJia2trYZGRlSqZQQ4uDgIBaLS0pKCCFsNjs8PLxv6qmurr5w4UJ6erpQKLSysvrkk0+8vb0JIVlZWfHx8Xw+nypp27ZtHh4efVMSvAmS9o6PQvY0i8TnTnxtO9xSyZEjba3PR0UEr9nzccieR/HHzE2N+6xIAAC6UCnKzJkz59mzZwkJCa880t/fPzg4+KVNIpEoNTW1srKy/0QZmUxGCFm5cuW8efO0tLQCAwOpKLN48WJ3d/e8vLzIyEjqmD5QWFj43XffWVtbDxw4sLa2tqio6NChQ0OGDHFxcXFzc3NzcyOEBAUFNTQ09E098OacOM99mPnsxKFNynMMZZCR4dG96wJXfR1x6PwPe7/og/IAAOhF1WtlDAwMVDnM3Ny8u6akpCSxWJyTk1NcXGxra6vi475RUqnU09NzwYIFL211cnJatWpVampqH1QikUj+9a9/HT58WF9fXy6X7969++HDh3K5/MmTJy4uLorDLCwsEGX6WHuH1OrDv6qxQ7lc3ljw9IOpEzwnuap4yojh5ks/9TsV88s/7xdp6+q+fg1/37QocMaE1+8HAKA/UPNlvy9OKincuHGD2oiPj1+7dq16H7d3pFKpl5eXkgOmTJnC4/H6ppjPP/9cX1+fEMJgMKZMmfLw4UPyQjRU8vTCm9MgFKmxN7lYJJdKlwT69eis/7dwFuenX5pqaxmGalieLW7veP1OAAD6iR5HmYqKipiYmPz8fCsrK39//2nTpqlyVl5enra2NrV9586dFStWdBnmKSgo4HA4MpnM2NhYV1e3ubk5IiLiypUrHA5HLpcTQtzc3NasWdPc3BwVFVVWVhYWFiYUCs+dOyeRSAghJ06cSExMvHv3rlgs9vHxCQoKqqioOHXqVG5uroWFRWhoqJOT04tVsdlsHR1lz4C+vr6vr29wcHB9fT0hZP369To6OhcuXCCErFu3buTIkYSQ9PT0uLg4Y2Pj2tpaT0/PefPmKU7vrkkikURGRj59+nTp0qX+/v6EED09PT09PcWJz58/J4S4uLiw2WxVnl5CyD//+c+kpCTFl42NjZaWr568AA0QtxqxmONd7Xt00pDBRi5OI38vqCDqiDIAAG+THkeZxMTElStXJicnx8XFHTlypK2tbe7cua886+bNm2FhYYcOHeLz+W1tbbdv3w4ICFC0NjY2RkREODk5bd++XS6XR0VFlZeXE0Lmz5+fnZ197949Qoi3t7e1tTUhZNasWTU1Nb6+voSQnJwcatTk6tWr5ubm9vb2ycnJFy9elEqllZWV7u7uhYWFhYWFe/fuPXXqVOesQFFl1szLy6upqSkqKooQkpaW1traWlpaKpfLExISQkJCUlJSDh065OPjs379eh6Pt3///tbW1kWLFhFClDQ9evQoPT2dEBIdHU1FGQW5XJ6amnr16lUGgxEQEKDivB4hpLS09MGDB533DBw4UMVz+5JI3H7xfo6mq+gBqVTdF0t1dAy3NdfW0e7pee+NsMx6WqqWcTne01LCHKCOntRmkIH+PHdHTVcBALTU4ygTGBhoa2sbFBSUkpJSW1sbGxs7ffr0AQOUvSy2tLTw+Xw7Ozs/Pz8Oh0MI4XK5/v7+iumSx48fC4XC58+ft7e36+rqLl++nFpYRAgJCAigokxqaqqfnx8hJCsra9myZVSrkZERteHl5TVmzJjy8vLk5GRCSFpa2vHjx3V0dAoKChISEhoaGkpKSuzte/Y+WIHFYlEbAoFg9+7dXC735s2bXl5eEonkzJkzcrncxsaGEEL9e+nSpfnz52tpaXXXpK+vb2Njo6ur297e/mJJ+/bt++233xTby5Yt+/TTT3tXdv8k6ZCm5JRouooekMnkau/ReBCrF+cZDzIkcvXkqqf8Oql+12SvWZaDWYgyANA7PY4y1NgGg8Fwdnb+9ddfm5qaCgsLnZ2dlZySnJxMLSqePn16TEyMRCLh8/mZmZljx46lDqDWDZWWloaFha1atWry5Mk7duygmtzc3GxtbYuLizMzM+vr6w0NDQUCwYtTJ1TaMDU1pb40NTWlZo6GDh1K7REK1XDLVFdXV0LI3LlzqYGoJ0+e1NXVEUIsLCwIIcbGxoQQkUhUVlbW0tLSXZO9vb2VldWxY8eKi4snTpzY5SG2bt2anZ19+vTpgoICQkhsbOyMGTMGDx78+sX3E3o62t6jbTRdRQ9IpbJz6u1RS6tR0Jv/jU2CFsJQz42gHK2HePWzn8IgA31NlwAAdNX7y34VIyKVlZXKo8zt27d37txJCGGxWFOnTk1MTCSExMfHK6IMm83+6aef6uvrKysr9+zZ4+zsvGnTJkUu8fPz++GHH+RyeUpKirm5uQZvqdJluqeiooLaiImJ4XK5hBBqyXRra6uSJmq/tbU1NV/2IhcXl++++27Xrl1Pnjzp6OgoKChwd3d/ZW0jRoyYPHmy4svGxkZqgKq/MRigu8x7rKar6IH2Duln6u1RR7ecXyOTyrS0e5ZLiksrGUov7VKdl+MIev0UAACU6P0ro+KCWSaTqeSwvLy88vLyyMhI6svGxkZq48GDB7W1tVReYTKZBw4cOH78eEZGBiEkNzd3586dx44do2agfHx8OBxOe3t7UlKSjY3NihUrel2zenV0/LEMZObMmR9//HHnpuLi4u6aVKGjo7N69WpqnRc1ovNKn376aeepKB6PFxsb29PHhb6gp99U/zz992cTx/ZgPqW+UZCVV0iYuOYXAKCr3keZlpYWamPYsGFKDrtx48aWLVs6z6SEhIRUVVXJZDIul7t8+XJCSGFhoZWV1e7du588eXLy5Mny8vKSkpKKigpq3MLQ0NDDwyM1NbWoqGjIkCH9Z7ZFMc9VVlamehOlqqqquLh4woQJL16MTBk2bJienp6WltZ7772npnqhlwazVL34WhVyQ/2GptrYf93qUZT574sJMpnMeIipdjf/YXpkgC4+fA0A3h49fkVT3P2Wz+cTQlxdXbubKCGENDQ0pKenf/HF/7lFKZvNjouLI4RwudyFCxcaGBjcv3/f0tLS19d33LhxR44cWb16dUNDQ+cFODNmzKBuVaf8NjCviVr1TTp9j1102e/s7MxkMltaWn777bdly5aZmJgQQu7du8disZQ0ubi4VFVVhYWFSSSScePG7dq1ixDy4MGD06dPOzo6fv7554aGhoSQ5uZmiUQSFhamWMQOGqGro11zc796+zx06tLmbzkPHudMnjBalePLK2rO/+PGioUzfvzuK/VWAgDwFujZbL2enl5OTg4hhM/n5+fns1iskJAQJcdfunTJ3t6+yx9jT09PakMkEl27do3ajo6Orq2tJYRoaWkRQkaPHk39+aeMHz/e1NRUV1dXcS5FIBAouiKdrkRRstGdlpYWsVhMbStmwSjNzc3URlVVVef9hoaGS5YsoR49PDz8p59+Onz4cGJi4ujRo5U0EUKKioqo2+Hk5eVRXT1+/LiqqurOnTvr16/Pz88XCASnT5/+9NNPZ82apbxsoKO1yz8c6/ze5h3/VcavfuXBzS2ir/56ZJChwe6Ny/ugNgAA2tFWrBVSLjc3d9q0aaGhoUlJSXfv3r1165aLi0t4eLiVlRV1AI/HKysrmzRpkp2dHbXnxIkT8fHx1dXVfD5/1KhR1CU1iYmJly9frq7+4xX86dOnenp61J3xsrKyHj58eP78eTs7uy+//LLzqAyDwRAIBCwW64MPPlDsjI+Pj4+Pp1Y/lZSUeHt7Hz58mLretrGxkcViiUSi2NjYtrY2Qgifzx87dqxiWXUXsbGxZ8+eVSSY3Nzc6upqGxsbJpOZk5Nz+vRpKirx+XyBQDB+/HjFiQ4ODsOHD6+rq6uqqiovL3d1dQ0NDaXSmJImMzOzvLw8oVC4ZMkS6t59jo6ODQ0NDQ0Nzc3NT58+LS4u/uSTT178sKqEhISamhpvb2/lk3qEkLKyMg6Hw2azEYb6IR0d7TkfuJ+58O+L1+6McbGzNB/S3ZGlZc8/X/8dv6LmlzM7RtuP6Msi4aXkcjn13obJZFK/zgCgcQzFrMproj4Ze82aNbNnz1ZLh11cv37dxMTkHf9E6PDw8JycHFU+GZvH402dOnXDhg0HDx7sm9qgp4rKni9YvTs7v+TDWeygxXPtRv6feFr5vPYfl27//K9bJsasuJPbPMbhniv9gkwmo0ZnzczMlN8rHAD6jJp/FRXXAqtFbm5udHT0hx9+OHXq1KysrI0bN6qxczpS79MLmvXecPO7cQcPnow7dCru2g3eMKuhI22thwwe1NgkLK+seVZQpquj89niWd+s+5PZEJVWsQEAvJvUHGXi4uLy8/N1dHQ2btz4+p99eOvWrezsbJlM1tHRYW9v/84O52ZnZ3O5XKFQWFpaqulaQJ2YBvrffPWnv6yY98vtBzfvPC4uf575e76JsaGb/fAvg/3nz5xibdHt3BMAAFDUFmW2bNmirq4UZs+eXVFRQd3qt598mLZGuLi4uLi4aLoKeFMGDzJc9vG0ZR+r9MmsAADQRb+e63VwcNi3b5+mqwAAAID+6x2dsgEAAIC3A6IMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQmI6mCwCAN6W9oyM9u5BfVVffKDQbYvzecHMXhxEMBkPTdQEAqBOiDMBbKL+Iv/f7f1xNuN/Q1Nx5/3BL00UB3ptDF5qaGGmqNgAA9UKUAXirSKWyvx44d/jHy0yDgbOne7zPHmttaTbIiFlfLygo5ienPj7KuXry5xuHtn+2atEsTRcLAKAGiDIAb4/WNsnHIXv+/Wv6koWzPl+xgGVooGgyMTayGznMb5rH2s9qDhyLWb3lWE5+6aHtn2mwWgAAtVApyly/fv3s2bOtra2EkIMHDzo4OLx4zP79+3k8np6enqOjo5GR0ebNm1+ckj969OiECRPef//9169bLVJTUw8cOEAIMTExsbW1zcjIkEqlhBAHBwexWFxSUkIIYbPZ4eHhfVlVYWHh1atXpVLphg0bqD1ZWVnx8fF8Pp8qadu2bR4eHn1ZEtCCXC5ftelIwt2Mb7eHzpkxpbvDhlkNPfLtV9+f+ufhHy9bmptsDPm4L4sEAFA7laLMnDlznj17lpCQ8Moj/f39g4ODX9okEolSU1MrKyv7T5SRyWSEkJUrV86bN09LSyswMJCKMosXL3Z3d8/Ly4uMjKSO6QNyufzRo0eXL1/OzMwkhHh5eSma3Nzc3NzcCCFBQUENDQ19Uw/Qzpl/Jvz3tZTwdcuU5BiKlhbjL58H1tY1bd1/dub748c6v9c3FQIAvAmqTjAZGBi8+iBCzM3Nu2tKSkoSi8U5OTnFxcW2trYqPu4bJZVKPT09FyxY8NJWJyenVatWpaam9k0xiYmJcrnc0NBQyTEWFhaIMm+HrP/wr9/LUWOH7e0dBw+f85g4evFHM1Q8Zeu65XfTfl/0xcGVwS//FeiRhdPGj7Qyff1+AAB6Ss3XyihZ53njxg1qIz4+fu3atep93N6RSqWdBz9eNGXKFB6P1zfFTJ8+nRDCZDLv3r3b3TFYRvvWSMst/et/XVVjh/LWFrmg+c+rejBbNHDggBV/+vC74z/99Xgc0X7dl4Ixo6wRZQBAI3r8+lVRURETE5Ofn29lZeXv7z9t2jRVzsrLy9PW1qa279y5s2LFii7DPAUFBRwORyaTGRsb6+rqNjc3R0REXLlyhcPhyOVyQoibm9uaNWuam5ujoqLKysrCwsKEQuG5c+ckEgkh5MSJE4mJiXfv3hWLxT4+PkFBQRUVFadOncrNzbWwsAgNDXVycnqxKjabraOj7BnQ19f39fUNDg6ur68nhKxfv15HR+fChQuEkHXr1o0cOZIQkp6eHhcXZ2xsXFtb6+npOW/ePMXp3TVJJJLIyMinT58uXbrU39+/8yMOGDBAlefzRSUlJc+fP1d8mZ+fP3DgwN51BbTUJrIwMxnrYtejk2Z9MPng9z/LW0UMQ6zNBgC66nGUSUxMXLlyZXJyclxc3JEjR9ra2ubOnfvKs27evBkWFnbo0CE+n9/W1nb79u2AgABFa2NjY0REhJOT0/bt2+VyeVRUVHl5OSFk/vz52dnZ9+7dI4R4e3tbW1sTQmbNmlVTU+Pr60sIycnJoUZNrl69am5ubm9vn5ycfPHiRalUWllZ6e7uXlhYWFhYuHfv3lOnTunp6XWpSpVZMy8vr6ampqioKEJIWlpaa2traWmpXC5PSEgICQlJSUk5dOiQj4/P+vXreTze/v37W1tbFy1aRAhR0vTo0aP09HRCSHR0dJco02v/+te/YmNjO++xtLRUS880UisUyeWaLkI1wlaJmnvskIwf49rTcbuhpsbDrIaWN7S+/uM3icQ1AtHr96N2hvq6A/V0NV0FALxBPY4ygYGBtra2QUFBKSkptbW1sbGx06dPVz6W0NLSwufz7ezs/Pz8OBwOIYTL5fr7+ytedh8/fiwUCp8/f97e3q6rq7t8+XJqYREhJCAggIoyqampfn5+hJCsrKxly5ZRrUZGf7yV9PLyGjNmTHl5eXJyMiEkLS3t+PHjOjo6BQUFCQkJDQ0NJSUl9vb2Pf1mKSwWi9oQCAS7d+/mcrk3b9708vKSSCRnzpyRy+U2NjaEEOrfS5cuzZ8/X0tLq7smfX19GxsbXV3d9vb2XpcEL/XNhSRJh1TTVagk/xlfvR3KpR2mQ4x7caKZ6eDyWuHrF3A2OT3hWdnr96N2Cz1d/MaO0nQVAPAG9fgzmKixDQaD4ezsTAhpamoqLCxUfkpycrK3tzchZPr06dTpfD6fWqdDodYNlZaWhoWFPXjwgMlk7tixg2pyc3OjrhHOzMysr6+XSCQCgeDF8QYqbZia/jFVb2pqSs0cDR06lNojFKrhxdrV1ZUQMnfu3KNHj7q4uOTk5NTV1RFCLCwsCCHGxsaEEJFIVFZWpqSJEGJlZXXs2LEtW7Zs37799asCIIQQwpDLejMkJZXJCS7BAgA66/21fooRkcrKSirWdOf27ds7d+4khLBYrKlTpyYmJhJC4uPjx44dSx3AZrN/+umn+vr6ysrKPXv2ODs7b9q0SZFL/Pz8fvjhB7lcnpKSYm5ursFbqnSZkKqoqKA2YmJiuFwuIYRaMt3a2qqkidpvbW1NzZeBGn23jDa3r43m3uf9ps4VTAxt7era3qxuq65pIFpquPw/dKb7TA9lrwOaoqujrekSAODN6v1LmOKCWSaTqeSwvLy88vLyyMhI6svGxkZq48GDB7W1tVReYTKZBw4cOH78eEZGBiEkNzd3586dx44do2agfHx8OBxOe3t7UlKSjY3NihUrel2zenV0dFAbM2fO/Pjj/7NypLi4uLumN2Tp0qWdL1rKyMjo7gY/bzGDAbS5JEJP7X9fdfUePcmTyeRaWj0YY+FX1lRU1TCM1bDyaICuDo2efwB4m/Q+yrS0tFAbw4YNU3LYjRs3tmzZMnHiRMWekJCQqqoqmUzG5XKXL19OCCksLLSystq9e/eTJ09OnjxZXl5eUlJSUVFBjVsYGhp6eHikpqYWFRUNGTJk8ODBva5ZvRTzXNS0kYpNlKqqquLi4gkTJrx4MXLvmJubd76pT319vVgsVkvPQA/6zLr65/cfZXtOclX9pOsJ9wghjAEq3TUKAKB/6nGUUdz9ls/nE0JcXV2VTJQ0NDSkp6d/8cUXnXey2ey4uDhCCJfLXbhwoYGBwf379y0tLX19fceNG3fkyJHVq1c3NDR0Xks8Y8YM6lZ1ym8D85rk/7P6pbs7/HbZ7+zszGQyW1pafvvtt2XLlpmYmBBC7t27x2KxlDS5uLhUVVWFhYVJJJJx48bt2rXrpQ/RZ3cZBo3wGjvqxJbFauxQKpPu+PbU309fnOLuouI6piZBc/R/X584zilkxUevX4DLyHduxRwA9BM9izJ6eno5OTlOTk58Pj8/P5/FYoWEhCg5/tKlS/b29oo7ylA8PT2pKCMSia5du0atT46OjnZ1dTU1NdXS0iKEjB49mvrzTxk/frypqWlTU5Onp2fnrgQCAbUhEolIpytRlGx0p6WlRTGMoZgFozQ3N1MbVVVVnfcbGhouWbLk5MmTIpEoPDzcx8enurq6tbV169atDAajuyZCSFFREXU7nLy8vC5lVFZWUhvUVcPwtnIcYeY4wky9fQ7V1/70z3ujzl7584pX371XJpNH7D0lFkt+Prze/j0r9VYCANCXtBVrhZTLzc2dNm1aaGhoUlLS3bt3b9265eLiEh4ebmX1x4sgj8crKyubNGmSnd0fN+k6ceJEfHx8dXU1n88fNWoUdUlNYmLi5cuXq6urqWOePn2qp6dH3RkvKyvr4cOH58+ft7Oz+/LLLzuPyjAYDIFAwGKxPvjgA8XO+Pj4+Ph4avVTSUmJt7f34cOHqettGxsbWSyWSCSKjY1ta2sjhPD5/LFjxyqWVXcRGxt79uxZRYLJzc2trq62sbFhMpk5OTmnT5+mohKfzxcIBOPHj1ec6ODgMHz48Lq6uqqqqvLycldX19DQUCqNKWkyMzPLy8sTCoVLlixR3LuPx+Ndu3bt2rVr1CU4dXV1xcXFQqGw84LthISEmpoab29v5ZN6hJCysjIOh8Nms2fNos2VsPCaRtuPqGsQRJ27xmTqK79XXntHx+7vztxMfHD6wJfT2GP7rMK3gFwup97bMJlM6tcZADSOIVfTPcWoT8Zes2bN7Nmz1dJhF9evXzcxMXnHPxE6PDw8JydHlU/G5vF4U6dO3bBhw8GDB/umNugPOjqkKzcdibmUNGuax7o/B1qav+R63uzcov1HY37PK9i7JXjT55/0fZG0JpPJqNFZMzMz5fcKB4A+o+ZfRcW1wGqRm5sbHR394YcfTp06NSsra+PGjWrsnI7U+/TC20dHRzv68IYJbnbfHIpJ/vWRN3u81xQ3m2GWhoYD6xsEBUX8O7z0tPTcYRZDrpyO+HDaJE3XCwCgBmqOMnFxcfn5+To6Ohs3bnz9zz68detWdna2TCbr6Oiwt7d/Z4dzs7OzuVyuUCgsLS3VdC1AA1+tnP//5vkcPn3p4vW7CXfSFPsZDMZEN7vvtq3689K5A/XVs3QOAEDj1BZltmzZoq6uFGbPnl1RUVFcXJyZmdlPPkxbI1xcXFxcXDRdBdCJuanxvvAV+8JXlPCryytrG5qahw4ZNHKExVCTQZouDQBAzfr1XK+Dg8O+ffs0XQUAjdlYm9lYq3mpFABAv/KOTtkAAADA2wFRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaHFekmAAACAASURBVAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGgMUQYAAABoDFEGAAAAaAxRBgAAAGhMR9MFAAD0nftPnl66cTf5tyz+87rmljazIYPeG27+4fTJn8zxsjI30XR1ANAbGJUBgHdCZl6x37Ltngs2HDt7TVdf/33P8Z/M8x3tPKq8uuGrnSftvD/bsveMoFmk6TIBoMcwKgMAb7+fryR/tvkYi2Xw9cYV/rPY+vp6nVsrq+vOxPxy+PSVX24/uHI6ws7WUlN1AkAvqBRlrl+/fvbs2dbWVkLIwYMHHRwcXjxm//79PB5PT0/P0dHRyMho8+bNDAajyzFHjx6dMGHC+++///p1q0VqauqBAwcIISYmJra2thkZGVKplBDi4OAgFotLSkoIIWw2Ozw8vG/qKSoq+vLLL7vsPHjwoFgsjo+P5/P5VEnbtm3z8PDom5IA3gLn4xKDNxz2mjxm345QQ6bBiwdYmg356/og/9le6/96zPvTzfevHh5uadr3dQJA76gUZebMmfPs2bOEhIRXHunv7x8cHPzSJpFIlJqaWllZ2X+ijEwmI4SsXLly3rx5WlpagYGBVJRZvHixu7t7Xl5eZGQkdUzfSE9P77Jn6NCh9vb2DAbDzc2NEBIUFNTQ0NBn9QC8BR7/XvD51u+nThlz9NuvtLSVTamPGW3H+f7r5X/e+XHInruXDurqYNAagB5U/V01MHjJW5kXmZubd9eUlJQkFotzcnKKi4ttbW1VfNw3SiqVenp6Lliw4KWtTk5Oq1atSk1N7bN6cnJyPvroo8577OzsOg9uWVhYIMrAW0woEndIpert84tvogYbG+375s/KcwxlxDCzyK8/X7vp0LGzv6xcNOvFA2QyWVNLGyFEVyjSUS3rMBgMY8OBPS0bAFSn5rcdL04qKdy4cYPaiI+PX7t2rXoft3ekUqmXl5eSA6ZMmcLj8fqmGLFYzGAwVqxYoeQYJU8vwFtgzld/v/d7kRo7lItb5XVVe7aFMA1UDRNeHmPenzJm096zW368TdTxCzfU2LCSu1cNHQFAN3ocZSoqKmJiYvLz862srPz9/adNm6bKWXl5edra2tT2nTt3VqxY0WWYp6CggMPhyGQyY2NjXV3d5ubmiIiIK1eucDgcuVxOCHFzc1uzZk1zc3NUVFRZWVlYWJhQKDx37pxEIiGEnDhxIjEx8e7du2Kx2MfHJygoqKKi4tSpU7m5uRYWFqGhoU5OTi9WxWazlb+v0tfX9/X1DQ4Orq+vJ4SsX79eR0fnwoULhJB169aNHDmSEJKenh4XF2dsbFxbW+vp6Tlv3jzF6d01SSSSyMjIp0+fLl261N/fn9qZnZ1dWFh45MgRc3PzESNGeHh4qPieDwC61SZiGRrMntGza8s+XTD913uZ8nYxQ2/AG6oLANSox38sExMTV65cmZycHBcXd+TIkba2trlz577yrJs3b4aFhR06dIjP57e1td2+fTsgIEDR2tjYGBER4eTktH37drlcHhUVVV5eTgiZP39+dnb2vXv3CCHe3t7W1taEkFmzZtXU1Pj6+hJCcnJyqFGTq1evmpub29vbJycnX7x4USqVVlZWuru7FxYWFhYW7t2799SpU3p6el2qUmXWzMvLq6mpKSoqihCSlpbW2tpaWloql8sTEhJCQkJSUlIOHTrk4+Ozfv16Ho+3f//+1tbWRYsWEUKUND169Ii6LCY6OloRZR4/flxTU5OYmEh9aWJism7durFjx76yQsr58+f//e9/K75sbm6mni6Ad5lc3OrpPUFHu2cvdB4TXfR0ddrFrQRRBoAOehxlAgMDbW1tg4KCUlJSamtrY2Njp0+fPmCAsl/4lpYWPp9vZ2fn5+fH4XAIIVwu19/fXzFd8vjxY6FQ+Pz58/b2dl1d3eXLl1MLiwghAQEBVJRJTU318/MjhGRlZS1btoxqNTIyoja8vLzGjBlTXl6enJxMCElLSzt+/LiOjk5BQUFCQkJDQ0NJSYm9vX1Pv1kKi8WiNgQCwe7du7lc7s2bN728vCQSyZkzZ+RyuY2NDSGE+vfSpUvz58/X0tLqrklfX9/GxkZXV7e9vb1zSV2u+a2vr9+5c+d33303atQoVYqsrq7Ozc3tvEf5DwXeEb/mlqQVVGi6ClVVNgjV3KNUOsxyaE9PGjBAd6jp4IqmNrWU0Cpp/9svv6mlK/WyszCZ5+6o6SoA1KDHUYYa22AwGM7Ozr/++mtTU1NhYaGzs7OSU5KTk729vQkh06dPj4mJkUgkfD4/MzNTMeRArRsqLS0NCwtbtWrV5MmTd+zYQTW5ubnZ2toWFxdnZmbW19cbGhoKBAJLy653faDShqnpH+snTU1NqdmZoUP/eBUTCtXwEunq6koImTt3LjUQ9eTJk7q6OkKIhYUFIcTY2JgQIhKJysrKWlpaumuyt7e3srI6duxYcXHxxIkTFZ1/88035eXlfD6/qKgoNTW1ra2to6MjNjb266+/fv3K4Z1VLWjJLa/RdBWqEkk61NmdnBC5zNBQpSULXRixDCoa1XO7vA6pvH/+CPR1MYUNb4ne/1dWjIhUVlYqjzK3b9/euXMnIYTFYk2dOpWaQ4mPj1dEGTab/dNPP9XX11dWVu7Zs8fZ2XnTpk2KXOLn5/fDDz/I5fKUlBRzc3MN3lKly4RURcUfb3ZjYmK4XC4hhFoy3draqqSJ2m9tbd1lAsjMzMzMzGzChAmEkBUrVhw+fPjhw4dPnz59o98RvPUcaXV/lNSUrBo1/tFnEKKl1dgo6MWp9Y1CoqWtlir0dLRnj7dTS1fqZTXYSNMlAKhH76OM4qJUJpOp5LC8vLzy8vLIyEjqy8bGRmrjwYMHtbW1VF5hMpkHDhw4fvx4RkYGISQ3N3fnzp3Hjh2jZqB8fHw4HE57e3tSUpKNjY3yNT59qaPjj3eQM2fO/Pjjjzs3FRcXd9ekChaLFR4evmbNGkX0eSU3NzexWKz4sqqqKjY2tqePC28f1xFmriPMNF2Fqg4b3nim1g4ZOroFxT2eX2tsaq6payRMY7XUMEBX+xOP0WrpCgBeqvdRpqWlhdoYNmyYksNu3LixZcuWzjMpISEhVVVVMpmMy+UuX76cEFJYWGhlZbV79+4nT56cPHmyvLy8pKSkoqKCGrcwNDT08PBITU0tKioaMmTI4MGDe12zeinmucrKylRvolRVVRUXF0+YMOHFi5Epenp648aNq62tVbGYWbNmzZr1v7fB4PF4x48fV/FcgLfWAIO0xzktolbVF2MTQn79LUMmlTEG9GZmCgD6Xo+jjOLut3w+nxDi6uqqZKVMQ0NDenr6F1980Xknm82Oi4sjhHC53IULFxoYGNy/f9/S0tLX13fcuHFHjhxZvXp1Q0PDwIH/+9IzY8YM6lZ1ym8D85qoVd+k0/fYRZf9zs7OTCazpaXlt99+W7ZsmYmJCSHk3r17LBZLSZOLi0tVVVVYWJhEIhk3btyuXbsIIfn5+X/7298sLS2DgoIU9w+sqqrq7vZ9AG8lK9NBI63UOSMmaWWW5vz+879urV4+79VHE0IIkUll0ReuDzAwGG7z8k9ioq7t09LSUvE+T4ONEIkA3qyeRRk9Pb2cnBwnJyc+n5+fn89isUJCQpQcf+nSJXt7e8UdZSienp5UlBGJRNeuXaPWJ0dHR7u6upqammppaRFCRo8eTf35p4wfP97U1LSpqcnT07NzVwLBH7PgIpGIdLoSRclGd1paWhQTNIpZMEpzczO1UVVV1Xm/oaHhkiVLTp48KRKJwsPDfXx8qqurW1tbt27dymAwumsihBQVFVG3w8nLy6O6yszMrKioqKioePLkSUBAwIQJE9LS0saMGePu7q68bIC3yT++XaX2Ppd9dfDsz9zZ06cMt1Zpoi02LuFZQdnVHyP8p09+sVUmk1GvA2ZmZrjzE0A/oa1YK6Rcbm7utGnTQkNDk5KS7t69e+vWLRcXl/DwcCsrK+oAHo9XVlY2adIkO7s/LnA7ceJEfHx8dXU1n88fNWoUdUlNYmLi5cuXq6urqWOePn2qp6dH3RkvKyvr4cOH58+ft7Oz+/LLLzuPyjAYDIFAwGKxPvjgA8XO+Pj4+Ph46h1SSUmJt7f34cOHqettGxsbWSyWSCSKjY1ta2sjhPD5/LFjxyqWVXcRGxt79uxZRYLJzc2trq62sbFhMpk5OTmnT5+mohKfzxcIBOPHj1ec6ODgMHz48Lq6uqqqqvLycldX19DQUCqNKWkyMzPLy8sTCoVLliyh7t03cuTIpqamuro6qVTa3Nzc2NgYEBDw4hBUQkJCTU2Nt7e38kk9QkhZWRmHw2Gz2Z1nnQDeQVPGO0VfvJ2Q8sjvg8n6+q+4Q8Fvab9HfHv60w+9/hq26KUHyOVy6r0Nk8mkfp0BQOMYilmV10R9MvaaNWtmz56tlg67uH79uomJyTv+idDh4eE5OTmqfDI2j8ebOnXqhg0bDh482De1AfRb99LzZi75erAx69CevzjajXjpMXK5/B+XE7879vME15GJ/73XYODLQw9GZQD6ITW/q1BcC6wWubm5W7dupa6SycrKmjRpkho7pyP1Pr0A74gp451SL36nzWD8v9Xf7Nj34++5hZ3fwrW2im8lpy39fOfew9EL/Dxux37bXY4BgP5Jze8q4uLi8vPzdXR0Nm7c+PqffXjr1q3s7GyZTNbR0WFvb//ODudmZ2dzuVyhUFhaWqrpWgBoaazzexk3v/8u6uIRzpXL3BRjY5a1hSnTYGBdQ1M5v1osaXd4z/qfJ7Z+PJuND20FoB21TTC9Cfn5+RwOp7i4mM1mr127tsvlw6AEJpgAXqpF1HbjzqNfH2SX8qtbWsXmpsbvDTefO23SpDEOWlqvDjGYYALoh/r1r6KDg8O+ffs0XQUAvD2YBvqfzPH6ZM4bvK0DAPSxd3TKBgAAAN4OiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAY4gyAAAAQGOIMgAAAEBjiDIAAABAYzqaLgAAAODlMnKLnuQUllfWymRya4shro42k8bYMxgMTdcF/QuiDAAA9C9iSfvfo+O/P3utuPw5IURLW0uLaHVIOwghVuYmIX+as/6zjwyZ+pouE/oLRBkAAOhHnuQULgz9trC0yps9LmTlR5PGOQ8ZbMTQYtQ1CJ5k5f87KW3X0Z+jYrix32/x8XDVdLHQLyDKAABAf3Ez5fHC0G+HmBidO7F9rItd5yZTk0EzfCbN8Jn0tMA/Ys/JWUu+PnPoqz/N99VQpdCPqBRlrl+/fvbs2dbWVkLIwYMHHRwcXjxm//79PB5PT0/P0dHRyMho8+bNL05nHj16dMKECe+///7r160WqampBw4cIISYmJjY2tpmZGRIpVJCiIODg1gsLikpIYSw2ezw8PC+qUcoFJ49e/bx48dSqdTV1TU4ONjMzIwQkpWVFR8fz+fzqZK2bdvm4eHRNyUBAPSZnGeli9bsG2U77PiBdYOMDLs7zHHUiLP/tX1jxPerNh21HWbOnujcl0VCP6RSlJkzZ86zZ88SEhJeeaS/v39wcPBLm0QiUWpqamVlZf+JMjKZjBCycuXKefPmaWlpBQYGUlFm8eLF7u7ueXl5kZGR1DF9oLm5ecOGDVVVVdSXqamp2dnZf//73w0NDd3c3Nzc3AghQUFBDQ0NfVMPAEBfksvlQev/xjQYeGTfl0pyDGXgwAEHdq1dHrpr+bpDObej9HQxw/BOU/XHb2BgoMph5ubm3TUlJSWJxeKcnJzi4mJbW1sVH/eNkkqlnp6eCxYseGmrk5PTqlWrUlNT+6aYH374YdSoUX/6059KS0svXboklUobGhq4XG5gYKDiGAsLC0QZAOgPMv/DP3HxVzV2WFBQ8ijrP4f2fGFibKTK8cyB+l9vCF4RFjkteI/bGCe11PDZfPZEpxFq6Qr6kpqTrJI1cjdu3KA24uPj165dq97H7R2pVOrl5aXkgClTpvB4vD6opLq6Wltbe8OGDdSXRkZGHA6HEFJUVNT5MCxBBIB+oriy7tQVdb48yuuq7EcNn/b+RNVPGT/GYYq7y2/3M+4V1amlhg8mOiDK0FGPo0xFRUVMTEx+fr6VlZW/v/+0adNUOSsvL09bW5vavnPnzooVK7oM8xQUFHA4HJlMZmxsrKur29zcHBERceXKFQ6HI5fLCSFubm5r1qxpbm6OiooqKysLCwsTCoXnzp2TSCSEkBMnTiQmJt69e1csFvv4+AQFBVVUVJw6dSo3N9fCwiI0NNTJ6SWZnc1m6+goewb09fV9fX2Dg4Pr6+sJIevXr9fR0blw4QIhZN26dSNHjiSEpKenx8XFGRsb19bWenp6zps3T3F6d00SiSQyMvLp06dLly719/cnhLBYrDVr1ihOnDp1KhVlBg4cqMrTSwjJzMz8z3/+0/n5ZLFYKp4LAKBhcplc3DbN272nb9im+0y69zCbIe0g2phjenf1+GefmJi4cuXK5OTkuLi4I0eOtLW1zZ0795Vn3bx5Myws7NChQ3w+v62t7fbt2wEBAYrWxsbGiIgIJyen7du3y+XyqKio8vJyQsj8+fOzs7Pv3btHCPH29ra2tiaEzJo1q6amxtfXlxCSk5NDjZpcvXrV3Nzc3t4+OTn54sWLUqm0srLS3d29sLCwsLBw7969p06d0tPT61KVKrNmXl5eTU1NUVFRhJC0tLTW1tbS0lK5XJ6QkBASEpKSknLo0CEfH5/169fzeLz9+/e3trYuWrSIEKKk6dGjR+np6YSQ6OhoKsp0F1kmTJjwygopt27dio2N7bxn6NChKp4LAP1HVWNzQ3Orpqt4tfI6gRp7k3d0ECJ3dR7Z0xOpU+Qd7Qx1RBl+vSC3vOb1+3lznIYNxeD8i3r8sw8MDLS1tQ0KCkpJSamtrY2NjZ0+ffqAAQOUnNLS0sLn8+3s7Pz8/KjBBi6X6+/vr0jfjx8/FgqFz58/b29v19XVXb58ObWwiBASEBBARZnU1FQ/Pz9CSFZW1rJly6hWI6M/plS9vLzGjBlTXl6enJxMCElLSzt+/LiOjk5BQUFCQkJDQ0NJSYm9vX1Pv1mKYnhDIBDs3r2by+XevHnTy8tLIpGcOXNGLpfb2NgQQqh/L126NH/+fC0tre6a9PX1bWxsdHV129vbuyvp999/J4SMHTtW+fwXALx9En8vSvq96NXHaVppWbU6u5NJCSGmJoN6et5QU2NCCJFK1VJF/KNnOfXqjGhqd/LzAIIrDV7Q4yhDjW0wGAxnZ+dff/21qampsLDQ2VnZWrjk5GRvb29CyPTp02NiYiQSCZ/Pz8zMHDt2LHUAtW6otLQ0LCxs1apVkydP3rFjB9Xk5uZma2tbXFycmZlZX19vaGgoEAgsLS27PASVNkxNTakvTU1NqZkjxciEUCjs6Xf6IldXV0LI3LlzqYGoJ0+e1NXVEUIsLCwIIcbGxoQQkUhUVlbW0tLSXZO9vb2VldWxY8eKi4snTnzJrLBcLr906ZKpqenGjRtxcQzAu0ZPW8tggK6mq3i1AepdNMRgEEIk7e09PY+6xoCo6ZVST1ebFk8+dNH7/4uKEZHKykrlUeb27ds7d+4khLBYrKlTpyYmJhJC4uPjFVGGzWb/9NNP9fX1lZWVe/bscXZ23rRpkyKX+Pn5/fDDD3K5PCUlxdzcXIO3VOkyIVVRUUFtxMTEcLlcQgi1ZLq1tVVJE7Xf2tqami970S+//FJbW7t79+5Bg3r8BgUA6G6hp8tCTxdNV/FqV3/N5P77kbp6Y2hrywl5/ryO/N/b4r1S5fN6Qoi6LpQJ9h0XOEPVaX3oP3r/41dcMMtkMpUclpeXV15eHhkZSX3Z2NhIbTx48KC2tpbKK0wm88CBA8ePH8/IyCCE5Obm7ty589ixY9SYhI+PD4fDaW9vT0pKsrGxWbFiRa9rVq+Ojg5qY+bMmR9//HHnpuLi4u6alCsqKrp8+fK3335LzUkRQqqrq6kb5Sm3du3a1atXK768f//+7NmzVX9cAABN0tYlWtr3H2XPmtazN6v30n4nDAZDt+ulkPBO0er1mS0tLdTGsGHDlBx248aNLVu27PsfUVFR1JyLTCajhisIIYWFhUZGRrt37961axfVW0lJiWJgw9DQkBqJKSoqam5uHjx4cK9rVi/FPFdZWZnqTZSqqqp79+79MTT6P5qamk6dOrVz505Fjrl//35mZqYqxejr6xt1YmBg0Gc39wMAeH0M/YG3f33c0tqm+ikd0o4bt+8x9AYSRu//lsFboMejMoo/kHw+nxDi6ura3UQJIaShoSE9Pf2LL77ovJPNZsfFxRFCuFzuwoULDQwM7t+/b2lp6evrO27cuCNHjqxevbqhoaHzop4ZM2ZQt6p7o5fBUqu+SafvsYsu+52dnZlMZktLy2+//bZs2TITExNCyL1791gslpImFxeXqqqqsLAwiUQybty4Xbt2Ub21tbXt2rVLS0uLWogkFotra2tLS0sPHz785r5lAIDe8XQbefNomBo7LC6t+nzToXM/cdd8pupI9j8vJ5Xxq7/b8edxPZyW6o7rKCu19AN9rGdRRk9PLycnx8nJic/n5+fns1iskJAQJcdfunTJ3t5ecUcZiqenJxVlRCLRtWvXqPXJ0dHRrq6upqamWlpahJDRo0dTf/4p48ePNzU1bWpq8vT07NyVQPDHpeYikYh0uhJFyUZ3WlpaxGIxta2YBaM0NzdTG4pPFaAYGhouWbLk5MmTIpEoPDzcx8enurq6tbV169atDAajuyZCSFFRETUek5eXR3UllUr379//7NkzQsjTp08VD6Gvrz98+HDllQMA9L2hxobTJzmqs8dJjmmPsn78+ZcJ4xynuL/6aqHc/JKjP/xjwawpG4I/VGcZQEPairVCyuXm5k6bNi00NDQpKenu3bu3bt1ycXEJDw+3svojw/J4vLKyskmTJtnZ/ZGOT5w4ER8fX11dzefzR40aRV1Sk5iYePny5erqP1bxPX36VE9Pj7ozXlZW1sOHD8+fP29nZ/fll192HpVhMBgCgYDFYn3wwQeKnfHx8fHx8dTqp5KSEm9v78OHD1PTUo2NjSwWSyQSxcbGtrW1EUL4fP7YsWO7u2tcbGzs2bNnFQkmNze3urraxsaGyWTm5OScPn2aikp8Pl8gEIwfP15xooODw/Dhw+vq6qqqqsrLy11dXUNDQ6k0pqTJzMwsLy9PKBQuWbKEundfdHQ0dTV0F46OjjNnzlR8mZCQUFNT4+3trXxSjxBSVlbG4XDYbPasWbOUHwkAqpPL5dR7GyaTSf06gxpN9xp3I/nRzxdvO9qPGDGs24/BIYRkZP/ny/DDlkMHX+V8M1AfF8q86xiKWZXXRH0y9po1a97Q1abXr183MTF5xz8ROjw8PCcnR5VPxubxeFOnTt2wYcPBgwf7pjaAd4FMJqNGZ83MzJTfKxx653lt40er9zzIzF/80YyQ5fOMjbu+/2wRtZ79mXsulmv/ntXVH795b7iyxAPvCDX/KiquBVaL3Nzc6OjoDz/8cOrUqVlZWRs3blRj53Sk3qcXAKC/MTc1Trqwd9uBc9+f+yXul2RPd9eJ453Mh5poa2s/r65Lz3rGu5fZJhavDJx58OvPWExVP9oF3m5qjjJxcXH5+fk6Ojpqub3brVu3srOzZTJZR0eHvb39Ozucm52dzeVyhUJhaWmppmsBAHizBujpHvz6szXL/U/+fP3yzXtJqY8VTbbDzD9bPPOzxbNdHPChj/C/1DbB9Cbk5+dzOJzi4mI2m7127doulw+DEphgAngTMMHU90St4vKqWqlUNszSFMMw8FL9+lfRwcFh3759mq4CAAA0xmDgAIf3ur3lBwB5nVvkAQAAAGgcogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMYDXTzwAAFApJREFUAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAAADQGKIMAAAA0BiiDAAAANAYogwAALy7pFJZQ1OzXC7XdCHQezqaLgAAAKBPyeXyG3ceXeTybv2aXlFdL5XKdHS0h1ua+nlPWDTPx8fDVdMFQs8gygAAwDskLePZlzt+uJeeN2Sw0fue4+ZbmxkPYtU3NJWUVcVe+zXqp+szpo478k3IaPsRmq4UVIUoAwAA74rzcYkh4ceHmhrv+2bNrA8ma2kxOrd2SDuucHknOHGeCzbEHN0UMGOypuqEHlEpyly/fv3s2bOtra2EkIMHDzo4OLx4zP79+3k8np6enqOjo5GR0ebNmxkMRpdjjh49OmHChPfff//161aL1NTUAwcOEEJMTExsbW0zMjKkUikhxMHBQSwWl5SUEELYbHZ4eHgfFyaXyyMiIjIyMvbv3+/s7JyVlRUfH8/n86mStm3b5uHh0cclAQDQ3cXrvOANh9mT3fbv+LMh0+DFA3S0dT4J8JnuPXHj9uMff77n6o8Rc3zd+75O6CmVosycOXOePXuWkJDwyiP9/f2Dg4Nf2iQSiVJTUysrK/tPlJHJZISQlStXzps3T0tLKzAwkIoyixcvdnd3z8vLi4yMpI7pYxcuXMjIyFB86ebm5ubmRggJCgpqaGjo+3oAAOgu9z9lQev/NsXd5djer7R1tJUcaTzI8O/fbVz5ReSyrw5l3vy7lblJnxUJvaPqBJOBwUsC7IvMzc27a0pKShKLxTk5OcXFxba2tio+7hsllUo9PT0XLFjw0lYnJ6dVq1alpqb2cVV5eXkXLlx4aZOFhQWiDAC8C8TtHdkFlWrscP2O/xowQG//jjXKcwxlwADdg3u++GhpeNg3J7/+aolaCtDT1XYdZaWWrqALNV8r8+KkksKNGzeojfj4+LVr16r3cXtHKpV6eXkpOWDKlCk8Hq/P6iGEtLa2/vzzzx4eHnfv3n2xVcnTCwDwNimtapi88oC6epO3S+Q1/M1/WWLEYqp4iqXZkKWf+p2OuXY1vYRoq+Fv5XDzwUWXdr1+P/CiHv94KioqYmJi8vPzrays/P39p02bpspZeXl52tp/BOE7d+6sWLGiyzBPQUEBh8ORyWTGxsa6urrNzc0RERFXrlzhcDjUcn83N7c1a9Y0NzdHRUWVlZWFhYUJhcJz585JJBJCyIkTJxITE+/evSsWi318fIKCgioqKk6dOpWbm2thYREaGurk5PRiVWw2W0dH2TOgr6/v6+sbHBxcX19PCFm/fr2Ojg41ZLJu3bqRI0cSQtLT0+Pi4oyNjWtraz09PefNm6c4vbsmiUQSGRn59OnTpUuX+vv7d37EH3/8ccmSJYrYBwAAatDWoqOt4++n7L3riz4O8Dl9/ppcLGIYGL2hukAtenyLvMTExJUrV86ePfs///nPkSNHuFyuKmfdvHkzLCzM2tqaENLW1nb79u3OrY2NjREREfr6+nv37t28efPAgQMbGxsJIfPnz1dc3+rt7W1tbe3o6Dhr1qx58+b5+voGBARMmjSJar169SqTybS3t6+trb148SKHwzlz5oy7u/uAAQMKCwv37t1LJZ4uDAwM9PT0lFfu5eUVGBhIbaelpSUmJpaWlpaUlFBXDqWkpOzYsWPw4MEbNmzw9/c/ffq0Ym5ISdOjR4/S09NFIlF0dHTnx+LxeEOGDHF0dFTlKe3i+++/n9bJ1q1bbWxsetEPAMBbSCJ2cxmp+pAMxcpi6Ihh5kQsfkNFgbr0eFQmMDDQ1tY2KCgoJSWltrY2NjZ2+vTpAwYMUHJKS0sLn8+3s7Pz8/PjcDiEEC6X6+/vr5guefz4sVAofP78eXt7u66u7vLly6mFRYSQgICAe/fuEUJSU1P9/PwIIVlZWcuWLaNajYz+SMpeXl5jxowpLy9PTk4mhKSlpR0/flxHR6egoCAhIaGhoaGkpMTe3r6n3yyFxWJRGwKBYPfu3Vwu9+bNm15eXhKJ5MyZM3K5nAoN1L+XLl2aP3++lpZWd036+vo2Nja6urrt7e2dS6qrq7t9+/a2bdt6V6RYLBYIBJ33KIbBAAAUou9k5PFrNF3Fq9U3tqixN7lUamU+pBcnWlsOLatuUksNTaK2v/786tUzmrLUe+zoYUM1XUUv9TjKUMMYDAbD2dn51//f3t0HNXXnexz/hYdiwoMIEVBQKxU0gleFrvWBoo51S+nD7tbeZahWRjuwnY4P3SlWdHtbrHXWaW1x4fZWrWUsq2u7tlTtDXVLRa8PpQ9zR8dewGEFA8qDBYEgCavE5P5xuimbaiAhEE54v/7KOfmd41d+Hvwkv9/5ndOn9Xp9bW2tRqOxc8jJkyeTk5OFEEuWLNm/f/+tW7caGhouXLgwc+ZMqYF031B9ff2aNWueffbZOXPm5ObmSm/NmDHj3nvv1el0Fy5caGtrCwgI6OzsHDdunM0fIaUNtVotbarVamnkaOzYHzvmxo0bjv5Nfy4+Pl4IkZqampqaKoQ4f/789evXhRARERFCiODgYCGE0Wi8cuWKwWC421sxMTHjx4/Pz8/X6XSJiYnSmS0Wy65du7KyssgfAAaV3viPlk6ju6vom76r25Wns5iVylFOHKdSjhLCNc80uG22DOef/M2e2+4uwXnOT2WyfiPS1NRkP8ocP358y5YtQojAwMCkpKSysjIhhFartUaZ+fPnHzhwoK2tramp6fXXX9doNBs2bLDmkocffnj37t0Wi+XUqVPh4eFuXFLFZn5PY2Oj9GL//v3SQJt0y3R3d7edt6T9kZGR0nCbRKvVqtXqzs5O6ZsVvf7HDwH19fVhYWGhoc58mACAn1scP3nmpAh3V9G3q9faiw+77K4LhZd3W7szX660Xu8QXq75hOnv5/tM8kyXnGowTFTLeD6Q81HGOmHW39/e6OPFixevXr26bds2aVOaBCOE+Pbbb1tbW6W84u/v/8YbbxQUFEiLqVRVVW3ZsiU/P18agVq4cGFhYWFPT8+JEycmTZq0atUqp2t2LZPJJL1YunTpk08+2fstnU53t7fuqKam5vjx41qt1mb/O++8k56enp6e3ucZHnzwQWv4E0LU1dXl5eX1eRSAkSZ+Qpi7S+iXvwf2awWQ/vLxvfj3K44e1GMyXbrcILz7mFLZT36+PsnTmcI4KJyPMgbDjwOZUVFRdpodO3Zs48aN1pEUIURWVlZzc7PZbC4pKVm5cqUQora2dvz48Vu3bj1//vyePXuuXr1aV1fX2NgofW8REBDwwAMPnDlz5vLly6GhoWPGjHG6ZteyjnNduWJ7hdh5S9Lc3KzT6RISEvqcd9xPc+bMmTPnpzW2z549++qrr7rkzAAge6OUjc0tVdV1mlgHwsQ3/1thMHYrQoMHry64hMNRxrr6bUNDgxAiPj6+90CJjfb29nPnzq1du7b3zvnz5xcXFwshSkpKnnrqKZVK9c0334wbN27RokWzZs3auXNnZmZme3u7Uqm0HvLQQw9JS9XZXwZmgKwPeb/bCr82+zUajb+/v8FgKC8vf+aZZ0JCQoQQX3/9dWBgoJ234uLimpub16xZc+vWrVmzZr322mtCiPXr169fv956ZukpENIL+4N3AOCRVKN8l/zCmds578hkun2m9OTeoqNvvb6279b/tLfov5VK5YIFMxQKh+/2/bmxwQEDPwnuyLEoc88991RWVk6bNq2hoaG6ujowMDArK8tO+08//TQmJsZmKuu8efOkKGM0Gj/77LO0tDQhRFFRUXx8vFqt9vLyEkJMnz5d+u9fMnv2bLVardfr582b1/tU1nt2jEaj6DUTxc6LuzEYDDf/eceddRRM0tXVJb1obm7uvT8gIGD58uV79uwxGo05OTkLFy784Ycfuru7N23apFAo7vaWEOLy5cvSzeEXL160XxUAjEyRY4P/9qc1Ljzh23snbdj2/skz5xYlze5P+08++5/z31cX5b244jeLXVgGBkN/o4xSqVy7du39999fWFi4Y8eOpqampKSkjIyM3vMzbLz77rvHjh3z8fF5++23V6xYERYWJoQoKyvrvajMoUOH/Pz8vLy8EhMTDx486OXlVVlZOX36dJuEpFAoFi9e3NDQ0HvirVar/e6776TX77///rZt2woKCqTNmpoarVYbGRkp3ZsthPjoo48mT548fvydF40+ePBgeXm5dfODDz6oqalZtmxZWFhYZWXloUOHpP0nTpzw9fVdvXq1teVjjz0WHBx89OhRnU735ZdfLlmyJD09XZriY+ethISEGTNmXLp06emnn+7j5w4AcIU1GY8d+aL8D1t37d65MV4Tbb9x+Xf/tz2v6Ne/nLv814uGpDoMiMI6qjJA0pjI888/n5KS4pIT2vj8889DQkJG+BOhc3JyKisr+/Nk7LNnzyYlJb344os7duwYmtqAkcBsNkvfzoaFhdlfKxzDUEubPvmpjbqGa5t+n/GrR5Lu+CiY26bbfyku3bnrrwlx0cf/8scAf2du4cYQc8H4X2/WucAuUVVVtWnTJmmWzPfff29d23fEcu2PFwBGlLEho8sPv7V47r/lbt+74ndbjnx+uqPjpyXHWq93fHz0xL+v/o+3/vPgb1OTTny4nRwjFy7+VFFcXFxdXe3j45OdnT3wZx+WlpZWVFSYzWaTyRQTEyNNoxmBKioqSkpKbty4UV9f7+5aAEDGgoP8tftyi499lZt34NU/7hVCjA4KCAxQdugNXQajEOIXM2P+689bf/lgv+bTYJhwWZTZuHGjq05llZKS0tjYKC31O0wepu0WcXFxcXFx7q4CADyBQqFY9siCZY8sqLp0pfT0uSuNLe36LnVI0MTIsJSFidETZbB+IGwM67He2NjY7du3u7sKAIAH0kyZoJkywd1VwAVG6JANAADwDEQZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgY0QZAAAgYz7uLgCD6PTp0zk5Oe6uAvAcFovFYDAIIVQqlZcXHwWBwZWcnJyamtp3Ows8UUtLywsvvDD4/8xwV76+vtHR0X5+fu4uBC42duzYiIgId1cBF/P29o6Ojh41apS7C8G/yM7O7s9/eQqLxeLuUjEo2tvba2tr3V3FyHXt2rVXXnnl5ZdfnjBhgrtrgSsVFRV1dHSsW7fO3YXAlfR6/UsvvZSdnR0TE+PuWvCT8PDwqKioPpsxwOSxxowZk5iY6O4qRq76+nohhEajmTp1qrtrgStptVqLxcLF5WFaW1uFEFOnTp09e7a7a4HDGOsFAAAyRpQBAAAyRpQBAAAyxlwZYFAEBQVlZGSEhIS4uxC42IIFC6T7seFJVCpVRkZGeHi4uwuBM7iDCQAAyBgDTAAAQMaIMgAAQMaIMgAAQMaIMgAAQMa4gwlwDb1ev2/fPoVCcfPmTW9v7+eee06lUvV5lMlkys7OfvTRR5cuXToERaKfHOpN57oeQ4+L1FPxrQzgAgaDITs7OzQ0dN26dRs2bFCr1bm5ud3d3X0eeODAAR6VNdw41JtOdz2GGBepByPKAC7w4YcfdnR0pKWlSZtpaWk1NTVHjhyxf1RFRcWlS5cGvzo4xqHedK7rMfS4SD0YUQYYqJ6entLS0mnTpvn6+kp7/Pz8YmNjv/jiC7PZfLejjEbjxx9/vHLlyqEqE/3iUG861/UYelykno0oAwxUdXW10WiMjo7uvTM6Orq1tdXO57nCwsLly5f7+fkNfoFwgEO96VzXY+hxkXo2ogwwUDqdTggREBDQe6e02dLScsdDvvrqq/Dw8ClTpgx+dXCMQ73pRNfDLbhIPRtRBhiozs5OIYTNrRD+/v7iLr8l29vby8rKli1bNjTlwSEO9aajXQ934SL1bEQZYKCUSqUQwsvrX64mhUIhhOjp6fl5+/feey8zM9OmPYYJh3rT0a6Hu3CRejbWlQEcYDQae9+WGRQUNHHixIiICCGEzdOSpc2goCCbM5SUlCQkJPAA3mHLod50qDHciIvUsxFlAAfU1dVt3rzZujl37tzNmzdLvyWNRmPvltLm5MmTbc5w6tSpysrK/Pz83jsLCgoKCgoyMzMff/zxwSod/eNQbzrUGG7ERerZiDKAAzQazdGjR212RkZGqlSqurq63jvr6upUKtV9991n03j16tVdXV3Wzba2tvz8/CeeeCIhISEqKmqQykb/OdSbDjWGG3GRejaiDDBQvr6+ixYtOnnypNlslgbXb9++XVVVlZKS4u3tLbW5fv16aGioECI2Nrb3sc3NzUKISZMmJSQkDHnhuAOHerM/jTEccJF6NuY0AS6Qnp6uUqkOHz4sbX7yySejR4+2riuq1WpXrVpVUlLivgLhAId6035jDB9cpB6Mb2UAFxg9enReXt6+fft2795tMpl6enrefPNN6VZPIURoaGhISIharXZvkegnh3rTfmMMH1ykHkxhsVjcXQMAAICTGGACAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAyRpQBAAAy9v92niGl8OHyIQAAAABJRU5ErkJggg==" width="60%" style="display: block; margin: auto;" /></p>
<p>As can be seen, the age of the tree has a non-linear effect on the predicted circumference of the tree (here for a out-of-sample tree):</p>
<div class="sourceCode" id="cb14"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb14-1" title="1">nd &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">age =</span> <span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">Tree =</span> <span class="kw">factor</span>(<span class="st">&quot;6&quot;</span>, <span class="dt">levels =</span> <span class="dv">1</span><span class="op">:</span><span class="dv">6</span>))</a>
<a class="sourceLine" id="cb14-2" title="2">PPD &lt;-<span class="st"> </span><span class="kw">posterior_predict</span>(post1, <span class="dt">newdata =</span> nd)</a>
<a class="sourceLine" id="cb14-3" title="3">PPD_df &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">age =</span> <span class="kw">as.factor</span>(<span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="kw">nrow</span>(PPD))),</a>
<a class="sourceLine" id="cb14-4" title="4">                     <span class="dt">circumference =</span> <span class="kw">c</span>(PPD))</a>
<a class="sourceLine" id="cb14-5" title="5"><span class="kw">ggplot</span>(PPD_df, <span class="kw">aes</span>(age, circumference)) <span class="op">+</span><span class="st"> </span><span class="kw">geom_boxplot</span>()</a></code></pre></div>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAu4AAAHPCAMAAAA726/2AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAABcRAAAXEQHKJvM/AAAgAElEQVR4nO2deYAU1bX/O4lJ3svyiy8vb0t+eb+X95K8mGd+7/1SM+yLiCBEatj3fVUQREACKohocEFAEES2YREEEYiyKo6I4IIILqhsAuICsorDMizD9NSva7nn1kyfc2u6p7uruut8/5iemnP6VtWdT1ffunXuORGDxQqNIn4fAIuVOTHurBCJcWeFSIw7K0Ri3FkhEuPOCpEYd1aIxLizQiTGnRUiMe6sEIlxZ4VIKcW9fNSTqWyOxUqxUot7pHYqm2OxUizGnRUiMe6sEIlxZ4VI1ca97A35O+POCrYSwv3A2M4DFl6R26P1mN6S24w7K9hKBPf3+hc+2U2/H7YP3Dlz5szCUunAuLOCrQRwL7/3vGGc7avvEX+YdLiyB+POCrQSwP3AFvPnWn2ds32q67ojFT0Yd1awlQDuJVHz5059s7O9IDZwH7jH7cG4s4KthGdmXmx5xvnt1NYFvfQCZ2Km5EhMX36PcWcFWQnjfv9s10bZioIuxdZvGzVTf2DcWT7q7HkPh0Rx3zv0SoXtpbp9eT+7J6bdf8u4s/zSladu1rSC1VGVT4K4XxryRcU/lBQslxs8dmelUd+8vWFPGWmN3m4NMLRxqiYSxH32jsp/6bFB/s64s9Kmq1NqxWDu+D5lX6052qZoJDHcN7xa+S/RNh/KDcadlbSuLr+j9YBZJaT9IRvm2nsI+xCB+wOKnSSE+9bV5s9j7xnGafGnnSPKpZ1xZyWrkj4Wq80OEfbP8hyaBxEO7QXuAxV7SQT39/sviWnO4EvGOn29sadP4Rnj2Dj3WJ5xZyWriQ6sHYhbzeWC5hpXcYfewuHPir0kgPun7XVLMw1jW8/tRvGEnn3nPldhnoZxZyWp0gaCVmJwPlfYtW9wh5nCvlKxG453ZwVBRzQPWuFOtB5x+f/mZtve7gput8S4s4Kg44D7C7jDqdqO/V6qiUMtTXO346rdMO6sICjaVOC+n/BYYpubkjjPsW5maz2v2g3jzgqEFnlMvMSGM000Le/Wo5R5h7iV/VSxF8adFQhF/2LB2u0M5XDMmmms+RRlv1t8PTyu2AvjzgqInulxc4eJpZQ12t2heTXh0EngPkSxD8adlRl9/lDvvo8eIc3R+yxWe54m7O8ImtsQDq2EQ3/FQTDurIxouRnxotWmrs1GoQPrrYR9NkzdXMQdWgp7L8VRMO6slGhP9wY1mzxKBizuznduJIkggWgTQes+3OEOwP0s7tBC2PsojpJxZ1VNJ3YcIh7fm1pvh7Q0PkfYHxAwPobbjwHNxERiG3D4Gne4Udhb00fJuLOqpI87m3PeaynzxZoOa70Jh54CRmKwchRoXoU7SNw/wx3EEWg3KE6DcWdVQfvr2SjNJ+xi5K3lE8MZCOAi4hWv1hAO23GHvl6DGbDnKc6DcWdVQQMdlGoTEyfDAbZPcIfxwj4Jt5dAAy/jDveCAxEzA3ZNcR6MO8tbl/IFSutxh5FeQ43Bwj4Kt++BBkbgDneCAzF2Z9xZKZK8kXwad1gp7DXLcQe4kWyG21+FPRDziP3A4STuwLizUiTPq3tU3ChSz3jgRpIgZC/AOhx3gHtd7RLuwLizUiVxo1iTGLt/LpbW9SAagDtRgpDTAOs03EFe3U/gDow7K1US657rEkunVwnU8ggHuLrXwu2fA6wLcQc5M0OsZmLcWSlSaR2B0ibcYRGwRlx7xdVfy8ftX0ADxEwlzGRqxAeKcWdVWS/2u6Fhv9cpq1xrRNyqvgxjFeLJK+BOzIq/AHsg7mVvBQci4h3sxAfKEuPOimmSTcpEwnwOWCJG1nCn2ZBoARogcF+oebQwBhyIW1UYLt1IHIIpxp3liq59C7d/DawRqyvmA83EUMNrMANfD9qfcIeHvMYqEGUwgHAwxbizDON+gcoY3C7nTZbgDjJe8QvcAXCvgdt3QgM67jAUHE7hDvCB2Y3bLTHuLNe0BzGPWAqsFeIOMoCLCPAFOzGYWQYON+EOt4MDcTMs4hRmE2ZLjDvLxRIxEJDxiuNwh47gsAt3ADsxFnkW7E1xhyeFPZ/MaP1mN0279SPKaolxZ7mWCj2J2/cDjP1wh9FeQw0v3D8F+1Dc4RvxYFcxNP+UDJcUYtzDofKju4pp6xknpKUREX51CGDsgju86DFW8Z4Vh8euVJ6ZZ9SHaIpxZ1laY2aUu43KrmsYu28xUWr+IWGWES0dcIe/gMNl3METdzG30408xjfa5Gm1RpKJOQzGnWXLuRGsTwSjx3TlladmFhGoGsbHXoMZyHqhfYA7eOJuTLMmb3oTs+rOXlYorIw7y9LZug5qdOrz6PZnlmwnbwIvAqxjcYfW4EBUG/DG3Th2l9ZSNYloGJ0Zd5a3igRq+VRxjEPtTHP3g4RZPmYiVu9BLQGNSOFYBdxj4xlVihiDcWdVSXKWj3gIdNZJFt2UuJ09BQ08iju0AwdiIpBxZ2VKGwVqeUTdUVhZPQe3fwWwtscdbgAHYmJFpKvW6tGHybizUqFD4hF+A8JB5F/UOuH2l4Dmm3EHmEakljtByMtM+jAZd1Yq5Hl1h0oxxLX3A6C5Be4AITHU4F7kb29BrGU1xbizUiHPsXstYa+J24uhASK9rsSdyrx0ZUTsG6DGaDLDr8G4s1KjFQCj19I6Il7xS2iACJmsBw5F5FEcJSMMHDHurFQI0mJoxMyLF+5yaV1X3KEbOBBpMQzGnZUhTfcazMC8CbFwegc0QNyqvuHxeTHFuLMyIpmB6zDuANmmiYVvcmaG+gf/ybETeQRMMe6sjGgE0LoDd+gi7B1xuxz8U1fvC1bxjLzJiqNg3FkZ0WSglVgJBKUAmuN2mVeDmLqJ6b2btWHUQiNLjDsrI9olYKWi0eEpEbFwejPgXp/eS0/tGeVRMO6s1OhKYc9GrceTBXanAq3EzAzYiYiWb8A+jD4Ixp2VEZX0sFCsR1QCkDEC2lbcwQt3A7KIKWr0Mu6sjGiKw+LNRFE6GY2+GXfwyvEF+agHK46CcWdlQuWNBa1EBkeZtGIv7gD5rMmEcy+Z1/e8seQCEINxZ2VGMqKFiM8aBw7EumbAna5rFJ2qtaMWh9hi3FmZkFxaR/DWFRw+xh08x+6GmQmgp/owGHdWRgRJjYhMAs2AZqKKI+OejBh3fyTi2amV1zrQ/BruwLgnI8bdJy2wBt99iLoWUCaSjJlh3JMR4+6Xvlj68Ow3yWkTmdLuc9yBccd0YGznAQuvwGbx1GkTJ7vv1hn3YEomPXofd7hF2InSGaZCh/t7/Quf7KbfLzYv9FtsGItGuh5tMO7BlLxVXY07fCLsirRGYcO9/N7zhnG2ry7yRM1rV2oYl9ssc3kw7j5p3YBbOk0hyiK5s69TS+uWW89V81S8hg33A1vMn2v1dfZmaUdr3eLo3nLIyLj7o6gTM0Ndm2UZLzIL8NfjNG2Mksaw4V5icb1T32xvfqxbz/Dm6DKRDuPuj+5zYK5H3KxCXg3tDbKNU2RJO0dhw93Wiy2dlMPr9OfMl6W67ELG3R9BuPoy3C7TYowk22DcUd0vat8stUc1q/XnzZdXG5n6PePug2QGR6K0kpxnpEtfMO6Y9g4VE5HP6xvMlzX2Rf7LVTGt/Cnj7oP2Ac3E2jvImqQtIhth3BFdGgKpG7bp1r6X6S+BlQczvkjWtCZqksrHTFdwB4NxRzVbLmU/rFuXikK+VfVbJ4HmVrjDFWGnUqIajDumDa/K30s7PmC+3NexDP7EuKdJW4c0azJwG2WVhcJuITxW2+Z8RRmvbMW9/O0i0DOatkZuvYe0kBDuW62Hcsdi7ZyOvT7VMWoYZe1dhWUZ9/TIqQP5GGGW2dfrUi18osdgv1116ctW3F/VSCElXhPB/f3+S2KaM/iSsU5fbxjFvVcZxvJ+F6QD454WfSQmEolZcyiPrsgCY/Snljo5CiruMwZK1dfayo3B9tV7pVarAFU+tjA3Adw/ba9bmhm7Te1pLnovnjprxuPuoFPGPUmVfFC0jwwBMB4RNI/C7e8A7uRS06zFvYS+eNvJiFdSJRiaVBN3bzHuSal8Yf3Yf68NlTZDhqsT/1iZsDT7ru7lR4+Apmo95cYRu+rlOU2bshBTb+1uy4FxzzY5lZFqULwPEzT3xO37shh3mfGpsuyCxTHc8eLHjzHu2alikW66M+EAhcKIhKOQEy8LcR9CD1YsO+Oea9oE/+DTuMM3Th6Z+l/h9rehAUV29cDiPu0sprcY9xzVX4FWqorvx81Na2NqsPMmNECnifEL9+N7QFs07S255SQLHqLNQhvbxbjnqF4HWKl58QvW6ruuVLD6huDivgsSNlVWDXuVEOMeNkGdL6JSjBFtaNvrXMbt8uoeuMHMGnpovtFyYNzDpnVAAFEN4EFhJxKSyhCxm+i9+IV7S3RofrY54x5SwVMkKk8AFAojJl7OQQPj6L34hXtr3PkWxj2kmgS0Ejnv5ACAaKGRsL9F7yU9uJfslTefHbQpcuOgXR2bcWdVkgx5IQIWPXFf6Jjb0oEI6cG9vIAcmk+3HDKBe91uqGoy7oHUWTF7QSzO8MY9erdlbYaDYSstuF+Jx1xohOWQCtzxIII+gDupzfHNMu6+64QIeCRK9FYlpd3m3lrDBcqlDenC/V30TvSxFOJOiXHPSm2E/88x3KEKuBsryBAER+nCHa8HMo1xZ+F6Hv4/xGAkt3EfiI5VHq467g1HoarNuAdSEPKSTxSD8Q/35eRDUTsyPRW4k7LsVcCdZ2aySjsBd+KpqX+4j6Bhs3IaZAR3nojMKa2GfzBRtLQKhcLShvt9ezBtqDruNxehasy4h1Tr4R98BHfwE/dpqPPequNOinEPp+Rgpgx3qEQAquzFnZ+qhksyfpdYvsG4I2Lcs1RTgAAiP7tM4Es34hvuC9Gh+XDGnYVrJBCwA3dIG+7RVXKi+0lNe0JurbCjb6qAOyWBu34EVTPGPaS6DRDZgjvcK+y96UaSwn0rTesrlkMKcOeISFYFtQdE1uIO5fVsc01FwGNSuG/Uag9EVVdbYzkw7iox7smoMSAyj/C40tu0drhAmE0liTuRQ7V1NuHeHl8v1ZhxD6TqASITSJ+zdJoCW77hvhNlbWIKcd+GDv7HcYhYduom+P/MIH3Oko9cHaUL984TMI1MYRBBNSMi38mj7DUOxDfLuGdAr43vPWzBRcraFf5BO8kWfMOdVIZwj7Yhj2CJ7XlUxja8omk74jLZVBDjnnZFx1r/neb7CTtkPFWUgcxh3CeiY5VXxKzrVdcQqYO2WG6cj280DZX3VGLcMT3jANKSKIzUBRBaQ7aRw7iTin+TD4UmVWLcMbUW/7/XcHsL+Af/lWwjOdw3uobbN2s95MbD9jdNFXBviJYKaJ4y3GeStHePf1O6cI+dy7GN6jdiYtwRXYZ/4BzcQU5EziQbSQ732hhGloZb9lTMzHjh3ugJVA0c3I2TRH73o0i8XHpw39Xo+7Gfc1p4vDVejDuiKzB1MBd3kLjjd22mksM9hjU6sdJVG2LZ04/7evTDZumVuDf5U3nvy2sj1l+na8TqGlKMO6YbxT+Y+L4Eu/Y02UayuCPVuGKaBbg32Y6qeapwP3uPfFTbR9P6ya2x8bea/uDe7XvdfmG+Xv2nKer3xolxx9RQ0EzciUIOMIItU+nCnVSKcHcroHVVf1Vk2MtmGvxR/d44Me6IZDGtp3CH5uCwkGwkuLhvQecRJ2QP7r80HNyv+5H6vXFi3BFdhbF7Ie4AAY/aHrKR4OJOKWtwv6nMxn1L5Dr1e+PEuGOqIQgg7kQPCHuNKNlGQHEvb082ED8NFVDcJ8+2cF/7s0g/9XvjxLgjklfAuwiP7o59Ot1IunBviiYa2HNLFXE3Sl0DmM7aTLmBPMEPKO5ldVv9/KFBv41Efn5S/d44Me6IzgDuPQmPEuu5at4D5XQj6cK9mhORFeRPkXi3kpt3Pz/imkhMrT9RvzVejDuiq564G+XbW2qDVOl7/cP9MTS+991cwj3Wu9ue2Uwsi1eJcUd0wXMwE9NgbbayERz3ixLAxVoHF47O90QKcCeVS7ibItLRqhRa3MuI6C9TUSCEKK1kKincl5Oh3p3tm15Na4aGvNxQVdynk7Tr8SOvbMV99PVm1z54H5Hlh1ZIcV/dpVZ+279S0yrnARHFyuqkcL+bvvaesxxoexVxjx6Ut69TtA6um1kke2uW4j43EnnZfB3VUXHvhCqcuE+2CRpG8C7X43Sg20gS97vReZVXUoa7W+mpmu2WP7jXqLPKer304wXq98YplLh/IBAi8ghcAsaa040kiftjqPMhxp0Sgvv/iF+uU6TxQRVK3B8VCBFD8xPAGFWMxkgf7mPR6Ns+jLtUM+e19PscRFAFDRU0EwsZLgLuDehG0oW718xMDfRWtqBmiHBve9h+nRf5b/V74xRK3OGGsStulzWtqSLwhl+4Q92QeMUnNMtV3Lf/0UxpUr74B5H4MB+1Qon7CkHIZNz+BSBUg27EH9zLt6wCLdK0Qrm1OX6WIldxN0Z+96bbOvwmEvktsthbqVDifqmlDXMjLNGD4b66K26F/MHdrfQUiXcrqLgbT/ydGUTQtUrPVctc2SJCibsxzWZ5GGHeB7Qrim8w7oaPT1VLX1tShD5V3T7cvTVaj+ktuR1K3PcLmtfj9iOMu63g4u6o8qOzvTMK2ro2D9w5c+bMwlL5h1Di3k3QfBNuP5TGW1VSjDsuJe53Vtq+bAxz4z7pcCV7KHGvKxDLx+1lwGAbuhHG3fAN96OPD73N1C1/E2e6x4X7qa7rKpWKCyXutbzGKpDt5VW6EQz3E60bgW7QtIZy6xY7WJhxr6SkcH/5byJCcTY37gtiA/eBFZZXhhJ3yFdNXN2NTY69haIRDPfXaJrtKI8q4N4drVagM+5Sf4g0GTkmptEFP4izVbi6b13QSy9wp/EMJe5NBWI1KY+11mrVroogYRT3zdoNeBaYAm2l5XC31mchpikC9/rkx2FU3N5Ci/u1sNP742xu3GMqW1HQpdj67e1upn4dQtwhbwZ97mXTtKbq/zSOexPcuRPg7jEz84HrM1CgDZYbi4/HvSm0uN8yVfwWn0WsEu6GsVS3L++fWLFH/xxm3OvQPqtVt6mm0oO7W0kVmnQrV3HfCf+a5+JscbiXFCyXG6EczECC3/q0D+NuBBb3stkLjlvaEw9vHO5Gjw3y91Di3k7grgh4ZNyNwOL+i0jVZmYsRdt8KDdCifstPJixlaW4960C7qfFH3aOcMXOhRJ3yHhKTUQajLulgOK+6vZTZaau7m0WZxvaynpZp6839vQpPGMcG/eFyxpK3OExkyLgkXE3Aov7WZhKr1waa++CAn22GRK2red2o3hCz75zn6swmRxK3PMF7YoIsCRxb4TPu7cE3Meh6Xe35Sjup5trjyizY3CxmvQrL31Xd1ICd1I5iPvVtuaJ1VxOOkSfrKNpN6ix5WI11ZV/uD9N2lvHJwHJetzF8w2qXlu0lW1XfkNwsZrqyj/cjVNyADNXa+cazlyN30MW4D5P60MbXxHnTk2ATREOdLlOLlZTfSWP+3mZ1PE27Qm54dTX9sbdpaTqqrqVdtyjmzTtTYV9hznFVZu8F+0JJ0+UH4c0+ornH1ysptqCW9VEcacr6NqlAGK3qnj29VbZiPtBq+Jaqy8p+1bn3McRdlnB6kPcQXaf4ii4WE11VZVuRnGnq54OsOyeE5FuBR334pr2qdVGkktagt4gBjxyvpc4zGRx52I1CSnpwUxtbSV68X40F3EfJDqJyPp9DHpxNO5QBxx24A7J4s7FahJSNa7u76DOhdmJ+9m22hR6Vhwu3kQg3XLoRSK1oKw+S2RiB7vi6TYXq/HWkfGt63R8irhBYtxtRSeZ9zBN36XseR4wLoBevAF3+BM4xD9TqLgHFYNcrMZLe+ygmJbUZ7+J6GZyNVMYcB9p90He+4QdbuiJVGqvA81xMbe2IDRJI+52YQ+KhA8Y7l8XlXKxGlBUBPiOJBy2iW6eRzeS87hDZsB6hAPASFzdZa61+CV0luTYnfhEwR7q0oeJr1VdpniDSrmI+0dwWSK+RMWz/G50VdTcx/12gJGYeQE7MeT7FOxjcAfIb6J9jjvA54HoNUsI7v/YSEQPUJNGlHIR9w3QzXspl7UtNO3GRQrak8XdY2m2Wz7jLheBb8UdvHB/B+yVkxs5koMZ4lYVlsgriqRguC9YKn6bGm9UKhdx3wTdfJj0uaxp+5SNJIW7Z+INt9KNe3TpLVqdhy9RZnikqT2LO3jhvhfsxHSgHMwQt5RQwPt2xWlgt6o71ttXquOKDM2ochF3GFTWpC/fyeL+Glq0dKaDO51WqUV8DdY0436lrdUHdXYTdvkQqAh38MJdppon4mYgnQ81dr9V2MfTp4Hh/q/XfOs718T0HWQ1k1q5iPt26GZq7J487pQGxPkmV0bYrerhfqdzZHWJmXVIx0B1hBfuW8DeEXdoBQ4XcAfAnbjXtYQQ3UOxeE+tXMQdAu2QghZCOY97FCa144dRlkbCsZ/GHeDiTEyry6s7MRYZBw5IuKepTsIenyFNCiF6cZ+9VhTplx/9VvFGTLmI+2PQza+QPtmPe3Rs7GhqT6LGawfh0IihxmRwIKY3xgj7I7i9GBogDnMuOHyDO3QVduJe1xKC+zGomIgvhaSVi7hDnLUiWDv7cW9r77gVUUn3DTg0oljmveBA1DDZLezEaZTB98d23OEFYa9NfCbhC0YVts6L9zwEhfMUeTWSxX1WEaZ7Mo87XDpn4vZ3gebGuEMvcCCeeU7wgPEoNEA8rSsS9obEOSSNOy/ec2kH/B++Jn1Q3E/OkoVM79AaucqabrIdvCYi3Uoz7jBtTjyRlCNrIoBLLr4gPvcw1CDmGXdBA3fgDovAgRgudRb2objdEi/e8xBclrRtpA+K+yMapXy7xFuAcIdDI7IpfAwO7XCHgeBADGZuFvYC3A6hGFpv3GE+OBAzZDcJe3vcbokX73lIToARX/QGgfv9WnM0uXofMXsRRNyJacJvwN4LdxgPDsQdDixGaorb5XBpEO4AM5V/Is4BogyIokGWePGeh+QNJXHZMUjcJ6DOR7IQ9w/BTlyc+4IDMREJMQAEjLL67HTc4Yru2AuJc4B/VCPCwRQv3vOQfDxOFIE3QoA73CdSY3dIlKl9gDvc6NFAKTyXfY84xtm2uR4ViA0z+/G576R48Z6HZPATFQGcFbjvb62NPkNaYRaQGLt/CZ1ADGYk7m/hDh2FvQdxCOJWpw0xz1gs/hF/IRpI9laVF++5JfMFJPiYKUi4X7YCdGuRgd0Qf0X8/8rzPTphGPQSUWl9sLD/mTiEkt6W+aaDhH29aKAh8WxAfAPlfUS0YKoC7uXTzc8mL95zayr8I4nQViN53Lui97KtUo/7AOcUiBAAo6Y4R2qdp3go2pC49h4WDVBJXmDqZgR1jKXPDi7oO42c7Z0H/4diwmOsbZ5BtWCqAu6rIg3NF16859Kr0M1EnLWRLO5ysXFlxU89Vw/3A6Jh6lGZ12AG0lHeQC2BGOo4EGkC5Lx7X/owlVoGn0iqpFv5mj71brhts7KVCrgP/d/m0Q6K9S4v3hMqEfdQijqRyeH+juvJk671kRtPHo57U/Vwh4AVLb4EmSUv3Eth2mMTsYutdhO3UJdeuLoPpw9Tqf3V/bxYqoD7IGvBzv/AtrqSc7xyEXeREqIGEcthisB9HBrOvg+brEuqarZbStzhGY+2AXdoIOzEU1W5knQO7lAiFqhTN5LwUJTO4Osh5wsmf2eyDZiqgPtiaxZJ4u6RXT5O2Yr7WSKm1NZK8xGJXjnXvVsE7qQyjbtc+LYId4B5EyJz62logMAZhnz1iYD4486Ebk0intFbl+41v0AaEctHqqgKuJ/5PzW7duv2k3a2WtcORbz7uYdv0mr2JBOkxHRlgPaAaiVqAHCPTtN0YtlDTI1hz5txB0jiQsyKyxxf03AHmXubCCKYI+yL6ZPw0mfrFr+e6PLpSqpI9L6bvhVxK8HGshH3s+28v2SHa08o2/Ad9+fNa2fecCqJl1xrRKTghUf8DXF7CUxErsUdVsAeiInIPsI+WHEa1dTlTXPmv6G8LsXNux9aM/cXU21NHvHjBPeXjbiLhQm1icuSqaDj/qLTLvXcVz4EIsbucHUnIlqM/o69BtFLMPfTiWgAbh8Ud/zV1AfNzPa7x6/jdQm5gMux+90J7jAbcYfpQMXlPUnc+63CVJh63CHQgbidlreqxNV9tLATCUuNj5ypGTJMTmQRe5Www3hK9Yi/Wjru3G83pwd1KO5z4bf4qtlqZSHupUCCYoYsSdw9JiLdqhbuJXAOOu4AKVi0l3AHsfKaTlqxo2XMWm8e8UQzdgzWc9XaJDEFYg+q8NxqCRYVL1E4KYfnoxLcYxbiLqOjFOVkAo67fBJGTJvLdKLrcAe4/NMZuKIHN+1S3ifufvbJdUQ4pOH6/lClxaiWenh9Q5nCcD/+xLDu3WJq+v0E95iFuMsrI5GK01TAcZfZc4mAxgFgJwa2ECtO//92zfjzoy+TF3dPQbg6HieUAnUQeyAC5i0huG/4XphmZuCBIrFozBSGe3SJfAz6uKbdL7fm2hfBzOEOz9cp3B8GO3H5hehbKotx+aOWuZe6lphKzgKQRFf7V13DxTkQuQ4sIUT//ntdxj0S08P3/CzBPWYj7nCbp8gQiOEul9/EyQ7EyhzuG2HPxGDmPnAg8lwC7lSy6KWOvV/S1/fyF7rXqddrfbJv9xZERO5SOCG4/xSGV6pREKYsxL0MSLiXdsJwf1OrOQpVI6dmS+ZwPwPnQISADQWH/bgDBBEQ+arLYSpTFV3rpah6TryaKncig5TjTgT3FvCkOQQZgKNAgiJxLI47EevaI+O4Q/QUlRZjFDgQ0eRthL0lbpcfqKRDXtKv8lXd69Tvk3DV7GfhwRfxCI3eYfgT6aQAACAASURBVPbhLiciiYRBppLEnVSKcd8HDRMZ6WRE5AHcQYxVqJiak9CAulKw34p6jbWwu9E7nSUrF4lpXFJZjbui6ENSuMv1CJXVOD5i2wv3E5pGJd81jEPQMhHy8gQ4EKUAjO62uQthjsLDOEVcaDYIwb38ie9ea+on3w7BzAwkCVMVbEsK9/JPZO3IWVoLVyVJZM2oGveSB/Njt2CPUGmuzsE59MQdYBawBpWgPTqlYV5eAzJHpDHNaaCNMnY0+EKIvitMIWIyFaeicF5SuLuFF4l3SYl71Elr0ZFYxyPPgXiGUyJCwOjV5YZxiaxVENNlO510Y+JWN2uEEP13v1zywb6Y9r793wk2Fkzcv3nnTXpl1tdZgDusSn4at1+GUxhGtPCyHdHYNIkVao7KVt3Zuv8TyU+7B0RYOQNISknnzcIVRNxLbjUfJHWj/tEnkse9Hlrzek/7lOMOMQDEB+wrOAVyicf2XjW1+qPp1bZhEYL7gERLMoECiHvUiU2qS62iAVQUh47jTirVuHt9ImX2dcVkaumJ5CMAckcI7ie7it9yoBTZM4IEKtkOBBE8STcScNzl7XaiEduhE4L7oUdmvWFpE5WGgVIAcYfcUvUJh+ccu6rWcsBxlzH7VMI5liMsJWouzczAujVynvFJ6/reiMoYYSrouK91zGmLJc8ZIUQPivzs16b+4+/jjdvdiyCKp06bONkdaBBA3FsLUuir95X7tAYfKhvx+VbVE3cxK95MNZfIMlDc1z3k/FJeOcXz3hkFrqjwC/0WG8aikRflXwKIu1hkqUqD/CK1CkjI54lIL9zLRZIXKtEASwjBvQQSn1XOinjZGObCfV670tif2rjybAYQ99tzAPc8D9w9q9axhBIcnt8jcS/tOMZ8Gd1bPnkOIO4Q66c4sqDjDlkviKkDuSKrGgFc0c8/p/J25JBUuM+IX+rlwv1j3bqWzNHlg+UA4t7We+weeNwhHTV1u91eOKhWNih1eXI9Tas3O9FiXFknBPcZjh78dXzuYBfu63Rr2flSVzq5AOIOmWeJhQumUNzv0Qjl22BmEHcoD3Uj4SCiDPon+yQp6mQsHZDWBRgBEIK7nIe8Jc7mwn2pbi1rX60/b77snmDqF4HDHRI+UKswDQL3GyjcnQIWqcT9JTOLS3MytTmkxaAqXxizrfFOd7oUpocgKmdNsi1kiRDcv92olSn9V33vibO5cH9et1JSrbEv8jutTPy/ChzuMBBQPDIjcF+GzjNOBtxTNhEpEplRQxFIwEgkaIzp0KLx0zcnf2n+s9hDsumos0UI7iLw4lzr0jibC/dt+grzZZkuU/X4MZg5NbFzo+5zyUGn16yGKQJ3POfvEsCd1NK4NylxhxgAamh+3nlW1ihtAYm9xSF0S9ceAiIEdyi/c1t8+LQL98O6tdKr0N9b1UP2nHMrqq5O2nA/QRbfqB1/lVbiPgjeSYVtHrLml1pUZ1m0WreJI6hWrYAskGpmZvxv4v5UYSLyAfPlvo5y/soH3Lt7fAvDJF6qcTfKZIWCk5r2rtxCFmEocYfCSBqZta3szYWFW6giLSnQrYy78dVvro37m8DdXFz8VMfYYLGsvausa+Zxh9jXPCLAt1bacHcJzSLmlhJ3+b1Apc9Nu+DZM1FHMmeE4P5rW//+3UirONtQ+0/r9PWGUdx7lWEs7+dKuJp53GVGIeI+r07wcZdfQGRVunQLcKcrg+eGVBORNb6sZNm7oECfbdai3NbTXJFePHXWjMfdV9XM474QUNmCO0CdlsRnZjKFuywlM0vZSBr1qDgCPDVO7gjB/XfP2uHuhxN+aJF53P8KqBCpNiFRpiKvBoF7Z7TqaeuU494dzoEoNpB+HXC+YfLVp5H9QnBfmHRjmcddZmokMqjAaiZq1bKR5GMmt6qH+yRo+QtlI+mUVe1Gq7HStwPIkOhb1STW8WYe9ysiuyEV9VLe0rm4Kyq++Yz7+ADgbnzyaP++j2R7Wg1vYbiPvt5M0PbgfQlHyGUe96gI9SYXJZdY4SCtK9+GuOUz7n2h5b8qG2FVX1ixmkjkZfN1VMdEB++Zx303oEI9Z7JSyWxWNuLz2B0KSpMlqVmpEoJ7jTp2gvJLP16QYGOZx/05QOVt0ieqaeq1eT7PzEAdFz8HMyGRqvLedYq5alSZx302oEIU2TKCgPuBu7RGReRXpSw2QEymslImBHdRCrD0+z9KsDE/r+50zgnfcZ9pBe50Ok6YH4NzSF9QDMsWgnvbw/brvEjwc0TOAFTeIn38xn21c4QtiQDdj8Up1OQ8X+kWgvv2P5o1HsoX/yCS6EPtzOMOWZO0caSP37jfJA6RKrFbz7HT0bfR91cs25HrK40yIWwicuR3b7qtw28ikd8S9e1JZR53SDmhDSB9UNy3tysA3aTVkBsFI22sUoa7LJhA5O2DWjM3UMnTv7A+1T2I0husqgt9zPTE35khM10TTo+cedxlvEl30gfF/SGNlD3GThnuMqU2UUpmLTgQT4ZLdNvcBKmDwEpICO5fF5WWvrakKCueqjYGVPqQPijuE7T+RZhWiUUWKXvMdBTeWAN3WOOFOwRCKLK2sqokBPc/RJbF/7FKyjzuTQGV20kfAvdxqPMJgbu8K6is+DS7Stw/hTcSQZl7hb1B/GJJSyO8P9KsqgnB/R8biSJAwS80KXEfTPokh/vZd7aD+mh3y41349fFKnGXZSCJtajlYqUolV8cHrsq0rezqiQE9wWwtDj4+d1rA0v0Y8vkcHcLS6vklhL3K3CIVM7tr+wwtsHU6jx4DkV/pFlVE3arumO9PTlxnBhskkoH7jsf7n/XfPJrBkq8UxV0Dd9xPw+HSK7NK1kysufYF8mJxkdEA76tdsoZYbWZrvnWd66J6TtByO9u52BpRj1vlMua/0Q24TPuMoNjf2UjtO7isXuqhBDdI0DlDER6q2bE9f1mYGkg2QaB+51oUqQ3Uo57GRxiP2UjtHjsnjIhRC/us/eIqS8/+m2CjaUe907iP70Ot/8FWFpPtkHgTirFuHsW+fUUrP9QrMhiVUkI7sfWit8eS7Cx1OMOq/RH4/aTwBJdztln3OVgJtlaMttEA76tZc0ZKccrxGMPUqnH3RMVsRToProNn3G/AA23VjaikHN5H85RM9WVEvdxCTaWDO4nXt94kG4QUGlCeJTZvKuqn6cd9/OTYgdAljI7Cw03VzWiUvTZ1vmaPp/+AmNVURVxPzl0tWFMfsTWQ8Pjs4iplTju50ebseD9D1MNAirNCA/D2N9O66VaiUrhfvt2TBsTx92uwJ5HFYqRpWTi84dXXRcvePuwPFUR9ymRnxtGrczNzJT3skG4gVhpKnEvoBsZrU1W7iTNE5HHRdbVbbhdDmZ4HtF3VST6YOPHDeOef2zZzVTnG7+dYGMJ475JkHAv4ZAFuEN9EKq2BtxuT1ceJisDwhJeQxDBkAQbSxh3mFEmKmFEAXf6KZLfuHvlGC6HJ7+LlIfJyoAQ3C9CdtFElwonjDtUHaIy+QNKN9GNoLifkU+OPta0NXLrqG1PHe6aB+5yMMPxu77L3yxiUDKOqJwkx+6K2ngY7kdl5t9KyrM/w5nDPQphbMnGVbNSJn+ziD0gSCAmLaoyM4PivoOiXdPsFJgTtEEpmpnxwt0YJr7BDtONsDIjf7OIfSZmNZ4lHAAlReE8HPfaZ1H1BtxJxX+tVQ/3g/Vt86N0G6wMyd8sYiVirekHuL3MkyWDxB137ufgvqkuRXuP+Kc51cPd+Mgcs9V5MgRFqQMvf7OIrRSk3IXbS71ZShJ3t5Kqmg2qwiGeWzFu+sfqXbAyIn+ziI0UpDTE7em7uruVbtxfNpd5549PdCkkK/XyN4uYLkghJiKzAXfPUpY7nZn5XC/Rmw3yN4tYAy9UwE5NzBu+4w71VamljlDni8czvsvfLGI1vXCHS6eiRhaOe42FqPSU4w6llYigzSg/VQ2O/M0i5jkQWOdx5TSV1Ly7W9XDHS7exKMBubxjpnovrPQLf6qaqSxi3iXcF1jWBicUjfiMO9xudyQcINPZavVeWOlXojG+SiWMexVKuF8cr9WlK3OYSj/uw1S4zxcNjyccRNHS+pzi0Xf5i7vn2N0wg4SbqhvBcc9rhKpm4rhvaRQ7vpbkl93XTjbJmp8QDmfbWvY8YnU5K4PyF3eYmVHUtE4S95TNzCx2DpEsD7LTys9e80WyhXOPNtfq9Nuu2AcrQ/IX95bBxx0e7JJToRNt+52qhdMlvKo6EEoz7ldXj+s//hXyDWkczKQK98lwiJ/iDu8L+xr1YbICoPTifraHBcJA6vm598wMgXvZEamh2ni5cdS+jKYO90ZwiH/GHcSdqCLlNisoSi/uf3ZIuJ94g5Ys7oPiZ1wcjbPsqcO9FrRMZDQdKuxJp5FhZUxpxf2kuHrnn8XfkPRgRtZkqqwulj11uMtlUUTKuruFna4kxgqK0or7m4DKu/gbYJEFsTTbFIH7enT1xtOAu4ZPRNZIFHeom0fNq0MQ8xS6EVZAlFbcZU1rarWSuHZ+TTdK4I5XbVzu4P65fIBVWfHltZW4y++RcbjD5da2uRGRK4cVIKUV9ymAyjziHdE2prWeKmAhKdyNz+Ty07djnza5hYQlKnGHGGWylsxRK0is5W7FObACorTivhRQWUu8wzDH7zuVjSaHu0vJlREWagvnQGVFKhkf+yrJG6v4gmIFRWnFXdZY/IZ+U61g434rnMPTuEPU8WjDq5WCr2rjXuaqtlsZ912Ayn66gSRxX4feqi5KOe6L4Bw+wx2KhL1QuRNWEJQQ7sVTp02c7LqIjdZjektuV8YdggW15+lGk8SdUqpxPycmU6mQ+/vFnnsrd8IKghLB/UK/xbGr3ciLYvvAnTNnzix01b6tjHs3gFCRbTLguL8rGs4jBud3CIeWyp2wgqBEcJ/XLob25TaQ+23S4UoOlXGXc9bKx5aBxn09tExMvUDlpL7KnbCCoARwL+04xnwZ3duJ7jvVdd2Rih6VcZeTeJ3pZn3H/cOOWi06F/UWaJmomQBJu3kpavCVAO4f61aBijm6c9+5IDZwH7jH7VEZ96mAyjN0s0niPrUI0wOJ426H39SmXM6KoJlbiBSC5YNte4dLytNgBUEJ4L5Of858Wao7czGnti7opRe4JmbicJerkqn654bvE5GiJHUekdMuKkZkvagWLj0a+0Tk381L87JACeC+VLeWn63W5TRL2YqCLnYNrsNmVosF/1ApMEtEC45VNIvhfnLyBNDtWl25MWGF7ZA63CEGeRxu3y3sNek0JJf3fMilk7JCCeD+vG7V9VxjX+QdiWv9q1YA1u8r4R616pvnjVGlEsZwn0YOzZ101CnDXdYHqYs7bACHvYqTYGWHEsB9m25dWpfp7iCrkoLlcgNZvHe4l9b/C2WzGO6TtD+NwjRcPOxJGe5ngGZi/SDciVKrmVhZpARwP6xbcw+FeoVHpD1cpZyxtar3aQ+pm8Vxvxv1LU457rLqKRFyf0KMdhrzctPsV0ITkQ+YL/d1dN/URdu4UMok7pQSxF3eTVMrTMRT08Xq02BlgxJ5zPRUx9gFrqy9GRtyWvxt5wjXuDxzuPckcY/PFK/EXRYKo+qDXB5uDXUmJ1rKhBVAJYJ7ce9VseFCvwvmnOR6Y0+fwjPGsXHugXnmcD+/Q0awD9BGyo2d8bPfStxlSm3iVjWm9+Y/soQuZM/KIiUYIjZrxuNmLO+2ntuN4gk9+859rsKEegpxvxUtFPYqFpeYVF1VoSuAe2P1YbJyQWlPq5Qk7qRSjLu8uitq+7FyRWHHXVZ/UtTlZuWKGHeh7urDZOWCwo67DCIgsiaxckm+4H51wkCpPK2r3Bhx1HKYpN2+B9M7Kcf9MnyOFDHKrFyRL7jvRq/bluwFn54TkW554H58lqbdQ8YxyJgZMuKRlTvyBfddmrYALRTWWpttOaQQ9zetePUaRBoBVwERIiceK5fkF+74M8pBKcf9tEhDSRS8uQpjd/WIiJUTynXc/yJobks4QMpTKtEZK4cUVNy7r8L0TMK4Q3p2qvjGEOHwPt0IK1cUTNxlyYw4fR73JiXungUTPnRGO3fQbbByRn7hvhK9erdzcH+/fQGoiZYvNwpGXI1rTYm75oW7scWqezoi0QrhrGyUX7hTmh3nnFRtJpA37kbJ60+/wAGP4RDjzgqR0oJ7ietB6FBtpGvLjkZn3Fn+KC24DyBptm8IM4h7FcrQs8KjtODegqS5nWWP4V6M5qsekHLcO4g9J1jgmJWTShPu8QWQTL0gcVdPRLpVPdwhPy8vVmLlPu43CtzJGu+sECnXcZdVgNWNsEKhrMf9pboxlG88QJmrUqiYFRr5hfuXRzD1TRj3eQ7MRYS9NuPOkkoT7oXoYqSnJO4pmoi8Kt5HZHhk3Flu+TIR+XUjyp6/Na5RJe6Pwzu34w71vT4PrDDJF9yNK6659praFrmB1CZV4n4DtNwVd/iTsFOF81hhkj+4u5VU9Q4heSfaBHcYKewtqnoSrBxWluMur+49cYcdwj6/iufAymX5cqtaQdXCvQ/g/iThMcY2d7lcxXNg5bJ8mYisoGrh3gtwn0h4lBU21rQ6E3j1BsvIetw7Ae6KamfHPufCGyxLWY77CMCdEwmwvJUm3JeiD00Xphz3OYD7NvKg3n124evxK1xZYVTgZ2ZWNtPybv+Gsq6Dlqm8GUe7m9b2u6p+FqzcVVpwH0nijoyw1bi3sd6W9wJhXgYtE1nCrrS1zQ2OJnQmrNxUWnCPfiUHMHdp98qNY8g9oxJ3UXZbI0qwr/W6uq8Wdo+026xQKJj53UFAc3z9SEvy6k7EzIgykVp79VGwQqFswZ2I8JIl3ImrO4yruBYNK+i4y7p4RPzuTrDHJ9OzNFXY+6uPghUKBRv3qBfuV0QocWuihT0iSeQq9VGwQqFg417uNZgxXrDN+W9STcywHYbwg1VW9uBOJhJY0yRmbbGFbv/F7rXz2yzg50wsI+i4Rz2v7oZxdf+WT8uUeyi7ojSzwqNg4361CrizWFVWsHGXg5la1T02FisAuJfU0F6hrYC7Xt1jY7F8x/1iM5Pl+vuoBqHWzO0pOUBWyOUz7k7du7z4AmO2IF/1iJQcICvk8hf3BwTNDYkGYTAzIDVHyAq3/MXdM8cX2Hum5ABZIZe/uHvFCMixOw9mWCmQv7h7Vj2tJ+yPp+YIWeGWv7jXqDLuU1NzhKxwy1/cm3g8NC2HmZlHU3SIrFArIdyLp06bOLmE2kwCd90rRgAcnk7kOFksXIngfqHfYsNYNPIivmlguO80ayM1/4JqESonUYOZaY655pcJHCeLRSgR3Oe1KzWMy22W4ZsGgrtYGL2JaLGBF+7n29vmhQkcJotFKQHcSzuOMV9G945im6Yq4w5PiaixSlMvB+PcQ420/A54UjIWK0ElgPvHupU0eo6+H9s0VRn3lXDxji/KYamLsNdU7Pf0paofI4ulUgK4r9OfM1+W6m9gm6Yq4w6lM7TOeJN/FvYbEzpoFis5JYD7Un2d+bJafx7Z/Hp7TG//qCLusqhpfbzJ8cLeLPEjZ7ESVgK4P69vMF/W2Ff1SpsbLWr/UBH3OoB7Q7zJMcLOJdxZmVACuG/TV5gvy/SXkM1Ss5BY8Xcq4g6z5loPvMkZwt498SNnsRJWArgf1heZL4XOvWmlTVOVx+5wJ6oNxpvcK4JmFiV00CxWckpoIvIB8+W+jmXYpqnKuLfxurqLx0i9OVcAKxNK5DHTUx2jhlHWvjD26+mKm44q4y4TX5OBBC+0zdeaTOOpRlZGlAjuxb1XGcbyfhfMScj17k2hyri/DLi/R7d6kaxVwGKlWAmGiM2a8bhJ57ae292bQpVxfx1w31H9I2Wxqq20BgBDLQHttVTuhsVKUmnFfTbgvjyVu2GxklRacZ8HuC8j3sFiZVJpxf09wJ3Km8RiZVJpxb1crEWNX+TEYvmgtOJ+VFzc80+ncjcsVpJKK+5vwGBGXRibxcqM0or724D7B6ncDYuVpNKK+zcib0btEuIdLFYmld48MyLl6aRU7oXFSlbpxf3SEIv2ERzwyAqE0pxFrPyNaaOe2Faeyp2wWEkr7UnzWKzgiHFnhUiMOytEYtxZIRLjzgqRGHdWiMS4s0Ikxp0VIjHurBCJcWeFSIw7K0RKMe4/1FisQKlCRd6U4m68EYnT31z//fg/uvVPv1bbI//5M7X9B9d/V+3wL//usYff/b3a/sPrr1E7/PxXHnu47qdq+4+u/47a4Rf/5rGH/7pWbf9f139L7fDLf/XYw/U/Udt/cr1HA//6S7X9W9f/L7XDtf/lsYd/+0X83+5KH+6I9mufqx2md/VoQfdID/yhdkztMLm3xx6aeuQFeVf7Wu3w8K0ee2i0Sm1/Szuvdhg/xGMPddaq7a9pHlHY99yltht5G9X2jWSBLUd33aO2X/HKvrW2jscehoz3cGDcDcbdEuOeCjHuphh3g3G3xLgbjLslxt0U424w7payAPdzRR5ZCA5s92jhjc/U9uIij2oI+9/x2MMWj0/kmSIPVPa+67GH146o7aeLrqoddisy5FvadFRtP1kUVTt8tMtjD0XH1fbjRR4N7PpIbY8WnVQ7HKWqrwu9t9vDIe24s1jBEePOCpEYd1aIFBLcy97w9mFVTdncl2nHfftwpfnA2M4DFqpuBA+N69R77mVlG1eHvqwyj9ZjekvZwmfzFpF3QdN0Sz3pdx+cMn/6BNX99MHphROXUneKooOKp06bOBm7rZc9SPSl+DPVl8JOdqWrXbwvwYHoS1cDeFc6DnRfihaovgQ73pXy1MlutJRm3PfOKGirsr/Xv/DJbvr9tMO+e956/1F9pnInC3UV7gfunDlzZmGpwqN44m37SWNJz8efXb58+cAZpMehbicM442uZ0iHXZ3OGGX3jEN5hw660G+xYSwaeZF0oPoS/kz0JdiprqzQLtaX0gHvS1cDeFcKB7IvoQWiL8FOdKU8daobHaUZ98vGMBXu5feeN4yzffU9pMeMSzGvQcqJ84/HKHGfdFh9iMaxASMU3x4vWTOA0a507cBpo2I/ytuR03DRQY/Efn6ob8aM0EHz2sUoutwm7gmA7EGiL8Wfqb6Et1Fd6W4X7UvpgPeltBNdKRzIvoQWiL4UdqIrXadOdaOjtA9m7lHhfmCL+XOtvo5yiFqf84fGKtoouf8TFe6nuq5Tz3lfHdL5hMJsz+l/1J4ecI3tFrOVtNpM2T/VlxhmveWhuNnuoNKOY8yX0b3jvwNkDxJ9af+Z7kvbruhKaJfqS8eB7EvHTnel7aDoS6cFsi9tO9GV8tQV3WjJX9xLrKPaiV/4QJcH0Vd/w5h+4HMV7gtiY8WBqvev0Bco925pzsO0bZX+WJlRNIT8PGzVV5gv/VvjqTLtDvpYn2/tR48fC1QRd7ovXW/DuxIcqL50HMi+dOx0V7qPHO1Lx4HsS9tOdKU8dUU3WvIXd1svtqTHvTFdnLBBYX3zOUOJ+6mtC3rpBYrJhP768+P7j/F4fN33VdpWPll/8PU59NX/bX2W+TJcx698dget058zX5bq8UdaRdxtYX0p7URXCgeyL8XVnepLx053pfsQ0b50HMi+tO3KrjRPXdGNloKA+/2zVdYV3fSC+aT1zINRNe4xla0o6FJMGY/qbd+6+MVYfaWqhUMtzyms0Uf0NoqZnxMF/cyL0TD9G9Rsd9BSexCyWn+ecKj4G+GA9yXYqa50HOi+lHvA+9K2K7rSdYh4XwoHqi9tu7IrzVNXdKOlAOC+d6g6IqX45bY6Wezm0eOGJ+6KT7thfKBPif281Km9KrRniTK46ezwsXrBC7R9mr6orHxzq7aqwczzunXdXWNfnRCHir8RDnhfSjvRlY4D3ZfuHdNfQIqudDWA96VwoPrSsSu60jp1RTda8h/3S0O+8GqjiAyKXG/ew3vjXlJAlu3+QC80Xx5XzA4ZxuA1CmPx4K3G6oICeuqmbFH/XuM2t3oQt9odtM0elS7TXyIcKv5GOOB96X4b2pW2g6Iv3S1gfSlwJ7vS1QDel44D2ZeOne5K+9QV3WjJf9xn05wIRTtRE++jdEcqHmPqQQ7/T+rW2t1n9G30u78iht22Jg+I/dik36s+gpXUDuwOOmxjWJj8raopvC/db0O70nZQ9GWFHSN9adsVXSkbIPrScSD70n0EaFfap67oRku+475BcRMI6k98WI3978ZUpM99V8Vj7J/c5kPKVN7R+nJ9Rld8x6y6Q9X4rdZMw8D+yiM42ZV6liYmIh8wX+7rWEY4VPwNdyD6ssLbsK60HRR9WeEDg/SlbVd0pWyA6EvHgexL1xGgXemcuqIbLfmN+9bV5s9j6mju4m7K5RXHPAczO0fQ9XKWtDJvvSYPIDrI1MhnVI2PHWH+nKB4NBz7qh12B36jCh30VMeoOaNcSDoYnrhTfel+G9qVLge8L90tYH3p2OmulA0QfSnm3am+dA3YsK6EU6e70VLacR/aSmV9v/+SmOYMphZoXB773Hnj8mT1eEeF+54+hWeMY+MUl+6SPrEbrK+7KAodnyk4qNr9++ZQ8+suqtKxu299nLwTdjqouPcqw1je7wLpYJB96fyZ7EvbruhKV7t4X9oOdF86DdBdCXug+lKcA9WX0ADalfLU6W60lO6YmQUF+mx6ku7T9vZokYyJKZvcu8/8FeQsoi0V7sUTevad+5xy6ufctIcWTKODZgzjRY/Ffx+OmbZwomItUFHhfPLjJjuoeOqsGY/HX7ekA9GX4s9UXwo72ZUV2sX6UjhQfSkbILrStQe8L6UD3pdgx7vSfepENzoKSQAwi2WKcWeFSIw7K0Ri3FkhEuPOCpEYd1aIxLizQiTGnRUiMe6sEIlxZ4VIjDsrRGLcWSES484KkRh3VojEuLNCJMadFSIx7qwQiXFnhUiMe1B0cECn3g0HZ9fwVAAAAilJREFU2gun93Zq83/1p74yf40+3uSX/9HjS18PLXfEuAdE+669zTAOf7ep+ftbP9hinP3vyLeuaWqcr/tUmbHxu9d6lcVjVUmMe0B0b8RMbf7byKnYz9/9KvZjTSQ/9nOEtQK/W6SRn8eWO2LcA6IDo8wSLb+P7DaMQ5HrYr9e/dtrzhtnr/l5Xky/iHzfo7Irq0pi3IOj4o8O/GfkA8PYGvmRmbnoum9fMLZFFInlWQmLcQ+KTt65xjCuN3E/9f3Ii4ZR/s/1DeOZyCC/jyunxLgHREU/NTOzW7gbkyK/evfc+N99YRibI7917Md8PLbcEeMeEP1LZG/s53UW7sazzQff8bRZ4O6rSGStZZ670cdjyx0x7sHQpUgkdnU/+sPIm7GNnVPh760j/7DFMMpmtPDv0HJJjHtApEX+fmybXg0jeu/LJ7/1/wbeOnzs5FdiN6yf/Esk8vuG//xHjzSZrKqJcQ+I9jX4cf1Fxuof/9erhnH7tyOWGl4xjBN9/vOH/zNJXTacVVUx7gHUY+0eGtbt5ut/Enna7yPJNTHuwdMdf3BStE+f5u+B5J4Y98DpvYiobD12n68HkoNi3AOng5Ffv2sYl754ocuzfh9KzolxD56e/afIj/42ck3/I34fSO6JcQ+grm6ev+S9i34fRS6KcWeFSIw7K0Ri3FkhEuPOCpEYd1aIxLizQiTGnRUiMe6sEIlxZ4VIjDsrRGLcWSES484Kkf4/aoxIRbPDZYAAAAAASUVORK5CYII=" width="60%" style="display: block; margin: auto;" /></p>
<p>If we were pharmacological, we could evaluate drug concentration using a first-order compartment model, such as</p>
<div class="sourceCode" id="cb15"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb15-1" title="1">post3 &lt;-<span class="st"> </span><span class="kw">stan_nlmer</span>(conc <span class="op">~</span><span class="st"> </span><span class="kw">SSfol</span>(Dose, Time, lKe, lKa, lCl) <span class="op">~</span><span class="st"> </span></a>
<a class="sourceLine" id="cb15-2" title="2"><span class="st">                    </span>(<span class="dv">0</span> <span class="op">+</span><span class="st"> </span>lKe <span class="op">+</span><span class="st"> </span>lKa <span class="op">+</span><span class="st"> </span>lCl <span class="op">|</span><span class="st"> </span>Subject), <span class="dt">data =</span> Theoph,</a>
<a class="sourceLine" id="cb15-3" title="3">                    <span class="dt">cores =</span> <span class="dv">2</span>, <span class="dt">seed =</span> <span class="dv">12345</span>, </a>
<a class="sourceLine" id="cb15-4" title="4">                    <span class="dt">QR =</span> <span class="ot">TRUE</span>, <span class="dt">init_r =</span> <span class="fl">0.25</span>, <span class="dt">adapt_delta =</span> <span class="fl">0.999</span>)</a>
<a class="sourceLine" id="cb15-5" title="5"><span class="kw">pairs</span>(post3, <span class="dt">regex_pars =</span> <span class="st">&quot;^l&quot;</span>)</a>
<a class="sourceLine" id="cb15-6" title="6"><span class="kw">pairs</span>(post3, <span class="dt">regex_pars =</span> <span class="st">&quot;igma&quot;</span>)</a></code></pre></div>
<p>However, in this case the posterior distribution is bimodal Thus, you should always be running many chains when using Stan, especially <code>stan_nlmer</code>.</p>
</div>
<div id="conclusion" class="section level1">
<h1>Conclusion</h1>
<p>There are model fitting functions in the <strong>rstanarm</strong> package that can do essentially all of what can be done in the <strong>lme4</strong> and <strong>gamm4</strong> packages — in the sense that they can fit models with multilevel structure and / or nonlinear relationships — and propagate the uncertainty in the parameter estimates to the predictions and other functions of interest. The documentation of <strong>lme4</strong> and <strong>gamm4</strong> has various warnings that acknowledge that the estimated standard errors, confidence intervals, etc. are not entirely correct, even from a frequentist perspective.</p>
<p>A frequentist point estimate would also completely miss the second mode in the last example with <code>stan_nlmer</code>. Thus, there is considerable reason to prefer the <strong>rstanarm</strong> variants of these functions for regression modeling. The only disadvantage is the execution time required to produce an answer that properly captures the uncertainty in the estimates of complicated models such as these.</p>
</div>



<!-- code folding -->


<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>