Загрузите образцы данных.

load flu

flu массив наборов данных имеет Date переменная и 10 переменных, содержащих оценочные показатели гриппа (в 9 различных регионах, оцененные по поискам Google ®, плюс общенациональная оценка Центра по контролю и профилактике заболеваний, CDC).

Для соответствия модели с линейными смешанными эффектами данные должны находиться в правильно отформатированном массиве наборов данных. Чтобы подогнать линейную модель смешанных эффектов с частотой гриппа в качестве ответов и областью в качестве предикторной переменной, объедините девять столбцов, соответствующих областям, в массив. Новый массив наборов данных, flu2, должна иметь переменную ответа, FluRate, номинальная переменная, Region, которая показывает, из какой области взята каждая оценка, и переменную группировки Date.

flu2 = stack(flu,2:10,'NewDataVarName','FluRate','IndVarName','Region');
flu2.Date = nominal(flu2.Date);

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

Регион - категориальная переменная. Можно указать контрасты для категориальных переменных с помощью DummyVarCoding аргумент пары имя-значение при подгонке модели. Если контрасты не заданы, fitlme использует 'reference' контраст по умолчанию. Потому что модель имеет перехват, fitlme принимает первый регион, NE, в качестве ссылки и создает восемь фиктивных переменных, представляющих другие восемь областей. Например, $$I[MidAtl$$] - фиктивная переменная, представляющая регион.MidAtl. Дополнительные сведения см. в разделе Фиктивные переменные.

Соответствующая модель:

где $${}_{\mathrm{yim}}$$ - наблюдение $$i$$ для уровня $$m$$ переменной группировки Date, $${}_{\mathrm{\beta j}}$$, $$j$$ = 0, 1,..., 8, являются коэффициентами с фиксированными эффектами, при этом $${\mathrm{\beta 0}}_{\mathrm{}}$$является коэффициентом для области NE. $${}_{\mathrm{}\mathrm{b0m}}$$ - случайный эффект для уровня $$m$$ переменной группировки Date, и $${}_{\mathrm{\alpha im}}$$ - погрешность наблюдения для наблюдения $$$$i. Случайный эффект имеет предыдущее распределение, $${}_{\mathrm{}}\mathrm{b0m\sim N}\mathrm{(}0{,}_{}^{\mathrm{}}$$startb2), а член ошибки имеет распределение$${,}_{}\mathrm{\epsilon im\sim N}{(}^0,$$start2).

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1|Date)')

lme = 
Linear mixed-effects model fit by ML

Model information:
    Number of observations             468
    Fixed effects coefficients           9
    Random effects coefficients         52
    Covariance parameters                2

Formula:
    FluRate ~ 1 + Region + (1 | Date)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    318.71    364.35    -148.36          296.71  

Fixed effects coefficients (95% CIs):
    Name                        Estimate    SE          tStat      DF 
    {'(Intercept)'     }          1.2233    0.096678     12.654    459
    {'Region_MidAtl'   }        0.010192    0.052221    0.19518    459
    {'Region_ENCentral'}        0.051923    0.052221     0.9943    459
    {'Region_WNCentral'}         0.23687    0.052221     4.5359    459
    {'Region_SAtl'     }        0.075481    0.052221     1.4454    459
    {'Region_ESCentral'}         0.33917    0.052221      6.495    459
    {'Region_WSCentral'}           0.069    0.052221     1.3213    459
    {'Region_Mtn'      }        0.046673    0.052221    0.89377    459
    {'Region_Pac'      }        -0.16013    0.052221    -3.0665    459


    pValue        Lower        Upper    
     1.085e-31       1.0334       1.4133
       0.84534    -0.092429      0.11281
        0.3206    -0.050698      0.15454
    7.3324e-06      0.13424      0.33949
       0.14902     -0.02714       0.1781
    2.1623e-10      0.23655      0.44179
       0.18705    -0.033621      0.17162
       0.37191    -0.055948      0.14929
     0.0022936     -0.26276    -0.057514

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.6443  


    Lower     Upper  
    0.5297    0.78368

Group: Error
    Name               Estimate    Lower      Upper
    {'Res Std'}        0.26627     0.24878    0.285

$$$$Значения p 7.3324e-06 и 2.1623e-10 соответственно показывают, что фиксированные эффекты заболеваемости гриппом в регионах WNCentral и ESCentral значительно отличаются по сравнению с показателями гриппа в регионе NE.

Доверительные пределы для среднеквадратического отклонения члена случайных эффектов, $${}_{\mathrm{startb}}$$, не включают 0 (0,5297, 0,78368), что указывает на то, что член случайных эффектов является значимым. Можно также проверить значимость терминов случайных эффектов с помощью compare способ.

Условный подогнанный ответ от модели при данном наблюдении включает в себя вклад фиксированных и случайных эффектов. Например, предполагаемый лучший линейный несмещенный предиктор (BLUP) частоты гриппа для региона WNCentral на неделе 10/9/2005

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

beta = fixedEffects(lme);
[~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS)
STATS.Level = nominal(STATS.Level);
y_hat = beta(1) + beta(4) + STATS.Estimate(STATS.Level=='10/9/2005')

y_hat = 1.2884

В предыдущем расчете beta(1) соответствует оценке для $${}_{\mathrm{\beta 0}}$$ и beta(4) соответствует оценке для $${}_{\mathrm{\beta 3}}$$. Вы можете просто отобразить подогнанное значение с помощью fitted способ.

F = fitted(lme);
F(flu2.Date == '10/9/2005' & flu2.Region == 'WNCentral')

ans = 1.2884

Предполагаемый предельный ответ для региона WNCentral на неделе 10/9/2005

Вычислите соответствующий краевой отклик.

F = fitted(lme,'Conditional',false);
F(flu2.Date == '10/9/2005' & flu2.Region == 'WNCentral')

ans = 1.4602

Загрузите образцы данных.

load('weight.mat');

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

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

tbl = table(InitialWeight,Program,Subject,Week,y);
tbl.Subject = nominal(tbl.Subject);
tbl.Program = nominal(tbl.Program);

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

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

Вычислите соответствующие значения и необработанные остатки.

F = fitted(lme);
R = residuals(lme);

Постройте график остаточных значений по сравнению с подходящими значениями.

plot(F,R,'bx')
xlabel('Fitted Values')
ylabel('Residuals')

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

figure()
gscatter(F,R,Program)