%% Fit Two-Compartment Model to PK Profiles of Multiple Individuals

%% Load Data
load('data10_32R.mat')
%% Convert to groupedData Format
gData = groupedData(data);
gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter'};
gData.Properties

%% Visualize Data
sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},...
            'Marker','+','LineStyle','none');

%% Set Up a Two-Compartment 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;

%% Define Dosing
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';

%% Map the Response Data to Corresponding Model Components
responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};

%% Specify Parameters to Estimate
paramsToEstimate   = {'log(Central)','log(Peripheral)','Q12','Cl_Central'};
estimatedParam     = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);

%% Perform a Pooled Fit
pooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',true)

%% Display Results
plot(pooledFit);

%% Perform an Unpooled Fit
unpooledFit =  sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',false);

%% Plot Unpooled Results
plot(unpooledFit);

%% Display Individual 1 Results
unpooledFit(1)

%% Pooled vs. Unpooled
t = [gData.Time;gData.Time];
res_pooled = vertcat(pooledFit.R);
res_pooled = res_pooled(:);
res_unpooled = vertcat(unpooledFit.R);
res_unpooled = res_unpooled(:);
plot(t,res_pooled,'o','MarkerFaceColor','w','markerEdgeColor','b')
hold on
plot(t,res_unpooled,'o','MarkerFaceColor','b','markerEdgeColor','b')
refl = refline(0,0); % A reference line representing a zero residual
title('Residuals versus Time');
xlabel('Time');
ylabel('Residuals');
legend({'Pooled','Unpooled'});


%% 
% Copyright 2012 The MathWorks, Inc.