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