diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 00000000..abcbaa10 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,353 @@ +# Contributing to Chi + +If you'd like to contribute to Chi (thanks!), please have a look at the [guidelines below](#workflow). + +If you're already familiar with our workflow, maybe have a quick look at the [pre-commit checks](#pre-commit-checks) directly below. + +## Pre-commit checks + +Before you commit any code, please perform the following checks: + +- [No style issues](#coding-style-guidelines): `$ flake8` +- [All tests pass](#testing): `$ python run-tests.py --unit` +- [The documentation builds](#building-the-documentation): `$ python run-tests.py --doctest` + +## Workflow + +We use [GIT](https://en.wikipedia.org/wiki/Git) and [GitHub](https://en.wikipedia.org/wiki/GitHub) to coordinate our work. When making any kind of update, we try to follow the procedure below. + +### A. Before you begin + +1. Create an [issue](https://guides.github.com/features/issues/) where new proposals can be discusssed before any coding is done. +2. Create a [branch](https://help.github.com/articles/creating-and-deleting-branches-within-your-repository/) of this repo (ideally on your own [fork](https://help.github.com/articles/fork-a-repo/)), where all changes will be made +3. Download the source code onto your local system, by [cloning](https://help.github.com/articles/cloning-a-repository/) the repository (or your fork of the repository). +4. [Install](#installation) Chi with the developer options. +5. [Test](#testing) if your installation worked, using the test script: `$ python run-tests.py --unit`. + +You now have everything you need to start making changes! + +### B. Writing your code + +5. Chi is developed in [Python](https://en.wikipedia.org/wiki/Python_(programming_language)), and makes heavy use of [NumPy](https://en.wikipedia.org/wiki/NumPy) (see also [NumPy for MatLab users](https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html) and [Python for R users](http://blog.hackerearth.com/how-can-r-users-learn-python-for-data-science)). +6. Make sure to follow our [coding style guidelines](#coding-style-guidelines). +7. Commit your changes to your branch with useful, descriptive commit messages: Remember these are publically visible and should still make sense a few months ahead in time. While developing, you can keep using the github issue you're working on as a place for discussion. [Refer to your commits](https://stackoverflow.com/questions/8910271/how-can-i-reference-a-commit-in-an-issue-comment-on-github) when discussing specific lines of code. +8. If you want to add a dependency on another library, or re-use code you found somewhere else, have a look at [these guidelines](#dependencies-and-reusing-code). + +### C. Merging your changes with Chi + +9. [Test your code!](#testing) +10. Chi has online documentation at http://chi.readthedocs.io/. To make sure any new methods or classes you added show up there, please read the [documentation](#documentation) section. +11. When you feel your code is finished, or at least warrants serious discussion, run the [pre-commit checks](#pre-commit-checks) and then create a [pull request](https://help.github.com/articles/about-pull-requests/) (PR) on [Chi's GitHub page](https://github.com/pints-team/chi). +12. Once a PR has been created, it will be reviewed by any member of the community. Changes might be suggested which you can make by simply adding new commits to the branch. When everything's finished, someone with the right GitHub permissions will merge your changes into Chi's main repository. + +## Installation + +To install Chi with all developer options, use: + +``` +$ git clone git@github.com:DavAug/chi.git +$ cd chi +$ pip install -e .[docs] +``` + +This will + +1. Install all the dependencies for Chi, including the ones for documentation (docs). +2. Tell Python to use your local Chi files when you use `import chi` anywhere on your system. + +At this point you will already have to have [Sundials](https://computing.llnl.gov/projects/sundials) installed, as explained in the non-dev [install instructions](https://chi.readthedocs.io/). + +You may also want to create a virtual environment first, using [virtualenv](https://docs.python.org/3/tutorial/venv.html) or [conda](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html). + +## Coding style guidelines + +Chi follows the [PEP8 recommendations](https://www.python.org/dev/peps/pep-0008/) for coding style. These are very common guidelines, and community tools have been developed to check how well projects implement them. + +We use [flake8](http://flake8.pycqa.org/en/latest/) to check our PEP8 adherence. To try this on your system, navigate to the Chi directory in a console and type + +``` +$ flake8 +``` + +When you commit your changes they will be checked against flake8 automatically (see [infrastructure](#infrastructure)). + +### Naming + +Naming is hard. In general, we aim for descriptive class, method, and argument names. Avoid abbreviations when possible without making names overly long, so `mean` is better than `mu`, but a class name like `AdaptiveMCMC` is fine. + +Class names are CamelCase, and start with an upper case letter, for example `SuperDuperMCMC`. Method and variable names are lower case, and use underscores for word separation, for example `x` or `iteration_count`. + +### Spelling + +To be consistent with the work so far, all Chi material in the repository (code, comments, docstrings, documentation, etc.) should be written in UK english, the only exception being when quoting other sources, e.g. titles of scientific articles in references. + +## Dependencies and reusing code + +While it's a bad idea for developers to "reinvent the wheel", it's important for users to get a _reasonably sized download and an easy install_. In addition, external libraries can sometimes cease to be supported, and when they contain bugs it might take a while before fixes become available as automatic downloads to Chi users. +For these reasons, all dependencies in Chi should be thought about carefully, and discussed on GitHub. + +Direct inclusion of code from other packages is possible, as long as their license permits it and is compatible with ours, but again should be considered carefully and discussed in the group. Snippets from blogs and stackoverflow can often be included without attribution, but if they solve a particularly nasty problem (or are very hard to read) it's often a good idea to attribute (and document) them, by making a comment with a link in the source code. + +## Testing + +All code requires testing. We use the [unittest](https://docs.python.org/3.3/library/unittest.html) package for our tests. (These tests typically just check that the code runs without error, and so, are more _debugging_ than _testing_ in a strict sense. Nevertheless, they are very useful to have!) + +To run quick tests, use + +``` +$ python run-tests.py --unit +``` + +### Writing tests + +Every new feature should have its own test. To create ones, have a look at the `test` directory and see if there's a test for a similar method. Copy-pasting this is a good way to start. + +Next, add some simple (and speedy!) tests of your main features. If these run without exceptions that's a good start! Next, check the output of your methods using any of these [assert methods](https://docs.python.org/3.3/library/unittest.html#assert-methods). + +Guidelines for writing unit tests: + +1. Unit tests should test a very small block of code (e.g. a single method) +1b. When writing tests, start from the simplest case, and then work up +2. Unit tests _test the public API_ (i.e. they never access private methods or variables). They test _if_ the code does what it promises to do, not _how_ it does it. For example, after running `my_object.set_x(4)`, you might check if `my_object.x()` returns 4, but you never check if `my_object._x == 4`: how the object stores its data is its own business. +3. There are hundreds of unit tests, and good developers run all of them several times a day. Therefore, unit tests should be _fast_. +4. If you're testing something stochastic, seed the number generator as part of the test, i.e. with `np.random.seed(1)`. + +## Documentation + +Chi is documented in several ways. + +First and foremost, every method and every class should have a [docstring](https://www.python.org/dev/peps/pep-0257/) that describes in plain terms what it does, and what the expected input and output is. + +These docstrings can be fairly simple, but can also make use of [reStructuredText](http://docutils.sourceforge.net/docs/user/rst/quickref.html), a markup language designed specifically for writing [technical documentation](https://en.wikipedia.org/wiki/ReStructuredText). For example, you can link to other classes and methods by writing ```:class:`chi.MechanisticModel` ``` and ```:meth:`simulate()` ```. + +In addition, we write a (very) small bit of documentation in separate reStructuredText files in the `docs` directory. Most of what these files do is simply import docstrings from the source code. But they also do things like add tables and indexes. If you've added a new class to a module, search the `docs` directory for that module's `.rst` file and add your class (in alphabetical order) to its index. If you've added a whole new module, copy-paste another module's file and add a link to your new file in the appropriate `index.rst` file. + +Using [Sphinx](http://www.sphinx-doc.org/en/stable/) the documentation in `docs` can be converted to HTML, PDF, and other formats. In particular, we use it to generate the documentation on http://pints.readthedocs.io/ + +### Docstring template + +1. Each docstring should start with a [single sentence](https://www.python.org/dev/peps/pep-0257/#one-line-docstrings) explaining what it does. + +2. If desired, [this can be followed by a blank line and one or several paragraphs](https://www.python.org/dev/peps/pep-0257/#multi-line-docstrings) providing a more detailed explanation. + These paragraphs can include code snippets or use LaTeX expressions for mathematics (see below). + +3. If the class is a subclass of some other Chi type, it may be good to + mention this here. For example: + + Extends :class:`MechanisticModel`. + +4. Simple arguments can be described textually. For example, a docstring could + be a single line "Sets the width parameter to `w`.". For complicated + functions or methods it may be good to include a parameters section: + + Parameters + ---------- + x : int + A variable `x` that should be an integer + y + A variable without a type hint + + This syntax can also be used for constructor arguments. + Note that default values for any arguments are already displayed + automatically in the function/method/constructor signature. + +5. Simple return types can be described textually, but complicated return types + (which are not encouraged) can use the syntax: + + Returns + ------- + samples + A list of samples. + likelihoods + A list of their corresponding log-likelihoods + +6. References to literature are highly encouraged, and go near the bottom of + the docstring: + + Adaptive covariance MCMC based on Haario et al. [1]_, [2]_. + + (rest of the docstring goes here) + + References + ---------- + .. [1] Johnstone, Chang, Bardenet, de Boer, Gavaghan, Pathmanathan, + Clayton, Mirams (2015) "Uncertainty and variability in models of + the cardiac action potential: Can we build trustworthy models?". + Journal of Molecular and Cellular Cardiology. + https://10.1016/j.yjmcc.2015.11.018 + + .. [2] Haario, Saksman, Tamminen (2001) "An adaptive Metropolis + algorithm". Bernoulli. + https://doi.org/10.2307/3318737 + + There is no standard format (e.g. APA style), but authors, titles, years, + and journals are recommended, as well as a link based on a + [DOI](https://www.doi.org/). + +7. Longer code snippets can go at the very end of a docstring + + Examples + -------- + :: + + errors = [ + pints.MeanSquaredError(problem1), + pints.MeanSquaredError(problem2), + ] + + # Equally weighted + e1 = pints.SumOfErrors(errors) + + # Differrent weights: + weights = [ + 1.0, + 2.6, + ] + e2 = pints.SumOfErrors(errors, weights) + +### Using code in documentation + +When referencing a variable in a docstring, please use the syntax ` ``x`` `. +Longer code snippets can be started using this form: + + """ + An example follows here:: + + print('Hello world') + + """ + +### Using Maths in documentation + +LaTeX expressions can be embedded in docstrings by using the syntax ```:math:`expression```` for inline mathematics, or a longer form for multi-line strings: + + r""" + Defines a :math:`\gamma` (log) prior with given shape parameter ``a`` + and rate parameter ``b``, with pdf + + .. math:: + f(x|a,b)=\frac{b^a x^{a-1} e^{-bx}}{\mathrm{\Gamma}(a)} + + """ + +Note that when using maths, it is best to define the docstring in a *raw string*, i.e. by writing ```r""" your stuff here """```. This will allow you to write e.g. `1 + \tau` instead of `1 + \\tau` and will stop flake8 from complaining about invalid escape sequences. + +### Building the documentation + +To test and debug the documentation, it's best to build it locally. To do this, make sure you have the relevant dependencies installed (see [installation](#installation)), navigate to your Chi directory in a console, and then type: + +``` +cd docs +make clean +make html +``` + +Next, open a browser, and navigate to your local PINTS directory (by typing the path, or part of the path into your location bar). Then have a look at `/docs/build/html/index.html`. + + +## Infrastructure + +### New releases on PyPI (for `pip`) + +Occasionally, we'll make a new release for Chi, and update the version on PyPI (which is where `pip` will download it from, for non-dev users). + +To do this: + +- Decide a new release is necessary, discuss it in the group. +- Make sure the version number has changed since the last release. +- Use the [GitHub releases page](https://github.com/DavAug/chi/releases/new) to create a new release. Each release will create a tag in the git repository, which should have the format `v1.2.3`. + - The first number is for big events, the second for regular releases (e.g. new features), the final for bugfixes and smaller improvements. This is subjective. + - Beyond that, there is no significance to these numbers (e.g. it doesn't matter if they're odd or even, `v0.9.9` is followed by `v0.9.10`). +- Check what has changed since the last release, and write some release notes to summarise what's new. +- Creating the new release in github **will automatically update PyPI**, so do this with care. + - Keep in mind that PyPI version numbers are eternal: You cannot modify a release, only create a new one with a new version number. +- Once the new release is done, create a PR to update the version number (final digit) to indicate that the code in the repo is no longer the version on PIP. + + +### Setuptools + +Installation of Chi _and dependencies_ is handled via [setuptools](http://setuptools.readthedocs.io/) + +Configuration files: + +``` +setup.py +MANIFEST.in +``` + +The `MANIFEST.in` file is used to list non-Python files to be included in the installation. + +### PIP + +It's always worth using an up-to-date version of pip. On older systems especially, having an up-to-date pip will prevent all kinds of version incompatibility issues: + +``` +$ pip install --upgrade pip +``` + +### GitHub Actions + +All committed code is tested using [GitHub Actions](https://help.github.com/en/actions), tests are published on https://github.com/pints-team/pints/actions. + +Configuration files: + +``` +.github/workflows/*.yml +``` + +Unit tests and flake8 testing is done for every commit. A nightly cronjob also tests the notebooks. + +### Codecov + +Code coverage (how much of our code is actually seen by the (linux) unit tests) is tested using [Codecov](https://docs.codecov.io/), a report is visible on https://codecov.io/gh/DavAug/chi. + +It is possible to measure code coverage locally using [coverage.py](https://coverage.readthedocs.io/en/coverage-5.5/). To run a particular test file (below `test_mechanistic_models.py`) and record coverage, run: + +``` +coverage run -m unittest test_mechanistic_models.py +``` + +To see the coverage for a particular file (below `_mechanistic_models.py`) triggered by this test, run: + +``` +coverage report -m _mechanistic_models.py +``` + +Configuration files: + +``` +.coveragerc +``` + +### Read the Docs + +Documentation is built using https://readthedocs.org/ and published on http://pints.readthedocs.io/. + +Configuration files: + +``` +.readthedocs.yaml +.requirements-docs.txt +``` + +### Flake8 + +[Style checking](#coding-style-guidelines) is performed using [flake8](http://flake8.pycqa.org/en/latest/). + +Configuration files: + +``` +.flake8 +``` + +### GitHub + +GitHub does some magic with particular filesnames. In particular: + +- The first page people see when they go to [our GitHub page](https://github.com/DavAug/chi) displays the contents of [README.md](README.md), which is written in the [Markdown](https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet) format. Some guidelines can be found [here](https://help.github.com/articles/about-readmes/). +- The license for using Chi is stored in [LICENSE.md](LICENSE.md), and [automatically](https://help.github.com/articles/adding-a-license-to-a-repository/) linked to by GitHub. +- This file, [CONTRIBUTING.md](CONTRIBUTING.md) is recognised as the contribution guidelines and a link is [automatically](https://github.com/blog/1184-contributing-guidelines) displayed when new issues or pull requests are created. + +## Acknowledgements + +We want to thank the [PINTS](https://github.com/pints-team/pints) developer team whose `CONTRIBUTING.md` helped us as guidance for this document. \ No newline at end of file diff --git a/README.md b/README.md index 71342535..47b82aef 100644 --- a/README.md +++ b/README.md @@ -43,11 +43,11 @@ pip install chi-drm ```python import chi ``` - + ### Tutorials Tutorials and more detailed explanations on how to use chi can be found in the [documentation's getting started](https://chi.readthedocs.io/en/latest/getting_started/index.html) section. - + ## Citation If you use this software in your work, please cite it using the following metadata: @@ -70,9 +70,11 @@ year = {2021} ``` ## Contributing -Pull requests are welcome. For major changes, please open an issue first to discuss what you would like to change. +There are lots of ways how you can contribute to Chi's development, and you are welcome to join in! +For example, you can report problems or make feature requests on the [issues](https://github.com/DavAug/chi/issues) pages. -Please make sure to update tests as appropriate. +Similarly, if you would like to contribute documentation or code you can create and issue, and then provide a pull request for review. +To facilitate contributions, we have written extensive [contribution guidelines](https://github.com/DavAug/chi/blob/main/CONTRIBUTING.md) to help you navigate the code. ## License [BSD-3-Clause](https://opensource.org/licenses/BSD-3-Clause) diff --git a/chi/_log_pdfs.py b/chi/_log_pdfs.py index 5323ce7d..c47cd8f5 100644 --- a/chi/_log_pdfs.py +++ b/chi/_log_pdfs.py @@ -25,14 +25,14 @@ class HierarchicalLogLikelihood(object): .. math:: \log p(\mathcal{D}, \Psi | \theta ) = - \sum _{ij} \log p(y_{ij} | \psi _{i} , t_{ij}) + \sum _{i} \log p(\mathcal{D}_i | \psi _{i}) + \sum _{i} \log p(\psi _{i}| \theta), where the first term is the sum over the log-likelihoods and the second term is the log-likelihood of the population model parameters. - :math:`\mathcal{D}=\{ (y_{ij}, t_{ij})\}` is the - data, where :math:`(y_{ij}, t_{ij})` is the :math:`j^{\mathrm{th}}` - measurement of log-likelihood :math:`i`. :math:`\Psi = \{ \psi_i\}` + :math:`\mathcal{D}=\{ \mathcal{D}_i\}` is the + data, where :math:`\mathcal{D}_i) = \{(y_{ij}, t_{ij})\}` denotes the + measurement from log-likelihood :math:`i`. :math:`\Psi = \{ \psi_i\}` denotes the parameters across the individual log-likelihoods. :param log_likelihoods: A list of log-likelihoods which are @@ -645,7 +645,7 @@ class LogLikelihood(pints.LogPDF): defined by ``observations`` and ``times`` .. math:: - p(\psi | \mathcal{D}) = \sum _{j=1} + p(\mathcal{D} | \psi) = \sum _{j=1} \log p(y_j | \psi, t_j), where :math:`\mathcal{D} = \{(y_j , t_j)\}` denotes the measurements. @@ -1157,9 +1157,9 @@ class LogPosterior(pints.LogPDF): an instance of a :class:`pints.LogPrior` of the same dimensionality as parameters in the log-likelihood. - Formally the log-posterior is given by the sum of the input log-likelihood - :math:`L(\psi | x^{\text{obs}})` and the input log-prior - :math:`\log p(\psi )` up to an additive constant + Formally the log-posterior is given by the sum of the log-likelihood, + :math:`\log p(\mathcal{D} | \psi)`, and the log-prior, + :math:`\log p(\psi )`, up to an additive constant .. math:: \log p(\psi | \mathcal{D}) = @@ -1167,7 +1167,7 @@ class LogPosterior(pints.LogPDF): where :math:`\psi` are the parameters of the log-likelihood and :math:`\mathcal{D}` are the observed data. The additive constant - is the normalisation of the log-posterior and is in general not known. + is the normalisation of the log-posterior and is in general unknown. Extends :class:`pints.LogPDF`. diff --git a/chi/_mechanistic_models.py b/chi/_mechanistic_models.py index e8378b13..972024a6 100644 --- a/chi/_mechanistic_models.py +++ b/chi/_mechanistic_models.py @@ -504,8 +504,9 @@ def simulate(self, parameters, times): values are returned. :type times: list, numpy.ndarray - :rtype: np.ndarray of shape (n_outputs, n_times) or - (n_times, n_outputs, n_parameters) + :rtype: np.ndarray of shape (n_outputs, n_times) or a Tuple of two + np.ndarray, the first of shape (n_outputs, n_times) and the second + of shape (n_times, n_outputs, n_parameters) """ # Reset simulation self._simulator.reset() @@ -622,7 +623,8 @@ def _add_dose_rate(self, compartment, drug_amount): # Set initial value to 0 and unit to unit of drug amount over unit of # time dose_rate.set_rhs(0) - dose_rate.set_unit(drug_amount.unit() / self.time_unit()) + if drug_amount.unit(): + dose_rate.set_unit(drug_amount.unit() / self.time_unit()) # Add the dose rate to the rhs of the drug amount variable rhs = drug_amount.rhs() diff --git a/chi/_predictive_models.py b/chi/_predictive_models.py index 67adaf1b..f2fe3803 100644 --- a/chi/_predictive_models.py +++ b/chi/_predictive_models.py @@ -263,7 +263,14 @@ def sample( individual = None # Check individual for individual predictive model else: - ids = self._posterior.individual + try: + ids = self._posterior.individual + except AttributeError: # pragma: no cover + # Posterior only contains information about one individual. + # Hack to make this work + individual = 'some value' + ids = ['some value'] + # Get default individual if individual is None: individual = str(ids.data[0]) diff --git a/docs/source/api/covariate_models.rst b/docs/source/api/covariate_models.rst index 326810df..f0eb1222 100644 --- a/docs/source/api/covariate_models.rst +++ b/docs/source/api/covariate_models.rst @@ -14,8 +14,8 @@ from are used to describe the variability within a subpopulation, while covariate models are used to describe the cross-subpopulation variability. -Functional classes ------------------- +Classes +------- - :class:`LinearCovariateModel` diff --git a/docs/source/api/error_models.rst b/docs/source/api/error_models.rst index acf92b35..d9b4de78 100644 --- a/docs/source/api/error_models.rst +++ b/docs/source/api/error_models.rst @@ -9,8 +9,8 @@ Error Models Error models in chi_ model the deviations between experimental observations and the :class:`MechanisticModel` predictions. -Functional classes ------------------- +Classes +------- - :class:`ConstantAndMultiplicativeGaussianErrorModel` - :class:`GaussianErrorModel` diff --git a/docs/source/api/inference.rst b/docs/source/api/inference.rst index b44d2ae1..b423d470 100644 --- a/docs/source/api/inference.rst +++ b/docs/source/api/inference.rst @@ -15,8 +15,8 @@ to easily explore different optimisation or sampling settings, e.g. using differ methods, fixing some parameters, or applying different transformations to the search space. -Functional classes ------------------- +Classes +------- - :class:`OptimisationController` - :class:`SamplingController` diff --git a/docs/source/api/log_pdfs.rst b/docs/source/api/log_pdfs.rst index 211bd1ab..75376c4d 100644 --- a/docs/source/api/log_pdfs.rst +++ b/docs/source/api/log_pdfs.rst @@ -10,8 +10,8 @@ Log-PDFs Log-PDFs in chi_ extend the log-pdfs provided in pints_ for the purpose of PKPD modelling. -Functional classes ------------------- +Classes +------- - :class:`HierarchicalLogLikelihood` - :class:`HierarchicalLogPosterior` diff --git a/docs/source/api/population_filters.rst b/docs/source/api/population_filters.rst index 40aac9a4..328c0b5f 100644 --- a/docs/source/api/population_filters.rst +++ b/docs/source/api/population_filters.rst @@ -8,8 +8,8 @@ Population Filters Population filters in chi_ can be used for population filter inference. -Functional classes ------------------- +Classes +------- - :class:`ComposedPopulationFilter` - :class:`GaussianFilter` diff --git a/docs/source/api/population_models.rst b/docs/source/api/population_models.rst index a3c7f886..0c7fde23 100644 --- a/docs/source/api/population_models.rst +++ b/docs/source/api/population_models.rst @@ -9,8 +9,8 @@ Population Models Population models in chi_ can be used to model the variation of mechanistic model parameters or error model parameters across individuals. -Functional classes ------------------- +Classes +------- - :class:`ComposedPopulationModel` - :class:`CovariatePopulationModel` diff --git a/docs/source/api/predictive_models.rst b/docs/source/api/predictive_models.rst index 4d5ae63e..1f5e0eb4 100644 --- a/docs/source/api/predictive_models.rst +++ b/docs/source/api/predictive_models.rst @@ -12,8 +12,8 @@ generate synthetic measurements. Each predictive model consists of a :class:`MechanisticModel` and an :class:`ErrorModel`. -Functional classes ------------------- +Classes +------- - :class:`PAMPredictiveModel` - :class:`PopulationPredictiveModel` diff --git a/docs/source/api/problems.rst b/docs/source/api/problems.rst index 6141517a..59da7d75 100644 --- a/docs/source/api/problems.rst +++ b/docs/source/api/problems.rst @@ -13,8 +13,8 @@ individuals separately or as a population, the :class:`ProblemModellingController` has been implemented. -Functional classes ------------------- +Classes +------- - :class:`ProblemModellingController` diff --git a/docs/source/conf.py b/docs/source/conf.py index 486b319e..5f26a865 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -37,6 +37,7 @@ 'sphinx.ext.mathjax', 'sphinx.ext.napoleon', 'sphinx.ext.viewcode', + 'sphinx_copybutton' ] # Autodoc defaults @@ -66,7 +67,7 @@ # General information about the project. project = u'Chi' -copyright = u'2022, David Augustin' +copyright = u'2023, David Augustin' author = u'David Augustin' # The version info for the project you're documenting, acts as replacement for diff --git a/docs/source/getting_started/code/2_mechanistic_model_1.py b/docs/source/getting_started/code/2_mechanistic_model_1.py new file mode 100644 index 00000000..3c3008e0 --- /dev/null +++ b/docs/source/getting_started/code/2_mechanistic_model_1.py @@ -0,0 +1,166 @@ +import argparse + + +# Set up argument parsing, so plotting and exports can be disabled for +# testing. +parser = argparse.ArgumentParser( + description='Run example scripts for chi docs.', +) +parser.add_argument( + '--test', + action='store_true', + help='Run testing version of script which ignores plotting.',) + +# Parse! +args = parser.parse_args() + +# Run testing verion of script +if args.test: + #1 + import chi + import numpy as np + + + class OneCompPKModel(chi.MechanisticModel): + def __init__(self): + super().__init__() + + def simulate(self, parameters, times): + a0, ke, v = parameters + times = np.array(times) + + c = a0 / v * np.exp(-ke * times) + + output = np.empty(shape=(self.n_outputs(), len(times))) + output[0] = c + + return output + + def has_sensitivities(self): + # We will not implement sensitivities in this part of the tutorial, so + # the output of this method is always False + return False + + def n_outputs(self): + return 1 + + def outputs(self): + return ['Drug concentration'] + + def n_parameters(self): + return 3 + + def parameters(self): + return ['Initial amount', 'Elimination rate', 'Volume of distribution'] + + + import os + # 2 + import plotly.graph_objects as go + + + # Define model + model = OneCompPKModel() + + # Run simulation + times = np.linspace(start=0, stop=10, num=200) + parameters = [ + 10, # Initial drug amount + 1, # Elimination rate + 2, # Volume of distribution + ] + simulation = model.simulate(parameters=parameters, times=times)[0] + + # Plot drug concentration + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/2_mechanistic_model_1.html') + + # Exit script + exit() + +#1 +import chi +import numpy as np + + +class OneCompPKModel(chi.MechanisticModel): + def __init__(self): + super().__init__() + + def simulate(self, parameters, times): + a0, ke, v = parameters + times = np.array(times) + + c = a0 / v * np.exp(-ke * times) + + output = np.empty(shape=(self.n_outputs(), len(times))) + output[0] = c + + return output + + def has_sensitivities(self): + # We will not implement sensitivities in this part of the tutorial, so + # the output of this method is always False + return False + + def n_outputs(self): + return 1 + + def outputs(self): + return ['Drug concentration'] + + def n_parameters(self): + return 3 + + def parameters(self): + return ['Initial amount', 'Elimination rate', 'Volume of distribution'] + + +import os +# 2 +import plotly.graph_objects as go + + +# Define model +model = OneCompPKModel() + +# Run simulation +times = np.linspace(start=0, stop=10, num=200) +parameters = [ + 10, # Initial drug amount + 1, # Elimination rate + 2, # Volume of distribution +] +simulation = model.simulate(parameters=parameters, times=times)[0] + +# Plot drug concentration +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/2_mechanistic_model_1.html') \ No newline at end of file diff --git a/docs/source/getting_started/code/2_mechanistic_model_2.py b/docs/source/getting_started/code/2_mechanistic_model_2.py new file mode 100644 index 00000000..3f816b4c --- /dev/null +++ b/docs/source/getting_started/code/2_mechanistic_model_2.py @@ -0,0 +1,496 @@ +import argparse + + +# Set up argument parsing, so plotting and exports can be disabled for +# testing. +parser = argparse.ArgumentParser( + description='Run example scripts for chi docs.', +) +parser.add_argument( + '--test', + action='store_true', + help='Run testing version of script which ignores plotting.',) + +# Parse! +args = parser.parse_args() + +# Run testing verion of script +if args.test: + # 1 + import os + + import chi + + + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'template.xml') + model = chi.PKPDModel(sbml_file=filename) + + # print(model._model.name()) + + + # 2 + import os + + import chi + + + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'one_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + + # print(model._model.code()) + # print(model.parameters()) + + # 3 + import os + + import chi + import numpy as np + import plotly.graph_objects as go + + + # Define model + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'one_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + model.set_outputs(['global.drug_concentration']) + + # Run simulation + times = np.linspace(start=0, stop=10, num=200) + parameters = [ + 10, # Initial drug amount + 1, # Elimination rate + 2, # Volume of distribution + ] + simulation = model.simulate(parameters=parameters, times=times)[0] + + # Plot drug concentation + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/2_mechanistic_model_2.html') + + + # 4 + model.set_administration(compartment='global') + + # print(model._model.code()) + + # 5 + times = np.linspace(start=0, stop=10, num=200) + parameters = [0, 1, 2] + + # 2 mg every day + model.set_dosing_regimen(dose=2, period=1) + sim1 = model.simulate(parameters=parameters, times=times)[0] + + # 4 mg every 2 days + model.set_dosing_regimen(dose=4, period=2) + sim2 = model.simulate(parameters=parameters, times=times)[0] + + # 2 mg every day with an infusion of 12 hours + model.set_dosing_regimen(dose=2, period=1, duration=0.5) + sim3 = model.simulate(parameters=parameters, times=times)[0] + + # Plot drug concentations + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=sim1, + mode='lines', + name='2mg every day' + )) + fig.add_trace(go.Scatter( + x=times, + y=sim2, + mode='lines', + name='4mg every 2 days' + )) + fig.add_trace(go.Scatter( + x=times, + y=sim3, + mode='lines', + name='2mg every day infused for 12h' + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/2_mechanistic_model_3.html') + + # 6 + # Define model + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'one_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + model.set_outputs(['global.drug_concentration']) + + # Set indirect administration + model.set_administration(compartment='global', direct=False) + + times = np.linspace(start=0, stop=10, num=200) + parameters = [0, 0, 10, 1, 2] + + # 2 mg every day + model.set_dosing_regimen(dose=2, period=1) + sim1 = model.simulate(parameters=parameters, times=times)[0] + + # 4 mg every 2 days + model.set_dosing_regimen(dose=4, period=2) + sim2 = model.simulate(parameters=parameters, times=times)[0] + + # 2 mg every day with an infusion of 12 hours + model.set_dosing_regimen(dose=2, period=1, duration=0.5) + sim3 = model.simulate(parameters=parameters, times=times)[0] + + # Plot drug concentations + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=sim1, + mode='lines', + name='2mg every day' + )) + fig.add_trace(go.Scatter( + x=times, + y=sim2, + mode='lines', + name='4mg every 2 days' + )) + fig.add_trace(go.Scatter( + x=times, + y=sim3, + mode='lines', + name='2mg every day infused for 12h' + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/2_mechanistic_model_4.html') + + # Print model + model = chi.PKPDModel(sbml_file=filename) + model.set_outputs(['global.drug_concentration']) + + # Set indirect administration + model.set_administration(compartment='global', direct=False) + + # print(model._model.code()) + + # 7 + from plotly.subplots import make_subplots + + + # Simulate concentration and sensitivities + model.set_dosing_regimen(dose=2, period=1) + model.enable_sensitivities(True) + simulation, sensitivities = model.simulate(parameters=parameters, times=times) + + # Plot results + fig = make_subplots(rows=3, cols=1, shared_xaxes=True) + fig.add_trace( + go.Scatter( + x=times, + y=simulation[0], + mode='lines', + name='Drug concentration' + ), + row=1, + col=1 + ) + fig.add_trace( + go.Scatter( + x=times, + y=sensitivities[:, 0, 1], + mode='lines', + name='Sensitivity w.r.t. initial amount' + ), + row=2, + col=1 + ) + fig.add_trace( + go.Scatter( + x=times, + y=sensitivities[:, 0, 3], + mode='lines', + name='Sensitivity w.r.t. elim. rate' + ), + row=3, + col=1 + ) + + fig.update_layout( + xaxis3_title='Time', + yaxis_title='c', + yaxis2_title='dc / da_0', + yaxis3_title='dc / dk_e', + template='plotly_white' + ) + # fig.show() + + # fig.write_html(directory + '/images/2_mechanistic_model_5.html') + + # Exit script + exit() + +# 1 +import os + +import chi + + +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'template.xml') +model = chi.PKPDModel(sbml_file=filename) + +print(model._model.name()) + + +# 2 +import os + +import chi + + +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'one_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) + +print(model._model.code()) +print(model.parameters()) + +# 3 +import os + +import chi +import numpy as np +import plotly.graph_objects as go + + +# Define model +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'one_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) +model.set_outputs(['global.drug_concentration']) + +# Run simulation +times = np.linspace(start=0, stop=10, num=200) +parameters = [ + 10, # Initial drug amount + 1, # Elimination rate + 2, # Volume of distribution +] +simulation = model.simulate(parameters=parameters, times=times)[0] + +# Plot drug concentation +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/2_mechanistic_model_2.html') + + +# 4 +model.set_administration(compartment='global') + +print(model._model.code()) + +# 5 +times = np.linspace(start=0, stop=10, num=200) +parameters = [0, 1, 2] + +# 2 mg every day +model.set_dosing_regimen(dose=2, period=1) +sim1 = model.simulate(parameters=parameters, times=times)[0] + +# 4 mg every 2 days +model.set_dosing_regimen(dose=4, period=2) +sim2 = model.simulate(parameters=parameters, times=times)[0] + +# 2 mg every day with an infusion of 12 hours +model.set_dosing_regimen(dose=2, period=1, duration=0.5) +sim3 = model.simulate(parameters=parameters, times=times)[0] + +# Plot drug concentations +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=sim1, + mode='lines', + name='2mg every day' +)) +fig.add_trace(go.Scatter( + x=times, + y=sim2, + mode='lines', + name='4mg every 2 days' +)) +fig.add_trace(go.Scatter( + x=times, + y=sim3, + mode='lines', + name='2mg every day infused for 12h' +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/2_mechanistic_model_3.html') + +# 6 +# Define model +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'one_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) +model.set_outputs(['global.drug_concentration']) + +# Set indirect administration +model.set_administration(compartment='global', direct=False) + +times = np.linspace(start=0, stop=10, num=200) +parameters = [0, 0, 10, 1, 2] + +# 2 mg every day +model.set_dosing_regimen(dose=2, period=1) +sim1 = model.simulate(parameters=parameters, times=times)[0] + +# 4 mg every 2 days +model.set_dosing_regimen(dose=4, period=2) +sim2 = model.simulate(parameters=parameters, times=times)[0] + +# 2 mg every day with an infusion of 12 hours +model.set_dosing_regimen(dose=2, period=1, duration=0.5) +sim3 = model.simulate(parameters=parameters, times=times)[0] + +# Plot drug concentations +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=sim1, + mode='lines', + name='2mg every day' +)) +fig.add_trace(go.Scatter( + x=times, + y=sim2, + mode='lines', + name='4mg every 2 days' +)) +fig.add_trace(go.Scatter( + x=times, + y=sim3, + mode='lines', + name='2mg every day infused for 12h' +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/2_mechanistic_model_4.html') + +# Print model +model = chi.PKPDModel(sbml_file=filename) +model.set_outputs(['global.drug_concentration']) + +# Set indirect administration +model.set_administration(compartment='global', direct=False) + +print(model._model.code()) + +# 7 +from plotly.subplots import make_subplots + + +# Simulate concentration and sensitivities +model.set_dosing_regimen(dose=2, period=1) +model.enable_sensitivities(True) +simulation, sensitivities = model.simulate(parameters=parameters, times=times) + +# Plot results +fig = make_subplots(rows=3, cols=1, shared_xaxes=True) +fig.add_trace( + go.Scatter( + x=times, + y=simulation[0], + mode='lines', + name='Drug concentration' + ), + row=1, + col=1 +) +fig.add_trace( + go.Scatter( + x=times, + y=sensitivities[:, 0, 1], + mode='lines', + name='Sensitivity w.r.t. initial amount' + ), + row=2, + col=1 +) +fig.add_trace( + go.Scatter( + x=times, + y=sensitivities[:, 0, 3], + mode='lines', + name='Sensitivity w.r.t. elim. rate' + ), + row=3, + col=1 +) + +fig.update_layout( + xaxis3_title='Time', + yaxis_title='c', + yaxis2_title='dc / da_0', + yaxis3_title='dc / dk_e', + template='plotly_white' +) +fig.show() + +fig.write_html(directory + '/images/2_mechanistic_model_5.html') \ No newline at end of file diff --git a/docs/source/getting_started/code/3_fitting_models_1.py b/docs/source/getting_started/code/3_fitting_models_1.py new file mode 100644 index 00000000..3adcae4e --- /dev/null +++ b/docs/source/getting_started/code/3_fitting_models_1.py @@ -0,0 +1,198 @@ +import argparse + + +# Set up argument parsing, so plotting and exports can be disabled for +# testing. +parser = argparse.ArgumentParser( + description='Run example scripts for chi docs.', +) +parser.add_argument( + '--test', + action='store_true', + help='Run testing version of script which ignores plotting.',) + +# Parse! +args = parser.parse_args() + +# Run testing verion of script +if args.test: + # 1 + import os + + import chi + import numpy as np + import plotly.graph_objects as go + + + # Implement the model + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'one_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + model.set_outputs(['global.drug_concentration']) + + # Set administration and dosing regimen + model.set_administration(compartment='global', direct=False) + model.set_dosing_regimen(dose=2, period=1) + + # Simulate treatment response + parameters = [0, 0, 10, 1, 2] + times = np.linspace(start=0, stop=3, num=200) + conc = model.simulate(parameters=parameters, times=times)[0] + + # Plot results + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=conc, + mode='lines', + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_1.html') + + + ### Generate synthetic data + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'two_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + model.set_administration(compartment='central', direct=False) + model.set_outputs(['central.drug_concentration']) + + # Synthethise data + error_model = chi.LogNormalErrorModel() + model = chi.PredictiveModel( + mechanistic_model=model, error_models=[error_model]) + model.set_dosing_regimen(dose=2, period=1) + + parameters = [0, 0, 0, 2, 10, 1, 5, 1, 10, 0.1] + times = [0.5, 1, 1.5, 2, 2.5, 3] + df = model.sample( + parameters=parameters, times=times, seed=1, include_regimen=True) + + # Save data (Only do this once, because I added some units changed names) + # directory = os.path.dirname(__file__) + # df.to_csv(directory + '/dataset_1.csv', index=False) + + + # 2 + import pandas as pd + + + # Import data + directory = os.path.dirname(__file__) + data = pd.read_csv(directory + '/dataset_1.csv') + + # Plot data + fig = go.Figure() + fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', + )) + fig.update_layout( + xaxis_title='Time in day', + yaxis_title='Drug concentration in ng/mL', + template='plotly_white' + ) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_2.html') + + # Exit script + exit() + +# 1 +import os + +import chi +import numpy as np +import plotly.graph_objects as go + + +# Implement the model +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'one_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) +model.set_outputs(['global.drug_concentration']) + +# Set administration and dosing regimen +model.set_administration(compartment='global', direct=False) +model.set_dosing_regimen(dose=2, period=1) + +# Simulate treatment response +parameters = [0, 0, 10, 1, 2] +times = np.linspace(start=0, stop=3, num=200) +conc = model.simulate(parameters=parameters, times=times)[0] + +# Plot results +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=conc, + mode='lines', +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_1.html') + + +### Generate synthetic data +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'two_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) +model.set_administration(compartment='central', direct=False) +model.set_outputs(['central.drug_concentration']) + +# Synthethise data +error_model = chi.LogNormalErrorModel() +model = chi.PredictiveModel( + mechanistic_model=model, error_models=[error_model]) +model.set_dosing_regimen(dose=2, period=1) + +parameters = [0, 0, 0, 2, 10, 1, 5, 1, 10, 0.1] +times = [0.5, 1, 1.5, 2, 2.5, 3] +df = model.sample( + parameters=parameters, times=times, seed=1, include_regimen=True) + +# Save data (Only do this once, because I added some units changed names) +# directory = os.path.dirname(__file__) +# df.to_csv(directory + '/dataset_1.csv', index=False) + + +# 2 +import pandas as pd + + +# Import data +directory = os.path.dirname(__file__) +data = pd.read_csv(directory + '/dataset_1.csv') + +# Plot data +fig = go.Figure() +fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', +)) +fig.update_layout( + xaxis_title='Time in day', + yaxis_title='Drug concentration in ng/mL', + template='plotly_white' +) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_2.html') \ No newline at end of file diff --git a/docs/source/getting_started/code/3_fitting_models_2.py b/docs/source/getting_started/code/3_fitting_models_2.py new file mode 100644 index 00000000..a0a4e151 --- /dev/null +++ b/docs/source/getting_started/code/3_fitting_models_2.py @@ -0,0 +1,1259 @@ +import argparse + + +# Set up argument parsing, so plotting and exports can be disabled for +# testing. +parser = argparse.ArgumentParser( + description='Run example scripts for chi docs.', +) +parser.add_argument( + '--test', + action='store_true', + help='Run testing version of script which ignores plotting.',) + +# Parse! +args = parser.parse_args() + +# Run testing verion of script +if args.test: + import os + + import chi + import numpy as np + import pandas as pd + import pints + import plotly.graph_objects as go + + # 3 + import pints + + + # Define mechanistic model + directory = os.path.dirname(__file__) + filename = os.path.join(directory, 'one_compartment_pk_model.xml') + model = chi.PKPDModel(sbml_file=filename) + model.set_administration(compartment='global', direct=False) + model.set_outputs(['global.drug_concentration']) + + # Define error model + error_model = chi.LogNormalErrorModel() + + # Define data + directory = os.path.dirname(__file__) + data = pd.read_csv(directory + '/dataset_1.csv') + + # Define prior distribution + log_prior = pints.ComposedLogPrior( + pints.GaussianLogPrior(mean=10, sd=2), # absorption rate + pints.GaussianLogPrior(mean=6, sd=2), # elimination rate + pints.LogNormalLogPrior(log_mean=0, scale=1), # volume of distribution + pints.LogNormalLogPrior(log_mean=-2, scale=0.5) # scale of meas. distrib. + ) + + # Define log-posterior using the ProblemModellingController + problem = chi.ProblemModellingController( + mechanistic_model=model, error_models=[error_model]) + problem.set_data( + data=data, + output_observable_dict={'global.drug_concentration': 'Drug concentration'} + ) + problem.fix_parameters(name_value_dict={ + 'dose.drug_amount': 0, + 'global.drug_amount': 0, + }) + problem.set_log_prior(log_prior=log_prior) + log_posterior = problem.get_log_posterior() + + + # 4 + # Run MCMC algorithm + n_iterations = 1000 + controller = chi.SamplingController(log_posterior=log_posterior, seed=1) + controller.set_n_runs(n_runs=3) + controller.set_parallel_evaluation(False) + controller.set_sampler(pints.HaarioBardenetACMC) + samples = controller.run(n_iterations=n_iterations, log_to_screen=False) + + # 5 + from plotly.colors import qualitative + from plotly.subplots import make_subplots + + + # Plot results + fig = make_subplots( + rows=4, cols=2, vertical_spacing=0.15, horizontal_spacing=0.1) + + # Plot traces and histogram of parameter + fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['dose.absorption_rate'].values[:, n_iterations//2::(n_iterations//2)//1].flatten(), + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[4], + ), + row=1, + col=1 + ) + fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + ), + row=1, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + ), + row=1, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + ), + row=1, + col=2 + ) + + # Plot traces and histogram of parameter + fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['global.elimination_rate'].values[:, n_iterations//2::(n_iterations//2)//1].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.2), + marker_color=qualitative.Plotly[4], + ), + row=2, + col=1 + ) + fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=2, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=2, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=2, + col=2 + ) + + # Plot traces and histogram of parameter + fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['global.volume'].values[:, n_iterations//2::(n_iterations//2)//1].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[4], + ), + row=3, + col=1 + ) + fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=3, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=3, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=3, + col=2 + ) + + # Plot traces and histogram of parameter + fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['Sigma log'].values[:, n_iterations//2::(n_iterations//2)//1].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.02), + marker_color=qualitative.Plotly[4], + ), + row=4, + col=1 + ) + fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=4, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=4, + col=2 + ) + fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=4, + col=2 + ) + + # Visualise prior distribution + parameter_values = np.linspace(4, 16, num=200) + pdf_values = np.exp([ + log_prior._priors[0]([parameter_value]) + for parameter_value in parameter_values]) + fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + ), + row=1, + col=1 + ) + + parameter_values = np.linspace(0, 12, num=200) + pdf_values = np.exp([ + log_prior._priors[1]([parameter_value]) + for parameter_value in parameter_values]) + fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=2, + col=1 + ) + + parameter_values = np.linspace(0, 12, num=200) + pdf_values = np.exp([ + log_prior._priors[2]([parameter_value]) + for parameter_value in parameter_values]) + fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=3, + col=1 + ) + + parameter_values = np.linspace(0, 0.6, num=200) + pdf_values = np.exp([ + log_prior._priors[3]([parameter_value]) + for parameter_value in parameter_values]) + fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=4, + col=1 + ) + + fig.update_layout( + xaxis_title='k_a', + yaxis_title='p', + yaxis2_title='k_a', + xaxis3_title='k_e', + yaxis3_title='p', + yaxis4_title='k_e', + xaxis5_title='v', + yaxis5_title='p', + yaxis6_title='v', + xaxis7_title='sigma', + yaxis7_title='p', + xaxis8_title='Number of iterations', + yaxis8_title='sigma', + template='plotly_white', + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1 + ), + margin=dict(t=0, r=0, l=0) + ) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_3.html') + + + # 6 + fig = go.Figure() + + # Plot histograms + fig.add_trace( + go.Histogram( + name='Run 1', + x=samples['dose.absorption_rate'].values[0, ::n_iterations//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[2], + ) + ) + fig.add_trace( + go.Histogram( + name='Run 2', + x=samples['dose.absorption_rate'].values[1, ::n_iterations//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[1], + ) + ) + fig.add_trace( + go.Histogram( + name='Run 3', + x=samples['dose.absorption_rate'].values[2, ::n_iterations//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[0], + ) + ) + + fig.update_layout( + template='plotly_white', + xaxis_title='Absorption rate in 1/day', + yaxis_title='Probability density', + barmode='overlay' + ) + fig.update_traces(opacity=0.75) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_4.html') + + + # 7 + fig = go.Figure() + + # Plot histograms + fig.add_trace( + go.Histogram( + name='Run 1', + x=samples['dose.absorption_rate'].values[0, n_iterations//2::(n_iterations//2)//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[2], + ) + ) + fig.add_trace( + go.Histogram( + name='Run 2', + x=samples['dose.absorption_rate'].values[1, n_iterations//2::n_iterations//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[1], + ) + ) + fig.add_trace( + go.Histogram( + name='Run 3', + x=samples['dose.absorption_rate'].values[2, n_iterations//2::n_iterations//1], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[0], + ) + ) + + fig.update_layout( + template='plotly_white', + xaxis_title='Absorption rate in 1/day', + yaxis_title='Probability density', + barmode='overlay' + ) + fig.update_traces(opacity=0.75) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_5.html') + + + # 8 + import arviz as az + + + # Summary for all samples + summary1 = az.summary(samples) + # print(summary1) + + # Summary for the samples post warmup + summary2 = az.summary(samples.sel(draw=slice(n_iterations//2, n_iterations))) + # print(summary2) + + # directory = os.path.dirname(__file__) + # summary1.to_csv(directory + '/3_fitting_models_summary_1.csv') + # summary2.to_csv(directory + '/3_fitting_models_summary_2.csv') + + + # 9 + # Fix model parameters + model = chi.ReducedMechanisticModel(mechanistic_model=model) + model.fix_parameters(name_value_dict={ + 'dose.drug_amount': 0, + 'global.drug_amount': 0, + }) + + # Set dosing regimen + dosing_regimen = problem.get_dosing_regimens()['1'] + model.set_dosing_regimen(dosing_regimen) + + # Extract model parameters from summary + parameters = [ + summary2['mean'].loc[parameter] for parameter in model.parameters() + ] + + # Simulate result + times = np.linspace(0.2, 3, 200) + simulation = model.simulate(parameters=parameters, times=times)[0] + + # Plot results + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + name='Fit' + )) + fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', + name='Measurements' + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' + ) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_6.html') + + + # 10 + # Define posterior predictive distribution + predictive_model = problem.get_predictive_model() + samples = samples.sel(draw=slice(n_iterations//2, n_iterations)) + posterior_model = chi.PosteriorPredictiveModel( + predictive_model=predictive_model, posterior_samples=samples) + + # Approximate distribution using sampling + n_samples = 10 + conc_samples = posterior_model.sample( + times=times, n_samples=n_samples, seed=1) + + # Reshape samples, so we can calculate mean and percentiles at the different + # time points + reshaped_samples = np.empty(shape=(n_samples, len(times))) + for sample_idx, sample_id in enumerate(conc_samples.ID.unique()): + reshaped_samples[ + sample_idx] = conc_samples[conc_samples.ID == sample_id].Value.values + + # Calculate mean, 5th and 95th percentile of the distribution at each time + # point + means = np.mean(reshaped_samples, axis=0) + lower = np.percentile(reshaped_samples, q=5, axis=0) + upper = np.percentile(reshaped_samples, q=95, axis=0) + + # Plot results + fig = go.Figure() + fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + name='Fit: mean parameters' + )) + fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', + name='Measurements' + )) + fig.add_trace(go.Scatter( + x=times, + y=means, + mode='lines', + name='Fit: mean posterior pred. distr.', + line=dict(color='black') + )) + fig.add_trace(go.Scatter( + x=times, + y=lower, + mode='lines', + name='Fit: 5th-95th perc. posterior pred. distr.', + line=dict(color='black', dash='dash') + )) + fig.add_trace(go.Scatter( + x=times, + y=upper, + mode='lines', + name='Fit: 5th-95th perc. posterior pred. distr.', + line=dict(color='black', dash='dash'), + showlegend=False + )) + fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white', + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1 + ), + ) + # fig.show() + + # directory = os.path.dirname(os.path.dirname(__file__)) + # fig.write_html(directory + '/images/3_fitting_models_7.html') + + # Exit script + exit() + + +import os + +import chi +import numpy as np +import pandas as pd +import pints +import plotly.graph_objects as go + +# 3 +import pints + + +# Define mechanistic model +directory = os.path.dirname(__file__) +filename = os.path.join(directory, 'one_compartment_pk_model.xml') +model = chi.PKPDModel(sbml_file=filename) +model.set_administration(compartment='global', direct=False) +model.set_outputs(['global.drug_concentration']) + +# Define error model +error_model = chi.LogNormalErrorModel() + +# Define data +directory = os.path.dirname(__file__) +data = pd.read_csv(directory + '/dataset_1.csv') + +# Define prior distribution +log_prior = pints.ComposedLogPrior( + pints.GaussianLogPrior(mean=10, sd=2), # absorption rate + pints.GaussianLogPrior(mean=6, sd=2), # elimination rate + pints.LogNormalLogPrior(log_mean=0, scale=1), # volume of distribution + pints.LogNormalLogPrior(log_mean=-2, scale=0.5) # scale of meas. distrib. +) + +# Define log-posterior using the ProblemModellingController +problem = chi.ProblemModellingController( + mechanistic_model=model, error_models=[error_model]) +problem.set_data( + data=data, + output_observable_dict={'global.drug_concentration': 'Drug concentration'} +) +problem.fix_parameters(name_value_dict={ + 'dose.drug_amount': 0, + 'global.drug_amount': 0, +}) +problem.set_log_prior(log_prior=log_prior) +log_posterior = problem.get_log_posterior() + + +# 4 +# Run MCMC algorithm +n_iterations = 20000 +controller = chi.SamplingController(log_posterior=log_posterior, seed=1) +controller.set_n_runs(n_runs=3) +controller.set_parallel_evaluation(False) +controller.set_sampler(pints.HaarioBardenetACMC) +samples = controller.run(n_iterations=n_iterations, log_to_screen=True) + +# 5 +from plotly.colors import qualitative +from plotly.subplots import make_subplots + + +# Plot results +fig = make_subplots( + rows=4, cols=2, vertical_spacing=0.15, horizontal_spacing=0.1) + +# Plot traces and histogram of parameter +fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['dose.absorption_rate'].values[:, n_iterations//2::(n_iterations//2)//1000].flatten(), + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[4], + ), + row=1, + col=1 +) +fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + ), + row=1, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + ), + row=1, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['dose.absorption_rate'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + ), + row=1, + col=2 +) + +# Plot traces and histogram of parameter +fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['global.elimination_rate'].values[:, n_iterations//2::(n_iterations//2)//1000].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.2), + marker_color=qualitative.Plotly[4], + ), + row=2, + col=1 +) +fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=2, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=2, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['global.elimination_rate'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=2, + col=2 +) + +# Plot traces and histogram of parameter +fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['global.volume'].values[:, n_iterations//2::(n_iterations//2)//1000].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[4], + ), + row=3, + col=1 +) +fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=3, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=3, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['global.volume'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=3, + col=2 +) + +# Plot traces and histogram of parameter +fig.add_trace( + go.Histogram( + name='Posterior samples', + legendgroup='Group 1', + x=samples['Sigma log'].values[:, n_iterations//2::(n_iterations//2)//1000].flatten(), + histnorm='probability density', + showlegend=False, + xbins=dict(size=0.02), + marker_color=qualitative.Plotly[4], + ), + row=4, + col=1 +) +fig.add_trace( + go.Scatter( + name='Run 1', + legendgroup='Group 2', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[0], + mode='lines', + line_color=qualitative.Plotly[2], + showlegend=False + ), + row=4, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 2', + legendgroup='Group 3', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[1], + mode='lines', + line_color=qualitative.Plotly[1], + showlegend=False + ), + row=4, + col=2 +) +fig.add_trace( + go.Scatter( + name='Run 3', + legendgroup='Group 4', + x=np.arange(1, n_iterations+1), + y=samples['Sigma log'].values[2], + mode='lines', + line_color=qualitative.Plotly[0], + showlegend=False + ), + row=4, + col=2 +) + +# Visualise prior distribution +parameter_values = np.linspace(4, 16, num=200) +pdf_values = np.exp([ + log_prior._priors[0]([parameter_value]) + for parameter_value in parameter_values]) +fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + ), + row=1, + col=1 +) + +parameter_values = np.linspace(0, 12, num=200) +pdf_values = np.exp([ + log_prior._priors[1]([parameter_value]) + for parameter_value in parameter_values]) +fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=2, + col=1 +) + +parameter_values = np.linspace(0, 12, num=200) +pdf_values = np.exp([ + log_prior._priors[2]([parameter_value]) + for parameter_value in parameter_values]) +fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=3, + col=1 +) + +parameter_values = np.linspace(0, 0.6, num=200) +pdf_values = np.exp([ + log_prior._priors[3]([parameter_value]) + for parameter_value in parameter_values]) +fig.add_trace( + go.Scatter( + name='Prior distribution', + legendgroup='Group 5', + x=parameter_values, + y=pdf_values, + mode='lines', + line_color='black', + showlegend=False + ), + row=4, + col=1 +) + +fig.update_layout( + xaxis_title='k_a', + yaxis_title='p', + yaxis2_title='k_a', + xaxis3_title='k_e', + yaxis3_title='p', + yaxis4_title='k_e', + xaxis5_title='v', + yaxis5_title='p', + yaxis6_title='v', + xaxis7_title='sigma', + yaxis7_title='p', + xaxis8_title='Number of iterations', + yaxis8_title='sigma', + template='plotly_white', + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1 + ), + margin=dict(t=0, r=0, l=0) +) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_3.html') + + +# 6 +fig = go.Figure() + +# Plot histograms +fig.add_trace( + go.Histogram( + name='Run 1', + x=samples['dose.absorption_rate'].values[0, ::n_iterations//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[2], + ) +) +fig.add_trace( + go.Histogram( + name='Run 2', + x=samples['dose.absorption_rate'].values[1, ::n_iterations//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[1], + ) +) +fig.add_trace( + go.Histogram( + name='Run 3', + x=samples['dose.absorption_rate'].values[2, ::n_iterations//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[0], + ) +) + +fig.update_layout( + template='plotly_white', + xaxis_title='Absorption rate in 1/day', + yaxis_title='Probability density', + barmode='overlay' +) +fig.update_traces(opacity=0.75) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_4.html') + + +# 7 +fig = go.Figure() + +# Plot histograms +fig.add_trace( + go.Histogram( + name='Run 1', + x=samples['dose.absorption_rate'].values[0, n_iterations//2::(n_iterations//2)//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[2], + ) +) +fig.add_trace( + go.Histogram( + name='Run 2', + x=samples['dose.absorption_rate'].values[1, n_iterations//2::n_iterations//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[1], + ) +) +fig.add_trace( + go.Histogram( + name='Run 3', + x=samples['dose.absorption_rate'].values[2, n_iterations//2::n_iterations//1000], + histnorm='probability density', + showlegend=True, + xbins=dict(size=0.5), + marker_color=qualitative.Plotly[0], + ) +) + +fig.update_layout( + template='plotly_white', + xaxis_title='Absorption rate in 1/day', + yaxis_title='Probability density', + barmode='overlay' +) +fig.update_traces(opacity=0.75) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_5.html') + + +# 8 +import arviz as az + + +# Summary for all samples +summary1 = az.summary(samples) +print(summary1) + +# Summary for the samples post warmup +summary2 = az.summary(samples.sel(draw=slice(n_iterations//2, n_iterations))) +print(summary2) + +directory = os.path.dirname(__file__) +summary1.to_csv(directory + '/3_fitting_models_summary_1.csv') +summary2.to_csv(directory + '/3_fitting_models_summary_2.csv') + + +# 9 +# Fix model parameters +model = chi.ReducedMechanisticModel(mechanistic_model=model) +model.fix_parameters(name_value_dict={ + 'dose.drug_amount': 0, + 'global.drug_amount': 0, +}) + +# Set dosing regimen +dosing_regimen = problem.get_dosing_regimens()['1'] +model.set_dosing_regimen(dosing_regimen) + +# Extract model parameters from summary +parameters = [ + summary2['mean'].loc[parameter] for parameter in model.parameters() +] + +# Simulate result +times = np.linspace(0.2, 3, 200) +simulation = model.simulate(parameters=parameters, times=times)[0] + +# Plot results +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + name='Fit' +)) +fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', + name='Measurements' +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_6.html') + + +# 10 +# Define posterior predictive distribution +predictive_model = problem.get_predictive_model() +samples = samples.sel(draw=slice(n_iterations//2, n_iterations)) +posterior_model = chi.PosteriorPredictiveModel( + predictive_model=predictive_model, posterior_samples=samples) + +# Approximate distribution using sampling +n_samples = 1000 +conc_samples = posterior_model.sample( + times=times, n_samples=n_samples, seed=1) + +# Reshape samples, so we can calculate mean and percentiles at the different +# time points +reshaped_samples = np.empty(shape=(n_samples, len(times))) +for sample_idx, sample_id in enumerate(conc_samples.ID.unique()): + reshaped_samples[ + sample_idx] = conc_samples[conc_samples.ID == sample_id].Value.values + +# Calculate mean, 5th and 95th percentile of the distribution at each time +# point +means = np.mean(reshaped_samples, axis=0) +lower = np.percentile(reshaped_samples, q=5, axis=0) +upper = np.percentile(reshaped_samples, q=95, axis=0) + +# Plot results +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=simulation, + mode='lines', + name='Fit: mean parameters' +)) +fig.add_trace(go.Scatter( + x=data.Time, + y=data.Value, + mode='markers', + name='Measurements' +)) +fig.add_trace(go.Scatter( + x=times, + y=means, + mode='lines', + name='Fit: mean posterior pred. distr.', + line=dict(color='black') +)) +fig.add_trace(go.Scatter( + x=times, + y=lower, + mode='lines', + name='Fit: 5th-95th perc. posterior pred. distr.', + line=dict(color='black', dash='dash') +)) +fig.add_trace(go.Scatter( + x=times, + y=upper, + mode='lines', + name='Fit: 5th-95th perc. posterior pred. distr.', + line=dict(color='black', dash='dash'), + showlegend=False +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white', + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1 + ), +) +fig.show() + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/3_fitting_models_7.html') \ No newline at end of file diff --git a/docs/source/getting_started/code/3_fitting_models_summary_1.csv b/docs/source/getting_started/code/3_fitting_models_summary_1.csv new file mode 100644 index 00000000..56bfea4e --- /dev/null +++ b/docs/source/getting_started/code/3_fitting_models_summary_1.csv @@ -0,0 +1,5 @@ +,mean,sd,hdi_3%,hdi_97%,mcse_mean,mcse_sd,ess_bulk,ess_tail,r_hat +dose.absorption_rate,8.898,3.516,0.399,13.007,1.503,1.125,10.0,11.0,1.21 +global.elimination_rate,1.525,1.959,0.207,6.662,0.938,0.714,9.0,11.0,1.25 +global.volume,6.018,2.529,0.47,9.233,1.05,0.783,8.0,11.0,1.26 +Sigma log,0.185,0.072,0.091,0.295,0.004,0.003,502.0,264.0,1.01 diff --git a/docs/source/getting_started/code/3_fitting_models_summary_2.csv b/docs/source/getting_started/code/3_fitting_models_summary_2.csv new file mode 100644 index 00000000..c2f39a47 --- /dev/null +++ b/docs/source/getting_started/code/3_fitting_models_summary_2.csv @@ -0,0 +1,5 @@ +,mean,sd,hdi_3%,hdi_97%,mcse_mean,mcse_sd,ess_bulk,ess_tail,r_hat +dose.absorption_rate,10.014,1.957,6.456,13.778,0.046,0.033,1802.0,3110.0,1.0 +global.elimination_rate,0.795,0.261,0.351,1.301,0.006,0.005,1850.0,2121.0,1.0 +global.volume,6.862,1.605,3.923,9.947,0.038,0.027,1757.0,2295.0,1.0 +Sigma log,0.178,0.053,0.095,0.278,0.001,0.001,1651.0,2381.0,1.0 diff --git a/docs/source/getting_started/code/dataset_1.csv b/docs/source/getting_started/code/dataset_1.csv new file mode 100644 index 00000000..fb59d2f3 --- /dev/null +++ b/docs/source/getting_started/code/dataset_1.csv @@ -0,0 +1,10 @@ +ID,Time,Time unit,Observable,Value,Observable unit,Duration,Dose,Dose unit +1,0.5,Day,Drug concentration,0.198,ng/mL,,, +1,1.0,Day,Drug concentration,0.123,ng/mL,,, +1,1.5,Day,Drug concentration,0.305,ng/mL,,, +1,2.0,Day,Drug concentration,0.184,ng/mL,,, +1,2.5,Day,Drug concentration,0.421,ng/mL,,, +1,3.0,Day,Drug concentration,0.306,ng/mL,,, +1,0.0,Day,,,,0.01,2.0,mg +1,1.0,Day,,,,0.01,2.0,mg +1,2.0,Day,,,,0.01,2.0,mg diff --git a/docs/source/getting_started/code/one_compartment_pk_model.xml b/docs/source/getting_started/code/one_compartment_pk_model.xml new file mode 100644 index 00000000..e36f1a18 --- /dev/null +++ b/docs/source/getting_started/code/one_compartment_pk_model.xml @@ -0,0 +1,36 @@ + + + + + + + + + + + + + + + + + -1 + elimination_rate + drug_amount + + + + + + + + + drug_amount + volume + + + + + + + diff --git a/docs/source/getting_started/code/pk_one_comp.xml b/docs/source/getting_started/code/pk_one_comp.xml new file mode 100644 index 00000000..097d194b --- /dev/null +++ b/docs/source/getting_started/code/pk_one_comp.xml @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + central + elimination_rate + drug + + + + + + + + diff --git a/docs/source/getting_started/code/template.xml b/docs/source/getting_started/code/template.xml new file mode 100644 index 00000000..a73b3c87 --- /dev/null +++ b/docs/source/getting_started/code/template.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/docs/source/getting_started/code/two_compartment_pk_model.xml b/docs/source/getting_started/code/two_compartment_pk_model.xml new file mode 100644 index 00000000..05d36cdc --- /dev/null +++ b/docs/source/getting_started/code/two_compartment_pk_model.xml @@ -0,0 +1,97 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + central + elimination_rate + drug + + + + + + + + + + + + + + + + + central + transition_rate_c2p + drug + + + + + + + + + + + + + + + + + peripheral + transition_rate_p2c + drug_peripheral + + + + + + + + + diff --git a/docs/source/getting_started/error_model.rst b/docs/source/getting_started/error_model.rst deleted file mode 100644 index 6e669c54..00000000 --- a/docs/source/getting_started/error_model.rst +++ /dev/null @@ -1,17 +0,0 @@ -.. currentmodule:: chi - -*************** -The error model -*************** - -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation - -.. autosummary:: - - chi.ErrorModel - chi.GaussianErrorModel - chi.LogNormalErrorModel - chi.MultiplicativeGaussianErrorModel - chi.ConstantAndMultiplicativeGaussianErrorModel - chi.ReducedErrorModel \ No newline at end of file diff --git a/docs/source/getting_started/fitting_models_to_data.rst b/docs/source/getting_started/fitting_models_to_data.rst new file mode 100644 index 00000000..74b6df4d --- /dev/null +++ b/docs/source/getting_started/fitting_models_to_data.rst @@ -0,0 +1,680 @@ +.. currentmodule:: chi + +.. _Dataset_1 : https://github.com/DavAug/chi/blob/main/docs/source/getting_started/data/dataset_1.csv +.. _Pints: https://github.com/pints-team/pints +.. _ArViz: https://python.arviz.org/en/stable/ +.. _Issue: https://github.com/DavAug/chi/issues + +********************************** +Fitting mechanistic models to data +********************************** + +In the previous tutorial, :doc:`mechanistic_model`, we have seen how we can +implement and simulate treatment response models in Chi. For example, using the same +1-compartment PK model from before, ``one_compartment_pk_model.xml``, we can simulate +the time course of drug concentration levels following repeated dose +adminstrations + +.. literalinclude:: code/3_fitting_models_1.py + :lines: 112-146 + +.. raw:: html + :file: images/3_fitting_models_1.html + +This ability to simulate treatment response is pretty cool in its own right, +but, at this point, our model has little to do with +real treatment responses. If our goal is to describe *real* treatment +responses, we need to somehow connect our model to reality. + +The most common approach to relate models to real treatment responses is to +compare the model predictions to measurements. Below, we have prepared an example +dataset with drug concentration measurements. These drug concentrations were +measured after repeatedly adminstering 2 mg of a drug every 24 hours. +You can download the dataset from the following link, if you would like to +follow this tutorial: Dataset_1_. + +.. csv-table:: Drug concentration measurements + :file: code/dataset_1.csv + :widths: 4, 12, 12, 12, 12, 12, 12, 12, 12 + :header-rows: 1 + +The dataset contains one column that identifies the measured individual (``ID``), +two columns that specify the time of a measurement or a dose administration +(``Time`` and ``Time unit``), three columns that specify the measured values +(``Observable``, ``Value``, ``Observable unit``), and three columns that specify +the dose administrations (``Dose``, ``Duration``, ``Dose unit``). + +If we download the file and save it in the same directory as the Python script, +we can visualise the measurements by executing the below script + +.. literalinclude:: code/3_fitting_models_1.py + :lines: 176-195 + +.. raw:: html + :file: images/3_fitting_models_2.html + +The figure shows that the treatment response dynamics indicated by the measurements +is similar to the treatment response simulated by the one-compartment PK model above. +But looking more closely at the magnitude of the values, +it appears that the measured values are much smaller +than the simulated ones. We can therefore conclude that, at this point, our +model does not provide an accurate description of the measured treatment +response. + +To find a better description of the treatment response, we have two options: +1. we can try to find parameters values that make the model output describe the +measurements more closely; or 2. we can define a new mechanistic model and see +whether this new model is able to describe the measurements better. This tutorial +will be about how we can find better model parameters. + +Estimating model parameters from data: Background +************************************************* + +Before we can try to find parameter values that most closely describe the treatment response +measurements, we first need to agree on what we mean by +"*most closely*" for the relationship between the mechanistic model output and the measurements. +The most naive way to define this notion of closeness is to use the distance +between the measurements and the model output, +i.e. the difference between the measured values and the +simulated values. If we used +the distance to define closeness, the model parameters that would most closely +describe the measurements would be those parameter values that make the mechanistic +model output perfectly match the measurements, leading to a distances of 0 ng/mL +at all measured time points. However, as outlined in Sections 1.3 and 1.4 of the +:doc:`quick_overview`, measurements of treatment responses are imprecise and +noisy, and will therefore not perfectly represent the treatment response dynamics. +Consequently, if we were to match the model outputs to measurements, +we would actually end up with an inaccurate description of the treatment response +because we would be paying too much attention to the measurement noise. + +One way to improve our notion of closeness is to incorporate the measurement +process into our computational model of the treatment response, thereby +explicitly stating that we do not expect the mechanistic model output to match +the measurements perfectly. In Chi, this can be done +using :class:`chi.ErrorModel` s. Error models promote the single value output +of mechanistic model simulations to a distribution of +values. This distribution characterises a +range of values around the mechanistic model output where measurements may be +expected. +For simulation, this distribution can be used to sample measurement values and +imitate the measurement process of real treatment responses, see +Section 1.3 in the :doc:`quick_overview` for an example. For parameter estimation, +the distribution can be used to quantify the likelihood with which the observed +measurements would have been produced by our model of the measurement process, +see Section 1.4 in the :doc:`quick_overview`. + +To account for measurement noise during the parameter estimation, we therefore +choose to quantify the closeness between the model output an the measurements +using likelihoods. Formally, the measurement process can be described by +distributions of measurement values of the form :math:`p(y | \psi, t, r)`, where :math:`y` +denotes the measurement value, :math:`\psi` denotes the model parameters, +:math:`t` denotes the time, and :math:`r` denotes the dosing +regimen. The likelihood of a single measurement :math:`((t_1, y_1), r)` is then +simply given by the value of the probability density evaluated at the measurement, +:math:`p(y_1 | \psi, t_1, r)`. This value depends on the chosen set of model +parameters, :math:`\psi`. The model parameters with the maximum likelihood are +the parameter values that most closely describe the measurements. + +.. note:: + The measurement distribution, :math:`p(y | \psi, t, r)`, is uniquely defined + by the mechanistic model output and the error model. For example for + the 1-compartment model, we may denote the simulated drug concentration + values by :math:`c(\psi, t, r)`, where the drug concentration values, :math:`c`, are + a function of the model parameters, :math:`\psi = (a_0, k_a, k_e, v)`, the time, + :math:`t`, and the dosing regimen, :math:`r`. + + 1. If we choose a :class:`chi.GaussianErrorModel` to describe the difference + between the model output and the measurements, we assume that measurements + are distributed according to a Normal distribution around the model output + + .. math:: + p(y | \psi, t, r) = \frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{e}^{-\big(y - c(\psi, t, r)\big)^2 / 2\sigma ^2}, + + where :math:`\sigma` defines the width of the distribution. For ease of notation, + we extend the definition of the model parameters to include :math:`\sigma`, + :math:`\psi = (a_0, k_a, k_e, v, \sigma)`. + + We can see that the model output + defines the mean or Expectation Value of the measurement distribution. + + 2. If we choose a :class:`chi.LogNormalErrorModel` to describe the difference + between the model output and the measurements, we assume that measurements + are distributed according to a lognormal distribution around the model output + + .. math:: + p(y | \psi, t, r) = \frac{1}{\sqrt{2\pi \sigma ^2}}\frac{1}{y}\mathrm{e}^{-\big(\log y - \log c(\psi, t, r) + \sigma / 2\big)^2 / 2\sigma ^2}. + + One can show that also for this distribution the model output defines the mean + or Expectation Value of the measurement distribution. + + The main difference between the two distributions is the shape. The + :class:`chi.GaussianErrorModel` is symmetrically distributed around the + model output, while :class:`chi.LogNormalErrorModel` is scewed in such a + way that measurements can never assume negative values. To visualise these + differences, we recommend simulating many measurements with different + error models, similar to Section 1.3 in :doc:`quick_overview`. But instead + of choosing different times, sample all measurements at the same time. You + can then histogram the samples, using for example ``go.Histogram``, as used + in Section 1.4.2 in :doc:`quick_overview`, to visualise the shape of + the probability density. + + +Assuming the independence of measurements, the likelihood for a dataset with :math:`n` measurements, +:math:`\mathcal{D}=((t_1, y_1), (t_2, y_2), \ldots, (t_n, y_n), r)`, is given +by the product of the individual likelihoods + +.. math:: + p(\mathcal{D}| \psi ) = \prod _{j=1}^n p(y_n | \psi, t_n, r), + +where *independence* refers to the assumption that measurements are +independently and identically distributed according to our model of the +measurement process (this does not have to be the case, and is especially unlikely +to be the case when our model fails to describe the measurement process accurately). + +While it is easy to numerically maximise the likelihood and +to derive maximum likelihood estimates in Chi, see Section 1.4.1 in the :doc:`quick_overview`, +the main objective of Chi is to support the Bayesian inference of model parameters. +Bayesian inference is conceptually different from maximum likelihood estimation +because it does not seek to find a single set of model parameters that "best" +describe the observed measurements. Instead, Bayesian inference acknowledges +the fact that noisy measurements cannot uniquely identify +the "best" model parameters, and therefore instead focuses on deriving +a distribution of parameter values which +are all consistent with the observed measurements, see Section 1.4.2 in :doc:`quick_overview` +for a more detailed discussion. This +distribution is derived from the likelihood using Bayes' rule + +.. math:: + p(\psi| \mathcal{D} ) = \frac{p(\mathcal{D}| \psi )\, p(\psi)}{p(\mathcal{D} )}, + +where :math:`p(\psi)` denotes the prior distribution of the model parameters. +The prior distribution can be used to quantify our belief of likely parameter +values for the model prior to the parameter estimation. + +Defining the log-posterior +************************** + +In Chi, you can estimate model parameters from measurements by inferring posterior +distributions of parameter values using Markov chain +Monte Carlo (MCMC) algorithms. In this tutorial, we will use MCMC algorithms +implemented in the open source Python package +Pints_, but in principle you can also use implementations from other libraries. +In Sections 1.4.1 and 1.4.2 in the :doc:`quick_overview`, +we showed in some detail how we can define (log-)posterior distributions, +:class:`chi.LogPosterior`, in Chi for this purpose. Here, we want to show +how we can use the :class:`chi.ProblemModellingController` to more easily +construct the :class:`chi.LogPosterior`, provided the measurements are available +in a CSV format of the form of Dataset_1_. + +The tricky bit when implementing log-posteriors for treatment response models +is often that those log-posteriors do not only depend on the treatment +response measurements, :math:`((t_1, y_1), (t_2, y_2), \ldots, (t_n, y_n))`, +but that they also depend on the administered dosing regimen, :math:`r`. +This can make it tedious to define log-posteriors, +especially when you are inferring parameters across measurements of multiple +individuals with difference dosing regimens. + +To simplify the process of constructing a :class:`LogPosterior`, we have +implemented the :class:`chi.ProblemModellingController`. +The :class:`chi.ProblemModellingController` facilitates the construction of +log-posteriors and reduces the workflow to a simple 4-step approach: + +- 1. Definition of the mechanistic model +- 2. Definition of the error model +- 3. Definition of the measurements +- 4. Definition of the prior distribution + +In the below code block, we illustrate this workflow for the above drug +concentration dataset, Dataset_1_. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 651-688 + +The first four blocks in the code define the individual components of the +log-posterior: the mechanistic model, the error model, the data, and the prior +distribution. Note that the administration of the dosing regimen is set +before passing the mechanistic model to the :class:`ProblemModellingController`. + +The prior distribution defines marginal distributions for the parameters, and +is implemented using Pints_. In Bayesian inference, we can use the prior +distribution to bias the inference +results towards feasible areas of the parameter space. In this +case, we have little prior knowledge about feasible parameter ranges, so we +choose relatively uninformative prior +distributions (below will be an illustration of the prior distribution). +Note that we seem to be specifying marginal prior distributions only for +4 of the 6 model parameters. This is because we fix the values of the initial drug +amounts to ``0`` prior to the inference +in the lines below, reflecting our knowledge that the subject had no prior +exposure to the drug before starting the trial. This reduces the number of model +parameters from 6 to 4. The fixing of model parmaters +is optional, but can sometimes improve the inference results when some model +parameters are already well understood. + +For the remaining 4 parameters, only positive values make biological sense, so +we choose prior distributions that focus on positive values. For +two model parameters, the volume of distribution and the scale parameter, +negative or zero values are particularly bad as they will break the simulation +(a volume of zero causes a division by zero error, and the lognormal distribution +is only defined for positive sigma). We therefore use ``pints.LogNormalLogPrior`` +to constrain those parameters to strictly positive values. + +In the final block of the code, we define the log-posterior. In the first line, +we specify the mechanistic model and the error model. In the next line, we set +the dataset. Note that we need to use the ``output_observable_dict`` to map the +output variable of the model, ``global.drug_concentration``, to the Observable name +in the dataset, ``Drug concentration``. Other specifications are not required, and +dosing regimens are automatically set, when the dosing regimen related columns, +``Dose`` and ``Duration``, are present in the dataset. In the following line, we +fix the initial drug amounts to ``0`` using :meth:`ProblemModellingController.fix_parameters`. +You can use this method to fix any parameters of the model. +In the last two lines, we set the log-prior and implement the log-posterior using +the :class:`ProblemModellingController.get_log_posterior` method. + +Inferring the posterior distribution +************************************ + +With this :class:`chi.LogPosterior` in place, we can infer the posterior +distribution using any MCMC algorithm of our choice. Recall that MCMC +algorithms infer distributions by sampling from them. +If we sample sufficiently many samples, +the histogram over the samples converges to the posterior +distribution. +To illustrate this, we run an MCMC algorithm to infer the above +defined posterior distribution of the 1-compartment PK +model. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 692-698 + +In the code block, we use an MCMC algorithm implemented in Pints_, called +``pints.HaarioBardenetACMC``. For technical reasons that we will discuss below, +we run the algorithm three times for 20,000 iterations. +Note that we use the :class:`chi.SamplingController` to set the number of +runs, the number of iterations, and to run the sampling algorithm. The +:class:`chi.SamplingController` can also do other things, such as running the +chains in parallel, but we will not go into this additional functionality in this +tutorial, and refer instead to the API reference. + +Executing the above code block will spawn a response of the below form. +The left most column indicates the current +iteration of the MCMC algorithm. The other columns show the total number of +log-posterior evaluations and the fraction of accepted MCMC proposals of each +chain. The target acceptance ratio depends on the details of the MCMC algorithm. +For the ``pints.HaarioBardenetACMC`` the target is ``0.23``. + +.. code-block:: bash + + Using Haario-Bardenet adaptive covariance MCMC + Generating 3 chains. + Running in sequential mode. + Iter. Eval. Accept. Accept. Accept. Time m:s + 0 3 0 0 0 0:00.0 + 1 6 0 0.5 0.5 0:00.0 + 2 9 0.333 0.667 0.333 0:00.0 + 3 12 0.25 0.75 0.5 0:00.0 + 20 63 0.714 0.571 0.476 0:00.0 + 40 123 0.756 0.61 0.561 0:00.0 + 60 183 0.738 0.475 0.59 0:00.0 + 80 243 0.716 0.407 0.642 0:00.1 + 100 303 0.693 0.406 0.653 0:00.1 + 120 363 0.736 0.421 0.686 0:00.1 + . + . + . + +When the three runs of the algorithm terminate, the generated samples are returned in +form of an ``xarray.Dataset``. +We can visualise the samples using the code +block documented at the end of this section (we move the code block to the end, +to avoid disrupting the flow of the tutorial with less relevant code snippets). + +.. raw:: html + :file: images/3_fitting_models_3.html + +The left column of the +figure shows the histogram over the samples across the three runs +in orange, as well as the probability density of the prior distribution in black. +The first row shows the +samples of the absorption rate, the second row shows the samples of the elimination +rate, the third row shows the samples of the volume of distribution, and the +fourth row shows the samples of the scale parameter of the error model, sigma. +The right column of the figure shows the +samples drawn at each iteration of the +MCMC algorithm runs. The samples from the different runs are illustrated in different +colours: run 1 (green); run 2 (red); run 3 (blue). + +The posterior distribution +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The orange distribution is the result of the inference -- the posterior distribution. +It contains all parameter values that are consistent with the drug concentration +measurements and our prior knowledge, assigning each set of parameter values +with a probability of being the data-generating parameter values. Notably, the +figure shows that our prior knowledge and Dataset_1_ are insufficient to +conclude on a single set of parameter values (see orange distribution). +Instead, the measurements only +allow us to refine our understanding of feasible parameter values. For example, +we can see in the second row of the figure that the marginal posterior distribution +substantially differs from the marginal prior distribution. This is because the +drug concentration measurements contain important information about the elimination rate, rendering +rates above 1.5 1/day or below 0.25 1/day as extremely unlikely for the +model of the treatment response. This in in stark contrast to the relatively wide +range of model parameters that we deemed feasible prior to the inference +(see black line). However, the measurements are not conclusive enough +to reduce the distribution of feasible parameters to a single value. Similarly, +for the volume of distribution (row 3) and the error scale parameter +(row 4), the measurements lead to substaintial updates relative to the +prior distribution. +In comparison, the measurements appear less informative about the absorption rate +(see row 1), given that the marginal posterior distribution of +the absorption rate is almost identical to its prior distribution. +We will have a closer look at an intuitive understanding of why the measurements +contain little information about +the absorption rate below. The take-away from this discussion is that inferring +distributions of parameter values consistent with the measurements is in many +ways more informative than just estimating a single set of model parameters, +as it can help to quantify the feasibility of different parameter values and the +overall remaining level of uncertainty. +In a Bayesian setting, inferring posterior distributions also +allows us to supplement the information contained in the measurements by prior +knowledge of feasible model parameters, facilitating the inference of parameter +values when treatment response measurements are limited. + +Assessing convergence +^^^^^^^^^^^^^^^^^^^^^ + +Similar to most numerical inference +algorithms, MCMC algorithms have practical +limitations, making a particular algorithm not equally well suited to all inference problems. +Some MCMC algorithms are, for example, only well suited to estimate model +parameters when the total number of parameters is small, while other MCMC algorithms only +work when the curvature of the posterior surface is not too extreme. While +understanding the limitations of each MCMC algorithm requires looking into the +algorithmic details, there is a simple way of testing the suitability of an +algorithm that does not require knowing about its details: the :math:`\hat{R}` metric. +We will motivate and provide an intuitive explanation of this test in this tutorial, +but will not go into the technical details. + +Let us begin this section by revisiting the right column in the figure above. The column +shows the samples from the three MCMC algorithm runs at each +iteration. For early iterations of the algorithm, +the samples from the MCMC runs look quite distinct -- each run appears to sample +from a different area of the parameter space. In contrast, for later iterations +the MCMC runs seem to converge and sample from the same area of the parameter space. Intuitively, +it does not really make sense for the samples from the MCMC runs to look different +-- after all, we use the same MCMC algorithm to sample from the same posterior distribution. +The histogram over the samples *should* therefore be identical within the limits of +sampling variation. However, if +we plot the samples from each of the three runs separately, +we find that the histograms actually look quite different. For illustrative purposes, +we focus on the absorption rate in the code block. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1025-1066 + +We select the samples of the absorption rate from the ``xarray.Dataset`` using +the name of the parameter in the mechanistic model, +``samples['dose.absorption_rate'].values``. The returned object is a numpy array +of shape ``(n_runs, n_iterations)``. We can select the samples from +run n using the index n-1, e.g. for run 1: ``samples['dose.absosption_rate'].values[0]``. +In addition, we subsample the MCMC samples in the above code block, which we can, +for now, regard as a cosmetic choice. + +.. raw:: html + :file: images/3_fitting_models_4.html + +We can see that the histograms over the samples from the three runs look very +different. This seems contradictory to the orange distribution which we presented +above as the histogram over the MCMC samples across the runs. But also, how can it +makes sense for the three histograms to differ when the +posterior distribution is the same? -- It does not, and in fact, it tells us that the three histograms in the figure above do not represent the posterior +distribution. Although, we did not disclose it until now, the orange +distribution is only based on the second half of each MCMC run! +All MCMC algorithms have a common limitation which we can see in the right column of the +figure presented earlier: they generate samples +from the posterior distribution using the latest sample as a starting point. +For some MCMC algorithms the internal sampling strategy is advanced enough to +sufficiently decorrelate the sample from this initial starting point, but for +many MCMC algorithms the initial starting point substantially influences the sampled value. That +means that if the last samples of two MCMC runs are far away from each other, the following +samples of the runs will also differ substantially from each other, see for example +the first 5000 iterations of the elimination rate in the second row of the figure +in the previous section. + +At the start of an MCMC run, there is no "last sampled" parameter value, +and we have to manually provide a starting point to the MCMC algorithm. +To reduce the number of manual steps in Chi, the :class:`chi.SamplingController` +automatically samples starting points from the prior distribution. +This means that the above runs of the MCMC algorithm start from three +different positions in parameter space. These starting points have little +to do with the posterior distribution, and are potentially far away from the area of +the parameter space that is interesting for the posterior distribution. +It therefore makes sense that during the early iterations of the runs, the sampled +values of the MCMC runs do not agree. Fortunately, MCMC algorithms are constructed +in such a way that their samples are, at least in theory, guaranteed to converge +to the posterior distribution, meaning that runs of MCMC algorithms starting from +different areas in parameter space can be expected to converge to the same area. +In practice, numerical limitations may nevertheless +make it impossible, or very time-consuming for samples from MCMC algorithms not +suited to the inference problem to +converge to the posterior distribution. Starting multiple runs from different +areas in parameter space therefore provides an effective way of testing whether +a given MCMC algorithm is suited for the inference of the problem at hand. If +the runs do not converge, the MCMC algorithm either needs more iterations to converge, or +is not suited for the inference problem, in which case the MCMC samples cannot +be expected to represent the posterior distribution. If the runs do converge, +we can be confident that the inference results are correct (although there is always +a chance for edge cases where this method may fail to identify the +limitations of the algorithm). + +So, starting multiple MCMC runs from different locations in parameter space is +a good idea to gain confidence in the inference results. We recommend running the +MCMC algorithm from 3-5 different starting points, randomly sampled from the prior distribution. +On the one hand, a large number of starting points will test the suitability +of the MCMC algorithm more thoroughly, but each run also comes at a computational +cost. 3-5 runs therefore provide a good tradeoff. In addition, starting points randomly sampled +from the prior usually make sure that the starting points are not too close together, +so that the suitability of the MCMC algorithm is indeed tested. At the same time, it also +ensures that the starting points are not extremely far away from each other, located in areas +of the parameter space that we do not even deem relevant for the inference result. +So, even if the MCMC algorithm would encounter problems in these extreme areas of parameter space, +it may not be very relevant for our inference problem. Randomly sampling from the +prior distribution therefore provides a good tradeoff between well dispersed +and not too extreme starting points. + +In conclusion, the initial iterations of the MCMC runs are very usefull to establish the +validity of the inference results. At the same time, they have little to do +with the posterior distribution, as outlined above. +It is therefore common practice to discard these early samples +as "warm-up" of the sampling algorithm. To this end, the earlier presented +"trace plot", visualising the samples of the different runs at each iteration +is extremely useful. We can see, for example, that somewhere around the 8000s iteration the +three runs of the MCMC algorithm converge. For this use case, we therefore choose the 10,000s iteration +as the cut-off for the warm-up (2000 extra iteration for good measure). +Below, we plot the histogram over the samples +from the three chains, this time using only the samples from the second half of +each run. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1073-1114 + +.. raw:: html + :file: images/3_fitting_models_5.html + +We can see that the histograms over the samples are now in much better agreement! + +This visual assessment of the convergence already gets us very far, but for those +that would prefer more quantitative metrics to assess the convergence of MCMC +runs, we recommend the open source library ArViz_. ArViz_ implements a number +of widely established metrics to assess the convergence of MCMC algorithms, +including the :math:`\hat{R}` metric and the effective sample size. We can estimate +both values (and many more) using the ``az.summary`` function. In the below code +block, we first estimate the values for all MCMC samples, and then estimate +the values only for the samples post warm-up. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1121-1130 + +The return of ``az.summary`` is a ``pandas.DataFrame``, containing the estimated +values, as illustrated below + +.. csv-table:: Summary: all samples + :file: code/3_fitting_models_summary_1.csv + :widths: 19, 9, 9, 9, 9, 9, 9, 9, 9, 9 + :header-rows: 1 + +.. csv-table:: Summary: samples post warm-up + :file: code/3_fitting_models_summary_2.csv + :widths: 19, 9, 9, 9, 9, 9, 9, 9, 9, 9 + :header-rows: 1 + +The summary output contains the estimates of the mean value, the standard deviation, +and many other values. For the assessment of the convergence, we would like to focus +on the :math:`\hat{R}` values in the right-most column of the dataframes, the ``r_hat`` columns. +We can see that the summary provides one estimate for each parameter. For the +samples post warm-up, the estimates are all ``1.00``, while the estimates for all +samples assume larger values. Loosely speaking, the :math:`\hat{R}` metric compares +the variance of samples across runs to the variance of samples within runs. If +samples from different runs are very different, the variance across runs is much larger +than the variance of samples within the same run, leading to an estimate :math:`>1`. +On the flip side, if the inter-run variance is the same as the intra-run variance, +the samples across runs look the same and the :math:`\hat{R}` metric returns a +value of :math:`1`. As a result, :math:`\hat{R}` estimates of :math:`1` indicate +convergence, while estimate :math:`>1` indicate that the runs have not yet fully converged. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 701-1018 + +Assessing convergence: Summary +****************************** + +1. Run the MCMC algorithm multiple times from different starting points. The :class:`chi.SamplingController` automatically samples starting points from the prior distribution. Recommended number of runs is 3-5. +2. Visually assess the convergence of the runs to the same area of parameter space using trace plots i.e. sample values on the y-axis and iteration at which each sample was obtained on the x-axis. +3. Discard the samples where the runs have not converged as "warm-up". If the runs do not converge, the MCMC algorithm has to run for more iterations until they converge, or the algorithm may not be suited for the inference problem. +4. Calculate the :math:`\hat{R}` value using the remaining samples. We recommend the ``az.rhat`` function implemented in ArViz_. Values :math:`<1.01` indicate good convergence. Larger values may indicate problems with the inference. +5. (Optional) For MCMC runs that require a large number of iterations for convergence, it can make sense to subsample the samples. Often, 1000-3000 samples are sufficient to represent a distribution. For ``n_iterations`` post warm-up, keeping every ``n_iterations // 1000`` th sample may therefore be sufficient. + +For more details, please refer to the previous section. + +Comparing model fits to data +**************************** + +Let us conclude this tutorial by comparing our model fit to the drug concentration +measurements. The simplest way to do this is to just focus on the means of the +posterior distributions, which we can extract from the summary dataframes presented +in Section 3.3.2. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1138-1177 + +.. raw:: html + :file: images/3_fitting_models_6.html + +In the above code block, we first fix the initial amount parameters of the mechanistic +model using the :class:`chi.ReducedMechanisticModel` class. This class is a wrapper +for mechanistic models that implements the :meth:`ReducedMechanisticModel.fix_parameters` +method. This method works analogously to the :meth:`ProblemModellingController.fix_parameters` +and fixes model parameters to constant values. In the next lines, we extract the +dosing regimen from the :class:`ProblemModellingController`. The controller formats +the dosing regimens from the dataset, Dataset_1_, into a useful format for the :class:`chi.PKPDModel`, returning the +dosing regimens in form of a dictionary, using the +IDs of the individuals as keys. This makes it possible to access the dosing regimen +of the individual with ``'ID 1'`` in Dataset_1_ using the key ``'1'``. +In the remaining lines leading up to the plotting, we extract the +estimates of the parameter means from the summary dataframe. + +We can see that our fit describes the measurements reasonably well. However, it is +a little bit unsatisfying that the uncertainty of the parameter estimates is not +reflected in the model fit at all. To resolve this shortcoming of just using the +parameter means to plot the model fit, Chi implements the :class:`chi.PosteriorPredictiveModel`, +which defines the posterior predictive distribution of the measurement +process + +.. math:: + + p(y | \mathcal{D}, r, t) = \int \mathrm{d} \psi \, p(y | \psi, r, t)\, p(\psi | \mathcal{D}). + +The posterior predictive distribution averages our model of the measurement +distribution over the parameter values that we found to be consistent with the data. +Before discussing this in more detail, let us take a look at how the posterior predictive +distribution compares to the model fit where we just used the means of the parameters. +To this end, we sample measurements repeatedly from :math:`p(y | \mathcal{D}, r, t)` +and estimate the mean and some of its percentiles at each time point. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1184-1206 + +.. raw:: html + :file: images/3_fitting_models_7.html + +In the first few lines of the code block, we implement the +:class:`chi.PosteriorPredictiveModel` by first retrieving the +:class:`chi.PredictiveModel` from the problem modelling controller, and then +combining it with the inferred posterior samples. +The :class:`chi.PredictiveModel` defines the measurement distribution, +:math:`p(y | \psi, r, t)`, for the purpose of simulating measurements. Note that +:math:`p(y | \psi, r, t)` is the same distribution as the one we used to construct the +likelihood in Section 3.1. +The important difference is that for the likelihood the values of :math:`(y, t, r)` +are fixed to the measurements, while for the :class:`chi.PredictiveModel`, the times +and the dosing regimen can be chosen at will to simulate any measurements of interest. + +In the following 2 lines of the code block, we sample measurements from the posterior predictive distribution -- +1000 samples at each simulated time point. The samples are returned as a ``pandas.DataFrame``, +which we reshape to a numpy array of shape ``(n_samples, n_times)`` in the next few lines. +This reshaping makes it very easy to calculate the mean value and the percentiles +of the samples at each time point, as demonstrated in the final lines of the code block. + +The figure visualises the estimated mean and percentiles of the posterior predictive +distribution on top of our earlier fit. The code to reproduce the figure is documented below. +The area between the percentile estimates marks the region of drug concentration values +which can be expected to be occupied by 90% of all measurements, +assuming, of course, that our model correctly describes the measurement process. +The estimated percentiles are more jittery than the plot with the mean parameter values +because we use sampling to estimate the percentiles. As a result, the estimates +are subject to sampling error, which can be reduced by increasing the number of +samples, ``n_samples``. However, for the purpose of this analysis, the jittering +is not substantial enough to warrant a larger number of samples. + +Plotting the posterior predictive distribution provides two important insights +that the fit with the means of the parameters cannot provide: 1. it can quantify +the uncertainty associated with the treatment response. This uncertainty is partially +due to measurement noise, but it also is due to the residual uncertainty in the +model parameters, as quantified by the posterior distribution. 2. It is +able to better reflect nonlinearities between model parameters and the model output +which are otherwise neglected when plotting the fit with the means of the model parameters. +We can see an indication of this in the figure by +noting that the fit with the parameter means is not the same as the mean of the +posterior predictive distribution. This difference is a result of the, in general, +nonlinear dependence of treatment response models on their parameters, meaning that the average +treatment response under parametric uncertainty is **not** the same as the treatment response +predicted with the mean parameters. In fact, these two predictions can differ substantially from +each other, especially when predictions are made for future times or for previously +unexplored dosing regimens. + +.. literalinclude:: code/3_fitting_models_2.py + :lines: 1208-1256 + +This concludes this tutorial. If you have any feedback or suggestions for +improvement, or would like to report any typos, mistakes or bugs, please do reach out to us, for example +by creating an Issue_. We are looking +forward to hearing from you! + +Reference to ErrorModel, LogPDF and PredictiveModel API +******************************************************* + +.. autosummary:: + + chi.ErrorModel + chi.GaussianErrorModel + chi.LogNormalErrorModel + chi.MultiplicativeGaussianErrorModel + chi.ConstantAndMultiplicativeGaussianErrorModel + chi.ReducedErrorModel + chi.LogLikelihood + chi.LogPosterior + chi.ProblemModellingController + chi.SamplingController + chi.PredictiveModel + chi.PosteriorPredictiveModel \ No newline at end of file diff --git a/docs/source/getting_started/images/2_mechanistic_model_1.html b/docs/source/getting_started/images/2_mechanistic_model_1.html new file mode 100644 index 00000000..fd10839b --- /dev/null +++ b/docs/source/getting_started/images/2_mechanistic_model_1.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/2_mechanistic_model_2.html b/docs/source/getting_started/images/2_mechanistic_model_2.html new file mode 100644 index 00000000..9d2987a1 --- /dev/null +++ b/docs/source/getting_started/images/2_mechanistic_model_2.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/2_mechanistic_model_3.html b/docs/source/getting_started/images/2_mechanistic_model_3.html new file mode 100644 index 00000000..383b98fa --- /dev/null +++ b/docs/source/getting_started/images/2_mechanistic_model_3.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/2_mechanistic_model_4.html b/docs/source/getting_started/images/2_mechanistic_model_4.html new file mode 100644 index 00000000..d806a7e0 --- /dev/null +++ b/docs/source/getting_started/images/2_mechanistic_model_4.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/2_mechanistic_model_5.html b/docs/source/getting_started/images/2_mechanistic_model_5.html new file mode 100644 index 00000000..c954704a --- /dev/null +++ b/docs/source/getting_started/images/2_mechanistic_model_5.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_1.html b/docs/source/getting_started/images/3_fitting_models_1.html new file mode 100644 index 00000000..26a3fcc3 --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_1.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_2.html b/docs/source/getting_started/images/3_fitting_models_2.html new file mode 100644 index 00000000..cade72bc --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_2.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_3.html b/docs/source/getting_started/images/3_fitting_models_3.html new file mode 100644 index 00000000..0f7917e0 --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_3.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_4.html b/docs/source/getting_started/images/3_fitting_models_4.html new file mode 100644 index 00000000..d1490722 --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_4.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_5.html b/docs/source/getting_started/images/3_fitting_models_5.html new file mode 100644 index 00000000..f9a3e93e --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_5.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_6.html b/docs/source/getting_started/images/3_fitting_models_6.html new file mode 100644 index 00000000..15874a64 --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_6.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/3_fitting_models_7.html b/docs/source/getting_started/images/3_fitting_models_7.html new file mode 100644 index 00000000..dfc47233 --- /dev/null +++ b/docs/source/getting_started/images/3_fitting_models_7.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/index.rst b/docs/source/getting_started/index.rst index 305a5de0..f98716f0 100644 --- a/docs/source/getting_started/index.rst +++ b/docs/source/getting_started/index.rst @@ -1,27 +1,20 @@ -.. _README: https://github.com/DavAug/chi/blob/main/README.md - *************** Getting started *************** .. currentmodule:: chi -This part of the documentation gets you started using chi. Each section will +This part of the documentation gets you started using Chi. Each section will give a brief introduction into the modelling framework and show examples of -how to implement models in chi. The covered topics include **simulation** and +how to implement models in Chi. The covered topics include **simulation** and **inference** of e.g. ODE models, PKPD models with drug administration and -non-linear mixed effects models / hierarchical models. For installation -instructions please refer to the README_. +non-linear mixed effects models / hierarchical models. .. toctree:: :numbered: :caption: Table of Contents - :maxdepth: 1 + :maxdepth: 2 quick_overview mechanistic_model - error_model - log_likelihood - log_posterior - mcmc_sampling - population_model \ No newline at end of file + fitting_models_to_data \ No newline at end of file diff --git a/docs/source/getting_started/log_likelihood.rst b/docs/source/getting_started/log_likelihood.rst deleted file mode 100644 index 8e5810bb..00000000 --- a/docs/source/getting_started/log_likelihood.rst +++ /dev/null @@ -1,12 +0,0 @@ -.. currentmodule:: chi - -****************** -The log-likelihood -****************** - -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation - -.. autosummary:: - - chi.LogLikelihood \ No newline at end of file diff --git a/docs/source/getting_started/log_posterior.rst b/docs/source/getting_started/log_posterior.rst deleted file mode 100644 index 2cd7c813..00000000 --- a/docs/source/getting_started/log_posterior.rst +++ /dev/null @@ -1,12 +0,0 @@ -.. currentmodule:: chi - -***************** -The log-posterior -***************** - -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation - -.. autosummary:: - - chi.LogPosterior \ No newline at end of file diff --git a/docs/source/getting_started/mcmc_sampling.rst b/docs/source/getting_started/mcmc_sampling.rst deleted file mode 100644 index 7309281f..00000000 --- a/docs/source/getting_started/mcmc_sampling.rst +++ /dev/null @@ -1,12 +0,0 @@ -**************************************************** -Inferring posterior distributions with MCMC sampling -**************************************************** - -.. currentmodule:: chi - -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation - -.. autosummary:: - - chi.SamplingController \ No newline at end of file diff --git a/docs/source/getting_started/mechanistic_model.rst b/docs/source/getting_started/mechanistic_model.rst index 11f70c20..50a1da8d 100644 --- a/docs/source/getting_started/mechanistic_model.rst +++ b/docs/source/getting_started/mechanistic_model.rst @@ -1,11 +1,942 @@ +.. _SBML: https://sbml.org/ +.. _Myokit: http://myokit.org/ +.. _Issue: https://github.com/DavAug/chi/issues + .. currentmodule:: chi ********************* The mechanistic model ********************* -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation +In Chi, mechanistic model is an umbrella term for time +series models describing the dynamics of treatment responses. In general, these models do +not have to be of a *mechanistic* nature, but in most cases they will have at least +a mechanistic element to them. +Popular mechanistic models include for example PKPD models, PBPK models and +QSP models which are based on mechanistic descriptions of the treatment response. + +Mechanistic models can be implemented in Chi in two ways: 1. using SBML_ files +(recommended); and 2. using the :class:`chi.MechanisticModel` interface. +The :class:`chi.MechanisticModel` interface has, perhaps, the +lower entry barrier, as it just involves implementing models +using any preferred Python code. This brings an enourmous amount of flexibility +and, for simple models, gets your model implemented fast. However, it has the +drawback that sensitivities of the model parameters and +mechanisms to administer dosing regimens have to implemented manually. It also +limits the ability to share models. + +In comparison, using SBML files to implement your models requires a small +amount of potentially unfamiliar SBML syntax, and therefore may initially take +a little longer to get started. However, for complex models, SBML +simplifies the model implementation and reduces the risk for implementation +errors. In Chi, SBML will also automate differentiation, the evaluation of parameter +sensitivities, and the implementation of dose administrations. Another benefit +of SBML files is that they are programming language-agnostic, meaning that SBML +files are supported by many simulation softwares facilitating the sharing and +reimplementation of models without forcing the continued use of Chi. We therefore +recommend using SBML files to implement your models. + +Below we will show how we can use either of those implementation strategies to +implement a 1-compartment PK model. + +**Use case: 1-compartment PK model:** + +A 1-compartment PK model is a semi-mechanistic description of the absorption, +distribution, metabolism and elimination of a drug using a simple differential +equation of the form + +.. math:: + \frac{\mathrm{d}a}{\mathrm{d}t} = -k_e\, a, + \quad c = \frac{a}{v}, + \quad a(t=0) = a_0, + +where :math:`a` is the drug amount in the compartment, :math:`t` is the time, +:math:`k_e` is the elimination rate of the drug, :math:`c` is the drug +concentration, :math:`v` is the volume of the compartment and :math:`a_0` is +the initial drug amount. + +Defining mechanistic models using the MechanisticModel interface +**************************************************************** + +The simplest way to implement mechanistic models in Chi is to use the +:class:`chi.MechanisticModel` interface. The :class:`MechanisticModel` is a +base class for mechanistic models with a small number of mandatory methods that +you need to implement in order to be able to use your implementation with all +of Chi's functionality. The base class takes the following form +(methods that don't have to be implemented are omitted): + +.. code-block:: python + + class MechanisticModel(object): + def simulate(self, parameters, times): + """ + Returns the numerical solution of the model outputs for the + specified parameters and times. + + The model outputs are returned as a 2 dimensional NumPy array of shape + ``(n_outputs, n_times)``. + + :param parameters: An array-like object with values for the model + parameters. + :type parameters: list, numpy.ndarray + :param times: An array-like object with time points at which the output + values are returned. + :type times: list, numpy.ndarray + + :rtype: np.ndarray of shape (n_outputs, n_times) + """ + raise NotImplementedError + + def has_sensitivities(self): + """ + Returns a boolean indicating whether sensitivities have been enabled. + """ + raise NotImplementedError + + def n_outputs(self): + """ + Returns the number of output dimensions. + """ + raise NotImplementedError + + def n_parameters(self): + """ + Returns the number of parameters in the model. + """ + raise NotImplementedError + + def outputs(self): + """ + Returns the output names of the model. + """ + raise NotImplementedError + + def parameters(self): + """ + Returns the parameter names of the model. + """ + raise NotImplementedError + +For full details of the base class, we refer to the API reference: +:class:`chi.MechanisticModel`. + +The main method of the mechanistic model is the :meth:`MechanisticModel.simulate` +method, implementing the simulation of the model. The other methods are +used for book-keeping and for enabling the computation of sensitivities. In +this example, we will not implement the sensitivities, so we can start +the model implementation by setting the return of the +:meth:`MechanisticModel.has_sensitivities` method to ``False`` + +.. code-block:: python + + class OneCompPKModel(chi.MechanisticModel): + def __init__(self): + super().__init__() + + def has_sensitivities(self): + # We will not implement sensitivities in this part of the tutorial, so + # the output of this method is always False + return False + +In the above code block, we define a new class called ``OneCompPKModel`` which will +implement the 1-compartment PK model. In the first line, the class +inherits from the :class:`chi.MechanisticModel` base class. In the next two +lines, we use ``super().__init__()`` to initialise properties of the the base +class. If you are not familiar with inheritance in Python, you can use these +three lines as a default setup of your model implementation. The remainder of +the code block implments the :meth:`MechanisticModel.has_sensitivities` method +by setting its return +to ``False``, indicating to other classes and functions in Chi +that the model does not implement sensitivities. This is a necessary method of +:class:`chi.MechanisticModel` in Chi, because some optimisation and inference +algorithms require sensitivites, and setting the output to ``False`` will +indicate to those algorithms that they cannot be used with this model. + +The other methods, besides the :meth:`MechanisticModel.simulate`, are used for +book-keeping and define the number of outputs of the model, the number of +parameters of the model and their names. In this case the 1-compartment PK +model has only 1 output, the concentration of the drug in the central +compartment. We can therefore continue the implementation of the model by setting +the return of :meth:`MechanisticModel.n_outputs` to ``1`` and the return of +:meth:`MechanisticModel.outputs` to ``['Drug concentration']`` + +.. code-block:: python + + class OneCompPKModel(chi.MechanisticModel): + def __init__(self): + super().__init__() + + def has_sensitivities(self): + # We will not implement sensitivities in this part of the tutorial, so + # the output of this method is always False + return False + + def n_outputs(self): + return 1 + + def outputs(self): + return ['Drug concentration'] + +Those two methods communicate to other classes and functions in Chi that the +model has only one output, the drug concentration. Note that the return of +:meth:`MechanisticModel.outputs` is a list of strings in order to facilitate +returning the names of mutliple outputs, in which case the return would take +the form ``['Output name 1', 'Output name 2', ..., 'Output name n']``. + +The remaining book-keeping methods return the number of parameters and +the names of the parameters, :meth:`MechanisticModel.n_parameters` and +:meth:`MechanisticModel.parameters`. A 1-compartment PK model has three model +parameters: 1. the initial drug amount, :math:`a_0`; 2. the elimination rate, +:math:`k_e`: and 3. the volume of distribution, :math:`v`. We can, thus, +implement the methods by setting the return of +:meth:`MechanisticModel.n_parameters` to ``3`` and the return of +:meth:`MechanisticModel.parameters` to +``['Initial amount', 'Elimination rate', 'Volume of distribution']``. + +.. code-block:: python + + class OneCompPKModel(chi.MechanisticModel): + def __init__(self): + super().__init__() + + def has_sensitivities(self): + # We will not implement sensitivities in this part of the tutorial, so + # the output of this method is always False + return False + + def n_outputs(self): + return 1 + + def outputs(self): + return ['Drug concentration'] + + def n_parameters(self): + return 3 + + def parameters(self): + return ['Initial amount', 'Elimination rate', 'Volume of distribution'] + +This leaves only :meth:`MechanisticModel.simulate` to implement. The +:meth:`MechanisticModel.simulate` is the heart of the mechanistic model, taking +the parameter values and the evaluation time points of the simulation as inputs +and returning the simulated model outputs. How you find the simulated +model outputs is in the :meth:`MechanisticModel` interface up to you! + +In thise case, the 1-compartment PK model is simple +enough so that we can find and implement the analytical solution to the +differential equation directly. Solving the differential equation shows that +the 1-compartment PK model yields an exponential decay of the drug +concentration in the central compartment + +.. math:: + c(t) = \frac{a_0}{v}\, \mathrm{e}^{-k_e t}. + +We can implement this exponential decay using numpy as follows + +.. literalinclude:: code/2_mechanistic_model_1.py + :lines: 96-130 + +The inputs to the :meth:`MechanisticModel.simulate` method are the values of +the model parameters and the time points for the model evaluation. As defined +by the base class, the parameter values, ``parameters``, have to take the form +of a list or an np.ndarray of length ``n_parameters``. These parameter values +are deconstructed in the first line of the method and assigned to the model +parameters. The time point input, ``times``, also has to take the form of a list or +an np.ndarray, but can be of any length, ``n_times``. In the second line of the method, +we use ``np.array(times)`` to convert the ``times`` input to an np.ndarray to enable the +computation of the drug concentration at all time points in parallel using numpy's vectorised operations. +The next line implements the model solution and calculates the drug +concentration for all time points in ``times``. Finally, the last two lines of +the method make sure that the output of the method adheres to the format defined +by the interface, :class:`chi.MechanisticModel`, reshaping the simulation result +to an np.ndarray of shape ``(n_outputs, n_times)`` by first instantiating an +empty array of shape ``(self.n_outputs(), len(times))``, i.e. ``(1, len(times))``, +and then populating the array with the simulated concentration values. + +This completes the implementation of the 1-compartment PK model and we can start +to model the pharmacokinetics of drugs. For example, we can simulate the time course +of the drug concentration following a bolus dose of 10 drug amount units + +.. literalinclude:: code/2_mechanistic_model_1.py + :lines: 135-162 + +.. raw:: html + :file: images/2_mechanistic_model_1.html + +As expected, the drug concentration decays exponentially with the rate defined +by the elimination rate, starting from an initial value of 5, defined by the +initial drug amount and the volume of distribution (:math:`10 / 2`). + +Strengths & limitations +^^^^^^^^^^^^^^^^^^^^^^^ + +The above implementation of the 1-compartment PK model can now be used together +with other classes and functions in Chi to simulate drug concentrations and to +infer parameter values from measurements. However, to get access to all of +Chi's features additional methods need to be implemented, including methods +mirroring :class:`PKPDModel.set_administration` and +:class:`PKPDModel.set_dosing_regimen` in order to facilitate dose +administrations that go beyond bolus doses at the start of the simulation. +Below, we provide an overview of the strengths and limitations of using the +:class:`chi.MechanisticModel` interface + +**Strengths** + +- Simple models can be implemented quickly. +- The interface is extremely flexible and allows for the use of any preferred Python code. +- Only knowledge of Python is required to implement the model. + +**Limitations** + +- No support for the implementation of complex models. +- No automated implementation of output names and parameter names. +- No automated implementation of changing model outputs. +- No automated support of renaming model outputs and model parameters. +- No automated implementation of dose administrations. +- No automated implementation of sensitivities. +- Sharing of models is limited to the use of Chi. + + +Defining mechanistic models using SBML files +******************************************** + +The systems biology markup language (SBML) can be used to define +mechanistic models in a programming language-agnostic way, facilitating the +sharing, the implementation, and the reproduction of computational models. In +chi, SBML also allow us to automatically implement sensitivities and dose +administrations of treatment response models. But before we go into more +detail of why SBML is a great choice to implement your models, let us implement +our first SBML model to see how it works in practice. + +Setting up a template +^^^^^^^^^^^^^^^^^^^^^ + +The first thing to note about SBML is that mechanistic models are not defined +in your Python scripts directly, but instead are defined in SBML files external +to your code (this makes the programming language-agnostic sharing of models possible). +These SBML files have a minimal biolerplate of five lines + +.. code-block:: xml + + + + + + + + + + +Create a new file called ``template.xml`` and copy-paste the above lines in there. +For all future models you will implement, you can use these lines as a starting point. +The first line specifies the XML version and the encoding, while the second and last line +in the file specify the XML namespace. This namespace is what makes an XML file +an SBML file. The remaining two lines begin the model definition. + +You can see +that the model tag has an ``id`` property with the value ``"template"``. +This id is not really used by chi, but Myokit_ uses it internally to name models, +and we can use this in this tutorial for debugging / making sure that the model +is implemented as expected. To this end, let us instantiate +a model from the SBML file using the code below + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 261-270 + +In this example, we instantiate the model using :class:`chi.PKPDModel` from +the SBML file. The first two lines define the absolute path to the SBML file +by first getting the absolute path of the Python script and then pointing to +the SBML file in the same directory as the script. The last line prints the +name of the ``_model`` property of the :class:`chi.PKPDModel` instance, which +is the compiled Myokit_ model. If everything works correctly, executing this +script should print the name of the model to the terminal, i.e. ``template``. + +.. note:: + We do not recommend accessing the ``_model`` property in your scripts + directly when you are modelling treatment responses (as indicated by the + leading ``_``, the Myokit_ model is supposed to be a private property of + :class:`chi.PKPDModel`), but for debugging SBML files it can be useful to + investigate ``_model`` directly. + +Implementing the model +^^^^^^^^^^^^^^^^^^^^^^ + +With this template in place, we can now implement our model. Let us +start by first sketching out the two blocks of definitions that we need for +this model: 1. parameter/variable definitions; and 2. rule definitions (including +assignment rules like :math:`c = \frac{a}{v}` and rate rules like +:math:`\mathrm{d}a/\mathrm{d}t = -k_e \,a`) + +.. code-block:: xml + + + + + + + + + + + + + + + + +All parameters and variables of the model will need to be defined inside the +``listOfParameters`` tags, while all assignment rules and rate rules will have +to be defined inside the ``listOfRules`` tags. + +Parameters are defined using parameter tags of the form ````. +The ``id`` property defines the name of the parameter and the ``value`` property defines +the value of the parameter to which it is initialised. Parameter tags can have more, optional +properties (see SBML_), but the ``id`` and the ``value`` properties are required. +Using the parameter tags, we can define the parameters and variables of the one-compartment +PK model :math:`(a, c, k_e, v)`. + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + +Let us check that our SBML file is implemented correctly by instantiating the +model and printing the compiled model back to us. To this end, let us +copy-paste the above model definition into an XML file with the name +``one_compartment_pk_model.xml`` and execute the below script + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 274-283 + +Here, we are using the ``code()`` method of the Myokit_ model to print +the model specification back to us. If the implementation is correct, you +should see a print out like this + +.. code-block:: python + + [[model]] + name: one_compartment_pk_model + + [global] + drug_amount = 1 + drug_concentration = 1 + elimination_rate = 1 + time = 0 bind time + in [1] + volume = 1 + +The model definition shows that the model contains 5 parameters, the 4 parameters +that we have defined, initialised to the value as specified in the SBML file, +and a fifth variable called ``time``. The time variable is automatically added +to the model and, as the name suggest, is the variable that will keep track of +time. + +At this point, the model is pretty useless because we haven't implemented the +rate rule that defines how the drug amount changes over time, and we also haven't +defined how the drug concentration is related to the drug amount. So let us +implement the rate rule, :math:`\mathrm{d}a/\mathrm{d}t = -k_e\, a` to bring this +PK model to life + +.. code-block:: xml + + + + + + + + + + + + + + + + + + -1 + elimination_rate + drug_amount + + + + + + + + +Before explaining the syntax in more detail, let us first see how the addition of +the rate rule has changed our model implementation. Executing the above Python +script again, now yields a model print out that looks like this + +.. code-block:: python + + [[model]] + name: one_compartment_pk_model + # Initial values + global.drug_amount = 1 + + [global] + dot(drug_amount) = -1 * elimination_rate * drug_amount + drug_concentration = 1 + elimination_rate = 1 + time = 0 bind time + in [1] + volume = 1 + +We can see that in comparison to our previous model print out two things have happened: +1. a new section with title ``# Initial values`` has appeared; and 2. the expression +of the ``drug_amount`` variable has changed. This is because rate rules promote (constant) +parameters to variables whose dynamics are now defined by ordinary differential equations. +In the model print out the ``dot(drug_amount)`` denotes the time derivative of the drug amount +variable, i.e. :math:`\mathrm{d}a/\mathrm{d}t`. The initial value of the drug amount at :math:`t=0` +is defined by the ``value`` property of the corresponding parameter, indicated in the +model print out by ``global.drug_amount = 1``. + +The SBML syntax for defining rate rules like this is to add ``rateRule`` tags to +the ``listOfRules`` that point to the relevant parameter in the ``listOfParameters`` + +.. code-block:: xml + + + + + + + + + + + + + + + + + + +This promotes the parameter to a variable that can change over time. The right +hand side of the differential equation is defined inside the ``math`` tags which +point to the MathML XML namespace which is used in SBML to define mathematical +expression. + +In MathML, mathematical expressions are encapsulated by the ```` +tags followed by a tag that indicates the mathetical operation and the relevant +parameters / variables. The most common operations are + +- addition: ```` +- subtraction: ```` +- multiplication: ```` +- division: ```` +- exponentiation: ```` + +Commutative operations, such as addition and multiplication, work with any +number of parameters / variables while other operators can only be applied +pairwise. For more information about the MathML syntax, we refer to +http://www.w3.org/1998/Math/MathML. + +Revisiting the ``rateRule`` in our model definition + +.. code-block:: XML + + + + + + -1 + elimination_rate + drug_amount + + + + +we can now see that the right hand side of the differential equation for the +drug amount is defined by the product of the constant ``-1``, +the ``elimination_rate`` parameter, and the ``drug_amount`` variable. Note that +constants are labelled by ```` tags, while parameters / variables +are labelled by ````. This labelling is used in SBML to indicate whether +parameters definitions need to be looked up in the list of parameters. + +We can now move on and complete the implementation of our model by adding the +assignment rule for the drug concentration, which relates the drug concentration +to the drug amount, :math:`c = \frac{a}{v}`. Fortunately, at this point you have +already learned almost all the necessary SBML syntax to do this. The only +difference is that assignment rules are indicated with ``assignmentRule`` tags +instead of the rate rule tag + +.. code-block:: xml + + + + + + + + + + + + + + + + + + -1 + elimination_rate + drug_amount + + + + + + + + + + drug_amount + volume + + + + + + + +Executing our Python script again, we obtain an updated model print out + +.. code-block:: python + + [[model]] + name: one_compartment_pk_model + # Initial values + global.drug_amount = 1 + + [global] + dot(drug_amount) = -1 * elimination_rate * drug_amount + drug_concentration = drug_amount / volume + elimination_rate = 1 + time = 0 bind time + in [1] + volume = 1 + +We can see that in contrast to the rate rule, the assignment rule did not promote +the expression for the drug concentration to a differential equation. Instead +it simply assigned the drug concentration to be equal to the drug amount divided +by the volume at all times. + +This completes the implementation of the one-compartment PK model and we can start +to model the pharmacokinetics of drugs. For example, we can simulate the time course +of the drug concentration following a bolus dose of 10 drug amount units + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 287-321 + +.. raw:: html + :file: images/2_mechanistic_model_2.html + +As expected, the drug concentration decays exponentially with the rate defined +by the elimination rate, starting from an initial value of 5, defined by the +initial drug amount and the volume of distribution (:math:`10 / 2`). + +.. note:: + The order of the parameter values when simulating the model follows a + simple pattern: the initial values of state variables come first, followed + by the values of constant parameters. If multiple initial values or parameters + exist, they are order alphabetically. However, to simplify keeping track of + the order of parameters, the model implements a ``parameters()`` method, + :meth:`chi.PKPDModel.parameters`, which returns the parameter names of the + model. + + Similarly, you can get the names of the simulation outputs of the model + with :meth:`chi.PKPDModel.outputs`. The outputs of the model can be changed + using :meth:`chi.PKPDModel.set_outputs` + +SBML in Chi has many more functionalities than those outlined in this tutorial, +including the definition of units for parameters and variables, the definition of +compartments, and the definition of reaction equations. While being clear about +units is of utmost importance across all applications of treatment response +modelling, the ability to define compartments and rate equations is +particularly useful when implementing large mechanistic models, +such as PBPK models and QSP models, because it allows you to +incrementally implement a model in modular blocks, focussing on one reaction at +a time. + +We +will not go into further detail about these elements of SBML in this tutorial. +However, for illustrative purposes we will show a version of the 1-compartment +PK model implementation which includes compartments, rate equations and units + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + central + elimination_rate + drug + + + + + + + + + +For a complete description of SBML, we refer to the SBML documentation: https://sbml.org/. + +Setting dosing regimens +^^^^^^^^^^^^^^^^^^^^^^^ + +In Chi, models implemented using SBML files have a simple two-step interface to +simulate treatment responses for different dosing strategies. The first +step is to specify the route of administration using +:meth:`PKPDModel.set_administration`. Once this step is completed, the dosing +regimen can be scheduled using :meth:`PKPDModel.set_dosing_regimen`. + +The :meth:`PKPDModel.set_administration` method defines the drug amount +variable of the model that gets elevated when dosages are administered. In our +1-compartment PK model from above, there is of course only one drug amount +variable in the model and this flexibility seems unnecessary, however for more +complicated models specifying drug amount +variables in different compartments can be used to emulate different +routes of administration. To facilitate this flexibility of the route of +administration the :meth:`PKPDModel.set_administration` method has one mandatory input parameter +and two optional input parameters. The mandatory input parameter, ``compartment``, +specifies the +compartment of the drug amount variable, the other two input parameters specify +the name of the drug amount +variable, ``amount_var``, (default is ``amount_var='drug_amount'``) and whether the administration is +direct, ``direct``, (default is ``direct=True``). If ``direct`` is set to ``False`` the drug is administered +to the drug amount variable through a dose copmartment. + +Let us first illustrate the direct aministration. In our above example, +``one_compartment_pk_model.xml``, we used rules to define our model and did not +specify any compartments. In this case, Myokit_ assigns all variables to a global +compartment called ``global``. We can therefore administer doses directly to +the drug amount variable by setting the compartment input to ``'global'`` + + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 329 + +To see what happens to the model when we set an administration, we can use +``print(model._model.code())`` to print the model. + +.. code-block:: python + + [[model]] + name: one_compartment_pk_model + # Initial values + global.drug_amount = 1 + + [global] + dose_rate = 0 bind pace + dot(drug_amount) = -1 * elimination_rate * drug_amount + dose_rate + drug_concentration = drug_amount / volume + elimination_rate = 1 + time = 0 bind time + in [1] + volume = 1 + +We can see that two things have happened: 1. a ``dose_rate`` variable has been added +to the model; and 2. the dose rate variable has been added to the right hand +side of the ``drug_amount`` rate equation. The ``dose_rate`` will be set by the +dosing regimen and implements the administration of dosages. + +We can now simulate the drug concentration for different dosing regimens. For +example, assuming the dose amount is in units of mg and the time is in units of +days, we may want to administer a dose of 2 mg every day, or a dose of 4 mg +every 2 days, or administer a total dose of 2 mg every day with an infusion +that lasts for 12 hours. + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 334-374 + +.. raw:: html + :file: images/2_mechanistic_model_3.html + +For full details of which dosing regimens are possible, we refer to API +documentation, :meth:`PKPD.set_dosing_regimen`. + +For completeness, let us administer the same dosing regimens indirectly to the +drug amount variable. + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 381-430 + +.. raw:: html + :file: images/2_mechanistic_model_4.html + +We can see that relative to the direct administration two things have happened: +1. Somehow the number of parameters have changed from 3 parameters to 5 parameters (see line 10 in the code block); +and 2. the peaks of the drug concentrations appear to be rounded off in the figure. + +To understand this, let us again investigate the Myokit_ model using +``print(model._model.code())``, and see how the route of administration has +changed the model + +.. code-block:: python + + [[model]] + name: one_compartment_pk_model + # Initial values + global.drug_amount = 1 + dose.drug_amount = 0 + + [dose] + absorption_rate = 1 + in [1] + dot(drug_amount) = -absorption_rate * drug_amount + global.dose_rate + + [global] + dose_rate = 0 bind pace + dot(drug_amount) = -1 * elimination_rate * drug_amount + dose.absorption_rate * dose.drug_amount + drug_concentration = drug_amount / volume + elimination_rate = 1 + time = 0 bind time + in [1] + volume = 1 + +We can see that the indirect administration makes changes to the model that are +conceptually similar to the direct administration, but instead of adding the +``dose_rate`` variable directly to the right hand side of the ``drug_amount`` +rate equation in the ``global`` compartment, it is added to the rate equation +of a ``drug_amount`` variable in a new ``dose`` compartment. From this ``dose`` +compartment the drug transitions to the ``global`` compartment at a constant rate -- +the ``absorption_rate``. + +These changes of the model explain the shaving off of the drug concentration +peaks: the drug is first adminstered to the dose compartment and then transitions +at a reduced rate to the observed compartment. It also explains why the number of +model parameters increases from 3 to 5: indirect administrations add the initial +value of the ``dose`` compartment ``drug_amount`` variable and the ``absorption_rate`` +parameter to the model. + +.. note:: + Whenever you are in doubt about the order of the parameters, use + :meth:`PKPDModel.parameters` to look up their order. + +Calculating sensitivities +^^^^^^^^^^^^^^^^^^^^^^^^^ + +In Chi, models implemented using SBML files also have a simple interface to +simulate the sensitivities of model parameters. The simulation of sensitivities +can be enabled using the :meth:`PKPDModel.enable_sensitivities` method which +modifies the output of :meth:`PKPDModel.simulate`. If sensitivities are enabled, +:meth:`PKPDModel.simulate` returns a tuple of numpy arrays instead of just a +single numpy array. The first numpy array is the simulated model output, as +before. The second numpy array contains the sensitivities of the +model outputs for each model parameter at all simulated time points and +is of shape ``(n_times, n_outputs, n_parameters)``. + +To see the simulation of sensitivities in action, we use again our 1-compartment +PK model example from above as a use case, ``'one_compartment_pk_model.xml'``, and simulate +the model output and its sensitivities for daily dose administrations of 2 mg. +For simplicity, we choose to visualise the sensitivities only for the initial +drug amount in the ``global`` compartment and the elimination rate. + +.. literalinclude:: code/2_mechanistic_model_2.py + :lines: 446-494 + +.. raw:: html + :file: images/2_mechanistic_model_5.html + +In the code block, you can see that we select the sensitivities by +1. indexing the relevant time +points in the first axis (in this case, we select all simulated time points with +``:``), 2. indexing the relevant output in the second axis, and 3. indexing +the relevant parameter in the third axis. The order of the parameters is +consistent across all methods of the :class:`PKPDModel` and can be looked up +using :meth:`PKPDModel.parameters`. + +In the figure, the sensitivities show that the influence of the initial drug +amount on the drug concentration decays exponentially over time, as indicated +by the decreasing magnitude of its sensitivity. The more intitial drug amount +is available, the larger the drug concentration is for early times of the +simulation. However, +drug concentration levels at later times will not be affected by elevated initial +drug amounts. This aligns with our +understanding of the model, as we know that the model clears drug amounts at an +exponential rate from the system. The initial drug amount can therefore only +influence the drug concentration at early times. + +We can similarly understand +the sensitivity with respect to the elimination rate. The figure shows that +the elimination rate has no influence on the drug concentration at :math:`t=0`. +Only for :math:`t>0`, the elimination rate begins to affect the drug +concentration level. This is because for early times, the drug concentration +value is dominiated by the initial drug amount. As the influence of the +initial drug amount on the drug concentration diminishes, the influence of the +elimination rate increases. Note that the sensitivity with respect to the +elimination rate is negative, in contrast to the positive sensitivity with +respect to the initial drug amount. As a result, larger elimination rates +negatively impact drug conentration levels. + +This concludes this tutorial. If you have any feedback or suggestions for +improvement, or would like to report any typos, mistakes or bugs, please do reach out to us, for example +by creating an Issue_. We are looking +forward to hearing from you! + +Reference to MechanisticModel API +********************************* .. autosummary:: diff --git a/docs/source/getting_started/population_model.rst b/docs/source/getting_started/population_model.rst deleted file mode 100644 index eaa2ac4c..00000000 --- a/docs/source/getting_started/population_model.rst +++ /dev/null @@ -1,18 +0,0 @@ -.. currentmodule:: chi - -******************** -The population model -******************** - -This site is still under construction. In the mean time, we encourage you to -have a look at the detailed API documentation - -.. autosummary:: - - chi.PopulationModel - chi.ComposedPopulationModel - chi.GaussianModel - chi.LogNormalModel - chi.PooledModel - chi.TruncatedGaussianModel - chi.ReducedPopulationModel \ No newline at end of file diff --git a/docs/source/getting_started/quick_overview.rst b/docs/source/getting_started/quick_overview.rst index c26d48bc..b72ca129 100644 --- a/docs/source/getting_started/quick_overview.rst +++ b/docs/source/getting_started/quick_overview.rst @@ -6,11 +6,11 @@ Quick overview ************** -The quick overview displays some of chi's main features to help you decide -whether chi is the right software for you. The running example for this +The quick overview displays some of Chi's main features to help you decide +whether Chi is the right software for you. The running example for this overview will be a 1-compartment pharmacokinetic model (any other deterministic time series model would have been an equally good choice, but the 1-compartment -pharmacokinetic model happens to be predefined in chi's model library; +pharmacokinetic model happens to be predefined in Chi's model library; :class:`chi.library.ModelLibrary`). The model is defined by an ordinary differential equation for the drug amount and an algebraic equation for the drug concentration @@ -29,7 +29,7 @@ units for simplicity). Simulating a mechanistic model ****************************** -Using the 1-compartment pharmacokinetic model from chi's model library, +Using the 1-compartment pharmacokinetic model from Chi's model library, simulation of the model simply involves specifying the set of model parameters :math:`\psi = (a_0, v, k_e)` and the time points of interest. @@ -49,7 +49,7 @@ The simulation results are of shape ``(n_outputs, n_times)`` which in this case is ``(1, 5)`` because we simuated the drug concentration for 5 different time points. For details on how to implement your own :class:`chi.MechanisticModel` and -many other details concerning mechanistic models in chi, we refer to +many other details concerning mechanistic models in Chi, we refer to section :doc:`mechanistic_model`. Visualisation of the simulation @@ -72,7 +72,7 @@ compartment in intervals of 0.5 Simulation of measurements ************************** -Measurements are in chi modelled using a :class:`chi.ErrorModel` on top of the +Measurements are in Chi modelled using a :class:`chi.ErrorModel` on top of the :class:`chi.MechanisticModel` . Formally this takes the form of .. math:: @@ -96,21 +96,21 @@ aestetic reasons) :file: images/1_simulation_2.html For details on how to implement a :class:`chi.ErrorModel` and -many other details concerning error models in chi, we refer to section -:doc:`error_model`. +many other details concerning error models in Chi, we refer to the API +reference :doc:`../api/error_models`. Inference of model parameters ***************************** While the simulation of mechanistic model outputs and measurements is an -interesting feature of chi in its own right, the inference of model parameters +interesting feature of Chi in its own right, the inference of model parameters from real-world measurements is arguably even more interesting. In this overview, we will use the simulated measurements from above to infer model parameters, but in practice the simulated measurements may be straightforwardly replaced by real-world measurements. -Inference in chi leverages the fact that the mechanistic model and the +Inference in Chi leverages the fact that the mechanistic model and the error model define a distribution for the modelled measurements :math:`y`. In the case of a Gaussian error model, as in the previous example, the distribution of the measurements is (yet again) a Gaussian distribution @@ -129,7 +129,7 @@ parameter values for the measurements \log p(\mathcal{D} | \psi , \sigma) = \sum _{j=1}^n \log p(y_j | \psi , \sigma , t_j). -In chi the log-likelihood of model parameters can be defined using +In Chi the log-likelihood of model parameters can be defined using :class:`chi.LogLikelihood`. A :class:`chi.LogLikelihood` defines the distribution of the measurements using a :class:`chi.MechanisticModel` and a :class:`chi.ErrorModel` , and couples the @@ -151,8 +151,8 @@ different parameter values can now be evaluated using the We can see that the data-generating parameter values have a larger log-likelihood than the made-up parameter values (which should intuitively make sense). For details on how to define a :class:`chi.LogLikelihood` and -many other details concerning log-likelihoods in chi, we refer to section -:doc:`log_likelihood`. +many other details concerning log-likelihoods in Chi, we refer to the API +reference :doc:`../api/log_pdfs`. Maximum likelihood estimation ----------------------------- @@ -166,7 +166,7 @@ a.k.a. maximum likelihood estimation \mathop{\text{arg}\,\text{max}}_{\psi, \sigma} \log p(\mathcal{D} | \psi , \sigma). -In chi the maximum likelihood estimates of the model parameters can be found +In Chi the maximum likelihood estimates of the model parameters can be found using any of the standard Python optimisation libraries such as scipy.optimize. We will use Pints_ and its implementation of the Covariance Matrix Adaption-Evolution Strategy (CMA-ES) optimisation algorithm, @@ -220,9 +220,9 @@ approximate the real data-generating processes, we can generalise the notion of data-generating parameters to being the effective set of parameter values that capture the most about the data-generating process within the limitations of the model approximation. Here, the maximum likelihood estimates can analogously -differ significantly from the sought after data-generating parameters. +differ substantially from the sought after data-generating parameters. -In chi the uncertainty of parameter estimates can be estimated using Bayesian +In Chi the uncertainty of parameter estimates can be estimated using Bayesian inference. In Bayesian inference Bayes' rule is used to define a distribution of likely parameter values conditional on the measurements, a.k.a. the posterior parameter distribution @@ -238,7 +238,7 @@ above and :math:`\log p(\psi , \sigma )` is the log-prior distribution of the model parameters. The prior distribution of the model parameters is a modelling choice and captures prior knowlegde about the model parameters. -In chi the log-posterior can be defined using :class:`chi.LogPosterior` which +In Chi the log-posterior can be defined using :class:`chi.LogPosterior` which is instantiated with a :class:`chi.LogLikelihood` and a :class:`pints.LogPrior`. For the sake of this overview, we will use uniform priors that constrain the parameters to values between 0 and 20. In practice, informative priors are @@ -256,8 +256,8 @@ to the :class:`chi.LogLikelihood` using :meth:`chi.LogPosterior.__call__` -8.096007012665059 For details on how to define a :class:`chi.LogPosterior` and -many other details concerning log-posteriors in chi, we refer to section -:doc:`log_posterior`. +many other details concerning log-posteriors in Chi, the API +reference :doc:`../api/log_pdfs`. While the :class:`chi.LogPosterior` allows us to evaluate the log-posterior up to the constant term for different parameter values it does not yet tell us @@ -280,14 +280,14 @@ parameters. .. raw:: html :file: images/1_simulation_5.html -For details on how to infer posterior distributions in chi and +For details on how to infer posterior distributions in Chi and many other details on MCMC sampling, we refer to section -:doc:`mcmc_sampling`. +:doc:`fitting_models_to_data`. Simulating a population model ***************************** -Above we have seen how chi may be used to infer model parameters from +Above we have seen how Chi may be used to infer model parameters from measurements. To this end, we assumed that measurements originate from a data-generating process that is approximated by a mechanistic model, an associated error model and a set of data-generating parameters. @@ -389,7 +389,7 @@ where :math:`\delta (x)` is Dirac's delta distribution which has non-zero probability only for :math:`x=0`. In this example, this ensures that all patients have the same initial drug amount :math:`a_0 = \theta _{a_0}`, the same compartment volume :math:`v = \theta _{v}` and the same measurement noise -:math:`\sigma = \theta _{\sigma}`. In chi delta distributions +:math:`\sigma = \theta _{\sigma}`. In Chi delta distributions can be implemented with a :class:`chi.PooledModel`. The distribution of the elimination rate in the population, :math:`p(k_e | \theta _{k_e})`, is a modelling choice. The simplest, sensible choice for @@ -412,7 +412,7 @@ over time) The resulting model has 5 population parameters :math:`\theta = (\theta _{a_0}, \theta _{v}, \mu _{k_e}, \sigma _{k_e}, \theta _{\sigma})`. -In chi we can define this population model from instances of +In Chi we can define this population model from instances of :class:`chi.PooledModel` and :class:`chi.LogNormalModel` using :class:`chi.ComposedPopulationModel`. From the population model we can then simulate the distribution of the individual parameters in the population via @@ -447,8 +447,8 @@ patients from above. :file: images/1_simulation_8.html For details on how to implement a :class:`chi.PopulationModel` and -many other details concerning population models in chi, we refer to -section :doc:`population_model`. +many other details concerning population models in Chi, the API +reference :doc:`../api/population_models`. Hierarchical inference ********************** @@ -483,7 +483,7 @@ log-likelihood of the bottom-level parameters to describe the observations; and 2. a contribution from the log-likelihood of the population parameters to describe the distribution of the bottom-level parameters. -In chi a hierarchical log-likelihood may be +In Chi a hierarchical log-likelihood may be defined using :class:`chi.HierarchicalLogLikelihood` with a list of bottom-level :class:`chi.LogLikelihood` instances and a :class:`chi.PopulationModel`. @@ -518,7 +518,7 @@ log-posterior .. math:: \log p(\Psi , \theta | \mathcal{D}) = - \log p(\mathcal{D}, \Psi | \theta) + \log p(\theta ) + \log p(\mathcal{D}, \Psi | \theta) + \log p(\theta ) + \text{constant}, where :math:`\log p(\theta )` is the log-prior of the population parameters. @@ -528,7 +528,7 @@ distribution of likely parameters. This will not remove the biases entirely, but make sure that the data-generating parameters are captured by the posterior distribution (as long as the prior choice is appropriate). -In chi we can define a hierarchical log-posterior +In Chi we can define a hierarchical log-posterior using :class:`chi.HierarchicalLogPosterior` which takes a :class:`chi.HierarchicalLogLikelihood` and a :class:`chi.LogPrior` as inputs. The posterior distribution can then be inferred using e.g. MCMC sampling, just @@ -557,19 +557,19 @@ parameters. Note that in this case the posteriors of :math:`\mu _{k_e}` and :math:`\sigma _{k_e}` are largely dominated by the prior distribution as 3 -patients are not informative enough to influence the posterior distribution -significantly. +patients are not informative enough to substantially influence the posterior +distribution. -This concludes the quick overview of chi. You now know how to use chi to +This concludes the quick overview of chi. You now know how to use Chi to simulate mechanistic model outputs, measurements and population models. You also know how to use chi to infer model parameters using maximum likelihood -estimation or Bayesian inference, and you know how to use chi to do +estimation or Bayesian inference, and you know how to use Chi to do hierarchical inference. There are a few things that the quick overview has not touched on. The two most interesting things being: 1. the :class:`chi.ProblemModellingController` which is a convenience -class that helps you to build your models more easiliy, especially when +class that helps you to build your models more easily, especially when measurement times and dosing regimens vary across individuals; 2. the :class:`CovariatePopulationModel` which allows you to define population models where some of the inter-individual is explained by covariates. -We hope you enjoyed this overview and have fun working with chi! \ No newline at end of file +We hope you enjoyed this overview and have fun working with Chi! \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst index 2b0ef0bf..a88412ac 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -2,12 +2,13 @@ .. Root of all chi docs .. _GitHub: https://github.com/DavAug/chi +.. _Myokit: https://myokit.com .. module:: chi .. toctree:: :hidden: - :maxdepth: 1 + :maxdepth: 2 getting_started/index.rst api/index.rst @@ -18,19 +19,79 @@ Welcome to Chi's documentation! **Chi** is an open source Python package hosted on GitHub_, which is designed for pharmacokinetic and pharmacodynamic (PKPD) modelling. -The main features of chi are +The main features of Chi are -- Simulation of mechanistic dose response models (differential equations) for arbitrary dosing regimens. -- Inference of mechanistic model parameters from data (classical or Bayesian). -- Simulation of the dose response variability in a population (hierarchical models/non-linear mixed effects models). -- Inference of population parameters from data (classical or Bayesian). -- Simulation of structured populations, where inter-individual variability can be partly explained by covariates. -- Inference of model parameters in a structured population from data (classical or Bayesian). +- Simulation of treatment responses to custom dosing regimens, using pharmacokinetic & pharmacodynamic (PKPD) models, physiology-based pharmacokinetic (PBPK) models, and/or quantitative systems pharmacology (QSP) models. +- Inference of model parameters from measurements, clinical factors, and/or genetic factors (classical or Bayesian). +- Simulation of treatment response variability across individuals (hierarchical models/non-linear mixed effects models). +- Inference of population parameters from measurements, clinical factors, and/or genetic factors (classical or Bayesian). +- Simulation of structured populations, where inter-individual variability can be partially explained by covariates, such as clinical or genetic factors. +- Inference of model parameters in a structured population from measurements, clinical factors, and/or genetic factors (classical or Bayesian). +- Dosing regimen optimisation and model-informed precision dosing (MIPD). -This page provides the API, or developer documentation for -chi. +This page provides tutorials to illustrate some of Chi's functionality, and a detailed API documentation as a complete reference to all of chi's functions and classes. .. note:: - This package is still in its infancy and is continuously being developed. - So if you find any bugs, please don't hesitate to reach out to us and share - your feedback. + Chi is being continuously developed and improved. + So if you find any bugs or have any suggestions for improvement, please + don't hesitate to reach out to us and share your feedback. + +Install instructions +-------------------- + +Chi can be installed in two steps (one step if you are using windows): +1. installation of a c-library called sundials; and 2. installation of chi. + +Step 1: Installation of sundials +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Chi uses the open source Python package Myokit_ to solve ordinary differential equations +and compute their sensitivities efficiently. Myokit_ does this using a c-library called sundials. +You can install sundials on your computer by entering the below commands in your terminal: + +- On Ubuntu, you can execute the below command to install sundials using ``apt-get``: + +.. code-block:: bash + + apt-get install libsundials-dev + + +- On MacOs, you can execute the below command to install sundials using ``brew``: + +.. code-block:: bash + + brew install sundials + +- On Windows, sundials does not need to be installed manually. Chi will install sundials automatically. + +Step 2: Installation of Chi +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Chi is distributed with PiPy which means that you can pip install Chi with + +.. code-block:: bash + + pip install chi-drm + +If you haven't installed sundials at this point, you will likely get some messages from ``myokit`` complaining +that it cannot find sundials on your machine. In that case, please go back to step 1. + +Note that you need to install ``chi-drm``, and not ``chi``, to install this package. +This has the simple reason that the name ``chi`` was already taken in PiPy when we wanted to +distribute our package. + +Now you are all done and have access to all of Chi's functionality. You can import chi in your python scripts with + +.. code-block:: python + + import chi + +We hope you enjoy using Chi. We are looking forward to seeing which insights you will generate for the pharmaceutical community. + +To get some idea what you can and what you cannot do with Chi, we recommend that you have a look at the tutorials on the following pages. +The API documentation can be found `here `_. + +.. note:: + Note that the package is distributed in PiPy under the name ``chi-drm`` + while in your python scripts you can import the package under the name + ``chi``. diff --git a/run-tests.py b/run-tests.py index 6b7f60db..a368a3ac 100644 --- a/run-tests.py +++ b/run-tests.py @@ -153,6 +153,9 @@ def doctest_example_code(): + '/docs/source/getting_started/code' scripts = os.listdir(script_dir) for script in scripts: + if script[-3:] != '.py': + # Makes sure that this only executes python scripts + continue script = script_dir + '/' + script p = subprocess.Popen([ 'python', diff --git a/setup.py b/setup.py index d690cf23..0c25ae64 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ setup( # Module name name='chi-drm', - version='0.2.3', + version='0.2.4', description='Package to model dose response dynamics', long_description=readme, long_description_content_type="text/markdown", @@ -41,8 +41,9 @@ ], extras_require={ 'docs': [ - 'furo', + 'sphinx-rtd-theme>=1.3', 'sphinx>=1.5, !=1.7.3', # For doc generation + 'sphinx-copybutton>=0.5.2' ], 'notebooks': [ 'jupyter==1.0.0',