nlmefit и nlmefitsaЗагрузите образцы данных.
load indomethacin
Данные в indomethacin.mat регистрирует концентрации препарата индометацина в кровотоке шести субъектов в течение восьми часов.
Постройте график рассеяния индометацина в кровотоке, сгруппированном по субъектам.
gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on
Specifying Mixed-Effects Models на странице рассматривается полезная модель для этого типа данных.
Создайте модель с помощью анонимной функции.
model = @(phi,t)(phi(1)*exp(-exp(phi(2))*t) + ...
phi(3)*exp(-exp(phi(4))*t));
Используйте nlinfit для соответствия модели всем данным, игнорируя эффекты, специфичные для субъекта.
phi0 = [1 2 1 1]; [phi,res] = nlinfit(time,concentration,model,phi0);
Вычислите среднеквадратичную ошибку.
numObs = length(time); numParams = 4; df = numObs-numParams; mse = (res'*res)/df
mse =
0.0304
Супер наложить модель на график рассеяния данных.
tplot = 0:0.01:8; plot(tplot,model(phi,tplot),'k','LineWidth',2) hold off
Нарисуйте рамочный график остатков по субъектам.
colors = 'rygcbm'; h = boxplot(res,subject,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(res,subject,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off
Прямоугольный график остатков по субъектам показывает, что прямоугольники в основном выше или ниже нуля, что указывает на то, что модель не учла специфичные для субъекта эффекты.
Чтобы учесть эффекты для конкретной темы, подгоняйте модель отдельно к данным для каждой темы.
phi0 = [1 2 1 1]; PHI = zeros(4,6); RES = zeros(11,6); for I = 1:6 tI = time(subject == I); cI = concentration(subject == I); [PHI(:,I),RES(:,I)] = nlinfit(tI,cI,model,phi0); end PHI
PHI =
2.0293 2.8277 5.4683 2.1981 3.5661 3.0023
0.5794 0.8013 1.7498 0.2423 1.0408 1.0882
0.1915 0.4989 1.6757 0.2545 0.2915 0.9685
-1.7878 -1.6354 -0.4122 -1.6026 -1.5069 -0.8731
Вычислите среднеквадратичную ошибку.
numParams = 24; df = numObs-numParams; mse = (RES(:)'*RES(:))/df
mse =
0.0057
Постройте график рассеяния данных и наложите модель на каждый объект.
gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on for I = 1:6 plot(tplot,model(PHI(:,I),tplot),'Color',colors(I)) end axis([0 8 0 3.5]) hold off
PHI дает оценки четырех параметров модели для каждого из шести предметов. Оценки значительно различаются, но взятая за 24-параметрическую модель данных, среднеквадратичная ошибка 0,0057 является значительным снижением от 0,0304 в исходной четырехпараметрической модели.
Нарисуйте рамочный график остатков по предметам.
h = boxplot(RES,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(RES,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off
Теперь рамочный график показывает, что большая модель составляет большую часть специфичных для субъекта эффектов. Разброс остатков (вертикальный масштаб бокса-графика) значительно меньше, чем на предыдущем боксе-графике, и боксы теперь в основном центрированы на ноль.
Хотя модель с 24 параметрами успешно учитывает вариации из-за конкретных субъектов в исследовании, она не рассматривает субъектов как представителей большей популяции. Распределение выборки, из которого взяты предметы, вероятно, более интересно, чем сама выборка. Цель моделей смешанных эффектов состоит в том, чтобы учитывать специфичные для субъекта вариации в более широком смысле, поскольку случайные эффекты, варьирующиеся вокруг популяционных средств.
Используйте nlmefit для подгонки модели смешанных эффектов к данным. Также можно использовать nlmefitsa вместо nlmefit .
Следующая анонимная функция, nlme_model , адаптирует четырехпараметрическую модель, используемую nlinfit к синтаксису вызова nlmefit путем предоставления отдельных параметров для каждого индивидуума. По умолчанию nlmefit присваивает случайные эффекты всем параметрам модели. Также по умолчанию, nlmefit предполагает диагональную ковариационную матрицу (без ковариации среди случайных эффектов), чтобы избежать гиперпараметризации и связанных с ней проблем сходимости.
nlme_model = @(PHI,t)(PHI(:,1).*exp(-exp(PHI(:,2)).*t) + ... PHI(:,3).*exp(-exp(PHI(:,4)).*t)); phi0 = [1 2 1 1]; [phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0)
phi =
2.8277
0.7729
0.4606
-1.3459
PSI =
0.3264 0 0 0
0 0.0250 0 0
0 0 0.0124 0
0 0 0 0.0000
stats =
struct with fields:
dfe: 57
logl: 54.5882
mse: 0.0066
rmse: 0.0787
errorparam: 0.0815
aic: -91.1765
bic: -93.0506
covb: [4x4 double]
sebeta: [0.2558 0.1066 0.1092 0.2244]
ires: [66x1 double]
pres: [66x1 double]
iwres: [66x1 double]
pwres: [66x1 double]
cwres: [66x1 double]
Среднеквадратичная ошибка 0,0066 сравнима с 0,0057 24-параметрической модели без случайных эффектов и значительно лучше, чем 0,0304 четырехпараметрической модели без случайных эффектов.
Оценочная ковариационная матрица PSI показывает, что дисперсия четвертого случайного эффекта по существу равна нулю, предполагая, что ее можно удалить, чтобы упростить модель. Для этого используйте 'REParamsSelect' пара «имя-значение» для указания индексов параметров, моделируемых случайными эффектами в nlmefit .
[phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3])
phi =
2.8277
0.7728
0.4605
-1.3460
PSI =
0.3270 0 0
0 0.0250 0
0 0 0.0124
stats =
struct with fields:
dfe: 58
logl: 54.5875
mse: 0.0066
rmse: 0.0780
errorparam: 0.0815
aic: -93.1750
bic: -94.8410
covb: [4x4 double]
sebeta: [0.2560 0.1066 0.1092 0.2244]
ires: [66x1 double]
pres: [66x1 double]
iwres: [66x1 double]
pwres: [66x1 double]
cwres: [66x1 double]
Логарифмическое правдоподобие logl почти идентичен тому, что было со случайными эффектами для всех параметров, информационный критерий Акайке aic сокращается с -91.1765 до -93.1750, а байесовский информационный критерий bic сокращается с -93.0506 до -94.8410. Эти меры поддерживают решение об отмене четвертого случайного эффекта.
Преобразование упрощенной модели с полной ковариационной матрицей позволяет идентифицировать корреляции между случайными эффектами. Для этого используйте CovPattern параметр для задания шаблона ненулевых элементов в ковариационной матрице.
[phi,PSI,stats] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3], ... 'CovPattern',ones(3))
phi =
2.8159
0.8263
0.5570
-1.1478
PSI =
0.4737 0.1145 0.0491
0.1145 0.0323 0.0029
0.0491 0.0029 0.0227
stats =
struct with fields:
dfe: 55
logl: 58.4419
mse: 0.0061
rmse: 0.0783
errorparam: 0.0781
aic: -94.8838
bic: -97.1744
covb: [4x4 double]
sebeta: [0.3018 0.1104 0.1170 0.1675]
ires: [66x1 double]
pres: [66x1 double]
iwres: [66x1 double]
pwres: [66x1 double]
cwres: [66x1 double]
Оценочная ковариационная матрица PSI показывает, что случайные эффекты на первых двух параметрах имеют относительно сильную корреляцию, и оба имеют относительно слабую корреляцию с последним случайным эффектом. Эта структура в ковариационной матрице более очевидна, если преобразовать PSI к корреляционной матрице с использованием corrcov .
RHO = corrcov(PSI) clf; imagesc(RHO) set(gca,'XTick',[1 2 3],'YTick',[1 2 3]) title('{\bf Random Effect Correlation}') h = colorbar; set(get(h,'YLabel'),'String','Correlation');
RHO =
1.0000 0.9264 0.4735
0.9264 1.0000 0.1070
0.4735 0.1070 1.0000
Включите эту структуру в модель, изменив спецификацию ковариационной картины на блок-диагональную.
P = [1 1 0;1 1 0;0 0 1] % Covariance pattern [phi,PSI,stats,b] = nlmefit(time,concentration,subject, ... [],nlme_model,phi0, ... 'REParamsSelect',[1 2 3], ... 'CovPattern',P)
P =
1 1 0
1 1 0
0 0 1
phi =
2.7830
0.8981
0.6581
-1.0000
PSI =
0.5180 0.1069 0
0.1069 0.0221 0
0 0 0.0454
stats =
struct with fields:
dfe: 57
logl: 58.0804
mse: 0.0061
rmse: 0.0768
errorparam: 0.0782
aic: -98.1608
bic: -100.0350
covb: [4x4 double]
sebeta: [0.3171 0.1073 0.1384 0.1453]
ires: [66x1 double]
pres: [66x1 double]
iwres: [66x1 double]
pwres: [66x1 double]
cwres: [66x1 double]
b =
-0.8507 -0.1563 1.0427 -0.7559 0.5652 0.1550
-0.1756 -0.0323 0.2152 -0.1560 0.1167 0.0320
-0.2756 0.0519 0.2620 0.1064 -0.2835 0.1389
Блок-диагональная ковариационная структура уменьшает aic от -94.9462 до -98.1608 и bic от -97.2368 до -100.0350 без существенного влияния на логарифмическое правдоподобие. Эти меры поддерживают ковариационную структуру, используемую в окончательной модели. Продукция b дает прогнозы трех случайных эффектов для каждого из шести субъектов. Они объединены с оценками фиксированных эффектов в phi для создания модели смешанных эффектов.
Постройте график модели смешанных эффектов для каждого из шести субъектов. Для сравнения также показана модель без случайных эффектов.
PHI = repmat(phi,1,6) + ... % Fixed effects [b(1,:);b(2,:);b(3,:);zeros(1,6)]; % Random effects RES = zeros(11,6); % Residuals colors = 'rygcbm'; for I = 1:6 fitted_model = @(t)(PHI(1,I)*exp(-exp(PHI(2,I))*t) + ... PHI(3,I)*exp(-exp(PHI(4,I))*t)); tI = time(subject == I); cI = concentration(subject == I); RES(:,I) = cI - fitted_model(tI); subplot(2,3,I) scatter(tI,cI,20,colors(I),'filled') hold on plot(tplot,fitted_model(tplot),'Color',colors(I)) plot(tplot,model(phi,tplot),'k') axis([0 8 0 3.5]) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') legend(num2str(I),'Subject','Fixed') end
Если очевидные отклонения в данных (видимые на предыдущих оконных графиках) игнорируются, обычный график вероятности остатков показывает разумное согласие с предположениями модели об ошибках.
clf; normplot(RES(:))
