乐趣区

关于数据挖掘:数据分享R语言回归模型诊断离群值分析学生考试成绩病人医护质量满意度婴儿死亡率和人均收入针叶树荫面积数据

一些规范的图形工具能够极大地帮忙了解数据集并评估所倡议模型的品质。

学生考试成绩

例如,数据集蕴含 600 个察看后果,用于国家统计教育核心对学生进行的一项十分大的钻研。数据集中的一些变量包含:

•性别:性别男性或女性。

•种族:种族或民族,具备西班牙裔,亚洲人,非洲裔美国人,白人的程度。

•学校类型,公立和私立。

•轨迹:管制位点,一个间断的协变量,批示受试者对影响他们的事件的自我感知管制水平(更高 = 更感知的管制)

•概念:自我概念,反映受试者自我感知的学术能力的间断协变量(更高 = 更多的自我感知能力)。

•动机,一个间断的协变量(更高 = 更多的自我感知的动机)。

•数学:标准化的数学考试问题。

• 浏览:标准化浏览测试问题。

要根本理解数据集中的变量,咱们能够应用摘要命令:

summary(hsb)

 

然而,因为变量太少,咱们实际上能够用成对命令把所有的连续变量绘制在一起。

 

图中显示,尽管动机是一个连续变量,但它只具备四个类别。这很可能是因为学生被要求在 4 个点的量表上评估他们的自我激励(0,0.33,0.66,和 1)。另外,数学和任何一个间断预测因素之间最强的关联仿佛是与浏览。

anova(M1, M2)

 

有两个标准图来查看 M1 是否仿佛违反了线性模型的假如。第一种是绘制残差 ei = yi – x ′ iβˆ作为预测值的函数ˆyi = x ′ iβˆ。除了通常的残差之外,还有两种类型的残差常常被思考。

1. 标准化残差:ri = ei σˆ。这些残差是单位的

2. 研究性残差:di = ei σˆ √ 1 – hi 

这样,如果模型是实在的,残差大概是 N(0,1)。

图 1 显示了残差与 M1 的拟合值之间的关系,应用的是通常的残差和 studentized 残差的一个版本,ei/ √ 1 – hi。其余统计学家更喜爱在标准化的尺度上绘制残差(除以 σˆ),这样 Y 轴上的单位对每张残差图都有雷同的解释。

plot(predict(M1), res, pch = 21, bg = "black", cex = cex, cex.axis = cex,

ex)

 因为 σ2(1-hi)是一个方差,咱们必须有(1-hi)>0=⇒hi<1=⇒(1-hi)<1。因而,studentized 残差 ei/ √ 1 – hi 在规模上大于原始的 ei。

通过绘制残差与实践上的正态分布的直方图,能够更容易地看到这一点,如图所示。

hist(res/sigma.hat, breaks = 50, freq = FALSE, cex.axis = cex,

  对于这个残差图,我集体更喜爱应用标准化的残差,即除以ˆσ。这是因为我更容易在 N(0,1)上发现大的残差值(规范正态在±2 之间有 95% 的机会,在±3 之间有 99.7% 的机会)。图 2 证实了残差稍微偏差负值,表明须要进一步建模。

线性预测模型

回顾一下满意度数据集,即病人自我报告对在医院承受的护理品质的满意度,由年龄、压力和病情的重大水平来预测。上面是潜在模型之间的一些比拟。

anova(lm(Satisfaction ~ Age

 

anova(lm(Satisfaction ~ Age, data = satis),

 

anova(lm(Satisfaction ~ Age + Severity, data = satis),

 

 

依据下面的测试后果,咱们能够决定思考应用年龄和压力的模型,因为它的 F 统计量比年龄和重大水平稍显重要。然而,这种办法会减少抉择谬误模型的机会,正如上面的模拟实验所显示的那样。

 

pval <- 2*pt(abs(beta.hat)/sqrt(diag(vcov(M))\[-1\]),

df = n-3, lower.tail = FALSE) # pvalues

sim.out\[ii,\] <- round(c(beta.hat\[1\], pval\[1\], beta.hat\[2\], pval\[2\]),2)

}

sim.out

 

在 10 次试验运行中,只有两次对两个协变量都有显著的 P 值,但对这些协变量的预计却有很大偏差。事实上,其中两个复制甚至把协变量的符号弄错了。尽管这个例子并不事实(协变量之间的相关性简直不可能达到 0.999),但它阐明了这样一个事实:当协变量高度相干时,回归很难弄清 y 的变动是由一个协变量还是另一个协变量引起的。换句话说,βˆ的估计值从一个随机样本到另一个随机样本会有很大的变动。这种景象被称为方差收缩

掂量预测因子之间的共线性的一个简略办法是所谓的方差收缩因子(VIF)。

其中 C 是协变量之间的样本相关矩阵。

上面的 R 代码阐明了 VIF 的两种模式之间的平等。

VIF <- diag(solve(cor(X)))

VIF

 

names(satis.star) <- c("Satisfaction", paste0(names(satis)\[-1\], ".star"))

head(satis.star)

 

diag(vcov(M.orig))\["Severity"\]/diag(vcov(M.trans))\["Severity.star"\]

 

VIF\["Severity"\]

 

重大水平和压力的方差膨胀系数为 2.00,这两个变量之间的相关度为 0.67。相比之下,模拟实验中的方差膨胀系数是 1 /(1-.9992)≈500。作为一条教训法令,VIF 大于 10 是值得关注的。

离群值

简略地说,离群值是指与其余观测值相比具备异样大的残差。然而,为了进步模型的拟合度而删除这些观测值的状况比拟少见,因为出错的简直总是模型,而不是数据。也就是说,异样值通常是模型中脱漏了重要协变量的后果。例如,数据集蕴含了 105 个世界国家的婴儿死亡率和人均收入的数据。数据集中的变量是

– 支出:以元计的人均收入。

– 婴儿:每 1000 个活产的婴儿死亡率。

– 地区:大陆,有非洲、美洲、亚洲、欧洲等级别。

– 石油:石油出口国,级别为无、有。

对数尺度的线性模型适宜于数据,残差与拟合值绘制在图 3 中。

coef(M) # 系数 
as.character(formula(M))\[c(2,1,3)\]), collapse = " "), # 模型

xlab = "Predicted Log-Infant Mortality",

这张图上的异样点是最下面的三角形,它对应的是沙特阿拉伯。目前还不分明为什么沙特阿拉伯的婴儿死亡率比模型预测的要高很多。一种可能性是,它的石油产量太高,以至于扭曲了平均收入。在这种状况下,咱们有以下抉择。

1. 将沙特阿拉伯排除在剖析之外,抵赖咱们并不齐全晓得起因(不举荐)。

2. 将所有 9 个产油国排除在剖析之外,并阐明咱们的预测只针对非产油国。

3. 将沙特阿拉伯排除在外,并评估这对咱们的剖析有何影响。

树木的针叶产生的阴凉面积

数据集蕴含以下对于 35 种针叶树的变量。- 针叶面积。树木的针叶产生的阴凉面积 – 高度:树木的高度 – 树干尺寸。树的直径 以下模型适宜于该数据。

X <- model.matrix(M)

H <- X %*% solve(crossprod(X), t(X)) # HAT matrix

head(diag(H))

 

range(h - diag(H))

 

在这四种类型的残差中,每一种都把大残差放大得越来越多。为了掂量每个观测点的总体影响,图中的库克间隔与杠杆率作了比照

n <- nrow(lforest)

hbar <- p/n # 均匀杠杆率

abline(v = 2*hbar, col = "gray60", lty = 2) # 2x 均匀杠杆率 

其中一个观测值的库克间隔简直是其余观测值的 3 倍以上(红色),而其中的 e 个观测值的均匀杠杆率是两倍(蓝色)。为了理解为什么这些观测值有很高的杠杆率,咱们能够看一下所有协变量的配对图(在这种状况下,这只是高度和树干尺寸)。

咱们能够看到,其中三个标记点(有影响力的点和两个高杠杆点)的 TrunkSize 值十分大(记得这些值是按对数比例绘制的)。剩下的那个高杠杆点的树干尺寸绝对于其高度来说十分小。至于为什么只有三棵树中树干尺寸最大的一棵被标记为离群点,咱们能够将 NeedleArea 与每个协变量(包含理论值和预测值)作比照。为了进行比拟,预测是在所有观测值和省略一个观测值的状况下进行的:要么是有影响力的观测值,要么是有最高杠杆的观测值。

# 用省略的值进行预测

omit.ind <- c(which(inf.ind), # 最有影响力的

which.max(h)) # 最高杠杆率

which(infl.ind) 中的谬误:"which" 的参数不合乎逻辑

names(omit.ind) <- c("inf", "lev")

omit.ind

 

yobs <- lforest\[, "NeedleArea"\] # 察看值

Mf <- lm(NeedleArea ~ Height + TrunkSize, data = lforest) # 所有观测值



for(jj in 1:length(opt.ind)) {

# 用省略的察看值建设模型

Mo <- lm(NeedleArea ~ Height + TrunkSize,

data = lforest, subset = -omit.ind\[jj\])

yfito <- predict(Mo, newdata = lforest) # 所有观测值的拟合值

for(ii in c("Height", "TrunkSize")) {

咱们能够看到,当移除有影响的观测值时,有无脱漏的预测值之间的差别更加显著。这在针孔面积与树干尺寸的图中最为显著,这是对具备高杠杆作用的三个点和具备高影响的点的预测。在这个非凡的案例中,咱们确定具备最大树干尺寸的三棵树的测量是不正确的,它们能够从剖析中移除。

退出移动版