Выполнить нелинейную регрессию методом наименьших квадратов
Для этой функции рекомендуется использовать Toolbox™ статистики и машинного обучения, Toolbox™ оптимизации и инструментарий глобальной оптимизации.
fitResults = sbiofit(sm,grpData,responseMap,estiminfo)sm с использованием нелинейной регрессии методом наименьших квадратов.
grpData является groupedData object указание данных для подгонки. responseMap определяет сопоставление между компонентами модели и данными ответа в grpData. estimatedInfo является EstimatedInfo object который определяет оценочные параметры в модели sm. fitResults является OptimResults object или NLINResults object или вектор этих объектов.
sbiofit использует первую доступную функцию оценки среди следующих: lsqnonlin(Панель инструментов оптимизации), nlinfit (Панель инструментов для статистики и машинного обучения) или fminsearch.
По умолчанию каждая группа в grpData подгоняется отдельно, что приводит к оценкам параметров, специфичных для группы. Если модель содержит активные дозы и варианты, они применяются до моделирования.
fitResults = sbiofit(sm,grpData,responseMap,estiminfo,dosing)dosing вместо использования активных доз модели sm при наличии таковых.
fitResults = sbiofit(sm,grpData,responseMap,estiminfo,dosing,functionName)functionName. Если указанная функция недоступна, выдается предупреждение и используется первая доступная функция по умолчанию.
fitResults = sbiofit(sm,grpData,responseMap,estiminfo,dosing,functionName,options)options для функции functionName.
fitResults = sbiofit(sm,grpData,responseMap,estiminfo,dosing,functionName,options,variants)variants вместо использования любых активных вариантов модели.
fitResults = sbiofit(_,Name,Value)Name,Value аргументы пары.
[ также возвращает вектор объектов SimData fitResults,simdata] = sbiofit(_)simdata с использованием любого из входных аргументов в предыдущих синтаксисах.
Примечание
sbiofit объединяет sbionlinfit и sbioparamestim функции оценки. Использовать sbiofit для выполнения нелинейной регрессии методом наименьших квадратов.
sbiofit моделирует модель с помощью SimFunction object, что автоматически ускоряет моделирование по умолчанию. Следовательно, нет необходимости запускать sbioaccelerate перед звонком sbiofit.
Фон
Этот пример показывает, как подогнать данные профиля ПК индивидуума к однокамерной модели и оценить фармакокинетические параметры.
Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от человека и вы хотите оценить объем центрального отделения и клиренс. Предположим, что концентрация лекарственного средства в сравнении с профилем времени следует за моноэкспоненциальным снижением C0e-ket, Ct - концентрация лекарственного средства в момент tC0 - начальная концентрация, ke - постоянная скорости удаления, которая зависит от зазора и объема центрального Cl/V.
Синтетические данные в этом примере были получены с использованием следующей модели, параметров и дозы:
Однокамерная модель с болюсным дозированием и устранением первого порядка
Объем центрального отсека (Central) = 1,70 литра
Параметр зазора (Cl_Central) = 0,55 л/час
Модель постоянной ошибки
Болюсная доза 10 мг
Загрузка данных и визуализация
Данные хранятся в виде таблицы с переменными. Time и Conc которые представляют собой временной курс концентрации в плазме человека после внутривенного введения болюса, измеренный в 13 различных временных точках. Переменные единицы для Time и Conc час и миллиграмм/литр соответственно.
clear all load('data15.mat')
plot(data.Time,data.Conc,'b+') xlabel('Time (hour)'); ylabel('Drug Concentration (milligram/liter)');
Преобразуйте в Формат groupedData
Преобразовать набор данных в groupedData объект, который является требуемым форматом данных для функции фитинга sbiofit для последующего использования. A groupedData объект также позволяет задавать независимые имена переменных и групп (если они существуют). Установка единиц измерения Time и Conc переменные. Установки являются необязательными и обязательными только для UnitConversion функция, которая автоматически преобразует соответствующие физические величины в одну согласованную систему единиц измерения.
gData = groupedData(data);
gData.Properties.VariableUnits = {'hour','milligram/liter'};
gData.Propertiesans = struct with fields:
Description: ''
UserData: []
DimensionNames: {'Row' 'Variables'}
VariableNames: {'Time' 'Conc'}
VariableDescriptions: {}
VariableUnits: {'hour' 'milligram/liter'}
VariableContinuity: []
RowNames: {}
CustomProperties: [1x1 matlab.tabular.CustomProperties]
GroupVariableName: ''
IndependentVariableName: 'Time'
groupedData автоматически установить имя IndependentVariableName к свойству Time переменная данных.
Построение однокамерной модели
Используйте встроенную библиотеку ПК для построения однокамерной модели с болюсным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отсека. Используйте configset объект для включения преобразования единиц измерения.
pkmd = PKModelDesign; pkc1 = addCompartment(pkmd,'Central'); pkc1.DosingType = 'Bolus'; pkc1.EliminationType = 'linear-clearance'; pkc1.HasResponseVariable = true; model = construct(pkmd); configset = getconfigset(model); configset.CompileOptions.UnitConversion = true;
Дополнительные сведения о создании моделей компартментальных ПК с использованием встроенной библиотеки SimBiology ® см. в разделе Создание фармакокинетических моделей.
Определение дозирования
Определите единичную болюсную дозу в 10 миллиграммов, заданную в момент времени = 0. Для получения подробной информации о настройке различных графиков дозирования см. Дозы в моделях SimBiology.
dose = sbiodose('dose'); dose.TargetName = 'Drug_Central'; dose.StartTime = 0; dose.Amount = 10; dose.AmountUnits = 'milligram'; dose.TimeUnits = 'hour';
Сопоставление данных ответа с соответствующим компонентом модели
Данные содержат данные о концентрации лекарственного средства, хранящиеся в Conc переменная. Эти данные соответствуют Drug_Central виды в модели. Поэтому сопоставьте данные с Drug_Central следующим образом.
responseMap = {'Drug_Central = Conc'};Укажите параметры для оценки
Параметры для соответствия в этой модели - объем центрального отсека (Central) и скорость клиренса (Cl_Central). В этом случае укажите логарифмическое преобразование для этих биологических параметров, поскольку они ограничены положительными. estimatedInfo объект позволяет при необходимости задавать преобразования параметров, начальные значения и границы параметров.
paramsToEstimate = {'log(Central)','log(Cl_Central)'};
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1],'Bounds',[1 5;0.5 2]);Параметры оценки
Теперь, когда вы определили однокамерную модель, данные для подгонки, сопоставленные данные ответа, параметры для оценки и дозирование, использование sbiofit для оценки параметров. Функция оценки по умолчанию, которая sbiofit использование будет изменяться в зависимости от доступных панелей инструментов. Чтобы увидеть, какая функция использовалась во время фитинга, проверьте EstimationFunction свойства соответствующего объекта результатов.
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose);
Отображение расчетных параметров и результатов графика
Обратите внимание, что оценки параметров не были далеки от истинных значений (1,70 и 0,55), которые использовались для генерации данных. Можно также попробовать различные модели ошибок, чтобы увидеть, могут ли они еще больше улучшить оценки параметров.
fitConst.ParameterEstimates
ans=2×4 table
Name Estimate StandardError Bounds
______________ ________ _____________ __________
{'Central' } 1.6993 0.034821 1 5
{'Cl_Central'} 0.53358 0.01968 0.5 2
s.Labels.XLabel = 'Time (hour)'; s.Labels.YLabel = 'Concentration (milligram/liter)'; plot(fitConst,'AxesStyle',s);
Использовать различные модели ошибок
Попробуйте использовать три другие поддерживаемые модели ошибок (пропорциональную, комбинацию постоянной и пропорциональной моделей ошибок и экспоненциальную).
fitProp = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','proportional'); fitExp = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','exponential'); fitComb = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','combined');
Использовать веса вместо модели ошибки
Можно указать веса как числовую матрицу, где количество столбцов соответствует количеству ответов. Установка для всех весов значения 1 эквивалентна модели постоянной ошибки.
weightsNumeric = ones(size(gData.Conc));
fitWeightsNumeric = sbiofit(model,gData,responseMap,estimatedParams,dose,'Weights',weightsNumeric);Можно также использовать дескриптор функции, который принимает вектор прогнозируемых значений отклика и возвращает вектор весов. В этом примере используйте дескриптор функции, эквивалентный модели пропорциональной ошибки.
weightsFunction = @(y) 1./y.^2;
fitWeightsFunction = sbiofit(model,gData,responseMap,estimatedParams,dose,'Weights',weightsFunction);Сравнение информационных критериев для выбора модели
Сравните логические данные, значения AIC и BIC каждой модели, чтобы увидеть, какая модель ошибок лучше всего подходит для данных. Большее значение правдоподобия указывает, что соответствующая модель лучше подходит для модели. Для AIC и BIC меньшие значения лучше.
allResults = [fitConst,fitWeightsNumeric,fitWeightsFunction,fitProp,fitExp,fitComb];
errorModelNames = {'constant error model','equal weights','proportional weights', ...
'proportional error model','exponential error model',...
'combined error model'};
LogLikelihood = [allResults.LogLikelihood]';
AIC = [allResults.AIC]';
BIC = [allResults.BIC]';
t = table(LogLikelihood,AIC,BIC);
t.Properties.RowNames = errorModelNames;
tt=6×3 table
LogLikelihood AIC BIC
_____________ _______ _______
constant error model 3.9866 -3.9732 -2.8433
equal weights 3.9866 -3.9732 -2.8433
proportional weights -3.8472 11.694 12.824
proportional error model -3.8257 11.651 12.781
exponential error model 1.1984 1.6032 2.7331
combined error model 3.9163 -3.8326 -2.7027
На основе информационных критериев модель постоянных ошибок (или равных весов) лучше всего подходит к данным, поскольку она имеет наибольшую стоимость средств к существованию и наименьшую AIC и BIC.
Отображение оценочных значений параметров
Отображение расчетных значений параметров каждой модели.
Estimated_Central = zeros(6,1); Estimated_Cl_Central = zeros(6,1); t2 = table(Estimated_Central,Estimated_Cl_Central); t2.Properties.RowNames = errorModelNames; for i = 1:height(t2) t2{i,1} = allResults(i).ParameterEstimates.Estimate(1); t2{i,2} = allResults(i).ParameterEstimates.Estimate(2); end t2
t2=6×2 table
Estimated_Central Estimated_Cl_Central
_________________ ____________________
constant error model 1.6993 0.53358
equal weights 1.6993 0.53358
proportional weights 1.9045 0.51734
proportional error model 1.8777 0.51147
exponential error model 1.7872 0.51701
combined error model 1.7008 0.53271
Заключение
Этот пример показал, как оценить параметры ПК, а именно объем центрального отсека и параметр зазора индивидуума, путем подгонки данных профиля ПК к однокамерной модели. Вы сравнили информационные критерии каждой модели и оценочные значения параметров различных моделей ошибок, чтобы увидеть, какая модель лучше всего объяснила данные. Окончательные подогнанные результаты предлагали как постоянные, так и комбинированные модели ошибок, предоставляющие самые близкие оценки к значениям параметров, используемым для генерации данных. Тем не менее, модель постоянной ошибки является лучшей моделью, о чем свидетельствуют критерии источника информации, AIC и BIC.
Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от трех человек, которые вы хотите использовать для оценки соответствующих фармакокинетических параметров, а именно объема центрального и периферического отделения (Central, Peripheral), скорость клиренса (Cl_Central) и межотраслевой клиренс (Q12). Предположим, что концентрация лекарственного средства по сравнению с профилем времени следует за бикспоненциальным снижением Be − bt, где Ct - концентрация лекарственного средства в момент t, а a и b - наклоны для соответствующих экспоненциальных понижений.
Набор синтетических данных содержит данные концентрации лекарственного средства в плазме, измеренные как в центральном, так и в периферическом отсеках. Данные были получены с использованием двухкамерной модели с инфузионной дозой и элиминацией первого порядка. Эти параметры использовались для каждого человека.
Central | Peripheral | Q12 | Cl_Central | |
|---|---|---|---|---|
| Индивидуальный 1 | 1.90 | 0.68 | 0.24 | 0.57 |
| Индивидуальный 2 | 2.10 | 6.05 | 0.36 | 0.95 |
| Индивидуальный 3 | 1.70 | 4.21 | 0.46 | 0.95 |
Данные хранятся в виде таблицы с переменными. ID, Time, CentralConc, и PeripheralConc. Он представляет собой временной ход концентраций в плазме, измеренных в восемь различных моментов времени как для центрального, так и для периферического отделения после инфузионной дозы.
load('data10_32R.mat')
Преобразовать набор данных в groupedData объект, который является требуемым форматом данных для функции фитинга sbiofit для последующего использования. A groupedData объект также позволяет задавать независимые имена переменных и групп (если они существуют). Установка единиц измерения ID, Time, CentralConc, и PeripheralConc переменные. Установки являются необязательными и обязательными только для UnitConversion функция, которая автоматически преобразует соответствующие физические величины в одну согласованную систему единиц измерения.
gData = groupedData(data);
gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter'};
gData.Properties
ans =
struct with fields:
Description: ''
UserData: []
DimensionNames: {'Row' 'Variables'}
VariableNames: {'ID' 'Time' 'CentralConc' 'PeripheralConc'}
VariableDescriptions: {}
VariableUnits: {1x4 cell}
VariableContinuity: []
RowNames: {}
CustomProperties: [1x1 matlab.tabular.CustomProperties]
GroupVariableName: 'ID'
IndependentVariableName: 'Time'
Создайте решетчатый график, показывающий профили ПК трех отдельных пользователей.
sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},... 'Marker','+','LineStyle','none');
Используйте встроенную библиотеку ПК для построения двухкамерной модели с инфузионным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отсека. Используйте объект configset для включения преобразования единиц измерения.
pkmd = PKModelDesign; pkc1 = addCompartment(pkmd,'Central'); pkc1.DosingType = 'Infusion'; pkc1.EliminationType = 'linear-clearance'; pkc1.HasResponseVariable = true; pkc2 = addCompartment(pkmd,'Peripheral'); model = construct(pkmd); configset = getconfigset(model); configset.CompileOptions.UnitConversion = true;
Предположим, что каждый человек получает инфузионную дозу в момент времени = 0, с общим инфузионным количеством 100 мг при скорости 50 мг/час. Для получения подробной информации о настройке различных стратегий дозирования см. Дозы в моделях SimBiology.
dose = sbiodose('dose','TargetName','Drug_Central'); dose.StartTime = 0; dose.Amount = 100; dose.Rate = 50; dose.AmountUnits = 'milligram'; dose.TimeUnits = 'hour'; dose.RateUnits = 'milligram/hour';
Данные содержат измеренные концентрации в плазме в центральном и периферийном отсеках. Сопоставьте эти переменные с соответствующими модельными видами, которые являются Drug_Central и Drug_Peripheral.
responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};
Параметрами для оценки в этой модели являются объемы центрального и периферийного отсеков (Central и Peripheral), межотраслевой клиренс Q12и степень очистки Cl_Central. В этом случае укажите преобразование журнала для Central и Peripheral поскольку они ограничены положительными. estimatedInfo объект позволяет задать преобразования параметров, начальные значения и границы параметров (необязательно).
paramsToEstimate = {'log(Central)','log(Peripheral)','Q12','Cl_Central'};
estimatedParam = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);
Подберите модель для всех данных, объединенных вместе, то есть оцените один набор параметров для всех пользователей. Метод оценки по умолчанию, который sbiofit использование будет изменяться в зависимости от доступных панелей инструментов. Просмотр функции оценки sbiofit используется для фитинга, проверьте EstimationFunction свойства соответствующего объекта результатов.
pooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',true)
pooledFit =
OptimResults with properties:
ExitFlag: 3
Output: [1x1 struct]
GroupName: []
Beta: [4x3 table]
ParameterEstimates: [4x3 table]
J: [24x4x2 double]
COVB: [4x4 double]
CovarianceMatrix: [4x4 double]
R: [24x2 double]
MSE: 6.6220
SSE: 291.3688
Weights: []
LogLikelihood: -111.3904
AIC: 230.7808
BIC: 238.2656
DFE: 44
DependentFiles: {1x3 cell}
EstimatedParameterNames: {'Central' 'Peripheral' 'Q12' 'Cl_Central'}
ErrorModelInfo: [1x3 table]
EstimationFunction: 'lsqnonlin'
Постройте график соответствия результатов исходным данным. Хотя были созданы три отдельных графика, данные были подстроены с использованием одного и того же набора параметров (то есть все три человека имели одинаковую подгоняемую линию).
plot(pooledFit);
Оцените один набор параметров для каждого индивидуума и проверьте, есть ли какие-либо улучшения в оценках параметров. В этом примере, поскольку существует три индивидуума, оцениваются три набора параметров.
unpooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',false);
Постройте график соответствия результатов исходным данным. Каждый индивидуум был подогнан по-разному (то есть каждая подогнанная линия уникальна для каждого индивидуума), и каждая строка, по-видимому, хорошо соответствовала отдельным данным.
plot(unpooledFit);
Просмотрите соответствующие результаты первого отдельного пользователя. MSE был ниже, чем у объединенной посадки. Это также справедливо для двух других лиц.
unpooledFit(1)
ans =
OptimResults with properties:
ExitFlag: 3
Output: [1x1 struct]
GroupName: 1
Beta: [4x3 table]
ParameterEstimates: [4x3 table]
J: [8x4x2 double]
COVB: [4x4 double]
CovarianceMatrix: [4x4 double]
R: [8x2 double]
MSE: 2.1380
SSE: 25.6559
Weights: []
LogLikelihood: -26.4805
AIC: 60.9610
BIC: 64.0514
DFE: 12
DependentFiles: {1x3 cell}
EstimatedParameterNames: {'Central' 'Peripheral' 'Q12' 'Cl_Central'}
ErrorModelInfo: [1x3 table]
EstimationFunction: 'lsqnonlin'
Создайте график остатков с течением времени для сравнения объединенных и незаполненных результатов аппроксимации. Рисунок показывает, что неохлажденные остатки посадки меньше, чем ожидаемые остатки объединенной посадки. Помимо сравнения остатков, для сравнения соответствующих результатов можно использовать другие строгие критерии.
t = [gData.Time;gData.Time]; res_pooled = vertcat(pooledFit.R); res_pooled = res_pooled(:); res_unpooled = vertcat(unpooledFit.R); res_unpooled = res_unpooled(:); plot(t,res_pooled,'o','MarkerFaceColor','w','markerEdgeColor','b') hold on plot(t,res_unpooled,'o','MarkerFaceColor','b','markerEdgeColor','b') refl = refline(0,0); % A reference line representing a zero residual title('Residuals versus Time'); xlabel('Time'); ylabel('Residuals'); legend({'Pooled','Unpooled'});
В этом примере показано, как выполнять объединенные и незаполненные оценки с использованием sbiofit. Как показано, незаполненная подгонка учитывает вариации из-за конкретных субъектов в исследовании, и в этом случае модель лучше подходит к данным. Однако объединенная подгонка возвращает параметры для всего населения. Если требуется оценить параметры для всего населения с учетом отдельных вариаций, используйте sbiofitmixed.
В этом примере показано, как оценить специфичные для категории (такие как молодой и старый, мужской и женский), индивидуальные и популяционные параметры с использованием данных профиля ПК от нескольких человек.
Фон
Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от 30 человек и вы хотите оценить фармакокинетические параметры, а именно объемы центрального и периферического отделения, клиренс и межкомпартментальный клиренс. Предположим, что концентрация лекарственного средства по сравнению с профилем времени следует за бикспоненциальным снижением 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.Propertiesans = 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];
Как показано на графике, незаполненная посадка обеспечивает наилучшее соответствие данным, поскольку она соответствует данным каждого отдельного пользователя. Это ожидалось, поскольку в нем использовалось наибольшее количество степеней свободы. Подгонка по категориям уменьшила число степеней свободы за счет подгонки данных к двум категориям (пол и возраст). В результате остатки были больше, чем незаполненная посадка, но все же меньше, чем популяция, которая оценивала только один набор параметров для всех индивидуумов. Соответствие категориям может быть хорошим компромиссом между незаполненным и объединенным фитингом при условии, что в данных существует любая иерархическая модель.
В этом примере используется модель дрожжевого гетеротримерного G-белка и экспериментальные данные, представленные в [1]. Дополнительные сведения о модели см. в разделе Фон в разделе Сканирование параметров, Оценка параметров и Анализ чувствительности в цикле дрожжевого гетеротримерного G-белка.
Загрузите модель G-белка.
sbioloadproject gproteinХранят экспериментальные данные, содержащие временной курс для фракции активного G-белка.
time = [0 10 30 60 110 210 300 450 600]'; GaFracExpt = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';
Создать groupedData объект на основе экспериментальных данных.
tbl = table(time,GaFracExpt); grpData = groupedData(tbl);
Сопоставьте соответствующий компонент модели с экспериментальными данными. Другими словами, укажите, какие виды в модели соответствуют какой переменной ответа в данных. В этом примере сопоставьте параметр модели GaFrac в экспериментальную переменную данных GaFracExpt от grpData.
responseMap = 'GaFrac = GaFracExpt';Использовать estimatedInfo объект для определения параметра модели kGd в качестве параметра, подлежащего оценке.
estimatedParam = estimatedInfo('kGd');Выполните оценку параметров.
fitResult = sbiofit(m1,grpData,responseMap,estimatedParam);
Просмотр расчетного значения параметра kGd.
fitResult.ParameterEstimates
ans=1×3 table
Name Estimate StandardError
_______ ________ _____________
{'kGd'} 0.11307 3.4439e-05
Предположим, что необходимо построить график результатов моделирования модели с использованием расчетного значения параметра. Можно либо повторно запустить sbiofit и укажите, чтобы вернуть дополнительный второй выходной аргумент, который содержит результаты моделирования, или используйте fitted метод извлечения результатов без повторного выполнения sbiofit.
[yfit,paramEstim] = fitted(fitResult);
Постройте график результатов моделирования.
sbioplot(yfit);
Этот пример показывает, как оценить временную задержку до введения болюсной дозы и продолжительность дозы с использованием однокамерной модели.
Загрузите образец набора данных.
load lagDurationData.matПостройте график данных.
plot(data.Time,data.Conc,'x') xlabel('Time (hour)') ylabel('Conc (milligram/liter)')
Преобразуйте в groupedData.
gData = groupedData(data);
gData.Properties.VariableUnits = {'hour','milligram/liter'};Создайте однокамерную модель.
pkmd = PKModelDesign; pkc1 = addCompartment(pkmd,'Central'); pkc1.DosingType = 'Bolus'; pkc1.EliminationType = 'linear-clearance'; pkc1.HasResponseVariable = true; model = construct(pkmd); configset = getconfigset(model); configset.CompileOptions.UnitConversion = true;
Добавьте два параметра, которые представляют временной лаг и продолжительность дозы. Параметр запаздывания определяет время запаздывания до введения дозы. Параметр длительности определяет продолжительность времени, которое требуется для введения дозы.
lagP = addparameter(model,'lagP'); lagP.ValueUnits = 'hour'; durP = addparameter(model,'durP'); durP.ValueUnits = 'hour';
Создайте объект дозы. Установите LagParameterName и DurationParameterName свойства дозы к названиям параметров запаздывания и длительности соответственно. Установите количество дозы в 10 миллиграмм, которое было количеством, использованным для генерации данных.
dose = sbiodose('dose'); dose.TargetName = 'Drug_Central'; dose.StartTime = 0; dose.Amount = 10; dose.AmountUnits = 'milligram'; dose.TimeUnits = 'hour'; dose.LagParameterName = 'lagP'; dose.DurationParameterName = 'durP';
Сопоставьте виды модели с соответствующими данными.
responseMap = {'Drug_Central = Conc'};Укажите параметры задержки и длительности в качестве параметров для оценки. Выполните преобразование параметров. Инициализируйте их как 2 и установите верхнюю и нижнюю границы.
paramsToEstimate = {'log(lagP)','log(durP)'};
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',2,'Bounds',[1 5]);Выполните оценку параметров.
fitResults = sbiofit(model,gData,responseMap,estimatedParams,dose,'fminsearch')fitResults =
OptimResults with properties:
ExitFlag: 1
Output: [1x1 struct]
GroupName: One group
Beta: [2x4 table]
ParameterEstimates: [2x4 table]
J: [11x2 double]
COVB: [2x2 double]
CovarianceMatrix: [2x2 double]
R: [11x1 double]
MSE: 0.0024
SSE: 0.0213
Weights: []
LogLikelihood: 18.7511
AIC: -33.5023
BIC: -32.7065
DFE: 9
DependentFiles: {1x2 cell}
EstimatedParameterNames: {'lagP' 'durP'}
ErrorModelInfo: [1x3 table]
EstimationFunction: 'fminsearch'
Просмотрите результат.
fitResults.ParameterEstimates
ans=2×4 table
Name Estimate StandardError Bounds
________ ________ _____________ ______
{'lagP'} 1.986 0.0051568 1 5
{'durP'} 1.527 0.012956 1 5
plot(fitResults)
Yi, T-M., Kitano, H. и Simon, M. (2003). Количественная характеристика дрожжевого гетеротримерного G-белкового цикла. PNAS. 100, 10764–10769.
[2] Gábor, A., and Banga, J.R. (2015). Надежная и эффективная оценка параметров в динамических моделях биологических систем. BMC системная биология. 9:74.
EstimatedInfo object | fminsearch | groupedData object | LeastSquaresResults object | NLINResults object | OptimResults object | sbiofitmixed | ga (инструментарий глобальной оптимизации) | particleswarm(Панель инструментов глобальной оптимизации) | patternsearch(Панель инструментов глобальной оптимизации) | fmincon(Панель инструментов оптимизации) | fminunc(Панель инструментов оптимизации) | lsqcurvefit(Панель инструментов оптимизации) | lsqnonlin(Панель инструментов оптимизации) | nlinfit(Набор инструментов для статистики и машинного обучения)