Рабочий процесс модели с линейными смешанными эффектами

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

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

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

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

Найдите два наблюдения в правом верхнем углу, которые выглядят как отклонения.

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

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

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

F(flu2.Region == 'ENCentral' & flu2.Date == '11/6/2005')

ans =

    1.4860

Случайная генерация значений ответа.

Случайная генерация значений ответа для национальной оценки 1,625, регион MidAtl, и дата 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

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