用广义线性模型拟合数据

打开脚本

这个例子展示了如何使用拟合和评估广义线性模型glmfit而且glmval．普通线性回归可用于拟合一条直线，或任何参数为线性的函数，以拟合具有正态分布误差的数据。这是最常用的回归模型;然而，它并不总是现实的。广义线性模型从两个方面扩展了线性模型。首先，通过引入连杆函数，放宽了参数线性的假设。其次，可以对非正态分布的误差分布进行建模

广义线性模型

回归模型根据一个或多个预测变量(通常表示为x1, x2等)定义响应变量(通常表示为y)的分布。最常用的回归模型，普通线性回归，将y建模为正态随机变量，其均值是预测因子的线性函数，b0 + b1*x1 +…，其方差为常数。在单个预测器x的最简单情况下，模型可以表示为每个点都有高斯分布的直线。

Mu = @(x) -1.9+.23*x;X = 5:.1:15;Yhat = (x);Dy = -3.5:.1:3.5;Sz =大小(dy);K =(长度(dy)+1)/2;X1 = 7*ones(sz);Y1 = mu(x1)+dy;Z1 = normpdf(y1,mu(x1)，1);X2 = 10*ones(sz); y2 = mu(x2)+dy; z2 = normpdf(y2,mu(x2),1); x3 = 13*ones(sz); y3 = mu(x3)+dy; z3 = normpdf(y3,mu(x3),1); plot3(x,yhat,zeros(size(x)),“b -”，.．.(x1, y1, z1,的r -， x1([k k])，y1([k k])，[0 z1(k)]，“:”，.．.x2, y2、z2的r -， x2([k k])，y2([k k])，[0 z2(k)]，“:”，.．.x3, y3、z3、的r -， x3([k k])，y3([k k])，[0 z3(k)]，“:”）;zlim ([0 1]);包含(“X”）;ylabel (“Y”）;zlabel (的概率密度）;网格在；视图(45 [-45]);

在广义线性模型中，响应的平均值被建模为预测因子线性函数g(b0 + b1*x1 +…)的单调非线性变换。变换g的逆函数称为“链接”函数。例如logit (sigmoid)链接和log链接。同样，y也可能是非正态分布，如二项分布或泊松分布。例如，具有对数链接和单个预测器x的泊松回归可以表示为每个点的泊松分布的指数曲线。

Mu = @(x) exp(-1.9+.23*x);X = 5:.1:15;Yhat = (x);X1 = 7*ones(1,5);Y1 = 0:4;Z1 = poisspdf(y1,mu(x1));X2 = 10*ones(1,7);Y2 = 0:6;Z2 = poisspdf(y2,mu(x2));X3 = 13*ones(1,9); y3 = 0:8; z3 = poisspdf(y3,mu(x3)); plot3(x,yhat,zeros(size(x)),“b -”，.．.[x1;x1],[日元;y1]、[z1;0(大小(y1))),的r -(x1, y1, z1,“r”。，.．.[x2;x2]、[y2;y2]、[z2;0(大小(y2))),的r -x2, y2, z2,“r”。，.．.[x3;x3]、[y3;y3]、[z3;0(大小(y3))),的r -z3、x3 y3,“r”。）;zlim ([0 1]);包含(“X”）;ylabel (“Y”）;zlabel (“概率”）;网格在；视图(45 [-45]);

Logistic回归拟合

这个例子涉及到一个实验，以帮助建模各种重量的汽车在里程测试中失败的比例。数据包括对重量的观察，测试的汽车数量和失败的数量。

一套汽车砝码重量= [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';每种重量测试的汽车数量。%测试= [48 42 31 34 31 21 23 23 21 16 17 21]';每个重量下未通过测试的汽车数量。Failed = [1 2 0 3 8 8 14 17 19 15 17 21]';每种重量下不合格的汽车的比例。比例=失败。/测试;情节(重量、比例、“年代”)包含(“重量”）;ylabel (“比例”）;

这张图是汽车故障的比例，作为重量的函数。可以合理地假设失败计数来自二项分布，其概率参数P随着权重的增加而增加。但是P到底是如何依赖于重量的呢?

我们可以试着把这些数据拟合成一条直线。

linearCoef = polyfit(权重，比例，1);linearFit = polyval(线性系数，权重);情节(重量、比例、“年代”、重量、linearFit的r -，[2000 4500]，[0 0]，凯西:”，[2000 4500]，[1 1]，凯西:”)包含(“重量”）;ylabel (“比例”）;

这种线性拟合存在两个问题:

1)这条线预测了小于0和大于1的比例。

2)比例不是正态分布的，因为它们必然是有界的。这违反了拟合简单线性回归模型所需的假设之一。

使用高阶多项式可能会有所帮助。

[cubicCoef,stats,ctr] = polyfit(权重，比例，3);cubicFit = polyval(cubicCoef,weight，[]，ctr);情节(重量、比例、“年代”、重量、cubicFit的r -，[2000 4500]，[0 0]，凯西:”，[2000 4500]，[1 1]，凯西:”)包含(“重量”）;ylabel (“比例”）;

然而，这种配合仍然存在类似的问题。从图中可以看出，当重量超过4000时，拟合比例开始下降;事实上，权重值越大，它就会变成负值。当然，正态分布假设仍然是不成立的。

相反，更好的方法是使用glmfit拟合逻辑回归模型。逻辑回归是广义线性模型的一种特殊情况，对于这些数据，它比线性回归更合适，原因有二。首先，它使用一种适合于二项分布的拟合方法。其次，逻辑链接将预测比例限制在[0,1]范围内。

对于逻辑回归，我们指定预测矩阵，矩阵的一列包含失败数，一列包含测试数。我们还指定了二项分布和logit链接。

[logitCoef,dev] = glmfit(重量，[失败的测试]，“二”，分对数的）;logitFit = glmval(logitCoef,weight，分对数的）;情节(重量、比例、“废话”、重量、logitFit的r -）;包含(“重量”）;ylabel (“比例”）;

如图所示，当重量变小或变大时，拟合比例渐近于0和1。

模型诊断

的glmfit函数提供了许多用于检查拟合和测试模型的输出。例如，我们可以比较两个模型的偏差值，以确定平方项是否会显著改善拟合。

[logitCoef2,dev2] = glmfit([weight重量。]^ 2],(未测试)“二”，分对数的）;Pval = 1 - chi2cdf(dev-dev2,1)

Pval = 0.4019

较大的p值表明，对于这些数据，二次项并没有显著改善拟合。两种拟合的曲线图显示，两种拟合的差异不大。

logitFit2 = glmval(logitCoef2，[weight weight.^2]，分对数的）;情节(重量、比例、“废话”、重量、logitFit的r -、重量、logitFit2“g -”）;传奇(“数据”，“线性条件”，“线性和二次项”，“位置”，“西北”）;

为了检查拟合的优度，我们也可以看看皮尔森残差的概率图。这些是标准化的，所以当模型与数据合理拟合时，它们大致具有标准的正态分布。(如果没有这种标准化，残差将有不同的方差。)

[logitCoef,dev,stats] = glmfit(weight，[失败的测试]，“二”，分对数的）;normplot (stats.residp);

残差图与正态分布很好地吻合。

评估模型预测

一旦我们对模型感到满意，我们就可以用它来进行预测，包括计算置信范围。在这里，我们预测在测试的100辆车中，在四个权重中每一个都无法通过里程测试的预期数量。

weightPred = 2500:500:4000;[failedPred,dlo,dhi] = glmval(logitCoef,weightPred，分对数的统计,.95,100);dlo errorbar (weightPred failedPred,济,“:”）;

二项式模型的链接函数

对于五个分布中的每一个glmfit万博1manbetx支持，有一个规范(默认)链接函数。对于二项分布，规范链接是logit。然而，对于二项式模型，还有其他三个联系是合理的。这四个都保持在区间[0,1]内的平均响应。

Eta = -5:.1:5;Plot (eta,1 ./ (1 + exp(-eta))，“- - -”埃塔,normcdf (eta),“- - -”，.．.1 -exp(-exp(Eta))“- - -”埃塔,exp (exp (eta)),“- - -”）;包含(“预测因子的线性函数”）;ylabel (“预测平均反应”）;传奇(分对数的，“probit”，互补的双对数的，“对数”，“位置”，“东”）;

例如，我们可以比较具有probit链接的拟合与具有logit链接的拟合。

probitCoef = glmfit(重量，[测试失败]，“二”，“probit”）;probitFit = glmval(probitCoef，权重，“probit”）;情节(重量、比例、“废话”、重量、logitFit的r -、重量、probitFit“g -”）;传奇(“数据”，“Logit模型”，“Probit模型”，“位置”，“西北”）;

数据通常很难区分这四个链接函数，通常要根据理论依据进行选择。