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

Открыть скрипт

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

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

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
    {'(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

Маленькое $p$- значения 0,0045885 и 3.0502e-76 указывают, что и прерывание и общенациональная оценка являются значительными. Кроме того, пределы достоверности для стандартного отклонения термина случайных эффектов$\sigma_{b}$, не включайте 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

$p$- значение 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

$p$- значение 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

$p$- значение 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

Документация Statistics and Machine Learning Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте