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

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 и Tomato являются переменными фиксированных эффектов, и среднее выражение изменяется в зависимости от блока (тип почвы), и графики внутри блоков (типы томатов в пределах типов почвы) независимо.

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

Спрогнозируйте значения отклика по исходным проектным значениям. Отобразите первые пять предсказаний с наблюдаемыми значениями отклика.

yhat = predict(lme);
[yhat(1:5) ds.Yield(1:5)]

ans = 5×2

  115.4788  104.0000
  135.1455  136.0000
  152.8121  158.0000
  160.4788  174.0000
   58.0839   57.0000

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

load carsmall

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

tbl = table(MPG,Weight,Model_Year);
lme = fitlme(tbl,'MPG ~ Weight + (1|Model_Year)');

Создайте предсказанные отклики на данные.

yhat = predict(lme,tbl);

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

figure()
gscatter(Weight,MPG,Model_Year)
hold on
gscatter(Weight,yhat,Model_Year,[],'o+x')
legend('70-data','76-data','82-data','70-pred','76-pred','82-pred')
hold off

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

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 и Tomato являются переменными фиксированных эффектов, и среднее выражение изменяется в зависимости от блока (тип почвы), и графики внутри блоков (типы томатов в пределах типов почвы) независимо.

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

Создайте новый массив набора данных с проектными значениями. Новый массив набора данных должен иметь те же переменные, что и исходный массив набора данных, который вы используете для подбора кривой модели lme.

dsnew = dataset();
dsnew.Soil = nominal({'Sandy';'Silty'});
dsnew.Tomato = nominal({'Cherry';'Vine'});
dsnew. Fertilizer = nominal([2;2]);

Спрогнозируйте условные и предельные отклики в исходных проектных точках.

yhatC = predict(lme,dsnew);
yhatM = predict(lme,dsnew,'Conditional',false);
[yhatC yhatM]

ans = 2×2

   92.7505  111.6667
   87.5891   82.6667

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

load carbig

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

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

X = [ones(406,1) Acceleration Horsepower];
Z = [ones(406,1) Acceleration];
Model_Year = nominal(Model_Year);
G = Model_Year;

Теперь подбирайте модель используя fitlmematrix с определенными матрицами проекта и сгруппированными переменными.

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

Создайте матрицы проекта, которые содержат данные, при которых можно предсказать значения отклика. Xnew должно иметь три столбца, как в X. Первый столбец должен быть столбцом 1с. А значения в последних двух столбцах должны соответствовать Acceleration и Horsepower, соответственно. Первый столбец Znew должен быть столбцом 1с, а второй столбец должен содержать то же Acceleration значения как в Xnew. Исходная сгруппированная переменная в G является модельным годом. Итак, Gnew должны содержать значения для модельного года. Обратите внимание, что Gnew должны содержать номинальные значения.

Xnew = [1,13.5,185; 1,17,205; 1,21.2,193];
Znew = [1,13.5; 1,17; 1,21.2]; % alternatively Znew = Xnew(:,1:2);
Gnew = nominal([73 77 82]);

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

yhat = predict(lme,Xnew,Znew,Gnew)

yhat = 3×1

    8.7063
    5.4423
   12.5384

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

Z = {ones(406,1),Acceleration};
G = {Model_Year,Model_Year};

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

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

Znew = {[1;1;1],[13.5;17;21.2]};
MY = nominal([73 77 82]);
Gnew = {MY,MY};

Спрогнозируйте ответы, используя новые матрицы проекта.

yhat = predict(lme,Xnew,Znew,Gnew)

yhat = 3×1

    8.6365
    5.9199
   12.1247

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

load carbig

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

tbl = table(MPG,Acceleration,Horsepower,Model_Year);

Теперь подбирайте модель используя fitlme с определенными матрицами проекта и сгруппированными переменными.

lme = fitlme(tbl,'MPG ~ Acceleration + Horsepower + (Acceleration|Model_Year)');

Создайте новые данные и сохраните их в новой таблице.

tblnew = table();
tblnew.Acceleration = linspace(8,25)';
tblnew.Horsepower = linspace(min(Horsepower),max(Horsepower))';
tblnew.Model_Year = repmat(70,100,1);

linspace создает 100 одинаково больших значений между нижним и верхним входными пределами. Model_Year фиксируется на 70. Можно повторить это для любого модельного года.

Вычислите и постройте график предсказанных значений и 95% доверительных пределов (несовпадающих).

[ypred,yCI,DF] = predict(lme,tblnew);
figure(); 
h1 = line(tblnew.Acceleration,ypred);
hold on;
h2 = plot(tblnew.Acceleration,yCI,'g-.');

Отображение степеней свободы.

DF(1)

ans = 389

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

[ypred,yCI,DF] = predict(lme,tblnew,'Simultaneous',true);
h3 = plot(tblnew.Acceleration,yCI,'r--');

Отображение степеней свободы.

DF

DF = 389

Вычислите одновременные доверительные границы с помощью метода Satterthwaite, чтобы вычислить степени свободы.

[ypred,yCI,DF] = predict(lme,tblnew,'Simultaneous',true,'DFMethod','satterthwaite');
h4 = plot(tblnew.Acceleration,yCI,'k:');
hold off
xlabel('Acceleration')
ylabel('Response')
ylim([-50,60])
xlim([8,25])
legend([h1,h2(1),h3(1),h4(1)],'Predicted response','95%','95% Sim',...
'95% Sim-Satt','Location','Best')

Отображение степеней свободы.

DF

DF = 3.6001