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(:))
