Модели смешанных эффектов Используя nlmefit и nlmefitsa

Этот пример подбирает модели смешанных эффектов, строит предсказания и остаточные значения, и интерпретирует результаты.

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

load indomethacin

Данные в indomethacin.mat концентрации записей метиндола препарата в кровотоке шести предметов более чем восемь часов.

Постройте график рассеивания метиндола в кровотоке, сгруппированном предметом.

gscatter(time,concentration,subject)
xlabel('Time (hours)')
ylabel('Concentration (mcg/ml)')
title('{\bf Indomethacin Elimination}')
hold on

Figure contains an axes object. The axes object with title blank I n d o m e t h a c i n blank E l i m i n a t i o n contains 6 objects of type line. These objects represent 1, 2, 3, 4, 5, 6.

Включая случайные эффекты в модели является эффективным, когда данные попадают в естественные группы. В этих данных группы являются просто индивидуумами при исследовании. Для получения дополнительной информации на моделях смешанных эффектов, которые составляют фиксированные эффекты и случайные эффекты, см. Модели Смешанных Эффектов.

Создайте модель с помощью анонимной функции.

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

Figure contains an axes object. The axes object with title blank I n d o m e t h a c i n blank E l i m i n a t i o n contains 7 objects of type line. These objects represent 1, 2, 3, 4, 5, 6.

Чертите диаграмму остаточных значений предметом.

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

Figure contains an axes object. The axes object contains 84 objects of type line.

Диаграмма остаточных значений предметом показывает, что поля в основном выше или ниже нуля, указывая, что модель перестала работать с учетом подчинено-специфичных эффектов.

С учетом подчинено-специфичных эффектов подбирайте модель отдельно к данным для каждого предмета.

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 = 4×6

    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

Figure contains an axes object. The axes object with title blank I n d o m e t h a c i n blank E l i m i n a t i o n contains 12 objects of type line. These objects represent 1, 2, 3, 4, 5, 6.

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

Figure contains an axes object. The axes object contains 84 objects of type line.

Теперь диаграмма показывает, что большая модель составляет большинство подчинено-специфичных эффектов. Распространение остаточных значений (вертикальная шкала диаграммы) намного меньше, чем в предыдущей диаграмме, и поля теперь в основном сосредоточены на нуле.

В то время как модель с 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 = 4×1

    2.8277
    0.7729
    0.4606
   -1.3459

PSI = 4×4

    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 = 4×1

    2.8277
    0.7728
    0.4605
   -1.3460

PSI = 3×3

    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 = 4×1

    2.8159
    0.8263
    0.5570
   -1.1478

PSI = 3×3

    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)
RHO = 3×3

    1.0000    0.9264    0.4735
    0.9264    1.0000    0.1070
    0.4735    0.1070    1.0000

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');

Figure contains an axes object. The axes object with title blank R a n d o m blank E f f e c t blank C o r r e l a t i o n contains an object of type image.

Включите эту структуру в модель путем изменения спецификации шаблона ковариации к диагонали блока.

P = [1 1 0;1 1 0;0 0 1] % Covariance pattern
P = 3×3

     1     1     0
     1     1     0
     0     0     1

[phi,PSI,stats,b] = nlmefit(time,concentration,subject, ...
                            [],nlme_model,phi0, ...
                            'REParamsSelect',[1 2 3], ...
                            'CovPattern',P)
phi = 4×1

    2.7830
    0.8981
    0.6581
   -1.0000

PSI = 3×3

    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 = 3×6

   -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

Figure contains 6 axes objects. Axes object 1 contains 3 objects of type scatter, line. These objects represent 1, Subject, Fixed. Axes object 2 contains 3 objects of type scatter, line. These objects represent 2, Subject, Fixed. Axes object 3 contains 3 objects of type scatter, line. These objects represent 3, Subject, Fixed. Axes object 4 contains 3 objects of type scatter, line. These objects represent 4, Subject, Fixed. Axes object 5 contains 3 objects of type scatter, line. These objects represent 5, Subject, Fixed. Axes object 6 contains 3 objects of type scatter, line. These objects represent 6, Subject, Fixed.

Если очевидные выбросы в данных (видимый в предыдущих диаграммах) проигнорированы, график нормального распределения остаточных значений показывает разумное соглашение с предположениями модели при ошибках.

figure
normplot(RES(:))

Figure contains an axes object. The axes object with title Normal Probability Plot contains 3 objects of type line.

Смотрите также

|

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте