Этот пример показывает, как построить простую нелинейную модель смешанных эффектов из клинических фармакокинетических данных.
Были собраны данные о 59 недоношенных детях, которым вводили фенобарбитал для профилактики судорог в течение первых 16 дней после рождения. Каждый индивидуум получал начальную дозу с последующей одной или несколькими поддерживающими дозами путем внутривенного введения болюса. Всего от 1 до 6 измерений концентрации были получены от каждого индивидуума в моменты времени, отличные от времени дозы, как часть обычного мониторинга, в общей сложности 155 измерений. Вес младенцев и баллы APGAR (показатель здоровья новорожденных) также были зарегистрированы.
В этом примере используются данные, описанные в [1], исследовании, финансируемом NIH/NIBIB grant P41-EB01975.
В этом примере требуется 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 ® выбирает преобразование журнала для всех расчетных параметров. Поэтому модель выглядит следующим образом:
и
θCl +ηCl, я,
где , и - фиксированные эффекты, случайные эффекты и оцененные значения параметров, соответственно, рассчитанные для каждой группы. Мы предоставляем некоторые произвольные начальные оценки для V и Cl в отсутствие лучших эмпирических данных.
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. Мы выбираем сфокусироваться на эффекте веса, моделируя две из этих ковариатных зависимостей: объем («Central»), изменяющийся с весом, и клиренс («Cl _ Central»), изменяющийся с весом. Наша модель:
и
+ηCl, я
% 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] Грасела ТХ-младший, Донн СМ. Фармакокинетика фенобарбитала неонатальной популяции, полученная из обычных клинических данных. Дев Фармакол Тер 1985:8 (6). 374-83.