%% Estimate Two-Compartment PK Parameters
%
% Load the sample data set.

load data10_32R.mat
gData = groupedData(data);
gData.Properties.VariableUnits = ["","hour","milligram/liter","milligram/liter"];

%
% Create a two-compartment PK model.

pkmd                 = PKModelDesign;
pkc1                 = addCompartment(pkmd,"Central");
pkc1.DosingType      = "Infusion";
pkc1.EliminationType = "linear-clearance";
pkc1.HasResponseVariable = true;
pkc2                 = addCompartment(pkmd,"Peripheral");
model                = construct(pkmd);
configset            = getconfigset(model);
configset.CompileOptions.UnitConversion = true;
responseMap = ["Drug_Central = CentralConc","Drug_Peripheral = PeripheralConc"];

%
% Provide model parameters to estimate.

paramsToEstimate   = ["log(Central)","log(Peripheral)","Q12","Cl_Central"];
estimatedParam     = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);

%
% Assume every individual receives an infusion dose at time = 0, with a total infusion amount of 100 mg at a rate of 50 mg/hour.

dose             = sbiodose("dose","TargetName","Drug_Central");
dose.StartTime   = 0;
dose.Amount      = 100;
dose.Rate        = 50;
dose.AmountUnits = "milligram";
dose.TimeUnits   = "hour";
dose.RateUnits   = "milligram/hour";

%
% Estimate model parameters. By default, the function estimates a set of parameter for each individual (unpooled fit).

fitResults =  sbiofit(model,gData,responseMap,estimatedParam,dose);

%
% Plot the results.

plot(fitResults);

%
% Plot all groups in one plot.

plot(fitResults,"PlotStyle","one axes");

%
% Change some axes properties.

s = struct;
s.Properties.XGrid = "on";
s.Properties.YGrid = "on";
plot(fitResults,"PlotStyle","one axes","AxesStyle",s);

%
% Compare the model predictions to the actual data.

plotActualVersusPredicted(fitResults)

%
% Use boxplot to show the variation of estimated model parameters.

boxplot(fitResults)

%
% Plot the distribution of residuals. This normal probability plot shows the deviation from normality and the skewness on the right tail of the distribution of residuals. The default (constant) error model might not be the correct assumption for the data being fitted.

plotResidualDistribution(fitResults)

%
% Plot residuals for each response using the model predictions on x-axis.

plotResiduals(fitResults,"Predictions")

%
% Get the summary of the fit results. stats.Name contains the name for each table from stats.Table, which contains a list of tables with estimated parameter values and fit quality statistics.

stats = summary(fitResults);
stats.Name
stats.Table

% Copyright 2020 The MathWorks, Inc.
%
%