Загрузите выборочные данные.
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(:))
