广义线性模型工作流

打开实时脚本

这个例子展示了如何拟合一个广义线性模型并分析结果。典型的工作流程包括以下步骤:导入数据，拟合广义线性模型，测试其质量，修改模型以提高其质量，并基于模型进行预测。在本例中，您使用Fisher虹膜数据来计算一朵花属于两个类别之一的概率。

加载数据

加载费雪虹膜数据。

负载fisheriris

提取第51到150行，它们的分类是versicololor virginica。

X = meas(51:end，:);

创建的逻辑响应变量真正的对于彩色和假virginica。

Y = strcmp(“多色的”、物种(51:结束));

拟合广义线性模型

对数据拟合二项式广义线性模型。

mdl = fitglm(X,y，“线性”，“分布”，“二”）

mdl =广义线性回归模型:logit(y) ~ 1 + x1 + x2 + x3 + x4分布=二项式估计系数:估计SE tStat pValue ________ ______ _______ ________(截距)42.638 25.708 1.6586 0.097204 x1 2.4652 2.3943 1.0296 0.30319 x2 6.6809 4.4796 1.4914 0.13585 x3 -9.4294 4.7372 -1.9905 0.046537 x4 -18.286 9.7426 -1.8769 0.060529 100个观测值，95个误差自由度离散度:1 Chi^2-统计量vs常数模型:127,p值= 1.95e-26

根据模型显示，有的p中的-valuespValue列都不小，这意味着可以简化模型。

检查和改进模型

确定系数的95%置信区间是否包括0。如果是这样，您可以删除具有这些间隔的模型项。

coefCI(mdl)

confint =5×2-8.3984 93.6740 -2.2881 7.2185 -2.2122 15.5739 -18.8339 -0.0248 -37.6277 1.0554

只有第四个预测器x3具有置信区间不包括0的系数。

的系数x1而且x2有大量p-values及其95%置信区间包括0。测试两个系数是否都可以为零。指定一个假设矩阵来选择的系数x1而且x2．

M = [0 1 0 0 0 0 0 0 1 0 0];p = coefTest(mdl,M)

P = 0.1442

的p-value约为0.14，这并不小。删除x1而且x2从模型中。

mdl1 = removeTerms(mdl，'x1 + x2'）

mdl1 =广义线性回归模型:logit(y) ~ 1 + x3 + x4分布=二项式估计系数:估计SE tStat pValue ________ ______ _______ __________(截距)45.272 13.612 3.326 0.00088103 x3 -5.7545 2.3059 -2.4956 0.012576 x4 -10.447 3.7557 -2.7816 0.0054092 100个观测值，97个误差自由度离散度:1 Chi^2-统计量vs常数模型:118,p值= 2.3e-26

或者，您可以使用stepwiseglm．

mdl2 = stepwiseglm(X,y，“不变”，“分布”，“二”，“上”，“线性”）

1.加上x4, Deviance = 33.4208, Chi2Stat = 105.2086, PValue = 1.099298e-24加上x3，偏差= 20.5635,Chi2Stat = 12.8573, PValue = 0.000336166加上x2，偏差= 13.2658,Chi2Stat = 7.29767, PValue = 0.00690441

mdl2 =广义线性回归模型:logit(y) ~ 1 + x2 + x3 + x4分布=二项式估计系数:估计SE tStat pValue ________ ______ _______ ________(截距)50.527 23.995 2.1057 0.035227 x2 8.3761 4.7612 1.7592 0.078536 x3 -7.8745 3.8407 -2.0503 0.040334 x4 -21.43 10.707 -2.0014 0.04535 100个观测值，96个误差自由度离散度:1 Chi^2统计量vs常数模型:125,p值= 5.4e-27

的p值(pValue)x2在系数表中均大于0.05，但stepwiseglm包括x2在模型中，因为p值(PValue)，以添加x2小于0.05。的stepwiseglm函数计算PValue使用with和without的配合x2，而函数是计算的pValue基于仅从最终模型计算出的近似标准误差。因此,PValue更可靠pValue．

识别异常值

检查杠杆图以寻找有影响的异常值。

plotDiagnostics (mdl2“杠杆”）

图中包含一个轴对象。标题为Case order plot的axis对象包含2个类型为line的对象。这些对象代表杠杆，参考线。

如果一个观察值的杠杆值大大超过p / n,在那里p系数的个数是和吗n是观测的数量。虚线参考线是建议阈值，由2 * p / n，在该图中对应0.08。一些观察结果的杠杆值大于10 * p / n(即0.40)。确定这些观察点。

idxOutliers = find(mdl2.Diagnostics. diagnostics .)杠杆> 10*mdl2.NumCoefficients/mdl2.NumObservations)

idxOutliers =4×119 21 57 85

当你拟合排除这些点的模型时，看看模型系数是否会发生变化。

oldCoeffs = mdl2.Coefficients.Estimate;mdl3 = fitglm(X,y，“线性”，“分布”，“二”，.．.“PredictorVars”2:4,“排除”, idxOutliers);newCoeffs = mdl3.Coefficients.Estimate;disp ([oldCoeffs newCoeffs])

50.5268 44.0085 8.3761 5.6361 -7.8745 -6.1145 -21.4296 -18.1236

模型系数mdl3和mdl2中的不同。这一结果表明，在高杠杆点的响应与简化模型的预测值不一致。

预测五彩缤纷的概率

使用mdl3预测一朵平均测量值的花是五彩缤纷的概率。生成预测的置信区间。

[newf,newc] = predict(mdl3,mean(X))

Newf = 0.4558

newc =1×20.1234 - 0.8329

该模型给出了近46%的概率，平均花是五颜六色的，有一个宽的置信区间。

另请参阅