Fitting
The model generators return a ModelFunctions object. The model field of a ModelFunctions object is a pure Julia function, allowing to implement the model fitting from scratch.
using AntibodyMethodsDoseResponseConvenience
model, params = accumulation_model(create_grid([1,2,3]))
typeof(model.model) <: FunctiontrueYet, since the model generators create a ModelFunctions object, it is convenient to construct the fitting objective with FittingObjectiveFunctions.jl. Then, only the minimization/maximization of the objective function remains to be implemented. Because these steps are always the same, the adaptive_dose_response_fit function summarizes the creation of the model function and the objective function, requiring only the implementation of a function minimizer.
If there is no reason to avoid the dependencies of AntibodyMethodsDoseResponseConvenience.jl, the workflow as described in the quick start guide should be used. FittingCondition and fit_condition expose the same options that are described here.
The setting
We consider the following measurement data (concentrations, responses, errors):
scatter(concentrations,responses, yerror = errors, xaxis = :log, legend = :none)
Simple model fitting
The data needs to be summarized in a FittingData object, as described in Models:
data = FittingData(concentrations, responses, errors)FittingData([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.00021410640885032888, 0.0007663741382536124, 0.0022904228427114396, 0.007480511858978876, 0.02234921802107462, 0.06800642522811402, 0.1999605075743781, 0.48627646991166545, 0.8640377791882656, 1.1066809652565421, 1.2718676867480627, 1.5954463405164325, 2.0827667975228743, 2.038420572824142, 2.129076978959971, 2.12713266152728], [9.038975931520719e-6, 7.48542293064205e-5, 0.0001715062199824497, 0.0006347118571287526, 0.0009298400467954249, 0.0034907971363670695, 0.002630265790858374, 0.04611199595837917, 0.014791696364573545, 0.06862939569552252, 0.03683894761707063, 0.09538214293549667, 0.15611494403288997, 0.03885184183872844, 0.146064853827133, 0.04829619497432603], Function[FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution, FittingObjectiveFunctions.normal_distribution])Next, a OneDimGrid needs to be created, ideally covering the concentration range:
grid = create_grid(LogRange(1e-10,1e-2,40))AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("ObmYcgCruP" => AdaptiveDensityApproximation.OneDimBlock("ObmYcgCruP", AdaptiveDensityApproximation.Interval{Float64, Float64}(8.376776400682924e-6, 1.3433993325989015e-5), ["PZ0hNQzpBe", "bjCfgdnRkW"], 1.0), "m2mnfvIIWx" => AdaptiveDensityApproximation.OneDimBlock("m2mnfvIIWx", AdaptiveDensityApproximation.Interval{Float64, Float64}(7.017038286703837e-9, 1.1253355826007646e-8), ["f7BJs6IbHA", "QCuKIQBoMW"], 1.0), "AXkRtlUZmM" => AdaptiveDensityApproximation.OneDimBlock("AXkRtlUZmM", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.6037187437513343e-10, 2.571913809059347e-10), ["1kIAUJawhi", "XdfXW9WC4a"], 1.0), "xKKjDLU8ns" => AdaptiveDensityApproximation.OneDimBlock("xKKjDLU8ns", AdaptiveDensityApproximation.Interval{Float64, Float64}(6.614740641230159e-10, 1.0608183551394483e-9), ["IWO6v2aSD7", "2o4jn6PXC0"], 1.0), "YAMJwlYch7" => AdaptiveDensityApproximation.OneDimBlock("YAMJwlYch7", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.00014251026703029993, 0.00022854638641349884), ["teOB8JlWy3", "d35bCQMLJw"], 1.0), "aO8JsUZzJx" => AdaptiveDensityApproximation.OneDimBlock("aO8JsUZzJx", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.000366524123707963, 0.0005878016072274918), ["d35bCQMLJw", "zfQs3TQaRn"], 1.0), "DNpbU89Gxw" => AdaptiveDensityApproximation.OneDimBlock("DNpbU89Gxw", AdaptiveDensityApproximation.Interval{Float64, Float64}(5.5410203300095034e-5, 8.886238162743407e-5), ["3Xsb2ZTZJf", "teOB8JlWy3"], 1.0), "5JqHgOLdHN" => AdaptiveDensityApproximation.OneDimBlock("5JqHgOLdHN", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.0015117750706156647, 0.0024244620170823334), ["3ojA4bSEiI", "FtUjCxyezj"], 1.0), "f7BJs6IbHA" => AdaptiveDensityApproximation.OneDimBlock("f7BJs6IbHA", AdaptiveDensityApproximation.Interval{Float64, Float64}(4.37547937507418e-9, 7.017038286703837e-9), ["M98LuOw7H3", "m2mnfvIIWx"], 1.0), "MnBvkqlPmu" => AdaptiveDensityApproximation.OneDimBlock("MnBvkqlPmu", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.7012542798525856e-9, 2.7283333764867695e-9), ["2o4jn6PXC0", "M98LuOw7H3"], 1.0)…))Finally, a function minimizer needs to be implemented. For this, we use Optim.jl, here:
function minimizer(f,∇f,init)
lower = zeros(length(init))
upper = [Inf for i in 1:length(init)]
return optimize(f,lower,upper, init, Fminbox(NelderMead()),
Optim.Options(g_tol = 1e-12, iterations =2000)).minimizer
endminimizer (generic function with 1 method)The implemented minimizer must take the objective function f, its gradient function ∇f, if applicable, and an initial parameter array init as arguments and return the minimizing parameters. Furthermore, since $K_\tau \geq 0$, the optimization domain should be limited.
Now, a model can be fitted to the data with adaptive_dose_response_fit:
result = adaptive_dose_response_fit(grid,data,
minimizer,
options = AdaptiveOptions(model = accumulation_model)
)AdaptiveResult(DoseResponseResult([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.00021579057077205242, 0.0007200106941876133, 0.00233851080138543, 0.007310962721306748, 0.022365018634090086, 0.06855254335691059, 0.19988219755490758, 0.4865341557734796, 0.8619866014673013, 1.116069451416069, 1.282043175222575, 1.5899059699965077, 1.9337471786888933, 2.0466840753162447, 2.085767901947363, 2.128150357841314]), AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("CFFr16mqVP" => AdaptiveDensityApproximation.OneDimBlock("CFFr16mqVP", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-10, 1.6037187437513343e-10), ["Z3KyK0Lpiy"], 4.832107825716488e-7), "42mX0upIT8" => AdaptiveDensityApproximation.OneDimBlock("42mX0upIT8", AdaptiveDensityApproximation.Interval{Float64, Float64}(4.37547937507418e-9, 7.017038286703837e-9), ["sMJgBlnMHk", "0GmhhkTb5Q"], 0.00021532439734844885), "iTsAt9PhSw" => AdaptiveDensityApproximation.OneDimBlock("iTsAt9PhSw", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.000366524123707963, 0.0005878016072274918), ["dvi92bYdAr", "Si8hzfA11e"], 1.0486616630043566e-5), "EPW04uCBRJ" => AdaptiveDensityApproximation.OneDimBlock("EPW04uCBRJ", AdaptiveDensityApproximation.Interval{Float64, Float64}(3.4551072945922186e-5, 5.5410203300095034e-5), ["dnCm8OkGAL", "U5fIpdFyyC"], 0.12881298051746065), "yJWkMirbGX" => AdaptiveDensityApproximation.OneDimBlock("yJWkMirbGX", AdaptiveDensityApproximation.Interval{Float64, Float64}(8.376776400682924e-6, 1.3433993325989015e-5), ["WEtKDDvB3M", "49y8BI1iOi"], 1.872669767651622e-6), "IUOe38I5vI" => AdaptiveDensityApproximation.OneDimBlock("IUOe38I5vI", AdaptiveDensityApproximation.Interval{Float64, Float64}(6.614740641230159e-10, 1.0608183551394483e-9), ["CcsayAsd7n", "8HVVqnKKp7"], 1.153937222274246e-5), "dvi92bYdAr" => AdaptiveDensityApproximation.OneDimBlock("dvi92bYdAr", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.00022854638641349884, 0.000366524123707963), ["5VNAmvdDYC", "iTsAt9PhSw"], 1.554851977125376e-5), "bvOVRPJ4sS" => AdaptiveDensityApproximation.OneDimBlock("bvOVRPJ4sS", AdaptiveDensityApproximation.Interval{Float64, Float64}(4.641588833612773e-8, 7.443803013251697e-8), ["NiHrS4W9P1", "zQaddlToMN"], 1.3017144684781726e-5), "fjXoq3p3tG" => AdaptiveDensityApproximation.OneDimBlock("fjXoq3p3tG", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.7012542798525856e-9, 2.7283333764867695e-9), ["8HVVqnKKp7", "sMJgBlnMHk"], 0.0007129551951835721), "aXK1bXt9fw" => AdaptiveDensityApproximation.OneDimBlock("aXK1bXt9fw", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.8047217668271702e-8, 2.8942661247167517e-8), ["62bZEsIJ0G", "NiHrS4W9P1"], 0.0004056230740643861)…)), [4.832107825716488e-7, 6.415183730561543e-6, 1.954877714344995e-5, 9.309041596690905e-10, 1.153937222274246e-5, 5.6546893420221556e-5, 0.0007129551951835721, 7.949876627202987e-5, 0.00021532439734844885, 0.0006877136582204332 … 0.320097502888623, 0.20583689191073756, 1.554851977125376e-5, 1.0486616630043566e-5, 9.65296207565328e-6, 1.313506765963266e-5, 0.020623421375504308, 0.03389682370230051, 0.03880926072454446, 0.02342825077704379], 1.7563840306558678, 1.8065428733825684)adaptive_dose_response_fit returns an AdaptiveResult object, that has the following fields:
optimizer: The estimated parameters (result of model fitting).objective_value: The objective function value for the estimated parameters.grid: A grid containing the estimated parameters as grid weights.result: TheDoseResponseResultobject corresponding to the grid (and the offset parameter).time: The elapsed time for the model fit (in seconds).
scatter(data, xaxis = :log, legend = :topleft, label = "data")
plot!(result.result, label = "fit result")
plot(DensityPlot(result.grid), xaxis = :log, color = 2, fill = 0,
fillalpha = 0.5, label = "fit result")
Adaptive model fitting
The model fit above can be improved in two areas. First, a regularization could be used. Second, the adaptive density approximation from AdaptiveDensityApproximation.jl could be used to reduce the number of parameters. Here, we recreate the default optimization from the AntibodyMethodsDoseResponseConvenience.jl package as described in the quick start guide to illustrate some of the available options.
Setting up the objective function properties
For the regularization, a log-posterior objective can be used, where a smoothing prior defines the regularization.
Log-posterior objectives differ from posterior objectives only by taking the logarithm of all functions. Mathematically, there is no difference, but for the computation of tiny probabilities, taking the logarithm upfront is numerically beneficial. Accordingly, the prior also needs to be defined as log-prior, i.e. as the logarithm of the prior.
While FittingObjectiveFunctions.jl expects standard functions as prior/log-prior, adaptive_dose_response_fit requires a prior-generating function (Here, the 500 is used to reproduce the scale = 500 from the quick start guide):
function log_prior_generator(centers, volumes, offset)
ℓV = log.(centers .+ volumes/2) .- log.(centers .- volumes/2)
if isnothing(offset)
return λ -> -500*(sum((λ[i]/ℓV[i]-λ[i+1]/ℓV[i+1])^2 for i in 1:length(λ)-1))/length(λ)^2
else
return λ -> -500*(sum((λ[i]/ℓV[i]-λ[i+1]/ℓV[i+1])^2 for i in 1:length(λ)-2) + λ[end]^2)/length(λ)^2
end
endlog_prior_generator (generic function with 1 method)The prior generator above will create a new log-prior function for each step of the adaptive fitting. The returned log-prior function reads:
\[\text{log-prior}(\lambda) = -\frac{500}{\text{length}(\lambda)^2} \left(\lambda_{\text{offset}} + \sum_{i=1}^{n-1} \left( \frac{\lambda_i}{\ell V_i} - \frac{\lambda_{i+1}}{\ell V_{i+1}} \right)^2 \right)\]
Since the default density-visualization of grids rescales the weights by using the visual interval lengths in a logarithmic scale ℓV (see Background: log-volume normalization), the smoothing should be applied to the rescaled parameters λ. Hence, the λ[i] are divided by ℓV[i].
Next, observe that the log-prior is just the logarithm of a normal distribution (up to a missing normalization) for the difference of the rescaled parameters. Essentially, the prior assumes that there is no difference between neighboring parameters where the scale $\frac{500}{\text{length}(\lambda)^2}$ expresses the strength/importance of this assumption. This is the aforementioned smoothing.
Using log-prior generating functions seems unnecessarily complicated, at first. However, during the adaptive fit, the underlying grid approximation changes, leading to different (visual) interval lengths. Defining a fixed function for the prior could not take the change of interval lengths into account, i.e. the parameters could not be rescaled properly. Hence, the prior needs to be recalculated after every change of the grid, which requires a function that generates the prior from the grid properties.
Without additional information about the measurement errors, a normal distribution is a sensible choice for the uncertainty distribution. Since the goal is a log-posterior objective, the logarithm of a normal distribution (y,m,Δy)-> -(y-m)^2/Δy^2 must be used:
data = FittingData(concentrations,responses, errors, distributions = (y,m,Δy)-> -(y-m)^2/Δy^2)FittingData([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.00021410640885032888, 0.0007663741382536124, 0.0022904228427114396, 0.007480511858978876, 0.02234921802107462, 0.06800642522811402, 0.1999605075743781, 0.48627646991166545, 0.8640377791882656, 1.1066809652565421, 1.2718676867480627, 1.5954463405164325, 2.0827667975228743, 2.038420572824142, 2.129076978959971, 2.12713266152728], [9.038975931520719e-6, 7.48542293064205e-5, 0.0001715062199824497, 0.0006347118571287526, 0.0009298400467954249, 0.0034907971363670695, 0.002630265790858374, 0.04611199595837917, 0.014791696364573545, 0.06862939569552252, 0.03683894761707063, 0.09538214293549667, 0.15611494403288997, 0.03885184183872844, 0.146064853827133, 0.04829619497432603], Function[Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"(), Main.var"#11#12"()])Here, y denotes the data point value, m denotes the model value calculated from the parameters and Δy denotes the measurement error. Uncertainty distributions must take the arguments in this order (y,m,Δy) and must return the distribution / log-distribution value.
Setting up the minimizers
As before, a minimizer needs to be defined:
function minimizer(f,∇f,init)
lower = zeros(length(init))
upper = [Inf for i in 1:length(init)]
return optimize(f,lower,upper, init, Fminbox(NelderMead()),
Optim.Options(g_tol = 1e-12, iterations =2000)).minimizer
endminimizer (generic function with 1 method)Likelihood and posterior objectives usually need to be maximized. However, Optim.jl only provides minimizers, as do some other optimization packages, expecting from the user to flip the sign of the function for a maximization. adaptive_dose_response_fit flips the sign of the posterior and log-posterior objectives automatically.
But, to recreate the default fitting from the quick start guide, a second minimizer is needed (using the LBFGS algorithm):
function post_minimizer(f,∇f,init)
lower = zeros(length(init))
upper = [Inf for i in 1:length(init)]
return optimize(f,lower,upper, init, Fminbox(LBFGS()),
Optim.Options(g_tol = 1e-12, iterations =2000)).minimizer
endpost_minimizer (generic function with 1 method)The second minimizer will be applied after the adaptive optimization to fine-tune the results with a gradient based minimizer.
Setting up the options
Before fitting the data, the fitting options (AdaptiveOptions) and the initial grid need to be defined. Among others, the objective and the prior-generator defined above need to be selected (the additional options in the example are needed to obtain the defaults from the quick start guide):
adaptive_options = AdaptiveOptions(objective = :log_posterior,
iterations = 30,
offset = eps(),
prior_generator = log_prior_generator
)AdaptiveOptions("Adaptive optimization", true, 30, AntibodyMethodsDoseResponse.accumulation_model, 2.220446049250313e-16, :log_posterior, Main.log_prior_generator, nothing, AntibodyMethodsDoseResponse.default_prior_gradient_generator, AntibodyMethodsDoseResponse.log_area_scaled_variation, maximum)The idea of the adaptive approximation is to start with a coarse grid, containing only 2 intervals (3 interval edges):
grid = create_grid(LogRange(1e-10,1e-2,3))AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("rORXhvax2O" => AdaptiveDensityApproximation.OneDimBlock("rORXhvax2O", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-10, 1.0e-6), ["mMf26k9w6y"], 1.0), "mMf26k9w6y" => AdaptiveDensityApproximation.OneDimBlock("mMf26k9w6y", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-6, 0.01000000000000001), ["rORXhvax2O"], 1.0)))With only 2 parameters, common optimizers find a good minimum even for suboptimal initial parameters. Then, the grid can be refined in regions of interest. Next, the previous result can be used as good initial point for the optimization. This process can be repeated several times, increasing the number intervals for regions of interest while not wasting computation time with small intervals for uninteresting regions.
Fitting
Now, the data can be fitted with adaptive_dose_response_fit:
adaptive_result = adaptive_dose_response_fit(grid,data,minimizer, options = fitting_options)To fine-tune the result with the gradient-based minimizer from above, adaptive_dose_response_fit needs to be called again with different options (i.e. no iterations to obtain a single fit and using the offset from the previous result):
post_options = AdaptiveOptions(objective = :log_posterior,
offset = result.optimizer[end],
prior_generator = log_prior_generator
)
adaptive_result = adaptive_dose_response_fit(adaptive_result.grid,data,
post_minimizer,
options = post_options
)AdaptiveResult(DoseResponseResult([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.0002151515291751624, 0.0006921587799722355, 0.0022905029655082786, 0.007443973653780975, 0.022853144957024545, 0.06797238519943972, 0.19966295081460653, 0.49571602988005475, 0.8642379125811639, 1.0893274365958179, 1.2773452268224166, 1.6085308539628609, 1.9318598992532872, 2.0464652408184993, 2.083818824385206, 2.1279853033212195]), AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("3JyfIUOX6B" => AdaptiveDensityApproximation.OneDimBlock("3JyfIUOX6B", AdaptiveDensityApproximation.Interval{Float64, Float64}(7.424145507812504e-5, 9.865527343750006e-5), ["qFUMw1vsMX", "E2XXh5c0Ng"], 0.16472499888424993), "rujLi5QtVb" => AdaptiveDensityApproximation.OneDimBlock("rujLi5QtVb", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-11, 1.563484375e-8), ["vquKmdZOGo"], 0.004341001429126158), "ZxffhPQV32" => AdaptiveDensityApproximation.OneDimBlock("ZxffhPQV32", AdaptiveDensityApproximation.Interval{Float64, Float64}(3.12596875e-8, 6.2509375e-8), ["trHP4iPjpW", "vquKmdZOGo"], 1.727235611365988e-77), "JHECggIBM8" => AdaptiveDensityApproximation.OneDimBlock("JHECggIBM8", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.0007822421875000004, 0.0015634843750000009), ["eBlfrZTAI9", "Hv0Ip73RC3"], 0.001501251526736913), "N4fL5MvOK2" => AdaptiveDensityApproximation.OneDimBlock("N4fL5MvOK2", AdaptiveDensityApproximation.Interval{Float64, Float64}(5.00005e-7, 7.500025e-7), ["poVc2Xv085", "JbWtsGfJBi"], 0.19869363240125143), "wdWYqyQ8h5" => AdaptiveDensityApproximation.OneDimBlock("wdWYqyQ8h5", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.3814659118652343e-6, 1.762931823730469e-6), ["0Nkvdq5e1X", "EgtC40YSEz"], 0.09853321569417299), "iMO3ai37cE" => AdaptiveDensityApproximation.OneDimBlock("iMO3ai37cE", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.3206909179687506e-5, 2.5413818359375012e-5), ["gMRPGbCd8X", "sGYplrh1wa"], 0.019068192982810157), "LKhlmOvMlS" => AdaptiveDensityApproximation.OneDimBlock("LKhlmOvMlS", AdaptiveDensityApproximation.Interval{Float64, Float64}(3.762072753906252e-5, 4.9827636718750027e-5), ["sGYplrh1wa", "qFUMw1vsMX"], 0.10975322269081782), "sGYplrh1wa" => AdaptiveDensityApproximation.OneDimBlock("sGYplrh1wa", AdaptiveDensityApproximation.Interval{Float64, Float64}(2.5413818359375012e-5, 3.762072753906252e-5), ["LKhlmOvMlS", "iMO3ai37cE"], 0.08048299314561375), "bDB0qNnb9H" => AdaptiveDensityApproximation.OneDimBlock("bDB0qNnb9H", AdaptiveDensityApproximation.Interval{Float64, Float64}(0.0031259687500000016, 0.006250937500000003), ["mvxhhAEn2R", "eBlfrZTAI9"], 0.0331808968368835)…)), [0.004341001429126158, 4.2168791772922093e-81, 1.727235611365988e-77, 3.074783068613537e-64, 0.01561989417632643, 0.10123508246723202, 0.11526152110184903, 0.19869363240125143, 0.14554036094467476, 0.15436608431246654 … 0.004093060419472806, 0.001501251526736913, 0.0249094597643404, 0.0331808968368835, 0.02981094997923824, 0.024666346416469635, 0.02136565728044595, 0.01171345131205015, 0.008112278353753133, 1.6356889642709812e-5], 3.320092090041977, 130.98534727096558)As before, the fitted dose-response curve describes the data well:
scatter(data, xaxis = :log, legend = :topleft, label = "data")
plot!(adaptive_result.result, label = "fit result", color = 3)
But comparing the estimated densities from the simple model fit and the adaptive model reveals the difference:
plot(DensityPlot(result.grid), xaxis = :log, color = 2,
fill = 0, fillalpha = 0.5, legend = :topleft, label = "simple fit")
plot!(DensityPlot(adaptive_result.grid), color = 3,
fill = 0, fillalpha = 0.5, label = "adaptive fit")
The same options in the convenience workflow
As mentioned above, the convenience workflow exposes the same options as adaptive_dose_response_fit. Although the adaptive model fitting section uses options that are already the defaults of the convenience workflow, this section creates these options explicitly to illustrate how to modify the options and to illustrate the convenience gain.
First, the log-prior-generator (with the scale 500) can be obtained with scaled_log_volume_prior:
log_prior_generator = scaled_log_volume_prior(500)#1 (generic function with 1 method)The FittingData object is created as before.
data = FittingData(concentrations,responses, errors, distributions = (y,m,Δy)-> -(y-m)^2/Δy^2)FittingData([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.00021410640885032888, 0.0007663741382536124, 0.0022904228427114396, 0.007480511858978876, 0.02234921802107462, 0.06800642522811402, 0.1999605075743781, 0.48627646991166545, 0.8640377791882656, 1.1066809652565421, 1.2718676867480627, 1.5954463405164325, 2.0827667975228743, 2.038420572824142, 2.129076978959971, 2.12713266152728], [9.038975931520719e-6, 7.48542293064205e-5, 0.0001715062199824497, 0.0006347118571287526, 0.0009298400467954249, 0.0034907971363670695, 0.002630265790858374, 0.04611199595837917, 0.014791696364573545, 0.06862939569552252, 0.03683894761707063, 0.09538214293549667, 0.15611494403288997, 0.03885184183872844, 0.146064853827133, 0.04829619497432603], Function[Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"()])The Optim.jl minimizers can be obtained with less boilerplate code, using minimizer_generator:
minimizer = minimizer_generator(NelderMead())
post_minimizer = minimizer_generator(LBFGS())#12 (generic function with 1 method)In general, the grid is created automatically, based on the concentration range of the dose-response data. However, if one concentration is 0, the automatic grid should not be used.
grid = create_grid(LogRange(1e-10,1e-2,3))AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("gCx8gXCWGR" => AdaptiveDensityApproximation.OneDimBlock("gCx8gXCWGR", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-10, 1.0e-6), ["7HYaQSnsls"], 1.0), "7HYaQSnsls" => AdaptiveDensityApproximation.OneDimBlock("7HYaQSnsls", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-6, 0.01000000000000001), ["gCx8gXCWGR"], 1.0)))The automatic grids are created with the LogRange function and are intended for logarithmic plots. A zero-concentration would lead to a DomainError, as LogRange demands positive numbers. To avoid an error, the automatic grid generator then substitutes eps() for 0 (only for the auto-generated grid) and raise a warning. However, too large grid domains lead to poor fit results. Hence, it is recommended to create the grid manually in those cases.
The same options as before are used:
adaptive_options = AdaptiveOptions(objective = :log_posterior, iterations = 30,
offset = eps(), prior_generator = log_prior_generator)
post_options = AdaptiveOptions(objective = :log_posterior, offset = eps(),
prior_generator = log_prior_generator)AdaptiveOptions("Adaptive optimization", true, 1, AntibodyMethodsDoseResponse.accumulation_model, 2.220446049250313e-16, :log_posterior, AntibodyMethodsDoseResponseConvenience.var"#1#6"{Int64}(500), nothing, AntibodyMethodsDoseResponse.default_prior_gradient_generator, AntibodyMethodsDoseResponse.log_area_scaled_variation, maximum)Finally, the objects, functions and options from above need to be summarized in a FittingCondition object:
condition = FittingCondition(data,
grid = grid,
options_1 = adaptive_options,
options_2 = post_options,
minimizer_1 = minimizer,
minimizer_2 = post_minimizer
)FittingCondition(FittingData([1.0e-10, 3.414548873833601e-10, 1.1659144011798312e-9, 3.981071705534969e-9, 1.3593563908785241e-8, 4.641588833612773e-8, 1.584893192461114e-7, 5.411695265464649e-7, 1.8478497974222906e-6, 6.309573444801943e-6, 2.1544346900318823e-5, 7.356422544596421e-5, 0.0002511886431509585, 0.0008576958985908955, 0.0029286445646252404, 0.01000000000000001], [0.00021410640885032888, 0.0007663741382536124, 0.0022904228427114396, 0.007480511858978876, 0.02234921802107462, 0.06800642522811402, 0.1999605075743781, 0.48627646991166545, 0.8640377791882656, 1.1066809652565421, 1.2718676867480627, 1.5954463405164325, 2.0827667975228743, 2.038420572824142, 2.129076978959971, 2.12713266152728], [9.038975931520719e-6, 7.48542293064205e-5, 0.0001715062199824497, 0.0006347118571287526, 0.0009298400467954249, 0.0034907971363670695, 0.002630265790858374, 0.04611199595837917, 0.014791696364573545, 0.06862939569552252, 0.03683894761707063, 0.09538214293549667, 0.15611494403288997, 0.03885184183872844, 0.146064853827133, 0.04829619497432603], Function[Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"(), Main.var"#1#2"()]), nothing, AdaptiveDensityApproximation.OneDimGrid{Dict{String, AdaptiveDensityApproximation.OneDimBlock}}(Dict{String, AdaptiveDensityApproximation.OneDimBlock}("gCx8gXCWGR" => AdaptiveDensityApproximation.OneDimBlock("gCx8gXCWGR", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-10, 1.0e-6), ["7HYaQSnsls"], 1.0), "7HYaQSnsls" => AdaptiveDensityApproximation.OneDimBlock("7HYaQSnsls", AdaptiveDensityApproximation.Interval{Float64, Float64}(1.0e-6, 0.01000000000000001), ["gCx8gXCWGR"], 1.0))), "", AdaptiveOptions("Adaptive optimization", true, 30, AntibodyMethodsDoseResponse.accumulation_model, 2.220446049250313e-16, :log_posterior, AntibodyMethodsDoseResponseConvenience.var"#1#6"{Int64}(500), nothing, AntibodyMethodsDoseResponse.default_prior_gradient_generator, AntibodyMethodsDoseResponse.log_area_scaled_variation, maximum), AdaptiveOptions("Adaptive optimization", true, 1, AntibodyMethodsDoseResponse.accumulation_model, 2.220446049250313e-16, :log_posterior, AntibodyMethodsDoseResponseConvenience.var"#1#6"{Int64}(500), nothing, AntibodyMethodsDoseResponse.default_prior_gradient_generator, AntibodyMethodsDoseResponse.log_area_scaled_variation, maximum), AntibodyMethodsDoseResponseConvenience.var"#12#14"{Optim.Options{Float64, Nothing}, NelderMead{Optim.AffineSimplexer, Optim.AdaptiveParameters}}(Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-12, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = true, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 2000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
, NelderMead{Optim.AffineSimplexer, Optim.AdaptiveParameters}(Optim.AffineSimplexer(0.025, 0.5), Optim.AdaptiveParameters(1.0, 1.0, 0.75, 1.0))), AntibodyMethodsDoseResponseConvenience.var"#12#14"{Optim.Options{Float64, Nothing}, LBFGS{Nothing, InitialStatic{Float64}, HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}}(Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-12, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = true, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 2000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
, LBFGS{Nothing, InitialStatic{Float64}, HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}(10, InitialStatic{Float64}
alpha: Float64 1.0
scaled: Bool false
, HagerZhang{Float64, Base.RefValue{Bool}}
delta: Float64 0.1
sigma: Float64 0.9
alphamax: Float64 Inf
rho: Float64 5.0
epsilon: Float64 1.0e-6
gamma: Float64 0.66
linesearchmax: Int64 50
psi3: Float64 0.1
display: Int64 0
mayterminate: Base.RefValue{Bool}
cache: Nothing nothing
, nothing, Optim.var"#19#21"(), Flat(), true)), nothing)Fitting the condition is done by calling fit_condition:
result = fit_condition(condition)
plot(DensityPlot(result.grid), xaxis = :log, fill = 0, fillalpha = 0.5)
Multi-threaded fitting
The fitting of a single condition cannot be parallelized. However, fitting more than one condition can be done in parallel. Consider an array of FittingCondition objects:
conditions = [condition_1,condition_2,condition_3]The conditions can be fitted in parallel with fit_conditions (observe the additional s):
fit_conditions(conditions)The number of threads cannot be changed after Julia is launched. The number of threads needs to be set in before, e.g. with julia -t 8 to run Julia with 8 threads. The number of threads can be checked with Threads.nthreads().
Changing the initial values
In each case discussed here, no initial parameter values were set for the model fitting. This is, because the wights of the gird are used as initial values. By default, create_grid sets all weights to 1:
grid = create_grid(LogRange(1e-10,1e-2,3))
export_weights(grid)2-element Vector{Float64}:
1.0
1.0Changing the weights of the gird before running the fitting-process allows to specify the initial values for the model fitting. E.g. for the manual definition of adaptive fitting:
import_weights!(grid,[eps(),1])
adaptive_result = adaptive_dose_response_fit(grid,data,minimizer, options = fitting_options)Or for the convenience workflow, e.g. after the FittingCondition has already been defined:
import_weights!(condition.grid,[eps(),1])
result = fit_condition(condition)