線形混合効果モデルのワクフロ

この例では，lme(線形混合効果モデル)の近似と解析の方法を説明します。

標本デ，タを読み込みます。

负载流感

デ，タセット配列流感には，変数日期と,インフルエンザ推定罹患率(谷歌®検索から推定される9地域の値とCDCによる全国の推定値)が格納されている10個の変数が含まれています。

デ，タを再編成してプロットします。

線形混合効果モデルを近似するには,データが適切な形式のデータセット配列になっていなければなりません。インフルエンザ羅患率を応答として使用することにより線形混合効果モデルをあてはめるには,地域に対応する9個の列を1つの配列にまとめます。新しいデ，タセット配列flu2には，応答変数FluRate，各推定の元になっている地域を示すノミナル変数地区，全国の推定値WtdILIおよびグル，プ化変数日期が含まれなければなりません。

Flu2 = stack(流感，2:10，“NewDataVarName”，“FluRate”，.．.“IndVarName”，“地区”）;flu2。日期= nominal(flu2.Date);

flu2を表格として定義します。

Flu2 = dataset2table(Flu2);

ンフルエンザ罹患率を全国の推定値に対してプロットします。

情节(flu2.WtdILI flu2.FluRate,“罗”)包含(“WtdILI”) ylabel (“流感率”）

各地域の@ @ンフルエンザ罹患率は全国の推定値と直接的な関係があることがわかります。

Lmeモデルを近似し，結果を解釈します。

予測子変数としての全国の推定値と，日期で変化するランダム切片により，線形混合効果モデルを近似します。

Lme = fitlme(flu2，' flate ~ 1 + WtdILI + (1|Date)'）

模型信息:观测数量468固定效应系数2随机效应系数52协方差参数2公式:波动~ 1 + WtdILI +(1 |日期)模型拟合统计量:AIC BIC log似然偏差286.24 302.83 -139.12 278.24固定效应系数(95% ci):名称估计SE tStat DF pValue{'(截取)'}0.16385 0.057525 2.8484 466 0.0045885 {'WtdILI'} 0.7236 0.032219 22.459 466 3.0502e-76下上0.050813 0.27689 0.66028 0.78691随机效应协方差参数(95% ci):组:日期(52级)Name1 Name2类型估计{'(截取)'}{'(截取)'}{'std'} 0.17146下上0.13227 0.22226组:错误名称估计下上{'Res std'} 0.30201 0.28217 0.32324

0.0045885および3.0502e-76という小さい値は，切片と全国の推定値の両方が有意であることを示しています。また，変量効果項 $\ sigma_ {b}$ の標準偏差の信頼限界(0.13227,0.22226)には,変量効果項が有意であることを示す0が含まれていません。

生の残差と近似値の対比をプロットします。

图();plotResiduals (lme,“安装”）

近似応答値の増加にともない残差の分散も増加します。この現象は不等分散性と呼ばれています。

外れ値のように右上に表示されている2件の観測値を検出します。

求(残差(lme) > 1.5)

Ans = 98 107

これらの観測値を削除して，モデルを再度近似します。

Lme = fitlme(flu2，' flate ~ 1 + WtdILI + (1|Date)'，“排除”[98107]);

モデルを改善します。

日期でグル，プ化された全国の推定値に独立したランダム項を含めると，モデルが改善されるか判定します。

Altlme = fitlme(flu2，' 1 + WtdILI +(1|日期)+ (WtdILI-1|日期)'，.．.“排除”[98107])

模型信息:观测数量466固定效应系数2随机效应系数104协方差参数3公式:波动~ 1 + WtdILI +(1 |日期)+ (WtdILI |日期)模型拟合统计量:AIC BIC log似然偏差179.39 200.11 -84.694 169.39固定效应系数(95% ci):名称估计SE tStat DF pValue{'(截取)'}0.17837 0.054585 3.2676 464 0.001165 {'WtdILI'} 0.70836 0.030594 23.153 464 2.123e-79下上0.0711 0.28563 0.64824 0.76849随机效应协方差参数(95% CIs):组:日期(52级)Name1 Name2类型估计{'(截取)'}{'(截取)'}{'std'} 0.16631下上0.12977 0.21313组:日期(52级)Name1 Name2类型估计下{'WtdILI'} {'WtdILI'} {'std'} 4.6788e-08 NaN上NaN组:错误名称估计下上{'Res Std'} 0.26691 0.24934 0.28572

WtdILI項の標準偏差の推定値はほぼ0なので，信頼区間を計算できません。これは，モデルのパラメ，タ，が多すぎ，(WtdILI-1 |日期)項が有意でないことを示しています。次のように比较メソッドを使用すると，これを正式に検定できます。比较(lme, altlme CheckNesting,真的)。

地域別にグル，プ化された切片の変量効果項を初期モデルlme三个月に追加します。

Lme2 = fitlme(flu2，'波动~ 1 + WtdILI +(1|日期)+(1|地区)'，.．.“排除”[98107]);

lme三个月モデルとlme2モデルを比較します。

比较(lme lme2,“CheckNesting”,真正的)

ans =理论似然比TEST模型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，' 1 + WtdILI +(1|日期)+ (1 + WtdILI|区域)'，.．.“排除”[98107])

模型信息:观测数量466固定效应系数2随机效应系数70协方差参数5公式:波动~ 1 + WtdILI +(1 |日期)+ (1 + WtdILI |区域)模型拟合统计量:AIC BIC log似然偏差13.338 42.348 0.33076 -0.66153固定效应系数(95% ci):名称估计SE tStat DF pValue{'(截距)'}0.1795 0.054953 3.2665 464 0.0011697 {'WtdILI'} 0.70719 0.04252 16.632 464 4.6451e-49上下0.071514 0.28749 0.62363 0.79074随机效应协方差参数(95% ci):组:日期(52级)Name1 Name2类型估计{'(截距)'}{'(截距)'}{'std'} 0.17634上下0.14093 0.22064组:区域(9个级别)Name1 Name2类型估计{'(拦截)'}{'(拦截)'}{'std'} 0.0077038 {'WtdILI'}{'(拦截)'}{'corr'} -0.059604 {'WtdILI'} {'WtdILI'} {'std'} 0.088069下上3.2352e-16 1.8344e+11 -0.99996 0.99995 0.051694 0.15004组:错误名称估计下上{'Res std'} 0.20976 0.19568 0.22486

地域(地区)別にグループ化された切片に関する変量効果項の標準偏差の推定値は0.0077037で信頼区間が非常に大きく,0が含まれています。これは，地域別にグル，プ化された切片の変量効果が有意でないことを示しています。切片とWtdILIの間で変量効果の相関は-0.059604です。この相関の信頼区間も非常に大きく，0が含まれています。これは，相関が有意でないことを示しています。

(1 + WtdILI |区域)変量効果項から切片を削除して，モデルを再度近似します。

Lme3 = fitlme(flu2，'波动~ 1 + WtdILI +(1|日期)+ (WtdILI - 1|区域)'，.．.“排除”[98107])

模型信息:观测数量466固定效应系数2随机效应系数61协方差参数3公式:波动~ 1 + WtdILI +(1 |日期)+ (WtdILI |区域)模型拟合统计量:AIC BIC log似然偏差9.3395 30.06 0.33023 -0.66046固定效应系数(95% ci):名称估计SE tStat DF pValue{'(截取)'}0.1795 0.054892 3.2702 464 0.0011549 {'WtdILI'} 0.70718 0.042486 16.645 464 4.0496e-49下上0.071637 0.28737 0.62369 0.79067随机效应协方差参数(95% ci):组:日期(52级)Name1 Name2类型估计{'(截取)'}{'(截取)'}{'std'} 0.17633下上0.14092 0.22062组:区域(9级)Name1 Name2类型估计下{'WtdILI'} {'WtdILI'} {'std'} 0.087925 0.054474上0.14192组:错误名称估计下上{'Res Std'} 0.20979 0.19585 0.22473

新しいモデルlme3の項はすべて有意です。

lme2とlme3を比較します。

比较(lme2 lme3,“CheckNesting”,真的,“NSim”, 100)

ans =模拟似然比测试: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下上0.00025064 0.053932

0.009901という値は，lme3がlme2より適切に当てはめられていることを示しています。

2次固定効果項をモデルlme3に追加します。

Lme4 = fitlme(flu2，'波动~ 1 + WtdILI^2 +(1|日期)+ (WtdILI - 1|区域)'，.．.“排除”[98107])

模型信息:观测数量466固定效应系数3随机效应系数61协方差参数3公式:波动~ 1 + WtdILI + WtdILI^2 +(1 |日期)+ (WtdILI |区域)模型拟合统计量:AIC BIC log似然偏差6.7234 31.588 2.6383 -5.2766固定效应系数(95% ci):名称估计SE tStat DF pValue{'(截取)'}-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上下-0.30385 0.17704 0.73406 1.3847 -0.18358 -0.010259随机效应协方差参数(95% ci):组:日期(52级)Name1 Name2类型估计{'(截取)'}{'(截取)'}{'std'} 0.16732上下0.13326 0.21009组:区域(9级)Name1 Name2类型估计下位{'WtdILI'} {'WtdILI'} {'std'} 0.087865 0.054443上位0.1418组:错误名称估计下位上位{'Res std'} 0.20979 0.19585 0.22473

0.028463という値は，2次項WtdILI ^ 2の係数が有意であることを示しています。

観測された応答と残差に対して，近似された応答をプロットします。

F =拟合(lme4);R =响应(lme4);图();情节(R、F“处方”)包含(“响应”) ylabel (“安装”）

観測された応答と比較すると，近似された応答の値の角度がほぼ45度になります。これは近似の精度が高いことを意味します。

残差と近似値の対比をプロットします。

图();plotResiduals (lme4“安装”）

モデルは改善されていますが，まだ不等分散性が認識されます。原因として，デ，タセットに存在しない(まり，)。

2005年11月6日のENCentral地域における近似されたインフルエンザ罹患率を特定します。

F (flu2。地区==“ENCentral”& flu2。日期= =“11/6/2005”）

Ans = 1.4860

応答値を無作為に生成します。

全国推定値が1.625でMidAtl地域および2006年4月23日に対する応答値を無作為に生成します。最初に新しいテ，ブルを定義します。日付と地域は元のテ，ブルでノミナルであるため，新しいテ，ブルでも同様に定義しなければなりません。

tblnew。日期= nominal(“4/23/2006”）;tblnew。WtdILI= 1.625; tblnew.Region = nominal(“MidAtl”）;Tblnew = struct2table(Tblnew);

次に，応答値を生成します。

随机(lme4 tblnew)

Ans = 1.2679