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

load('fertilizer.mat');

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

Сохраните данные в массиве набора данных под названием ds, в практических целях и определить Tomato, Soil, и Fertilizer как категориальные переменные.

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

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

где $$i$$ = 1, 2,..., 60 соответствует наблюдениям, $$j$$ = 2,..., 5 соответствует типам томатов, и $$k$$ = 1, 2, 3 соответствует блокам (грунту). $$S_k$$ представляет $$k$$ I тип грунта, и $$(S*T{)}_{jk}$$ представляет $$j$$ th томатный тип, вложенный в $$k$$ Тип грунта. $$I[T{]}_{ij}$$ - фиктивная переменная, представляющая уровень $$j$$ томатного типа.

Случайные эффекты и ошибка наблюдения имеют следующие предыдущие распределения: $$b_{0k}\sim N(0,{\sigma }_S^2)$$, $$b_{0jk}\sim N(0,{\sigma }_{S*T}^2)$$, и $${ϵ}_{ijk}\sim N(0,{\sigma }^2)$$.

lme = fitlme(ds,'Yield ~ Fertilizer + (1|Soil) + (1|Soil:Tomato)');

Вычислите оценки ковариационных параметров (оценки $${\sigma }_S^2$$ и $${\sigma }_{S*T}^2$$) членов случайных эффектов.

psi = covarianceParameters(lme)

psi=2×1 cell array
    {[3.8000e-17]}
    {[  352.8481]}

Вычислите остаточное отклонение ($${\sigma }^2$$).

[~,mse] = covarianceParameters(lme)

mse = 151.9007

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

load('weight.mat');

weight содержит данные продольного исследования, где 20 субъектов случайным образом назначены для 4 программ упражнений, и их потеря веса регистрируется в течение шести 2-недельных периодов времени. Это моделируемые данные.

Сохраните данные в массиве набора данных. Определите Subject и Program как категориальные переменные.

ds = dataset(InitialWeight,Program,Subject,Week,y);
ds.Subject = nominal(ds.Subject);
ds.Program = nominal(ds.Program);

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

Для 'reference' фиктивное переменное кодирование, fitlme использует Программу A как ссылку и создает необходимые переменные фиктива $$I[.]$$. Эта модель соответствует

где $$i$$ соответствует номеру наблюдения, $$i=1,2,...,120$$, и $$m$$ соответствует номеру субъекта, $$m=1,2,...,20$$. $${\beta }_j$$ являются коэффициентами с фиксированными эффектами, $$j=0,1,...,8$$, и $$b_{0m}$$ и $$b_{1m}$$ являются случайными эффектами. $$IW$$ означает начальный вес и $$I[.]$$ является фиктивной переменной, представляющей тип программы. Для примера, $$I[PB{]}_i$$ - фиктивная переменная, представляющая программу B.

Случайные эффекты и ошибка наблюдения имеют следующие предыдущие распределения:

и

lme = fitlme(ds,'y ~ InitialWeight + Program + (Week|Subject)');

Вычислите оценки ковариационных параметров для случайных эффектов.

[psi,mse,stats] = covarianceParameters(lme)

psi = 1x1 cell array
    {2x2 double}

mse = 0.0105

stats=2×1 cell array
    {3x7 classreg.regr.lmeutils.titleddataset}
    {1x5 classreg.regr.lmeutils.titleddataset}

mse - расчетное остаточное отклонение. Это оценка для $${\sigma }^2$$.

Чтобы увидеть оценки ковариационных параметров для членов со случайными эффектами ($${\sigma }_0^2$$, $${\sigma }_1^2$$, и $${\sigma }_{0,1}$$), индексируйте в psi.

psi{1}

ans = 2×2

    0.0572    0.0490
    0.0490    0.0624

Оценка отклонения термина случайных эффектов для точки пересечения, $${\sigma }_0^2$$, составляет 0,0572. Оценка отклонения термина случайных эффектов для недели, $${\sigma }_1^2$$, составляет 0,0624. Оценка ковариации терминов случайных эффектов для точки пересечения и недели, $${\sigma }_{0,1}$$, составляет 0,0490.

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

stats{1}

ans = 
    Covariance Type: FullCholesky

    Group      Name1                  Name2                  Type        
    Subject    {'(Intercept)'}        {'(Intercept)'}        {'std' }    
    Subject    {'Week'       }        {'(Intercept)'}        {'corr'}    
    Subject    {'Week'       }        {'Week'       }        {'std' }    


    Estimate    Lower      Upper  
    0.23927     0.14365    0.39854
    0.81971     0.38662    0.95658
     0.2497     0.18303    0.34067

На отображении показано имя параметра группирования (Group), переменные случайных эффектов (Name1, Name2), тип ковариационных параметров (Type), оценка (Estimate) для каждого параметра и 95% доверительных интервалов для параметров (Lower, Upper). Приведенные в этой таблице оценки связаны с оценками в psi следующим образом.

Стандартное отклонение термина случайных эффектов для точки пересечения составляет 0,23927 = sqrt (0,0527). Аналогично, стандартное отклонение термина случайных эффектов для недели составляет 0,2497 = sqrt (0,0624). Наконец, корреляция между терминами случайных эффектов точки пересечения и недели 0,81971 = 0 .0490/( 0,23927 * 0,2497).

Обратите внимание, что это отображение также показывает, какой ковариационный шаблон вы используете при подборе модели. В этом случае ковариационный шаблон FullCholesky. Чтобы изменить ковариационный шаблон для членов со случайными эффектами, вы должны использовать 'CovariancePattern' аргумент пары "имя-значение" при подборе модели.

Вторая камера stats включает аналогичную статистику для остаточного стандартного отклонения. Отображение содержимого второй камеры.

stats{2}

ans = 
    Group    Name               Estimate    Lower       Upper  
    Error    {'Res Std'}        0.10261     0.087882    0.11981

Оценка остаточного стандартного отклонения является квадратным корнем mse, 0,10261 = sqrt (0,0105 ).

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

load carbig

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

где $$m=1,2,...,13$$ представляет уровни для переменной Model_Year, и $$k=1,2,...,8$$ представляет уровни для переменной Origin. $$MP{G}_{imk}$$ - мили на галлон для i observation,|m| го модельного года, and'k|-го источники, которые соответствуют i-му наблюдению. Условия случайных эффектов и ошибка наблюдения имеют следующие предыдущие распределения:

Здесь, термин случайных эффектов $$b_{1m}$$ представляет первый случайный эффект на уровне $$m$$ первой сгруппированной переменной. Термин случайных эффектов $$b_{10m}$$ соответствует первому члену случайных эффектов (1), для точки пересечения (0), в $$m$$ 1й уровень ( $$m$$ ) первой сгруппированной переменной. Аналогично $$b_{11m}$$ является уровнем $$m$$ для первого предиктора (1) в первом члене (1) случайных эффектов.

Точно так же, $$b_{2k}$$ означает второй термин случайных эффектов на уровне $$k$$ второй сгруппированной переменной.

$${\sigma }_{10}^2$$ - отклонение термина случайных эффектов для точки пересечения, $${\sigma }_{11}^2$$ является отклонением термина случайных эффектов для ускорения предиктора, и $${\sigma }_{10,11}$$ - ковариация членов случайных эффектов для точки пересечения и ускорения предиктора. $${\sigma }_2^2$$ является отклонением второго термина случайных эффектов, и $${\sigma }^2$$ - остаточное отклонение.

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

X = [ones(406,1) Acceleration Weight];
Z = {[ones(406,1) Acceleration],[Weight]};
Model_Year = nominal(Model_Year);
Origin = nominal(Origin);
G = {Model_Year,Origin};

Подгонка модели с помощью матриц проекта.

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Weight'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'},{'Weight'}},'RandomEffectGroups',{'Model_Year','Origin'});

Вычислите оценки ковариационных параметров для случайных эффектов.

[psi,mse,stats] = covarianceParameters(lme)

psi=2×1 cell array
    {2x2 double  }
    {[6.6778e-08]}

mse = 9.0750

stats=3×1 cell array
    {3x7 classreg.regr.lmeutils.titleddataset}
    {1x7 classreg.regr.lmeutils.titleddataset}
    {1x5 classreg.regr.lmeutils.titleddataset}

Остаточное отклонение mse - 9.0755. psi является массивом ячеек 2 на 1 и stats - массив ячеек 3 на 1. Чтобы увидеть содержимое, вы должны индексировать в эти массивы ячеек.

Во-первых, индекс в первую камеру psi.

psi{1}

ans = 2×2

    8.2648   -0.8699
   -0.8699    0.1158

Первая камера psi содержит ковариационные параметры для коррелированных случайных эффектов для точки пересечения $${\sigma }_{10}^2$$ как 8.5160, и для ускорения $${\sigma }_{11}^2$$ как 0.1087. Оценка ковариации членов случайных эффектов для точки пересечения и ускорения $${\sigma }_{10,11}$$ -0,8387.

Теперь индексируйте во вторую камеру psi.

psi{2}

ans = 6.6778e-08

Вторая камера psi содержит оценку отклонения термина случайных эффектов для веса $${\sigma }_2^2$$.

Индекс в первую камеру stats.

stats{1}

ans = 
    Covariance Type: FullCholesky

    Group         Name1                   Name2                   Type        
    Model_Year    {'Intercept'   }        {'Intercept'   }        {'std' }    
    Model_Year    {'Acceleration'}        {'Intercept'   }        {'corr'}    
    Model_Year    {'Acceleration'}        {'Acceleration'}        {'std' }    


    Estimate    Lower       Upper   
     2.8749       1.0481      7.8853
    -0.8894     -0.98663    -0.32528
    0.34023      0.19356     0.59803

В этой таблице показаны стандартные оценки отклонений для терминов случайных эффектов для точки пересечения и ускорения. Обратите внимание, что стандартные оценки отклонений являются квадратными корнями диагональных элементов в первой камере psi. В частности, 2.9182 = sqrt (8.5160) и 0.32968 = sqrt (0.1087). Корреляция является функцией ковариации точки пересечения и ускорения, и стандартных отклонений точки пересечения и ускорения. Ковариация точки пересечения и ускорения является недиагональным значением в первой камере psi, -0,8387. Итак, корреляция составляет -.8387/( 0,32968 * 2,92182) = -0.87.

Для точки пересечения и ускорения сгруппированной переменной Model_Year.

Индекс во вторую камеру stats.

stats{2}

ans = 
    Covariance Type: FullCholesky

    Group     Name1             Name2             Type           Estimate  
    Origin    {'Weight'}        {'Weight'}        {'std'}        0.00025842


    Lower         Upper    
    9.0892e-05    0.0007347

Вторая камера stats имеет стандартную оценку отклонения и 95% доверительные пределы для стандартного отклонения термина случайных эффектов для Weight. Это сгруппированная переменная Origin.

Индекс в третью камеру stats.

stats{3}

ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0125      2.8024    3.2383

Третья камера stats содержит оценку остаточного стандартного отклонения и 95% доверительных пределов. Оценка остаточного стандартного отклонения является квадратным корнем mse, sqrt (9.0755) = 3.0126.

Создайте 99% доверительные интервалы для ковариационных параметров.

[~,~,stats] = covarianceParameters(lme,'Alpha',0.01);
stats{1}

ans = 
    Covariance Type: FullCholesky

    Group         Name1                   Name2                   Type        
    Model_Year    {'Intercept'   }        {'Intercept'   }        {'std' }    
    Model_Year    {'Acceleration'}        {'Intercept'   }        {'corr'}    
    Model_Year    {'Acceleration'}        {'Acceleration'}        {'std' }    


    Estimate    Lower      Upper    
     2.8749     0.76334       10.827
    -0.8894     -0.9932    0.0022801
    0.34023     0.16213      0.71399

stats{2}

ans = 
    Covariance Type: FullCholesky

    Group     Name1             Name2             Type           Estimate  
    Origin    {'Weight'}        {'Weight'}        {'std'}        0.00025842


    Lower         Upper    
    6.5453e-05    0.0010202

stats{3}

ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0125      2.7395    3.3127