Линейный рабочий процесс модели Смешанных Эффектов

Этот пример показывает, как соответствовать и анализировать линейную модель смешанных эффектов (ЛБМ).

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

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')

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

Соответствуйте модели LME и интерпретируйте результаты.

Соответствуйте линейной модели смешанных эффектов общенациональной оценкой как переменная прогноза и случайное прерывание, которое отличается 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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.17146     0.13227


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

Переоборудуйте модель путем удаления этих наблюдений.

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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.16631     0.12977


    Upper  
    0.21313

Group: Date (52 Levels)
    Name1           Name2           Type         Estimate      Lower    Upper
    'WtdILI'        'WtdILI'        'std'        4.6824e-08    NaN      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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.17634     0.14093


    Upper  
    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.441e-16    1.7247e+11
     -0.99996       0.99995
     0.051696       0.15003

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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.17633     0.14092


    Upper  
    0.22062

Group: Region (9 Levels)
    Name1           Name2           Type         Estimate    Lower       Upper  
    'WtdILI'        'WtdILI'        'std'        0.087925    0.054474    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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.16732     0.13326


    Upper  
    0.21009

Group: Region (9 Levels)
    Name1           Name2           Type         Estimate    Lower       Upper 
    'WtdILI'        'WtdILI'        'std'        0.087865    0.054443    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 в модели. Эта сила произойти из-за другого предиктора, который не существует в наборе данных, следовательно не в модели.

Найдите подходящее значение уровня гриппа для области ENCentral, дата 11/6/2005.

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

Документация