В этом примере показано, как создать простую нелинейную модель смешанных эффектов из клинических фармакокинетических данных.
Данные были собраны по 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.