.. _uncertainty:

Uncertainty and Standard Errors
######################################

Parameter estimation (also referred to as "model fitting") aims to find a solution in parameter space that has a likelihood at least as good as all other nearby solutions. We know that a locally optimal solution has been found when changing any parameter in any direction decreases the likelihood (or, equivalently, increases the cost).

.. _likelihood_hessian:

Likelihood Hessian 
=====================

We can also quantify the confidence in our estimates from the shape of the likelihood surface at the local minimum: changing some parameters will lower the likelihood by a lot, suggesting that we have high confidence (or low uncertainty) in their estimated value; changing others, however, will lower the likelihood by only a little or possibly not at all, suggesting that we have low confidence (high uncertainty) in our estimate.

.. note::
    See :ref:`sum_link_d1cmp_ka_stderr_tut` for the :ref:`tut_script` used here.

First we return to the example used in the last chapter -- a one compartment model with absorption -- only here we reduce it to a one dimensional problem by fixing *CL* and *V* to the optimal values and allowing only *KA* to vary. We can then plot |objv| for a range of values of *KA* as a line (:numref:`fig_d1cmp_ka_objv`).

.. figure:: /_autogen/surfaces/d1cmp_ka_stderr_copy.pyml_output/d1cmp_ka_stderr_fit.pyml_output/objv_line.*
    :name: fig_d1cmp_ka_objv
    :align: center
    :width: 50%

    |objv| |vs| *KA* with other parameters fixed at their true values

Although |objv| as a function of *KA* is not symmetric, it is smooth and has a clearly defined minimum. We can therefore compute properties of the curve (|eg| its gradient) at any point and approximate it with a polynomial using a `Taylor series <https://en.wikipedia.org/wiki/Taylor_series>`_ expansion.

More specifically, using a second order Taylor series at the minimum requires only the minimum value, |objv|, and its second derivative because the gradient is zero by definition. This approximates the |objv| function with a quadratic (:numref:`fig_d1cmp_ka_objv_overlay`).

.. figure:: /_autogen/surfaces/d1cmp_ka_stderr_copy.pyml_output/d1cmp_ka_stderr_fit.pyml_output/objv_line_overlay.*
    :name: fig_d1cmp_ka_objv_overlay
    :align: center
    :width: 50%

    |objv| |vs| *KA* with other parameters fixed at their true values, plus the quadratic approximation at the minimum

In effect, this approximates the likelihood surface with a Gaussian whose peak coincides with that of the true likelihood (:numref:`fig_d1cmp_ka_lhood_overlay`).

.. figure:: /_autogen/surfaces/d1cmp_ka_stderr_copy.pyml_output/d1cmp_ka_stderr_fit.pyml_output/lhood_line_overlay.*
    :name: fig_d1cmp_ka_lhood_overlay
    :align: center
    :width: 50%

    Likelihood |vs| *KA* with other parameters fixed at their true values, plus the gaussian approximation at the maximum

.. note::
    See :ref:`sum_link_d1cmp_ka_v_stderr_tut`, :ref:`sum_link_d1cmp_ka_cl_stderr_tut` and :ref:`sum_link_d1cmp_v_cl_stderr_tut`,  for the :ref:`tutorial scripts <tut_script>` used here.

In more than one dimension, the second derivative is a matrix and is known as the :term:`hessian`. In two dimensions, for example, the quadratic approximation to |objv| resembles a basin whose shape reflects the true function at its minimum (:numref:`figtable_2d_approx`).

.. comment
    The following table is pretty big so \clearpage makes sure it starts on a fresh page in the PDF

.. only:: latex

    .. raw:: latex

        \clearpage

.. list-table:: Comparison of |objv| surface (left) with quadratic approximation (right)
    :name: figtable_2d_approx

    * - .. thumbnail:: /_autogen/surfaces/d1cmp_ka_v_stderr_copy.pyml_output/d1cmp_ka_v_stderr_fit.pyml_output/objv_contour.*
      - .. thumbnail:: /_autogen/surfaces/d1cmp_ka_v_stderr_copy.pyml_output/d1cmp_ka_v_stderr_fit.pyml_output/se_contour.*
    * - .. thumbnail:: /_autogen/surfaces/d1cmp_ka_cl_stderr_copy.pyml_output/d1cmp_ka_cl_stderr_fit.pyml_output/objv_contour.*
      - .. thumbnail:: /_autogen/surfaces/d1cmp_ka_cl_stderr_copy.pyml_output/d1cmp_ka_cl_stderr_fit.pyml_output/se_contour.*
    * - .. thumbnail:: /_autogen/surfaces/d1cmp_v_cl_stderr_copy.pyml_output/d1cmp_v_cl_stderr_fit.pyml_output/objv_contour.*
      - .. thumbnail:: /_autogen/surfaces/d1cmp_v_cl_stderr_copy.pyml_output/d1cmp_v_cl_stderr_fit.pyml_output/se_contour.*

.. comment
    David does not really understand this section, so it is commented out.

    .. _standard_errors:

    Standard Errors
    =================
    .. note::
        See :ref:`sum_link_d1cmp_ka_v_cl_stderr_tut` for the :ref:`tut_script` used here.

    The hessian matrix (or its inverse) that is used to approximate the |objv| surface,

    .. literalinclude:: /_autogen/indiv_examples/uncertainty/d1cmp_ka_v_cl_stderr/fe_hess_inv.txt

    contains useful information about the shape of the surface, and therefore about our certainty in the parameter estimates.

    Specifically, the hessian in *n* dimensions has *n* axes whose direction indicates the level of correlation between population parameters (|ie| :term:`fixed effects`) and length indicates the amount of variation in those directions.

    Using an eigenvalue decomposition of the hessian, we can compute the directions of variation (from the unit-length eigenvectors) and also the magnitudes (the square root of the eigenvalues):

    .. literalinclude:: /_autogen/indiv_examples/uncertainty/d1cmp_ka_v_cl_stderr/fe_hess_inv_eigs.txt

    The greater the magnitude, the greater the uncertainty in the parameter estimate. In this example, the smallest uncertainty (highest confidence) is associated with the third parameter, *CL*, while the greatest uncertainty (lowest confidence) is associated with the second parameter, *V*.


.. _confidence_intervals:

Confidence Intervals
======================
.. note::
    See :ref:`sum_link_d1cmp_ka_cl_stderr_tut` for the :ref:`tut_script` used here.

With a gaussian approximation to the likelihood, we can sample new values for the population parameters from the distribution. By sampling many possible solutions from the distribution, we can compute a range for each parameter in which a given percentage of the sampled values lie. These ranges are called :term:`confidence intervals` and typically encompass 90% or 95% of the sampled values (:numref:`fig_d1cmp_ka_cl_confint`).

.. figure:: /_autogen/surfaces/d1cmp_ka_cl_stderr_copy.pyml_output/d1cmp_ka_cl_stderr_fit.pyml_output/objv_contour_samples.*
    :name: fig_d1cmp_ka_cl_confint
    :align: center
    :width: 50%

    |objv| surface with population parameters sampled from the approximating gaussian and 90% confidence intervals in red for *KA* and *CL*

.. literalinclude:: /_autogen/surfaces/d1cmp_ka_cl_stderr_copy.pyml_output/d1cmp_ka_cl_stderr_fit.pyml_output/conf_ints.txt

The purpose of the confidence interval is to compare one estimate with another, for example by confirming that our estimated confidence interval contains a value published in research literature by a third party.

When computing confidence intervals, we can get different estimates depending on how we parameterized the model. When population parameters are required to be positive, for example, we might optimise over values for their logarithm (which spans the whole of the real line). This can change the |objv| surface to be closer to a quadratic, and effectively samples from a lognormal distribution over the population parameters (:numref:`fig_d1cmp_ka_cl_log_confint`).

.. figure:: /_autogen/surfaces/d1cmp_ka_cl_log_stderr_copy.pyml_output/d1cmp_ka_cl_log_stderr_fit.pyml_output/objv_contour_samples.*
    :name: fig_d1cmp_ka_cl_log_confint
    :align: center
    :width: 50%

    |objv| surface with population parameters sampled from the approximating gaussian and 90% confidence intervals in red for log(*KA*) and log(*CL*)

One outcome of this is that (in the limit) the confidence interval is less likely to be centred on the estimated value:

.. literalinclude:: /_autogen/surfaces/d1cmp_ka_cl_log_stderr_copy.pyml_output/d1cmp_ka_cl_log_stderr_fit.pyml_output/conf_ints_log.txt

.. _standard_errors:

Standard Errors
=================

The :ref:`likelihood_hessian` can be used to compute standard errors for parameter estimates, which are similar to the :ref:`confidence_intervals` described above and computed in a similar manner.


