Выполните регрессию нелинейного метода наименьших квадратов
Statistics and Machine Learning Toolbox™, Optimization Toolbox™ и Global Optimization Toolbox рекомендуются для этой функции.
fitResults = sbiofit(sm,grpData,ResponseMap,estiminfo)sm использование регрессии нелинейного метода наименьших квадратов.
grpData isa groupedData object определение данных, чтобы соответствовать. ResponseMap задает отображение между компонентами модели и данными об ответе в grpData. estimatedInfo EstimatedInfo object это задает предполагаемые параметры в модели sm. fitResults isa OptimResults object или NLINResults object или вектор из этих объектов.
sbiofit использует первую доступную функцию оценки среди следующего: lsqnonlin (Optimization Toolbox), nlinfit (Statistics and Machine Learning Toolbox), или 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 симулирует модель с помощью a SimFunction object, который автоматически ускоряет симуляции по умолчанию. Следовательно не необходимо запуститься sbioaccelerate прежде чем вы вызовете sbiofit.
Фон
В этом примере показано, как соответствовать PK индивидуума профилируют данные к модели с одним отсеком и оценить фармакокинетические параметры.
Предположим, что вы имеете данные о концентрации плазмы препарата от индивидуума и хотите оценить объем центрального отсека и разрешения. Примите, что концентрация препарата по сравнению с профилем времени следует за моноэкспоненциальным снижением , где концентрация препарата во время t, начальная концентрация, и константа скорости устранения, которая зависит от разрешения и объема центрального отсека .
Синтетические данные в этом примере были сгенерированы с помощью следующей модели, параметров и дозы:
Модель с одним отсеком с болюсным введением и устранением первого порядка
Объем центрального отсека (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 для дальнейшего использования. 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 переменная данных.
Создайте модель с одним отсеком
Пользуйтесь встроенной библиотекой PK, чтобы создать модель с одним отсеком с болюсным введением и устранением первого порядка, где уровень устранения зависит от разрешения и объема центрального отсека. Используйте 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;
Для получения дополнительной информации при создании разделенных на отсеки моделей PK с помощью 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'};Задайте параметры, чтобы оценить
Параметры, чтобы поместиться в эту модель являются объемом центрального (Центрального) отсека и уровень раскрываемости преступлений (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
Заключение
Этот пример показал, как оценить параметры PK, а именно, объем центрального отсека и параметр разрешения индивидуума, путем подбора кривой данным о профиле PK к модели с одним отсеком. Вы сравнили информационные критерии каждой модели и оценили, что значения параметров различных ошибочных моделей видели, какая модель лучше всего объяснила данные. Финал соответствовал результатам, предложенным обоих, которых постоянные и объединенные ошибочные модели предоставили, самые близкие оценки к значениям параметров раньше генерировали данные. Однако постоянная ошибочная модель является лучшей моделью, как обозначено логарифмической правдоподобностью, AIC и критериями информации о BIC.
Предположим, что у вас есть данные о концентрации плазмы препарата от трех индивидуумов, которых вы хотите использовать, чтобы оценить соответствующие фармакокинетические параметры, а именно, объем центрального и периферийного отсека (Central, Peripheral), уровень раскрываемости преступлений (Cl_Central), и межразделенное на отсеки разрешение (Q12). Примите, что концентрация препарата по сравнению с профилем времени следует за снижением biexponential , где Ct является концентрацией препарата во время t, и a и b являются наклонами для соответствующих экспоненциальных снижений.
Синтетический набор данных содержит данные о концентрации плазмы препарата, измеренные и в центральных и в периферийных отсеках. Данные были сгенерированы с помощью модели 2D отсека с капельным внутривенным введением и устранением первого порядка. Эти параметры использовались для каждого индивидуума.
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Время, CentralConc, и PeripheralConc. Это представляет ход времени плазменных концентраций, измеренных в восьми различных моментах времени и для центральных и для периферийных отсеков после капельного внутривенного введения.
load('data10_32R.mat')
Преобразуйте набор данных в groupedData объект, который является необходимым форматом данных для подходящей функции sbiofit для дальнейшего использования. groupedData объект также позволяет вам установить независимую переменную и имена переменных группы (если они существуют). Установите модули IDВремя, 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'
Создайте график решетки, который показывает профили PK трех индивидуумов.
sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},... 'Marker','+','LineStyle','none');
Пользуйтесь встроенной библиотекой PK, чтобы создать модель 2D отсека с дозированием вливания и устранением первого порядка, где уровень устранения зависит от разрешения и объема центрального отсека. Используйте объект 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.
В этом примере показано, как оценить специфичный для категории (такой как молодой по сравнению со старым, штекерным по сравнению с розеткой), отдельно-специфичные, и параметры всего населения с помощью PK профилируют данные от нескольких индивидуумов.
Фон
Предположим, что вы имеете данные о концентрации плазмы препарата от 30 индивидуумов и хотите оценить фармакокинетические параметры, а именно, объемы центрального и периферийного отсека, разрешения и межразделенного на отсеки разрешения. Примите, что концентрация препарата по сравнению с профилем времени следует за снижением biexponential , где концентрация препарата во время t, и и наклоны для соответствующих экспоненциальных снижений.
Загрузка данных
Эти синтетические данные содержат ход времени плазменных концентраций 30 индивидуумов после дозы шарика (100 мг), измеренных в разное время и для центральных и для периферийных отсеков. Это также содержит категориальные переменные, а именно, Sex и Age.
clear
load('sd5_302RAgeSex.mat')Преобразуйте в groupedData Формат
Преобразуйте набор данных в groupedData объект, который является необходимым форматом данных для подходящего функционального sbiofit. groupedData объект также позволяет вам независимую переменную набора и имена переменных группы (если они существуют). Установите модули IDВремя, 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];
Настройте модель 2D отсека
Пользуйтесь встроенной библиотекой PK, чтобы создать модель 2D отсека с дозированием вливания и устранением первого порядка, где уровень устранения зависит от разрешения и объема центрального отсека. Используйте 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;
Для получения дополнительной информации при создании разделенных на отсеки моделей PK с помощью 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);Отображение результатов
Постройте подходящие результаты по сравнению с исходными данными для каждого индивидуума (группа).
plot(unpooledFit);
Для необъединенной подгонки, sbiofit всегда возвращается, каждый заканчивается объект для каждого индивидуума.
Исследуйте оценки параметра на зависимости от категории
Исследуйте необъединенные оценки, чтобы видеть, существуют ли какие-либо специфичные для категории параметры, то есть, если некоторые параметры связаны с одной или несколькими категориями. Если существуют какие-либо зависимости от категории, может быть возможно уменьшать количество степеней свободы путем оценки только специфичных для категории значений для тех параметров.
Сначала извлеките ID и значения категории для каждого ID
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 требует Statistics and Machine Learning 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 всегда возвращается, каждый заканчивается объект, не один для каждого уровня категории. Это вызвано тем, что и штекерные и розеточные индивидуумы считаются частью той же оптимизации с помощью тех же ошибочных параметров модели и параметров ошибок, так же для молодых и старых индивидуумов.
Постройте результаты
Постройте специфичные для категории предполагаемые результаты.
plot(categoryFit);
Для 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.
Постройте результаты
Постройте подходящие результаты по сравнению с исходными данными. Несмотря на то, что отдельный график был сгенерирован для каждого индивидуума, данные были адаптированы с помощью того же набора параметров (то есть, у всех индивидуумов была та же подходящая линия).
plot(pooledFit);
Сравните остаточные значения
Сравните остаточные значения 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]. Для получения дополнительной информации о модели, смотрите раздел Background в Сканировании Параметра, Оценке Параметра и Анализе чувствительности в Дрожжах Гетеротримерный Цикл Белка 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)
[1] И, T-M., Kitano, H. и Саймон, M. (2003). Количественная характеристика дрожжей гетеротримерный цикл белка G. PNAS. 100, 10764–10769.
[2] Габор, A., и Банга, J.R. (2015). Устойчивая и эффективная оценка параметра динамические модели биологических систем. Системная биология BMC. 9:74.
EstimatedInfo object | groupedData object | LeastSquaresResults object | NLINResults object | OptimResults object | sbiofitmixed | nlinfit (Statistics and Machine Learning Toolbox) | fmincon (Optimization Toolbox) | fminunc (Optimization Toolbox) | fminsearch | lsqcurvefit (Optimization Toolbox) | lsqnonlin (Optimization Toolbox) | patternsearch (Global Optimization Toolbox) | ga (Global Optimization Toolbox) | particleswarm (Global Optimization Toolbox)