www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/phenobarbdemo.m

    %% Modeling the Population Pharmacokinetics of Phenobarbital in Neonates
%
% This example shows how to build a simple nonlinear mixed-effects model
% from clinical pharmacokinetic data.
%
% Data were collected on 59 pre-term infants given phenobarbital for
% prevention of seizures during the first 16 days after birth. Each
% individual received an initial dose followed by one or more sustaining
% doses by intravenous bolus administration. A total of between 1 and 6
% concentration measurements were obtained from each individual at times
% other than dose times, as part of routine monitoring, for a total of 155
% measurements. Infant weights and APGAR scores (a measure of newborn
% health) were also recorded.
%
% This example uses data described in [1], a study funded by NIH/NIBIB
% grant P41-EB01975.
%
% This example requires Statistics and Machine Learning Toolbox(TM).

% Copyright 2004-2014 The MathWorks, Inc.

%% Load the Data
% These data were downloaded from the website
% |http://depts.washington.edu/rfpk/| (no longer active) of the Resource
% Facility for Population Pharmacokinetics as a text file, and saved as a
% dataset array for ease of use.
load pheno.mat ds

%% Visualize the Data in a Trellis Plot
t = sbiotrellis(ds, 'ID', 'TIME', 'CONC', 'marker', 'o',...
       'markerfacecolor', [.7 .7 .7], 'markeredgecolor', 'r', ...
       'linestyle', 'none');

% Format the plot.
t.plottitle = 'States versus Time';
t.updateFigureForPrinting();


%% Describe the Data
% In order to perform nonlinear mixed-effects modeling on this dataset, we
% need to convert the data to a |groupedData| object, a container for
% holding tabular data that is divided into groups. The |groupedData|
% constructor automatically identifies commonly use variable names as the
% grouping variable or the independent (time) variable. Display the
% properties of the data and confirm that |GroupVariableName| and
% |IndependentVariableName| are correctly identified as |'ID'| and
% |'TIME'|, respectively.
data = groupedData(ds);
data.Properties

%% Define the Model
% We will fit a simple one-compartment model, with bolus dose
% administration and linear clearance elimination, to the data.
pkmd = PKModelDesign;
pkmd.addCompartment('Central', 'DosingType', 'Bolus', 'EliminationType', ...
    'Linear-Clearance', 'HasResponseVariable', true);
model = pkmd.construct;

% The model species |Drug_Central| corresponds to the data variable |CONC|.
responseMap = 'Drug_Central = CONC';

%% Specify Initial Estimates for the Parameters
% The parameters fit in this model are the volume of the central
% compartment (|Central|) and the clearance rate (|Cl_Central|). NLMEFIT
% calculates fixed and random effects for each parameter. The underlying
% algorithm assumes parameters are normally distributed. This assumption
% may not hold for biological parameters that are constrained to be
% positive, such as volume and clearance. We need to specify a transform
% for the estimated parameters, such that the transformed parameters do
% follow a normal distribution. By default, SimBiology(R) chooses a log
% transform for all estimated parameters. Therefore, the model is:
% 
% $$log(V_i)=log(\phi_{V,i})=\theta_V+\eta_{V,i}$$
% 
% and
%
% $$log(Cl_i)=log(\phi_{Cl,i})=\theta_{Cl}+\eta_{Cl,i},$$
%
% where $\theta$, $eta$, and $\phi$ are the fixed effects, random effects,
% and estimated parameter values respectively, calculated for each group
% $i$. We provide some arbitrary initial estimates for V and Cl in the
% absence of better empirical data.
estimatedParams = estimatedInfo({'log(Central)', 'log(Cl_Central)'}, ...
    'InitialValue', [1 1]);

%% Extract the Dosing Information from the Data
% Create a sample dose that targets species |Drug_Central| and use the
% |createDoses| method to generate doses for each infant based on the
% dosing amount listed in variable |DOSE|.
sampleDose = sbiodose('sample', 'TargetName', 'Drug_Central');
doses = createDoses(data, 'DOSE', '', sampleDose);

%% Fit the Model
% Slightly loosen the tolerances to speed up the fit.
fitOptions.Options = statset('TolFun', 1e-3, 'TolX', 1e-3);
[nlmeResults, simI, simP] = sbiofitmixed(model, data, responseMap, ...
    estimatedParams, doses, 'nlmefit', fitOptions); 

%% Visualize the Fitted Model with the Data
% Overlay the fitted simulation results on top of the observed data already
% displayed on the trellis plot. For the population results, simulations
% are run using the estimated fixed effects as the parameter values. For
% the individual results, simulations are run using the sum of the fixed
% and random effects as the parameter values.
t.plot(simP, [], '', 'Drug_Central');
t.plot(simI, [], '', 'Drug_Central',...
    'legend',{'Observed', 'Fit-Pop.', 'Fit-Indiv.'});

%% Examine Fitted Parameters and Covariances
disp('Summary of initial results');
disp('Parameter Estimates Without Random Effects:');
disp(nlmeResults.PopulationParameterEstimates(1:2,:));
disp('Estimated Fixed Effects:');
disp(nlmeResults.FixedEffects);
disp('Estimated Covariance Matrix of Random Effects:');
disp(nlmeResults.RandomEffectCovarianceMatrix);

%% Generate a Box Plot of the Estimated Parameters
% This example uses MATLAB(R) plotting commands to visualize the results.
% Several commonly used plots are also available as built-in options when
% performing parameter fits through the SimBiology(R) desktop interface.

% Create a box plot of the calculated random effects.
boxplot(nlmeResults);

%% Generate a Plot of the Residuals over Time

% The vector of observation data also includes NaN's at the time points at
% which doses were given. We need to remove these NaN's to be able to plot
% the residuals over time.
timeVec = data.(data.Properties.IndependentVariableName);
obsData = data.CONC;
indicesToKeep = ~isnan(obsData);
timeVec = timeVec(indicesToKeep);

% Get the residuals from the fitting results.
indRes = nlmeResults.stats.ires(indicesToKeep);
popRes = nlmeResults.stats.pres(indicesToKeep);

% Plot the residuals. Get a handle to the plot to be able to add more data
% to it later.
resplot = figure;
plot(timeVec,indRes,'d','MarkerFaceColor','b','markerEdgeColor','b');
hold on;
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','b');
hold off;

% Create a reference line representing a zero residual, and set its
% properties to omit this line from the plot legend.
refl = refline(0,0);
refl.Annotation.LegendInformation.IconDisplayStyle = 'off';

% Label the plot.
title('Residuals versus Time');
xlabel('Time');
ylabel('Residuals');
legend({'Individual','Population'});

%% Extract Group-dependent Covariate Values from the Dataset
% Get the mean value of each covariate for each group.
covariateLabels = {'WEIGHT', 'APGAR'};
covariates = grpstats(ds, data.Properties.GroupVariableName, 'mean',...
    'DataVars', covariateLabels);

%% Examine NLME Parameter Fit Results for Possible Covariate Dependencies

% Get the parameter values for each group (empirical Bayes estimates). 
paramValues = nlmeResults.IndividualParameterEstimates.Estimate;
isCentral = strcmp('Central', nlmeResults.IndividualParameterEstimates.Name);
isCl = strcmp('Cl_Central', nlmeResults.IndividualParameterEstimates.Name);

% Plot the parameter values vs. covariates for each group.
figure;
subplot(2,2,1);
plot(covariates.mean_WEIGHT,paramValues(isCentral), '.');
ylabel('Volume');

subplot(2,2,3);
plot(covariates.mean_WEIGHT,paramValues(isCl), '.');
ylabel('Clearance');
xlabel('weight');

subplot(2,2,2);
plot(covariates.mean_APGAR, paramValues(isCentral), '.');

subplot(2,2,4);
plot(covariates.mean_APGAR, paramValues(isCl), '.');
xlabel('APGAR');

%% Create a CovariateModel to Model the Covariate Dependencies
% Based on the plots, there appears to be a correlation between volume and
% weight, clearance and weight, and possibly volume and APGAR score. We
% choose to focus on the effect of weight by modeling two of these
% covariate dependencies: volume ("Central") varying with weight and
% clearance ("Cl_Central") varying with weight. Our model is:
% 
% $$log(V_i)=log(\phi_{V,i})=\theta_V+\theta_{V/weight}*weight_i+\eta_{V,i}$$
% 
% and
%
% $$log(Cl_i)=log(\phi_{Cl,i})=\theta_{Cl}+\theta_{Cl/weight}*weight_i+\eta_{Cl,i}$$

% Define the covariate model.
covmodel            = CovariateModel;
covmodel.Expression = ({'Central = exp(theta1 + theta2*WEIGHT + eta1)','Cl_Central = exp(theta3 + theta4*WEIGHT + eta2)'});

% Use constructDefaultInitialEstimate to create a initialEstimates struct.
initialEstimates = covmodel.constructDefaultFixedEffectValues;
disp('Fixed Effects Description:');
disp(covmodel.FixedEffectDescription);

%%
% Update the initial estimates using the values estimated from fitting the base model.
initialEstimates.theta1 = nlmeResults.FixedEffects.Estimate(1);
initialEstimates.theta3 = nlmeResults.FixedEffects.Estimate(2);
covmodel.FixedEffectValues = initialEstimates;


%% Fit the Model

[nlmeResults_cov, simI_cov, simP_cov] = sbiofitmixed(model, data, responseMap, ...
    covmodel, doses, 'nlmefit', fitOptions); 

%% Examine Fitted Parameters and Covariances
disp('Summary of results when modeling covariate dependencies');
disp('Parameter Estimates Without Random Effects:');
disp(nlmeResults_cov.PopulationParameterEstimates);
disp('Estimated Fixed Effects:');
disp(nlmeResults_cov.FixedEffects);
disp('Estimated Covariance Matrix:');
disp(nlmeResults_cov.RandomEffectCovarianceMatrix);

%% Visualize the Fitted Model with the Data
t.plot(simP_cov, [], '', 'Drug_Central');
t.plot(simI_cov, [], '', 'Drug_Central',...
    'legend',{'Observed', 'Fit-Pop..', 'Fit-Indiv.', 'Cov. Fit-Pop.', 'Cov. Fit-Indiv.'});


%% Compare the Residuals to Those from a Model Without Covariate Dependencies
% The following plot shows that the population residuals are smaller in the
% covariate model fit compared to the original fit.

% Get the residuals from the fitting results.
indRes = nlmeResults_cov.stats.ires(indicesToKeep);
popRes = nlmeResults_cov.stats.pres(indicesToKeep);

% Return to the original residual plot figure and add the new data.
figure(resplot);
hold on;
plot(timeVec,indRes,'d','MarkerFaceColor','r','markerEdgeColor','r');
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','r');
hold off;

% Update the legend.
legend('off');
legend({'Individual','Population','Individual(Cov.)','Population(Cov.)'});

%% References
% [1] Grasela TH Jr, Donn SM. Neonatal population pharmacokinetics of
% phenobarbital derived from routine clinical data. _Dev Pharmacol Ther_
% 1985:8(6). 374-83.
% <http://www.ncbi.nlm.nih.gov/pubmed/4075936?dopt=Abstract PubMed
% Abstract>