В этом примере показано, как создать простую нелинейную модель смешанных эффектов из клинических фармакокинетических данных.
Данные были собраны по 59 недоношенным детям, данным фенобарбитал для предотвращения занятости в течение первых 16 дней после рождения. Каждый индивидуум получил начальную дозу, сопровождаемую одной или несколькими дозами поддержки внутривенным администрированием шарика. В общей сложности между 1 и 6 концентрациями измерения получались от каждого индивидуума время от времени кроме времен дозы, как часть стандартного контроля, для в общей сложности 155 измерений. Младенческие веса и баллы APGAR (мера новорожденного здоровья) были также зарегистрированы.
Этот пример использует данные, описанные в [1], исследование, финансируемое NIH/NIBIB, предоставляют P41-EB01975.
Этот пример требует Statistics and Machine Learning Toolbox™.
Эти данные были загружены с веб-сайта http://depts.washington.edu/rfpk/
(более не активный) из Средства Ресурса для Фармакокинетики Населения как текстовый файл и сохраненный как массив набора данных для простоты использования.
load pheno.mat ds
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();
Для того, чтобы выполнить нелинейное моделирование смешанных эффектов на этом наборе данных, мы должны преобразовать данные в groupedData
объект, контейнер для содержания табличных данных, который разделен на группы. groupedData
конструктор автоматически обычно идентифицирует имена переменных использования как сгруппированную переменную или независимого политика (время) переменная. Отобразите свойства данных и подтвердите тот GroupVariableName
и IndependentVariableName
правильно идентифицированы как 'ID'
и 'TIME'
, соответственно.
data = groupedData(ds); data.Properties
ans = struct with fields:
Description: ''
UserData: []
DimensionNames: {'Observations' 'Variables'}
VariableNames: {'ID' 'TIME' 'DOSE' 'WEIGHT' 'APGAR' 'CONC'}
VariableDescriptions: {}
VariableUnits: {}
VariableContinuity: []
RowNames: {}
CustomProperties: [1×1 matlab.tabular.CustomProperties]
GroupVariableName: 'ID'
IndependentVariableName: 'TIME'
Мы будем подбирать простую модель с одним отсеком, с администрированием дозы шарика и линейным устранением разрешения, к данным.
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';
Параметры помещаются в эту модель, объем центрального отсека (Central
) и уровень раскрываемости преступлений (Cl_Central
). NLMEFIT вычисляет зафиксированные и случайные эффекты для каждого параметра. Базовый алгоритм принимает, что параметры нормально распределены. Это предположение не может содержать для биологических параметров, которые ограничиваются быть положительными, такими как объем и разрешение. Мы должны задать преобразование для предполагаемых параметров, таких, что преобразованные параметры действительно следуют за нормальным распределением. По умолчанию SimBiology® выбирает, журнал преобразовывают для всех предполагаемых параметров. Поэтому модель:
и
где , , и фиксированные эффекты, случайные эффекты и оцененные значения параметров соответственно, вычисленный для каждой группы . Мы обеспечиваем некоторые произвольные первоначальные оценки для V и Статья в отсутствие лучших эмпирических данных.
estimatedParams = estimatedInfo({'log(Central)', 'log(Cl_Central)'}, ... 'InitialValue', [1 1]);
Создайте демонстрационную дозу, которая предназначается для разновидностей Drug_Central
и используйте createDoses
метод, чтобы сгенерировать дозы для каждого младенца на основе объема дозирования, перечисленного в переменной DOSE
.
sampleDose = sbiodose('sample', 'TargetName', 'Drug_Central'); doses = createDoses(data, 'DOSE', '', sampleDose);
Немного ослабьте допуски, чтобы ускорить подгонку.
fitOptions.Options = statset('TolFun', 1e-3, 'TolX', 1e-3); [nlmeResults, simI, simP] = sbiofitmixed(model, data, responseMap, ... estimatedParams, doses, 'nlmefit', fitOptions);
Наложите подходящие результаты симуляции сверху наблюдаемых данных, уже отображенных на графике решетки. Для результатов населения симуляции запущены с помощью предполагаемых фиксированных эффектов в качестве значений параметров. Для отдельных результатов симуляции запущены с помощью суммы фиксированных и случайных эффектов как значения параметров.
t.plot(simP, [], '', 'Drug_Central'); t.plot(simI, [], '', 'Drug_Central',... 'legend',{'Observed', 'Fit-Pop.', 'Fit-Indiv.'});
disp('Summary of initial results');
Summary of initial results
disp('Parameter Estimates Without Random Effects:');
Parameter Estimates Without Random Effects:
disp(nlmeResults.PopulationParameterEstimates(1:2,:));
Group Name Estimate _____ ______________ ________ 1 {'Central' } 1.4086 1 {'Cl_Central'} 0.006137
disp('Estimated Fixed Effects:');
Estimated Fixed Effects:
disp(nlmeResults.FixedEffects);
Name Description Estimate StandardError __________ ______________ ________ _____________ {'theta1'} {'Central' } 0.34257 0.061246 {'theta2'} {'Cl_Central'} -5.0934 0.079501
disp('Estimated Covariance Matrix of Random Effects:');
Estimated Covariance Matrix of Random Effects:
disp(nlmeResults.RandomEffectCovarianceMatrix);
eta1 eta2 _______ _______ eta1 0.20323 0 eta2 0 0.19338
Этот пример использует команды графического вывода MATLAB®, чтобы визуализировать результаты. Несколько обычно используемых графиков также доступны как встроенные опции при выполнении подгонок параметра через интерфейс рабочего стола SimBiology®.
% Create a box plot of the calculated random effects.
boxplot(nlmeResults);
% 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'});
Получите среднее значение каждого коварианта для каждой группы.
covariateLabels = {'WEIGHT', 'APGAR'}; covariates = grpstats(ds, data.Properties.GroupVariableName, 'mean',... 'DataVars', covariateLabels);
% 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');
На основе графиков, кажется, существует корреляция между объемом и весом, разрешением и весом, и возможно счетом APGAR и объемом. Мы принимаем решение фокусироваться на эффекте веса путем моделирования двух из этих ковариационных зависимостей: объем ("Центральное") меняться в зависимости от веса и разрешения ("Cl_Central"), меняющийся в зависимости от веса. Наша модель:
и
% 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:');
Fixed Effects Description:
disp(covmodel.FixedEffectDescription);
{'Central' } {'Cl_Central' } {'Central/WEIGHT' } {'Cl_Central/WEIGHT'}
Обновите первоначальные оценки с помощью значений, оцененных от того, чтобы подбирать базовую модель.
initialEstimates.theta1 = nlmeResults.FixedEffects.Estimate(1); initialEstimates.theta3 = nlmeResults.FixedEffects.Estimate(2); covmodel.FixedEffectValues = initialEstimates;
[nlmeResults_cov, simI_cov, simP_cov] = sbiofitmixed(model, data, responseMap, ... covmodel, doses, 'nlmefit', fitOptions);
disp('Summary of results when modeling covariate dependencies');
Summary of results when modeling covariate dependencies
disp('Parameter Estimates Without Random Effects:');
Parameter Estimates Without Random Effects:
disp(nlmeResults_cov.PopulationParameterEstimates);
Group Name Estimate _____ ______________ _________ 1 {'Central' } 1.2973 1 {'Cl_Central'} 0.0061576 2 {'Central' } 1.3682 2 {'Cl_Central'} 0.0065512 3 {'Central' } 1.3682 3 {'Cl_Central'} 0.0065512 4 {'Central' } 0.99431 4 {'Cl_Central'} 0.0045173 5 {'Central' } 1.2973 5 {'Cl_Central'} 0.0061576 6 {'Central' } 1.1664 6 {'Cl_Central'} 0.00544 7 {'Central' } 1.0486 7 {'Cl_Central'} 0.004806 8 {'Central' } 1.1664 8 {'Cl_Central'} 0.00544 9 {'Central' } 1.2973 9 {'Cl_Central'} 0.0061576 10 {'Central' } 1.2973 10 {'Cl_Central'} 0.0061576 11 {'Central' } 1.1664 11 {'Cl_Central'} 0.00544 12 {'Central' } 1.2301 12 {'Cl_Central'} 0.0057877 13 {'Central' } 1.1059 13 {'Cl_Central'} 0.0051132 14 {'Central' } 1.1059 14 {'Cl_Central'} 0.0051132 15 {'Central' } 1.2301 15 {'Cl_Central'} 0.0057877 16 {'Central' } 1.1664 16 {'Cl_Central'} 0.00544 17 {'Central' } 1.1059 17 {'Cl_Central'} 0.0051132 18 {'Central' } 1.0486 18 {'Cl_Central'} 0.004806 19 {'Central' } 1.0486 19 {'Cl_Central'} 0.004806 20 {'Central' } 1.1664 20 {'Cl_Central'} 0.00544 21 {'Central' } 1.605 21 {'Cl_Central'} 0.0078894 22 {'Central' } 1.3682 22 {'Cl_Central'} 0.0065512 23 {'Central' } 3.2052 23 {'Cl_Central'} 0.017654 24 {'Central' } 3.3803 24 {'Cl_Central'} 0.018782 25 {'Central' } 0.89394 25 {'Cl_Central'} 0.0039908 26 {'Central' } 3.9653 26 {'Cl_Central'} 0.022619 27 {'Central' } 1.6927 27 {'Cl_Central'} 0.0083936 28 {'Central' } 3.3803 28 {'Cl_Central'} 0.018782 29 {'Central' } 1.0486 29 {'Cl_Central'} 0.004806 30 {'Central' } 1.605 30 {'Cl_Central'} 0.0078894 31 {'Central' } 1.2973 31 {'Cl_Central'} 0.0061576 32 {'Central' } 4.1819 32 {'Cl_Central'} 0.024064 33 {'Central' } 1.5218 33 {'Cl_Central'} 0.0074154 34 {'Central' } 1.5218 34 {'Cl_Central'} 0.0074154 35 {'Central' } 2.3292 35 {'Cl_Central'} 0.012173 36 {'Central' } 1.3682 36 {'Cl_Central'} 0.0065512 37 {'Central' } 1.1664 37 {'Cl_Central'} 0.00544 38 {'Central' } 1.2301 38 {'Cl_Central'} 0.0057877 39 {'Central' } 1.6927 39 {'Cl_Central'} 0.0083936 40 {'Central' } 1.1059 40 {'Cl_Central'} 0.0051132 41 {'Central' } 1.5218 41 {'Cl_Central'} 0.0074154 42 {'Central' } 2.7323 42 {'Cl_Central'} 0.014659 43 {'Central' } 0.99431 43 {'Cl_Central'} 0.0045173 44 {'Central' } 1.2973 44 {'Cl_Central'} 0.0061576 45 {'Central' } 0.94279 45 {'Cl_Central'} 0.0042459 46 {'Central' } 1.1059 46 {'Cl_Central'} 0.0051132 47 {'Central' } 2.4565 47 {'Cl_Central'} 0.012951 48 {'Central' } 0.89394 48 {'Cl_Central'} 0.0039908 49 {'Central' } 1.2301 49 {'Cl_Central'} 0.0057877 50 {'Central' } 1.1059 50 {'Cl_Central'} 0.0051132 51 {'Central' } 0.99431 51 {'Cl_Central'} 0.0045173 52 {'Central' } 0.99431 52 {'Cl_Central'} 0.0045173 53 {'Central' } 1.5218 53 {'Cl_Central'} 0.0074154 54 {'Central' } 1.605 54 {'Cl_Central'} 0.0078894 55 {'Central' } 1.1059 55 {'Cl_Central'} 0.0051132 56 {'Central' } 0.84763 56 {'Cl_Central'} 0.003751 57 {'Central' } 1.8827 57 {'Cl_Central'} 0.0095009 58 {'Central' } 1.2973 58 {'Cl_Central'} 0.0061576 59 {'Central' } 1.1059 59 {'Cl_Central'} 0.0051132
disp('Estimated Fixed Effects:');
Estimated Fixed Effects:
disp(nlmeResults_cov.FixedEffects);
Name Description Estimate StandardError __________ _____________________ ________ _____________ {'theta1'} {'Central' } -0.48453 0.067952 {'theta3'} {'Cl_Central' } -5.9575 0.12199 {'theta2'} {'Central/WEIGHT' } 0.53203 0.040788 {'theta4'} {'Cl_Central/WEIGHT'} 0.61957 0.074264
disp('Estimated Covariance Matrix:');
Estimated Covariance Matrix:
disp(nlmeResults_cov.RandomEffectCovarianceMatrix);
eta1 eta2 ________ _______ eta1 0.029793 0 eta2 0 0.04644
t.plot(simP_cov, [], '', 'Drug_Central'); t.plot(simI_cov, [], '', 'Drug_Central',... 'legend',{'Observed', 'Fit-Pop..', 'Fit-Indiv.', 'Cov. Fit-Pop.', 'Cov. Fit-Indiv.'});
Следующий график показывает, что остаточные значения населения меньше в ковариационной подгонке модели по сравнению с исходной подгонкой.
% 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.)'});
[1] Грэзела Т младший, Донн СМ. Неонатальная фармакокинетика населения фенобарбитала выведена из стандартных клинических данных. Фармакол Dev Там 1985:8 (6). 374-83.