Дисперсионный Анализ со случайными эффектами

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

В обычной модели ANOVA каждая сгруппированная переменная представляет фиксированный коэффициент. Уровни этого фактора являются фиксированным множеством значений. Цель состоит в том, чтобы определить, приводят ли различные уровни факторов к различным значениям отклика.

Настройте модель

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

load mileage

The anova2 функция работает только с сбалансированными данными, и она выводит значения сгруппированных переменных из номеров строка и столбец входа матрицы. The anovan функция, с другой стороны, требует, чтобы вы явным образом создали векторы с сгруппированной переменной значениями. Создайте эти векторы следующим образом.

Создайте массив, указывающий фабрику для каждого значения в километре. Этот массив равен 1 для первого столбца, 2 для второго и 3 для третьего.

factory  = repmat(1:3,6,1);

Создайте массив, указывающий модель автомобиля для каждого значения пробега. Этот массив равен 1 для первых трех строк пробега и 2 для остальных трех строк.

carmod = [ones(3,3); 2*ones(3,3)];

Превратите эти матрицы в векторы и отобразите их.

mileage = mileage(:);
factory = factory(:);
carmod = carmod(:);
[mileage factory carmod]
ans = 18×3

   33.3000    1.0000    1.0000
   33.4000    1.0000    1.0000
   32.9000    1.0000    1.0000
   32.6000    1.0000    2.0000
   32.5000    1.0000    2.0000
   33.0000    1.0000    2.0000
   34.5000    2.0000    1.0000
   34.8000    2.0000    1.0000
   33.8000    2.0000    1.0000
   33.4000    2.0000    2.0000
      ⋮

Подбор модели случайных эффектов

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

[pvals,tbl,stats] = anovan(mileage, {factory carmod}, ... 
'model',2, 'random',1,'varnames',{'Factory' 'Car Model'});

Figure N-Way ANOVA contains objects of type uicontrol.

В версии фиксированных эффектов этой подгонки, которую вы получаете, опуская входы 'random',1 в предыдущем коде эффект модели автомобиля значителен, с p-значением 0,0039. Но в этом примере, который учитывает случайное изменение эффекта переменной 'Car Model' от одной фабрики к другой, эффект все еще значителен, но с более высоким значением p 0,0136.

F-статистика для моделей со случайными эффектами

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

В примере, описанном в Set Up the Model, эффект переменной 'Factory' может варьироваться в зависимости от моделей автомобилей. В этом случае средний квадрат взаимодействия занимает место среднего квадрата ошибки в F-статистике.

Найдите F-статистику.

F = 26.6756 / 0.02
F = 1.3338e+03

Степени свободы для статистической функции являются степенями свободы для числителя (2) и знаменателя (2) средних квадратов.

Найдите значение p.

pval = 1 - fcdf(F,2,2)
pval = 7.4919e-04

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

Найдите коэффициенты этих линейных комбинаций.

stats.ems
ans = 4×4

    6.0000    0.0000    3.0000    1.0000
    0.0000    9.0000    3.0000    1.0000
    0.0000    0.0000    3.0000    1.0000
         0         0         0    1.0000

Это возвращает ems поле stats структура.

Отображение текстовых представлений линейных комбинаций.

stats.txtems
ans = 4x1 cell
    {'6*V(Factory)+3*V(Factory*Car Model)+V(Error)'  }
    {'9*Q(Car Model)+3*V(Factory*Car Model)+V(Error)'}
    {'3*V(Factory*Car Model)+V(Error)'               }
    {'V(Error)'                                      }

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

В некоторых случаях не существует ни одного термина, ожидаемое значение которого совпадало бы с тем, которое требуется для знаменателя F-статистики. В этом случае знаменатель является линейной комбинацией средних квадратов. Структура stats содержит поля, дающие определения знаменателей для каждой F-статистики. The txtdenom поле, stats.txtdenom, содержит текстовое представление и denom поле содержит матрицу, которая задает линейную комбинацию отклонений членов в модели. Для сбалансированных моделей, таких как эта, denom матрица, stats.denom, содержит нули и таковые, потому что знаменатель - это всего один средний квадрат термина.

Отобразите txtdenom поле.

stats.txtdenom
ans = 3x1 cell
    {'MS(Factory*Car Model)'}
    {'MS(Factory*Car Model)'}
    {'MS(Error)'            }

Отобразите denom поле.

stats.denom
ans = 3×3

    0.0000    1.0000   -0.0000
    0.0000    1.0000   -0.0000
    0.0000   -0.0000    1.0000

Отклонение компоненты

Для модели, описанной в Set Up the Model, рассмотрите пробег для конкретного автомобиля конкретной модели, изготовленной на случайном заводе. Отклонение этого автомобиля является суммой компонентов или вкладов, по одному от каждого из случайных членов.

Отображение имен случайных членов.

stats.rtnames
ans = 3x1 cell
    {'Factory'          }
    {'Factory*Car Model'}
    {'Error'            }

Вы не знаете отклонений, но их можно оценить из данных. Напомним, что ems поле stats структура выражает ожидаемое значение среднего квадрата каждого члена как линейной комбинации неизвестных отклонений для случайных членов и неизвестных квадратичных форм для фиксированных членов. Если вы возьмете ожидаемые средние квадратные выражения для случайных членов и приравняете эти ожидаемые значения к вычисленным средним квадратам, то получите систему уравнений, которую можно решить для неизвестных отклонений. Эти решения являются отклонением оценками компонентов.

Отображение оценки отклонения компонента для каждого термина.

stats.varest
ans = 3×1

    4.4426
   -0.0313
    0.1139

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

Создайте штриховой график компонентов.

bar(max(0,stats.varest))

Figure contains an axes. The axes contains an object of type bar.

gca.xtick = 1:3;
gca.xticklabel = stats.rtnames;

Можно также вычислить доверительные границы для оценки отклонения. The anovan функция делает это путем вычисления доверительных границ для ожидаемых средних квадратов отклонения и нахождения более низких и верхних пределов для каждого компонента отклонения, содержащего все эти ограничения. Эта процедура приводит к набору границ, который консервативен для сбалансированных данных. (То есть 95% доверительных границ будут иметь вероятность не менее 95% содержащих истинные отклонения, если количество наблюдений для каждой комбинации сгруппированных переменных одинаковое.) Для несбалансированных данных это приближения, которые не гарантированы как консервативные.

Просмотрите оценки отклонений и доверительные пределы для оценок отклонений каждого компонента.

[{'Term' 'Estimate' 'Lower' 'Upper'};
 stats.rtnames, num2cell([stats.varest stats.varci])]
ans=4×4 cell array
    {'Term'             }    {'Estimate'}    {'Lower' }    {'Upper'   }
    {'Factory'          }    {[  4.4426]}    {[1.0736]}    {[175.6038]}
    {'Factory*Car Model'}    {[ -0.0313]}    {[   NaN]}    {[     NaN]}
    {'Error'            }    {[  0.1139]}    {[0.0586]}    {[  0.3103]}