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>