R语言数据分析与挖掘(第四章):回归分析(4)

您所在的位置:网站首页 logistic模型的检验评价 R语言数据分析与挖掘(第四章):回归分析(4)

R语言数据分析与挖掘(第四章):回归分析(4)

2024-04-28 23:49| 来源: 网络整理| 查看: 265

前面我们介绍的回归方法,一般适用于数值型数据对象,对于分类数据类型就不再适用。对于分类数据对象,我们需要引入广义线性回归方法,比如logistic回归和poisson回归模型。这里我们介绍logistic回归。

logistic回归又称logistic回归分析,是一种广义的线性回归分析模型,常用于数据挖掘,疾病自动诊断,经济预测等领域。例如,探讨引发疾病的危险因素,并根据危险因素预测疾病发生的概率等。以胃癌病情分析为例,选择两组人群,一组是胃癌组,一组是非胃癌组,两组人群必定具有不同的体征与生活方式等。因此因变量就为是否胃癌,值为“是”或“否”,自变量就可以包括很多了,如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的,也可以是分类的。然后通过logistic回归分析,可以得到自变量的权重,从而可以大致了解到底哪些因素是胃癌的危险因素。同时根据该权值可以根据危险因素预测一个人患癌症的可能性。

R语言中用于实现logistic回归的函数是glm(),其基本书写格式为:

代码语言:javascript复制glm(formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = list(...), model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)

参数介绍:

Formula:指定用于拟合的模型公式,类似于Im中的用法:

Family: 指定描述干扰项的概率分 布和模型的连接函数, 默认值为gaussian, 若需进行logistic同归,则需设置为binomial(link = "logit");

Data:指定用于回归的数据对象,可以是数据框、列表或能被强制转换为数据框的数据对象:

Weights:一个向量,用于指定每个观测值的权重:

Subset:一个向量,指定数据中需要包含在模型中的观测值;

Na.ction:一个函数,指定当数据中存在缺失值时的处理办法,用法与Im中的一致;

Start:一个数值型向量,用于指定现行预测器中参数的初始值;

Etastart:一个数值型向量,用于指定现行预测器的初始值;

Mustart:一个数值型向量,用于指定均值向量的初始值:

Offset:指定用于添加到线性项中的一组系数恒为1的项:

Contol:指定控制拟合过程的参数列表,其中epsilon 表示收敛的容忍度,maxit表示迭代的最大次数,trace 表示每次迭代是否打印具体信息;

Model: 逻辑值,指定是否返回“模型框架”,默认值为TRUE:

Method;指定用于拟合的方法,“glm.ft”表示用于拟合,“model.frame"表示可以返回模型框架;

X:逻辑值,指定是否返回“横型矩阵”,默认值为FALSE:

Y:逻辑值,制度是否能够返回响应变量,默认值为TRUE;

Contrasts:模型中因子对照的列表。

  下面利用iris 数据集进行操作演练,由于iris数据集中的分类变量Specics中有三种元素:setosa、versicolor 和virginica,即鸢尾花的有三个不同的种类,在建模之前,先对数据集进行处理,将数据集中Species属于setosa类的数据剔除,然后利用剩余的数据进行建模分析,具体操作如下:

代码语言:javascript复制> iris iris$Species log1 summary(log1) Call: glm(formula = Species ~ ., family = binomial(link = "logit"), data = iris) Deviance Residuals: Min 1Q Median 3Q Max -2.01105 -0.00541 -0.00001 0.00677 1.78065 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -42.638 25.707 -1.659 0.0972 . Sepal.Length -2.465 2.394 -1.030 0.3032 Sepal.Width -6.681 4.480 -1.491 0.1359 Petal.Length 9.429 4.737 1.991 0.0465 * Petal.Width 18.286 9.743 1.877 0.0605 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 138.629 on 99 degrees of freedom Residual deviance: 11.899 on 95 degrees of freedom AIC: 21.899 Number of Fisher Scoring iterations: 10

上述代码表示:选择iris数据集中第51行到150行的数据,将该数据集中变量

Species列中记录为virginica 的替换为1,否则替换为0,然后利用清洗好的数据进行logistic回归;模型的输出结果显示:解释变量Sepal.Length和Sepal.Width没能通过显著性水平为0.05的检验。

下面基于前面介绍的AIC准则(R语言数据分析与挖掘(第四章):回归分析(3)——变量的选择)进行逐步回归:

代码语言:javascript复制> log2 summary(log2) Call: glm(formula = Species ~ Sepal.Width + Petal.Length + Petal.Width, family = binomial(link = "logit"), data = iris) Deviance Residuals: Min 1Q Median 3Q Max -1.75795 -0.00412 0.00000 0.00290 1.92193 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -50.527 23.995 -2.106 0.0352 * Sepal.Width -8.376 4.761 -1.759 0.0785 . Petal.Length 7.875 3.841 2.050 0.0403 * Petal.Width 21.430 10.707 2.001 0.0453 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 138.629 on 99 degrees of freedom Residual deviance: 13.266 on 96 degrees of freedom AIC: 21.266 Number of Fisher Scoring iterations: 10

不难发现逐步回归剔除了变量:Sepal.Length,对逐步回归的结果进行详细展示,可以看到剩余变量的参数估计均通过了显著性水平的0.05的检验,说明构建模型得到了数据的支持。

下面根据WMCR对模型的预测能力进行评估:

代码语言:javascript复制> pred prob yhat0.5) > table(iris$Species,yhat) yhat 0 1 0 48 2 1 1 49

上述代码表示:首先将模型的预测结果存储到变量pred中,再根据前面介绍的模型进行logit变换的逆变换,输出结果存储到变量prob,此时该变量中的值即为响应变量取值为1的概率值,即变量Species=virginica的概率值,然后分别计算变量prob中大于0.5和小于等于0.5的记录总数,其中0.5即为阈值,由于原始的鸢尾花数据中,种类为versicolor和virginica的记录各有50条,故阈值a取值为0.5。最后利用函数table( )统计原始数据中的记录和预测结果的记录情况(“0”表示versicolor,“1”表示virginica), 不难发现,输出的表格中,数字“48”和“49”均表示预测正确的总数,数字“2”表示真实种类为versicolor, 而预测结果为virginica 的记录总数,

类似地,数字“1”表示真实种类为virginica,而预测结果为versicolor 的记录总数。除此之外,还可以利用图形展示模型的预测效果,业界一般采用ROC曲线对logistic 回归模型的效果进行刻画,R语言的RORC包中有专门的函数用于刻画ROC曲线,具体操作如下:

代码语言:javascript复制> library(ROCR) > pred2 performance(pred2,'auc')@y.values [[1]] [1] 0.9972 > perf plot(perf,col=2,type="l",lwd=2) > f=function(x){ y=x;return(y)} > curve(f(x),0,1,col=4,lwd=2,lty=2,add=T)


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3