Фон

Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от 30 человек и вы хотите оценить фармакокинетические параметры, а именно объемы центрального и периферического отделения, клиренс и межкомпартментальный клиренс. Предположим, что концентрация лекарственного средства по сравнению с профилем времени следует за бикспоненциальным снижением $${}_{\mathrm{Ct}}={}^{}\mathrm{Ae-at}{}^{+}$$Be-bt, $${}_{}$$где Ct - концентрация лекарственного средства в момент t, $$$$а a $$$$и b - наклоны для соответствующих экспоненциальных понижений.

Загрузить данные

Эти синтетические данные содержат временной курс концентраций в плазме 30 особей после болюсной дозы (100 мг), измеренной в разное время как для центрального, так и для периферического отделения. Он также содержит категориальные переменные, а именно: Sex и Age.

clear
load('sd5_302RAgeSex.mat')

Преобразуйте в Формат groupedData

Преобразовать набор данных в groupedData объект, который является требуемым форматом данных для функции фитинга sbiofit. A groupedData объект также позволяет задавать независимые имена переменных и групп (если они существуют). Установка единиц измерения ID, Time, CentralConc, PeripheralConc, Age, и Sex переменные. Установки являются необязательными и обязательными только для UnitConversion функция, которая автоматически преобразует соответствующие физические величины в одну согласованную систему единиц измерения.

gData = groupedData(data);
gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter','',''};
gData.Properties

ans = struct with fields:
                Description: ''
                   UserData: []
             DimensionNames: {'Row'  'Variables'}
              VariableNames: {1x6 cell}
       VariableDescriptions: {}
              VariableUnits: {1x6 cell}
         VariableContinuity: []
                   RowNames: {}
           CustomProperties: [1x1 matlab.tabular.CustomProperties]
          GroupVariableName: 'ID'
    IndependentVariableName: 'Time'

IndependentVariableName и GroupVariableName свойства были автоматически установлены в значение Time и ID переменные данных.

Визуализация данных

Просмотрите данные ответа для каждого отдельного пользователя.

t = sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},...
                'Marker','+','LineStyle','none');
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Настройка двухкамерной модели

Используйте встроенную библиотеку ПК для построения двухкамерной модели с инфузионным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отсека. Используйте configset объект для включения преобразования единиц измерения.

pkmd                                    = PKModelDesign;
pkc1                                    = addCompartment(pkmd,'Central');
pkc1.DosingType                         = 'Bolus';
pkc1.EliminationType                    = 'linear-clearance';
pkc1.HasResponseVariable                = true;
pkc2                                    = addCompartment(pkmd,'Peripheral');
model                                   = construct(pkmd);
configset                               = getconfigset(model);
configset.CompileOptions.UnitConversion = true;

Дополнительные сведения о создании моделей компартментальных ПК с использованием встроенной библиотеки SimBiology ® см. в разделе Создание фармакокинетических моделей.

Определение дозирования

Предположим, что каждый человек получает болюсную дозу 100 мг в момент времени = 0. Для получения подробной информации о настройке различных стратегий дозирования см. Дозы в моделях SimBiology.

dose             = sbiodose('dose','TargetName','Drug_Central');
dose.StartTime   = 0;
dose.Amount      = 100;
dose.AmountUnits = 'milligram';
dose.TimeUnits   = 'hour';

Сопоставление данных ответа с соответствующими компонентами модели

Данные содержат измеренную концентрацию плазмы в центральном и периферийном отсеках. Сопоставьте эти переменные с соответствующими компонентами модели, которые являются Drug_Central и Drug_Peripheral.

responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};

Укажите параметры для оценки

Указать объемы центрального и периферийного отсеков Central и Peripheral, межотраслевой клиренс Q12и клиренс Cl_Central в качестве параметров для оценки. estimatedInfo объект позволяет дополнительно указать преобразования параметров, начальные значения и границы параметров. Поскольку оба Central и Peripheral должны быть положительными, укажите логарифмическое преобразование для каждого параметра.

paramsToEstimate    = {'log(Central)', 'log(Peripheral)', 'Q12', 'Cl_Central'};
estimatedParam      = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);

Оценка индивидуальных параметров

Оцените один набор параметров для каждого человека, установив 'Pooled' аргумент пары имя-значение для false.

unpooledFit =  sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',false);

Показать результаты

Постройте график соответствия результатов исходным данным для каждого человека (группы).

t = plot(unpooledFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Для незаполненной посадки, sbiofit всегда возвращает один объект результатов для каждого отдельного объекта.

Анализ оценок параметров для зависимостей категорий

Изучите незаполненные оценки, чтобы выяснить, существуют ли какие-либо специфичные для категории параметры, то есть если некоторые параметры связаны с одной или несколькими категориями. При наличии каких-либо зависимостей категорий можно уменьшить число степеней свободы путем оценки только специфичных для категории значений для этих параметров.

Сначала извлеките значения идентификаторов и категорий для каждого идентификатора.

catParamValues = unique(gData(:,{'ID','Sex','Age'}));

Добавьте переменные в таблицу, содержащую оценку каждого параметра.

allParamValues            = vertcat(unpooledFit.ParameterEstimates);
catParamValues.Central    = allParamValues.Estimate(strcmp(allParamValues.Name, 'Central'));
catParamValues.Peripheral = allParamValues.Estimate(strcmp(allParamValues.Name, 'Peripheral'));
catParamValues.Q12        = allParamValues.Estimate(strcmp(allParamValues.Name, 'Q12'));
catParamValues.Cl_Central = allParamValues.Estimate(strcmp(allParamValues.Name, 'Cl_Central'));

Печать оценок каждого параметра для каждой категории. gscatter требуются Toolbox™ статистики и машинного обучения. Если он отсутствует, используйте другие альтернативные функции печати, такие как plot.

h           = figure;
ylabels     = {'Central','Peripheral','Cl\_Central','Q12'};
plotNumber  = 1;
for i = 1:4
    thisParam = estimatedParam(i).Name;
    
    % Plot for Sex category
    subplot(4,2,plotNumber);
    plotNumber  = plotNumber + 1;
    gscatter(double(catParamValues.Sex), catParamValues.(thisParam), catParamValues.Sex);
    ax          = gca;
    ax.XTick    = [];
    ylabel(ylabels(i));
    legend('Location','bestoutside')
    % Plot for Age category
    subplot(4,2,plotNumber);
    plotNumber  = plotNumber + 1;
    gscatter(double(catParamValues.Age), catParamValues.(thisParam), catParamValues.Age);
    ax          = gca;
    ax.XTick    = [];
    ylabel(ylabels(i));
    legend('Location','bestoutside')
end
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Исходя из сюжета, кажется, что молодые люди, как правило, имеют более высокие объемы центрального и периферического отделений (Central, Peripheral), чем старые люди (то есть объемы, по-видимому, зависят от возраста). Кроме того, мужчины, как правило, имеют более высокие показатели клиренса (Cl_Central), чем женщины, но наоборот для параметра Q12 (то есть клиренс и Q12, по-видимому, зависят от пола).

Параметры, специфичные для категории оценки

Используйте 'CategoryVariableName' имущества estimatedInfo для указания категории, используемой во время фитинга. Использовать 'Sex' как группа, подходящая для зазора Cl_Central и Q12 параметры. Использовать 'Age' как группа, подходящая для Central и Peripheral параметры.

estimatedParam(1).CategoryVariableName = 'Age';
estimatedParam(2).CategoryVariableName = 'Age';
estimatedParam(3).CategoryVariableName = 'Sex';
estimatedParam(4).CategoryVariableName = 'Sex';
categoryFit = sbiofit(model,gData,responseMap,estimatedParam,dose)

categoryFit = 
  OptimResults with properties:

                   ExitFlag: 3
                     Output: [1x1 struct]
                  GroupName: []
                       Beta: [8x5 table]
         ParameterEstimates: [120x6 table]
                          J: [240x8x2 double]
                       COVB: [8x8 double]
           CovarianceMatrix: [8x8 double]
                          R: [240x2 double]
                        MSE: 0.4362
                        SSE: 205.8690
                    Weights: []
              LogLikelihood: -477.9195
                        AIC: 971.8390
                        BIC: 1.0052e+03
                        DFE: 472
             DependentFiles: {1x3 cell}
    EstimatedParameterNames: {'Central'  'Peripheral'  'Q12'  'Cl_Central'}
             ErrorModelInfo: [1x3 table]
         EstimationFunction: 'lsqnonlin'

При подборе по категориям (или группам), sbiofit всегда возвращает один объект результатов, а не один для каждого уровня категории. Это связано с тем, что как мужчины, так и женщины считаются частью одной и той же оптимизации, используя одну и ту же модель ошибок и параметры ошибок, аналогично для молодых и старых людей.

Результаты графика

Постройте график оценочных результатов для конкретной категории.

t = plot(categoryFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Для Cl_Central и Q12 параметры, все мужчины имели одинаковые оценки, и аналогично для женщин. Для Central и Peripheral параметры, все молодые люди имели одинаковые оценки и аналогично для старых людей.

Оценка параметров для всего населения

Чтобы лучше сравнить результаты, подберите модель для всех данных, объединенных вместе, то есть оцените один набор параметров для всех пользователей, установив 'Pooled' аргумент пары имя-значение для true. Предупреждающее сообщение указывает на то, что эта опция игнорирует любую информацию, относящуюся к конкретной категории (если она существует).

pooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',true);

Warning: CategoryVariableName property of the estimatedInfo object is ignored when using the 'Pooled' option.

Результаты графика

Постройте график соответствия результатов исходным данным. Хотя для каждого индивидуума был сформирован отдельный график, данные подгонялись с использованием одного и того же набора параметров (то есть все индивидуумы имели одинаковую подгоняемую линию).

t = plot(pooledFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Сравнить остатки

Сравнить остатки CentralConc и PeripheralConc отклики для каждой посадки.

t = gData.Time;
allResid(:,:,1) = pooledFit.R;
allResid(:,:,2) = categoryFit.R;
allResid(:,:,3) = vertcat(unpooledFit.R);

h = figure;
responseList = {'CentralConc', 'PeripheralConc'};
for i = 1:2
    subplot(2,1,i);
    oneResid = squeeze(allResid(:,i,:));
    plot(t,oneResid,'o');
    refline(0,0); % A reference line representing a zero residual
    title(sprintf('Residuals (%s)', responseList{i}));
    xlabel('Time');
    ylabel('Residuals');
    legend({'Pooled','Category-Specific','Unpooled'});
end
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Как показано на графике, незаполненная посадка обеспечивает наилучшее соответствие данным, поскольку она соответствует данным каждого отдельного пользователя. Это ожидалось, поскольку в нем использовалось наибольшее количество степеней свободы. Подгонка по категориям уменьшила число степеней свободы за счет подгонки данных к двум категориям (пол и возраст). В результате остатки были больше, чем незаполненная посадка, но все же меньше, чем популяция, которая оценивала только один набор параметров для всех индивидуумов. Соответствие категориям может быть хорошим компромиссом между незаполненным и объединенным фитингом при условии, что в данных существует любая иерархическая модель.