线性混合效应模型工作流

打开脚本

这个例子展示了如何拟合和分析线性混合效应模型(LME)。

加载示例数据。

负载流感

的流感数据集数组具有日期变量，以及包含估计流感发病率的10个变量（在9个不同地区，通过Google®搜索估计，再加上CDC的全国性估计）。

重新组织和绘制数据。

为了适应线性混合效果模型，数据必须在正确格式化的数据集数组中。要拟合一个以流感率作为响应的线性混合效应模型，请将与区域相对应的九列组合成一个数组。新的数据集数组，flu2，必须具有响应变量FluRate，名义变量地区这显示了每个估计来自哪个地区，全国范围的估计WtdILI，以及分组变量日期．

flu2=堆叠（flu，2:10，“NewDataVarName”，“FluRate”，．..“因德瓦名”，“地区”)；flu2.日期=标称（flu2.日期）；

定义flu2作为一个桌子．

flu2 = dataset2table (flu2);

绘制流感发病率与全国估计值的对比图。

地块（flu2.WtdILI、flu2.FluRate、，“罗”)包含(“WtdILI”) ylabel (“流感率”）

你可以看到地区的流感率与全国的估计有直接关系。

拟合LME模型并解释结果。

拟合线性混合效应模型，以全国估计值作为预测变量，随机截距随时间变化日期．

lme=fitlme（flu2，'FluRate ~ 1 + WtdILI + (1|Date)'）

lme=线性混合效应模型通过ML模型拟合信息：观测次数468固定效应系数2随机效应系数52协方差参数2公式：FluRate~1+WtdILI+（1 |日期）模型拟合统计：AIC BIC对数似然偏差286.24 302.83-139.12 278.24固定效应系数（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上下0.050813 0.27689 0.66028 0.78691随机效应协方差参数（95%CI）：组：日期（52个级别）名称1 Name类型估计{（Intercept）}（Intercept）}{（Intercept）}std）}0.17146上下0.13227 0.22226组：错误名称估计上下{'Res Std'}0.30201 0.28217 0.32324

小 $p$ -值0.0045885和3.0502e-76表示截距和全国估计值都是显著的。同时，随机效应项标准差的置信限， $\sigma_{b}$ ，不包括0（0.13227，0.22226），这表明随机效应项是显著的。

绘制原始残差与拟合值的对比图。

图（）；绘制残差（lme，“合适的”）

残差方差随拟合响应值的增大而增大，称为异方差。

找到右上方的两个异常值。

查找（残差（lme）>1.5）

Ans = 98 107

通过删除这些观测值来重新构建模型。

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

改进模型。

确定是否包括一个独立的随机项为分组的全国估计日期提高了模型。

altlme = fitlme (flu2,“FluRate~1+WtdILI+（1 |日期）+（WtdILI-1 |日期）”，．..“排除”,[98,107])

altlme=线性混合效应模型通过ML模型拟合信息：观测次数466固定效应系数2随机效应系数104协方差参数3公式：FluRate~1+WtdILI+（1 |日期）+（WtdILI |日期）模型拟合统计：AIC BIC对数似然偏差179.39 200.11-84.694 169.39固定效应系数（95%CIs）：名称估计值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级）名称1名称2类型估计值{（截距）{（截距）}std0.16631上下0.12977 0.21313组：日期（52级）名称1名称2类型估计下{'WtdILI'}{'WtdILI'}{'std'}4.6788e-08南上南组：错误名称估计下上{'Res std'}0.26691 0.24934 0.28572

的估计标准差WtdILI项接近0，其置信区间无法计算。这表明模型被过度参数化了(WtdILI-1 |日期)术语不重要。您可以使用比较方法如下:比较(lme, altlme CheckNesting,真的)．

将按区域分组的截距的随机效果项添加到初始模型中伦敦金属交易所．

lme2=fitlme（flu2，'FluRate ~ 1 + WtdILI +(1|日期)+(1|区域)'，．..“排除”[98107]);

比较模型伦敦金属交易所和lme2．

比较(lme lme2,“检查嵌套”,真正的)

ans =理论似然比检验模型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 10 0

的 $p$ -值为0表示lme2是一个更好的适合比伦敦金属交易所．

现在，检查为截距和全国平均值添加一个潜在相关的随机效应项是否可以改进模型lme2．

lme3=fitlme（flu2，'FluRate ~ 1 + WtdILI +(1|日期)+ (1 + WtdILI|区域)'，．..“排除”,[98,107])

公式:FluRate ~ 1 + WtdILI + (1 | Date) + (1 + WtdILI | Region)模型拟合统计量:AIC BIC LogLikelihood Deviance 13.338 42.348 0.33076 -0.66153固定效应系数(95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.1795 0.054953 3.2665 0.0011697 {'WtdILI'} 0.70719 0.04252 16.632 464 4.6451e-49 Lower Upper 0.071514 0.28749 0.62363 0.79074随机效应协方差参数(95% CIs):组:Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17634 Lower Upper 0.14093 0.22064 Group:区域(9个级别)Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.0077038 {'WtdILI'} {'corr'} -0.059604 {'WtdILI'} {'WtdILI'} {'std'} 0.088069 Lower Upper 3.235e -16 1.8344e+11 -0.99996 0.99995 0.051694 0.15004 Group: Error Name Estimate Lower Upper {'Res std'} 0.20976 0.19568 0.22486 .

按Region分组的截距随机效应项的标准差估计为0.0077037，其置信区间非常大，包括0。这说明按区域分组的拦截随机效应不显著。截距随机效应和WtdILI为-0.059604。其置信区间也非常大，包括零。这表明相关性不显著。

通过消除模型中的截距来重新安装模型(1 + WtdILI |区域)随机项。

lme3=fitlme（flu2，'FluRate ~ 1 + WtdILI + (1|Date) + (WtdILI - 1|Region)'，．..“排除”,[98,107])

公式:FluRate ~ 1 + WtdILI + (1 | Date) + (WtdILI | Region)模型拟合统计量: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个级别)Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17633 Lower Upper 0.14092 0.22062 Group: Region(9个级别)Name1 Name2 Type Estimate Lower {'WtdILI'} {'WtdILI'} 0.087925 0.054474 Upper 0.14192 Group: Error Name Estimate Lower Upper {'Res std'} 0.20979 0.19585 0.22473 . Date(52个级别)Name1 Name2 Type Estimate {'(Intercept)'} {' Intercept)} 0.17633 Lower Upper 0.14092 0.22062 Group: Region(9个级别)Name1 Name2 Type Estimate Lower {'WtdILI'} {'WtdILI'} 0.087925 0.054474 Upper 0.14192

新模型中的所有术语lme3是重要的。

比较lme2和lme3．

比较(lme2 lme3,“检查嵌套”,真的,“NSim”, 100)

ans =模拟似然比测试:NSIM = 100, ALPHA = 0.05模型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$ -value 0.009901表示指定lme3是一个更好的适合比lme2．

向模型中添加一个二次固定效应项lme3．

lme4 = fitlme (flu2,'FluRate ~ 1 + WtdILI^2 +(1|日期)+ (WtdILI - 1|区域)'，．..“排除”,[98,107])

公式:FluRate ~ 1 + WtdILI + WtdILI^2 + (1 | Date) + (WtdILI | Region)模型拟合统计量:AIC BIC对数似然偏差6.7234 31.588 2.6383 -5.2766固定效应系数(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随机效应协方差参数(95% CIs):组:Date(52个级别)Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.16732 Lower Upper 0.13326 0.21009 Group: Region(9个级别)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 . Date(52个级别)Name1 Name2 Type Estimate {'(Intercept)'} {' Intercept)'} 0.16732 Lower Upper 0.13326 0.21009

的 $p$ -value为0.028463表示二次项系数WtdILI ^ 2是显著的。

绘制拟合响应与观测响应和残差的曲线。

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

拟合与观测的响应值形成几乎45度角，表明拟合良好。

绘制残差与拟合值的关系。

图();plotResiduals (lme4“合适的”）

虽然它有所改进，但您仍然可以在模型中看到一些异方差。这可能是由于数据集中不存在另一个预测值，因此不在模型中。

请于11/6/2005找到符合中环地区的流感率数值。

F（flu2.0）区域==“ENCentral”& flu2。日期= ='11/6/2005'）

ans = 1.4860

随机生成响应值。

对2006年4月23日的全国估计值1.625随机生成响应值。首先，定义新表。因为Date和Region在原始表中是标称的，所以必须在新表中以类似的方式定义它们。

tblnew。日期= nominal(“4/23/2006”); tblnew.WtdILI=1.625；tblnew.区域=标称(“Midtl”）;tblnew = struct2table (tblnew);

现在，生成响应值。

随机（lme4，tblnew）

ans = 1.2679