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

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. Для получения дополнительной информации смотрите Фиктивные Переменные Индикатора.

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

где $$y_{im}$$ наблюдение $$i$$ для уровня $$m$$ из группирующей переменной Date, $${\beta }_j$$, $$j$$ = 0, 1..., 8, коэффициенты фиксированных эффектов, с $${\beta }_0$$ будучи коэффициентом для области NE. $$b_{0m}$$ случайный эффект для уровня $$m$$ из группирующей переменной Date, и $${\epsilon }_{im}$$ ошибка наблюдения для наблюдения $$i$$. Случайный эффект имеет предшествующее распределение, $$b_{0m}\sim N(0,{\sigma }_b^2)$$ и остаточный член имеет распределение, $${\epsilon }_{im}\sim N(0,{\sigma }^2)$$.

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    Lower 
    '(Intercept)'        '(Intercept)'        'std'        0.6443      0.5297


    Upper  
    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.

Пределы достоверности для стандартного отклонения термина случайных эффектов, $${\sigma }_b$$, не включайте 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) соответствует оценке для $${\beta }_0$$ и beta(4) соответствует оценке для $${\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(fullfile(matlabroot,'examples','stats','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)