.. _builtin_tut_example:


Generate data and Fit using a Two Compartment Model
#####################################################

In this example we will demonstrate a :ref:`tut_script` using the same two compartment model with absorption and bolus dosing, as used in the :ref:`builtin_fit_example` and :ref:`builtin_gen_example`, see :numref:`fig_two_comp_depot_diagram_tut`:-

.. _fig_two_comp_depot_diagram_tut:

.. figure:: /_autogen/quick_start/builtin_tut_example/compartment_diagram.*
    :width: 80%
    :align: center
    
    Two compartment model with depot dosing for :ref:`tut_script`.
    
This |popy| model is also used when :ref:`creating_a_tut_script`.

.. note::
    
   See the :ref:`sum_link_builtin_tut_example_tut` obtained by the |popy| developers for this example, including input script and output data file.
 
A :ref:`tut_script` can be used as a theoretical tool to quickly investigate identifiability of |pkpd| models, because the true |fx| parameters and the structure of the data are known.
 
The |popy| Manual makes extensive use of :ref:`tut_scripts<tut_script>` to create examples to illustrate different :ref:`indiv_pkpd_index`.
 
.. _running_a_tut_script:
 
Running the Tutorial Script
============================

This tutorial example requires a single input file:-

.. code-block::

    c:\PoPy\examples\builtin_tut_example.pyml
                            
:ref:`open_a_popy_command_prompt` to setup the |popy| environment in this folder:-

.. code-block::

    c:\PoPy\examples\

With the |popy| environment enabled, call :ref:`popy_run` on the :ref:`tut_script` from the command line:-

.. code-block:: console

    popy run builtin_tut_example.pyml

When the tut script has terminated, you can view the output of the fit using :ref:`popy_view`, by typing the following command:-

.. code-block:: console

    popy view builtin_tut_example.pyml.html

Note the extra '.html' extension in the above command. This command opens a local .html file in your web browser to summarise the result of the generating process.

You can compare your local html output with the pre-computed documentation output, see :ref:`sum_link_builtin_tut_example_tut`. You should expect some minor numerical differences when comparing results with the documentation.

.. _syntax_tut_script_builtin:

Syntax of Tut Script
======================

The major structural difference between a :ref:`fit_script` or :ref:`gen_script` and a :ref:`tut_script` is that the generating |effects| are encoded in |gen_effects| and the fitting |effects| are encoded in |fit_effects|. For example the |gen_effects| section for this tutorial example is as follows:-

.. literalinclude:: 
    /_autogen/quick_start/builtin_tut_example/tut_sections/GEN_EFFECTS.pyml
    :language: pyml      
    
And the |fit_effects| section is as follows:-

.. literalinclude:: 
    /_autogen/quick_start/builtin_tut_example/tut_sections/FIT_EFFECTS.pyml
    :language: pyml     
       
The |gen_effects| get passed to the :ref:`gen_script` and the |fit_effects| get passed to the :ref:`fit_script`. From the examples above you can see that the GEN_EFFECTS->POP section has:-

.. code-block:: pyml

    f[X] = true_value
    
Whereas the FIT_EFFECTS->POP section has:-

.. code-block:: pyml

    f[X] ~ P starting_value
    
Reflecting the fact that the |fx| are known constants for a :ref:`gen_script`, but are unknown values to be estimated in a :ref:`fit_script`.

.. _summary_of_builtin_tut_results:

Summary of Tut Results
=============================

The :ref:`tut_script` generates an output folder containing four new scripts:-

.. code-block::

    builtin_tut_example.pyml_output/
        builtin_tut_example_gen.pyml
        builtin_tut_example_fit.pyml
        builtin_tut_example_comp.pyml
        builtin_tut_example_tutsum.pyml
                                              
See :ref:`files_generated_by_tut_script` for more info. These four scripts are run in order. The :ref:`gen_script` is described in :ref:`builtin_gen_example` and the :ref:`fit_script` is described in :ref:`builtin_fit_example`. Therefore here we will focus on the :ref:`comp_script` outputs, which are fitted |fx| and generated |fx| plots and their objective values. The simplest comparison output is a visual comparison of the true and fitted |fx| |pk| curves and the synthetic measurement data, see :numref:`table_true_vs_fit_fx_plots_builtin`.

.. _table_true_vs_fit_fx_plots_builtin:

.. list-table:: Fitted model |pk| curves vs true model |pk| curves for first three individuals 
    :class: image_table
    
    * - .. thumbnail:: /_autogen/quick_start/builtin_tut_example/images/fit_sim_grph_outputs/indOBS_vs_TIME/000001.*
      - .. thumbnail:: /_autogen/quick_start/builtin_tut_example/images/fit_sim_grph_outputs/indOBS_vs_TIME/000002.*
      - .. thumbnail:: /_autogen/quick_start/builtin_tut_example/images/fit_sim_grph_outputs/indOBS_vs_TIME/000003.* 

The solid blue lines in :numref:`table_true_vs_fit_fx_plots_builtin` show the predicted |pk| curves for the fitted model |fx| values. The dotted blue lines show the |pk| curves for the true |fx| values that were used to generated the data set (in the :ref:`gen_script`). The blue dots are the target |cdvcen| values from the data file.

The target |cdvcen| values have measurement noise added, so blue dot data points do |not| lie on the true |fx| curves. The graphs show that the |pk| curves for the fitted |fx| are very similar to the true |fx| curves. Most divergence occurs away from the data points and most agreement close to the data points, for example for Individual 2 the fitted and true curves differ over time period [0,30], but are very similar in the period [30,50]. In the period [0,30], the model tends to impute the curve for the median individual, in the absence of actual data. 

.. _fit_vs_true_fx_builtin:

Fit vs True f[X] values
------------------------------

If the :ref:`comp_script` has been run, the :ref:`tutsum_script` output also contains a convenient table to compare the initial, fitted and true |fx| values, see :numref:`table_true_vs_fit_fx_main_values_builtin`.

.. _table_true_vs_fit_fx_main_values_builtin:

.. csv-table:: Comparison of main initial, fitted and true |fx| values
    :file: ../../../_autogen/quick_start/builtin_tut_example/fx_comp_main.csv
    :header-rows: 1

:numref:`table_true_vs_fit_fx_main_values_builtin` shows that the :pyml:`f[KA], f[CL], f[V1], f[Q]` parameters are recovered reasonably well, in the sense that the fitted values are much closer to the true values, compared to the initial starting values. However the :pyml:`f[V2]` fitted value is quite different from the true value.

You can also compare the proportional noise estimate:-

.. _table_true_vs_fit_fx_noise_values_builtin:

.. csv-table:: Comparison of fitted and true proportional noise |fx| values
    :file: ../../../_autogen/quick_start/builtin_tut_example/fx_comp_noise.csv
    :header-rows: 1
    
Here the :pyml:`f[PNOISE]` starts at 0.1, is estimated at 0.141, which can be compared with the true value 0.15. Generally the proportional noise is identifiable in a |pkpd| model. This is because every row of the data set and the current :pyml:`f[PNOISE]` parameter has an influence on the overall likelihood. 
    
The comparison of the covariance matrix estimates is quite long, as it displays a row for 5*5 |fx| comparisons. However the diagonal elements (only) are shown here:-

.. _table_true_vs_fit_fx_diag_isv_values_builtin:

.. csv-table:: Comparison of fitted and true isv variance diagonal |fx| values
    :file: fx_comp_variance_diag.csv
    :header-rows: 1
    
This shows that :pyml:`f[CL_isv]` and :pyml:`f[Q_isv]`, the inter-subject variances of the clearances are estimated reasonably well. The :pyml:`f[KA_isv]` estimate is far too high. The :pyml:`f[V1_isv]` and :pyml:`f[V2_isv]` estimates are worse than the starting value estimates.

The fact that :pyml:`f[V1_isv]` and :pyml:`f[V2_isv]` are badly estimated compared to the true values is not that surprising, since |vols| are harder to estimate in |pk| models, compared to |clears| which are rates. The difficulty in estimating :pyml:`f[KA_isv]` may be due to the relative lack of data before the |cmax| peak of the |pk| curve, there being only 5 data points sampled per individual. 

There can be multiple reasons for the fitted values not agreeing with the true parameters, for example:-

* Too few observations in the data set.
* Too few individuals in the data set.
* Too much noise added to the measurement data.
* False minima on the likelihood surface.
* A fundamental difficulty in identifying some |pk| parameters, for example if the model is over-parametrised relative to the data set or even unidentifiable.

Given the relatively small amount of data generated (250 data points), the results in :numref:`table_true_vs_fit_fx_main_values_builtin` are adequate. As shown in the next section the objective function for the fitted |fx| solution is actually lower than for the true |fx| solution here, see below for discussion.

.. _fit_vs_true_obj_builtin:

Fit vs True Objective value
------------------------------

It is difficult to know by just comparing the fitted |fx| and true |fx| in section :ref:`fit_vs_true_fx_builtin`, if the fitting method has done well or badly. We can run a form of sanity check by computing the objective function of the true |fx| and comparing this with the objective function of the fitted |fx|. This is done by optimising the |rx| for both solutions and using the synthetic data file to compute the likelihood.

The rational is that the fitted |fx| objective value should always be **lower** than the true |fx| objective value, because the fitted model estimates can take advantage of correlations in the random measurement noise to get a better fit to the synthetic data. If the fitted objective function is **higher** then the |popy| fitting method has ended up in a sub optimal local minima, because the known true values are a better minima.
        
In this example, the true model objective function is:-

.. code-block:: pyml

    -881.0027

Compared with the fitted model objective function:-
    
.. code-block:: pyml

    -896.8752
      
This indicates that the fitted |fx| pass the sanity check and perhaps justifies the lack of agreement with parameters such as :pyml:`f[V2]`. However it does |not| say if the global optimum solution was found, |ie| whether or not |fx| is a true global minima of the likelihood surface.

The difference between the two objective values above is partially dependent on the amount of noise applied to the measurements |ie| :pyml:`f[PNOISE]`, the number of individuals simulated, the number of observations per individual and the number of parameters in the model. More random noise in the `synthetic data` creates more likelihood minima that are away from the true |fx| solution.

If the model is practically identifiable then increasing the number of data points and individuals should lead to a convergence between the true and fitted |fx| and the objective functions above.
      
.. _rerun_tut_builtin:

Re-run the tutorial
=====================

You can familiarise yourself with |popy|'s various features by tweaking the tutorial example and re-running. 
A simple way of avoiding overwriting previous results is to do:-

.. code-block:: console
    
    copy builtin_tut_example.pyml builtin_tut_example_v2.pyml

and modifying the script with your chosen text editor.    

When you are happy with the edited file, do:-

.. code-block:: console

    popy run builtin_tut_example_v2.pyml

For example, you can adjust the amount of data generated, in this section:-

.. code-block:: pyml

    GEN_EFFECTS:
        ID: |
            c[ID] = sequential(50)
            t[DOSE] = 2.0
            t[OBS] ~ unif(1.0, 50.0; 5)
            # t[OBS] = range(1.0, 50.0; 5)
                
Note the 'range' function can sample time points evenly (instead of randomly). This usually makes the model fitting easier as the data points span the time range better. You can experiment with the number of doses or make a random sample of dose times for each individual.

You can also edit the underlying model or compartment structure, see the |model_params| or |derivatives| sections.

Another possibility is changing the fitted model only, i.e take a copy of this file:-

.. code-block::

    builtin_tut_example.pyml_output/
        builtin_tut_example_fit.pyml

Name it 'builtin_tut_example_fit_v2.pyml', change the model structure or compartment. Then do:-

.. code-block:: console

    popy run builtin_tut_example_fit_v2.pyml
    
Using a different model from the underlying generative model should result in a worse fit (using say an |akaike_inf| to compare fits).


