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

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