Загрузите выборочные данные.

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