В этом примере показано, как соответствовать и анализировать линейную модель смешанных эффектов (ЛБМ).
load flu
flu
массив набора данных имеет Date
переменная и 10 переменных, содержащих оцененные уровни гриппа (в 9 различных областях, оцененных от поисковых запросов Google®, плюс общенациональная оценка от CDC).
Чтобы подбирать линейно смешанную модель эффектов, ваши данные должны быть в правильно отформатированном массиве набора данных. Чтобы подбирать линейную модель смешанных эффектов с уровнями гриппа как ответы, объедините эти девять столбцов, соответствующих областям в массив. Новый массив набора данных, flu2
, должен иметь переменную отклика FluRate
, номинальная переменная Region
это показывает, какая область каждая оценка от, общенациональная оценка WtdILI
, и сгруппированная переменная Date
.
flu2 = stack(flu,2:10,'NewDataVarName','FluRate',... 'IndVarName','Region'); flu2.Date = nominal(flu2.Date);
Задайте flu2
как table
.
flu2 = dataset2table(flu2);
Постройте уровни гриппа по сравнению с общенациональной оценкой.
plot(flu2.WtdILI,flu2.FluRate,'ro') xlabel('WtdILI') ylabel('Flu Rate')
Вы видите, что уровни гриппа в областях имеют непосредственную связь с общенациональной оценкой.
Подбирайте линейную модель смешанных эффектов с общенациональной оценкой как переменный предиктор и случайная точка пересечения, которая варьируется Date
.
lme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date)')
lme = Linear mixed-effects model fit by ML Model information: Number of observations 468 Fixed effects coefficients 2 Random effects coefficients 52 Covariance parameters 2 Formula: FluRate ~ 1 + WtdILI + (1 | Date) Model fit statistics: AIC BIC LogLikelihood Deviance 286.24 302.83 -139.12 278.24 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.16385 0.057525 2.8484 466 0.0045885 {'WtdILI' } 0.7236 0.032219 22.459 466 3.0502e-76 Lower Upper 0.050813 0.27689 0.66028 0.78691 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17146 Lower Upper 0.13227 0.22226 Group: Error Name Estimate Lower Upper {'Res Std'} 0.30201 0.28217 0.32324
Маленькое - значения 0,0045885 и 3.0502e-76 указывают, что и точка пересечения и общенациональная оценка являются значительными. Кроме того, пределы достоверности для стандартного отклонения термина случайных эффектов, не включайте 0 (0.13227, 0.22226), который указывает, что термин случайных эффектов является значительным.
Постройте необработанные остаточные значения по сравнению с подходящими значениями.
figure();
plotResiduals(lme,'fitted')
Отклонение остаточных значений увеличивается с увеличением подходящих значений отклика, который известен heteroscedasticity.
Найдите эти два наблюдения относительно правых верхних, которые появляются как выбросы.
find(residuals(lme) > 1.5)
ans = 98 107
Переоборудуйте модель путем удаления этих наблюдений.
lme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date)','Exclude',[98,107]);
Определите если включая независимый случайный термин для общенациональной оценки, сгруппированной Date
улучшает модель.
altlme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date) + (WtdILI-1|Date)',... 'Exclude',[98,107])
altlme = Linear mixed-effects model fit by ML Model information: Number of observations 466 Fixed effects coefficients 2 Random effects coefficients 104 Covariance parameters 3 Formula: FluRate ~ 1 + WtdILI + (1 | Date) + (WtdILI | Date) Model fit statistics: AIC BIC LogLikelihood Deviance 179.39 200.11 -84.694 169.39 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.17837 0.054585 3.2676 464 0.001165 {'WtdILI' } 0.70836 0.030594 23.153 464 2.123e-79 Lower Upper 0.0711 0.28563 0.64824 0.76849 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.16631 Lower Upper 0.12977 0.21313 Group: Date (52 Levels) Name1 Name2 Type Estimate Lower {'WtdILI'} {'WtdILI'} {'std'} 4.6939e-08 NaN Upper NaN Group: Error Name Estimate Lower Upper {'Res Std'} 0.26691 0.24934 0.28572
Предполагаемое стандартное отклонение WtdILI
термин - почти 0, и его доверительный интервал не может быть вычислен. Это - индикация, что модель сверхпараметрируется и (WtdILI-1|Date)
термин не является значительным. Можно официально протестировать это использование compare
метод можно следующим образом: compare(lme,altlme,'CheckNesting',true)
.
Добавьте случайный термин эффектов для точки пересечения, сгруппированной областью к первоначальной модели lme
.
lme2 = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date) + (1|Region)',... 'Exclude',[98,107]);
Сравните модели lme
и lme2
.
compare(lme,lme2,'CheckNesting',true)
ans = THEORETICAL LIKELIHOOD RATIO TEST Model DF AIC BIC LogLik LRStat deltaDF pValue lme 4 177.39 193.97 -84.694 lme2 5 62.265 82.986 -26.133 117.12 1 0
- значение 0 указывает на тот lme2
лучшая подгонка, чем lme
.
Теперь проверяйте при добавлении потенциально коррелированого термина случайных эффектов для точки пересечения, и средний национальный показатель улучшает модель lme2
.
lme3 = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date) + (1 + WtdILI|Region)',... 'Exclude',[98,107])
lme3 = Linear mixed-effects model fit by ML Model information: Number of observations 466 Fixed effects coefficients 2 Random effects coefficients 70 Covariance parameters 5 Formula: FluRate ~ 1 + WtdILI + (1 | Date) + (1 + WtdILI | Region) Model fit statistics: AIC BIC LogLikelihood Deviance 13.338 42.348 0.33076 -0.66153 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.1795 0.054953 3.2665 464 0.0011697 {'WtdILI' } 0.70719 0.04252 16.632 464 4.6451e-49 Lower Upper 0.071514 0.28749 0.62363 0.79074 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17634 Lower Upper 0.14093 0.22064 Group: Region (9 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.0077038 {'WtdILI' } {'(Intercept)'} {'corr'} -0.059604 {'WtdILI' } {'WtdILI' } {'std' } 0.088069 Lower Upper 3.2081e-16 1.8499e+11 -0.99996 0.99995 0.051693 0.15004 Group: Error Name Estimate Lower Upper {'Res Std'} 0.20976 0.19568 0.22486
Оценка для стандартного отклонения термина случайных эффектов для точки пересечения, сгруппированной областью, 0.0077037, ее доверительный интервал является очень большим и включает нуль. Это указывает, что случайные эффекты для точки пересечения, сгруппированной областью, незначительны. Корреляция между случайными эффектами для точки пересечения и WtdILI
-0.059604. Его доверительный интервал является также очень большим и включает нуль. Это - индикация, что корреляция не является значительной.
Переоборудуйте модель путем устранения точки пересечения из (1 + WtdILI | Region)
термин случайных эффектов.
lme3 = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date) + (WtdILI - 1|Region)',... 'Exclude',[98,107])
lme3 = Linear mixed-effects model fit by ML Model information: Number of observations 466 Fixed effects coefficients 2 Random effects coefficients 61 Covariance parameters 3 Formula: FluRate ~ 1 + WtdILI + (1 | Date) + (WtdILI | Region) Model fit statistics: AIC BIC LogLikelihood Deviance 9.3395 30.06 0.33023 -0.66046 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.1795 0.054892 3.2702 464 0.0011549 {'WtdILI' } 0.70718 0.042486 16.645 464 4.0496e-49 Lower Upper 0.071637 0.28737 0.62369 0.79067 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17633 Lower Upper 0.14092 0.22062 Group: Region (9 Levels) Name1 Name2 Type Estimate Lower {'WtdILI'} {'WtdILI'} {'std'} 0.087925 0.054474 Upper 0.14192 Group: Error Name Estimate Lower Upper {'Res Std'} 0.20979 0.19585 0.22473
Все условия в новой модели lme3
являются значительными.
Сравните lme2
и lme3
.
compare(lme2,lme3,'CheckNesting',true,'NSim',100)
ans = SIMULATED LIKELIHOOD RATIO TEST: NSIM = 100, ALPHA = 0.05 Model DF AIC BIC LogLik LRStat pValue lme2 5 62.265 82.986 -26.133 lme3 5 9.3395 30.06 0.33023 52.926 0.009901 Lower Upper 0.00025064 0.053932
- значение 0,009901 указывает на тот lme3
лучшая подгонка, чем lme2
.
Добавьте квадратичный термин фиксированных эффектов в модель lme3
.
lme4 = fitlme(flu2,'FluRate ~ 1 + WtdILI^2 + (1|Date) + (WtdILI - 1|Region)',... 'Exclude',[98,107])
lme4 = Linear mixed-effects model fit by ML Model information: Number of observations 466 Fixed effects coefficients 3 Random effects coefficients 61 Covariance parameters 3 Formula: FluRate ~ 1 + WtdILI + WtdILI^2 + (1 | Date) + (WtdILI | Region) Model fit statistics: AIC BIC LogLikelihood Deviance 6.7234 31.588 2.6383 -5.2766 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} -0.063406 0.12236 -0.51821 463 0.60456 {'WtdILI' } 1.0594 0.16554 6.3996 463 3.8232e-10 {'WtdILI^2' } -0.096919 0.0441 -2.1977 463 0.028463 Lower Upper -0.30385 0.17704 0.73406 1.3847 -0.18358 -0.010259 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.16732 Lower Upper 0.13326 0.21009 Group: Region (9 Levels) Name1 Name2 Type Estimate Lower {'WtdILI'} {'WtdILI'} {'std'} 0.087865 0.054443 Upper 0.1418 Group: Error Name Estimate Lower Upper {'Res Std'} 0.20979 0.19585 0.22473
- значение 0,028463 указывает что коэффициент квадратичного термина WtdILI^2
является значительным.
F = fitted(lme4); R = response(lme4); figure(); plot(R,F,'rx') xlabel('Response') ylabel('Fitted')
Подходящие по сравнению с наблюдаемыми значениями отклика формируют почти угол в 45 градусов, указывающий на хорошую подгонку.
Постройте остаточные значения по сравнению с подходящими значениями.
figure();
plotResiduals(lme4,'fitted')
Несмотря на то, что это улучшилось, можно все еще видеть некоторый heteroscedasticity в модели. Эта сила произойти из-за другого предиктора, который не существует в наборе данных, следовательно не в модели.
F(flu2.Region == 'ENCentral' & flu2.Date == '11/6/2005')
ans = 1.4860
Случайным образом сгенерируйте значения отклика для национальной оценки 1,625, область Мидэтл и дата 4/23/2006. Во-первых, задайте новую таблицу. Поскольку Дата и область номинальны в исходной таблице, необходимо задать их так же в новой таблице.
tblnew.Date = nominal('4/23/2006'); tblnew.WtdILI = 1.625; tblnew.Region = nominal('MidAtl'); tblnew = struct2table(tblnew);
Теперь сгенерируйте значение отклика.
random(lme4,tblnew)
ans = 1.2679