Продольный анализ

В этом примере показано, как выполнить продольный анализ с помощью mvregress.

Загрузите выборочные данные.
Постройте графики данных.
Задайте матрицы проекта.
Подгонка модели.
Постройте график подобранной модели.
Задайте сокращенную модель.
Подгонка под уменьшенную модель.
Проведите тест коэффициента вероятности.

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

Загрузите продольные данные выборки.

load longitudinalData

Матрица Y содержит данные отклика для 16 индивидуумов. Ответ представляет собой уровень препарата в крови, измеренный в пяти временных точках (t = 0, 2, 4, 6 и 8). Каждая строка Y соответствует индивидууму, и каждый столбец соответствует временной точке. Первые восемь субъектов - женщины, а вторые восемь - мужчины. Это моделируемые данные.

Постройте графики данных.

Постройте график данных для всех 16 субъектов.

figure()
t = [0,2,4,6,8];
plot(t,Y)
hold on

hf = plot(t,Y(1:8,:),'^');
hm = plot(t,Y(9:16,:),'o');
legend([hf(1),hm(1)],'Female','Male','Location','NorthEast')

title('Longitudinal Response')
ylabel('Blood Drug Level')
xlabel('Time')
hold off

Задайте матрицы проекта.

Позвольте _yij обозначить ответ для индивидуума i = 1,..., n измеренный в моменты времени _tij, j = 1,..., d. В этом примере n = 16 и d = 5. Позвольте Gi обозначить пол отдельных i, где Gi = 1 для самцов и 0 для самок.

Рассмотрите подбор кривой квадратичной продольной модели с отдельным наклоном и точкой пересечения для каждого пола,

$y_{i j} = β_{0} + β_{1} G_{i} + β_{2} t_{i j} + β_{3} t_{i j}^{2} + β_{4} G_{i} \times t_{i j} + β_{5} G_{i} \times t_{i j}^{2} + ε_{i j},$

где $ε_{i} = (ε_{i 1}, \dots, ε_{i d})^{'} \sim M V N (0, Σ)$ . Корреляция ошибок учитывает кластеризацию внутри индивидуума.

Чтобы подогнать эту модель используя mvregress, данные отклика должны быть в n -by d матрице. Y уже в соответствующем формате.

Затем создайте массив ячеек n длины из матриц d -by K design. Для этой модели существует K = 6 параметров.

Для отдельных i матрица проекта 5 на 6 является

$X {i} = (\begin{matrix} 1 & G_{i} & t_{i 1} & t_{i 1}^{2} & G_{i} \times t_{i 1} & G_{i} \times t_{i 1}^{2} \\ 1 & G_{i} & t_{i 2} & t_{i 2}^{2} & G_{i} \times t_{i 2} & G_{i} \times t_{i 2}^{2} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & G_{i} & t_{i 5} & t_{i 5}^{2} & G_{i} \times t_{i 5} & G_{i} \times t_{i 5}^{2} \end{matrix}),$

соответствующий вектору параметра

$β = (\begin{array}{l} β_{0} \\ β_{1} \\ ⋮ \\ β_{5} \end{array}) .$

Матрица X1 имеет матрицу проекта для самки, и X2 имеет матрицу проекта для самца.

Создайте массив ячеек из матриц проекта. Первые восемь индивидуумов - самки, а вторые восемь - самцы.

X = cell(8,1);
X(1:8) = {X1};
X(9:16) = {X2};

Подгонка модели.

Подгонка модели с помощью максимальной оценки правдоподобия. Отобразите оцененные коэффициенты и стандартные ошибки.

[b,sig,E,V,loglikF] = mvregress(X,Y);
[b sqrt(diag(V))]

ans =

   18.8619    0.7432
   13.0942    1.0511
    2.5968    0.2845
   -0.3771    0.0398
   -0.5929    0.4023
    0.0290    0.0563

Коэффициенты на условиях взаимодействия (в последних двух строках b) не появляются значительными. Можно использовать значение целевой функции логарифмической правдоподобности для этой подгонки, loglikF, чтобы сравнить эту модель с моделью без условий взаимодействия с помощью теста коэффициента правдоподобия.

Постройте график подобранной модели.

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

Yhatf = X1*b;
Yhatm = X2*b;

figure()
plot(t,Y)
hold on

plot(t,Y(1:8,:),'^',t,Y(9:16,:),'o')

hf = plot(t,Yhatf,'k--','LineWidth',3);
hm = plot(t,Yhatm,'k','LineWidth',3);
legend([hf,hm],'Females','Males','Location','NorthEast')

title('Longitudinal Response')
ylabel('Blood Drug Level')
xlabel('Time')
hold off

Задайте сокращенную модель.

Подгонка модели без условий взаимодействия,

$y_{i j} = β_{0} + β_{1} G_{i} + β_{2} t_{i j} + β_{3} t_{i j}^{2} + ε_{i j},$

где $ε_{i} = (ε_{i 1}, \dots, ε_{i d})^{'} \sim M V N (0, Σ)$ .

Эта модель имеет четыре коэффициента, которые соответствуют первым четырем столбцам матриц проекта X1 и X2 (для женщин и мужчин, соответственно).

X1R = X1(:,1:4);
X2R = X2(:,1:4);

XR = cell(8,1);
XR(1:8) = {X1R};
XR(9:16) = {X2R};

Подгонка под уменьшенную модель.

Подгонка этой модели с помощью максимальной оценки правдоподобия. Отобразите оцененные коэффициенты и их стандартные ошибки.

[bR,sigR,ER,VR,loglikR] = mvregress(XR,Y);
[bR,sqrt(diag(VR))]

ans =

   19.3765    0.6898
   12.0936    0.8591
    2.2919    0.2139
   -0.3623    0.0283

Проведите тест коэффициента вероятности.

Сравните две модели с помощью теста коэффициента правдоподобия. Нулевая гипотеза заключается в том, что редуцированная модель достаточна. Альтернатива состоит в том, что уменьшенная модель является неадекватной (по сравнению с полной моделью с терминами взаимодействия).

Тестовая статистическая величина коэффициента правдоподобия сравнивается с хи-квадратным распределением с двумя степенями свободы (для двух понижаемых коэффициентов).

LR = 2*(loglikF-loglikR);
pval = 1 - chi2cdf(LR,2)

pval =

    0.0803

Значение p 0.0803 указывает, что гипотеза null не отклоняется на уровне 5% значимости. Поэтому нет достаточных доказательств того, что дополнительные условия улучшают подгонку.

См. также

mvregress | mvregresslike

Документация