# BYOM, byom_doseresp.m

- Author: Tjalling Jager
- Date: April 2017
- Web support: http://www.debtox.info/byom.html
- Back to index walkthrough_doseresp.html

BYOM is a General framework for simulating model systems in terms of ordinary differential equations (ODEs). The model itself needs to be specified in derivatives.m, and call_deri.m may need to be modified to the particular problem as well. The files in the engine directory are needed for fitting and plotting. Results are shown on screen but also saved to a log file (results.out).

**The model:** Log-logistic dose response fitting.

**This script:** Example for reproduction data, using the likelihood function based on a normal distribution of the residuals (various transformation are optional).

## Contents

## Initial things

Make sure that this script is in a directory somewhere **below** the BYOM folder.

clear, clear global % clear the workspace and globals global DATA W X0mat % make the data set and initial states global variables global glo % allow for global parameters in structure glo global pri zvd % global structures for optional priors and zero-variate data diary off % turn of the diary function (if it is accidentaly on) % set(0,'DefaultFigureWindowStyle','docked'); % collect all figure into one window with tab controls set(0,'DefaultFigureWindowStyle','normal'); % separate figure windows pathdefine % set path to the BYOM/engine directory glo.basenm = mfilename; % remember the filename for THIS file for the plots glo.saveplt = 0; % save all plots as (1) Matlab figures, (2) JPEG file or (3) PDF (see all_options.txt)

## The data set

Data are entered in matrix form, time in rows, scenarios (exposure concentrations) in columns. First column are the exposure times, first row are the concentrations or scenario numbers. The number in the top left of the matrix indicates how to calculate the likelihood:

- -2 for binomial likelihood (dose-response fitting of survival data)
- -1 for multinomial likelihood (for survival data)
- 0 for log-transform the data, then normal likelihood
- 0.5 for square-root transform the data, then normal likelihood
- 1 for no transformation of the data, then normal likelihood

% Cumulative reproduction Daphnia after 15 days exposure to cadmium. % Exposure in mg/L. DATA{1} = [1 15 0 93.2 0.085 84.3 0.14 34.4 0.23 13.6]; % Weight factors (number of individuals per observation in data set) W{1} = [9 8 6 1];

## Initial values for the state variables

Initial states, scenarios in columns, states in rows. First row are the 'names' of all scenarios.

X0mat(1,:) = DATA{1}(1,2); % assume second element of first row is the scenario (exposure duration) X0mat(2,:) = 0; % initial values state 1 (not used)

## Initial values for the model parameters

Model parameters are part of a 'structure' for easy reference.

glo.x_EC = 50; % the x in the ECx glo.logsc = 1; % plot the dose-response curve on log-scale % syntax: par.name = [startvalue fit(0/1) minval maxval]; par.ECx = [0.1 1 0 1e6]; % ECx, with x in glo.x_EC par.Y0 = [95 1 0 1e6]; % control response par.beta = [4 1 0 20]; % slope factor of the log-logistic curve

## Time vector and labels for plots

Specify what to plot. If time vector glo.t is not specified, a default is used, based on the data set. Note that t is now used for the exposure concentrations!

% specify the y-axis labels for each state variable glo.ylab{1} = 'reproduction'; % specify the x-axis label (same for all states) glo.xlab = 'concentration (mg/L)'; glo.leglab1 = 'time '; % legend label before the 'scenario' number glo.leglab2 = '(d)'; % legend label after the 'scenario' number prelim_checks % script to perform some preliminary checks and set things up % Note: prelim_checks also fills all the options (opt_...) with defauls, so % modify options after this call, if needed.

## Calculations and plotting

Here, the functions are called that will do the calculation and the plotting. Profile likelihoods are used to make robust confidence intervals. Note that the dose-response curve is plotted on a log scale for the exposure concentration. If the control is truly zero, it is plotted as an open symbol, at a low concentration on the x-axis. However, for the fitting, it is truly zero.

glo.useode = 0; % use the function in simplefun.m opt_optim.fit = 1; % fit the parameters (1), or don't (0) opt_optim.it = 0; % show iterations of the optimisation (1, default) or not (0) opt_plot.annot = 1; % extra subplot in multiplot for fits: 1) box with parameter estimates, 2) overall legend % optimise and plot (fitted parameters in par_out) par_out = calc_optim(par,opt_optim); % start the optimisation calc_and_plot(par_out,opt_plot); % calculate model lines and plot them disp(['x in ECx is: ',num2str(glo.x_EC)])

Provisional goodness-of-fit measures for each data set 0.9944 ================================================================================= Results of the parameter estimation with BYOM version 4.0 ================================================================================= Filename : byom_doseresp Analysis date : 28-Apr-2017 (11:11) Data entered : data state 1: 4x1, continuous data, no transf. Search method: Nelder-Mead simplex direct search, 2 rounds. The optimisation routine has converged to a solution Total 113 simplex iterations used to optimise. Minus log-likelihood has reached the value 15.5298 (AIC=37.0596). ================================================================================= ECx 0.1262 (fit: 1, initial: 0.1) Y0 94.7 (fit: 1, initial: 95) beta 4.864 (fit: 1, initial: 4) ================================================================================= Time required: 0.2 secs Plots result from the optimised parameter values. x in ECx is: 50

## Profiling the likelihood

By profiling you make robust confidence intervals for one or more of your parameters. Use the name of the parameter as it occurs in your parameter structure *par* above. You do not need to run the entire script before you can make a profile.

Options for the profiling can be set using opt_prof (see prelim_checks)

calc_proflik(par_out,'ECx',opt_prof); % calculate a profile

95% confidence interval from the profile ================================================================================= ECx interval: 0.1193 - 0.1349 Time required: 1.9 secs

## Likelihood region

Another way to make intervals on model predictions is to use a sample of parameter sets taken from the joint likelihood-based conf. region. This is done by the function calc_likregion.m. It first does profiling of all fitted parameters to find the edges of the region. Then, Latine-Hypercube sampling, keeping only those parameter combinations that are not rejected at the 95% level in a lik.-rat. test.

Options for the likelihood region can be set using opt_likreg (see prelim_checks.m).

% % UNCOMMENT LINE(S) TO CALCULATE % opt_likreg.detail = 2; % detailed (1) or a course (2) calculation % opt_likreg.subopt = 0; % number of sub-optimisations to perform to increase robustness % calc_likregion(par_out,2000,opt_likreg); % second argument number of samples (-1 to re-use saved sample from previous runs)

## Plot results with confidence intervals

The following code can be used to make a standard plot (the same as for the fits), but with confidence intervals (and sampling error for survival data, if needed). Options for confidence bounds on model curves and calculation of error ellipse can be set using opt_conf (see prelim_checks).

Use opt_conf.type to tell calc_conf which sample to use, and how to use it: 1) Bayesian MCMC sample (default); CI on predictions are 95% ranges on the model curves from the sample 2) parameter sets from the 95% joint likelihood region; CIs on predictions are min-max of the curves from sample of the conf. region 3) as 1 but only 95% of the sample with the highest posterior density used; CIs on predictions are min-max of the curves from sample selection

% % UNCOMMENT LINE(S) TO CALCULATE % opt_conf.type = 2; % use the values from the slice sampler (1) or likelihood region (2) to make intervals % out_conf = calc_conf(par_out,opt_conf); % calculate confidence intervals on model curves % calc_and_plot(par_out,opt_plot,out_conf); % call the plotting routine again to plot fits with CIs