Модели смешанных эффектов Используя 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 почти идентично тому, чем это было со случайными эффектами для всех параметров, критерий информации о Akaike 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.8149
0.8293
0.5613
-1.1407
PSI =
0.4768 0.1152 0.0500
0.1152 0.0321 0.0032
0.0500 0.0032 0.0236
stats =
struct with fields:
dfe: 55
logl: 58.4732
mse: 0.0061
rmse: 0.0782
errorparam: 0.0781
aic: -94.9464
bic: -97.2371
covb: [4x4 double]
sebeta: [0.3028 0.1104 0.1179 0.1662]
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.9316 0.4708
0.9316 1.0000 0.1180
0.4708 0.1180 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(:))
