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

Этот пример показывает аппроксимацию и анализ линейной модели смешанных эффектов (LME).

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

load flu

The 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, date 11/6/2005.

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

ans =

    1.4860

Случайным образом сгенерируйте значения отклика.

Случайным образом сгенерируйте значения отклика для национальной оценки 1,625, область MidAtl и 4/23/2006 даты. Сначала определите новую таблицу. Поскольку Дата (Date) и Область (Region) являются номинальными в исходной таблице, необходимо определить их аналогично в новой таблице.

tblnew.Date = nominal('4/23/2006');
tblnew.WtdILI = 1.625;
tblnew.Region = nominal('MidAtl');
tblnew = struct2table(tblnew);

Теперь сгенерируйте значение отклика.

random(lme4,tblnew)

ans =

    1.2679

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