R语言CH7-统计回归分析课件

1、R语言基础与数据科学应用沈刚 主编人民邮电出版社第七章 统计回归分析内容导航CONTENTS数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4定性数据如果一个样本的取值属于一组已知的且互不重叠的类型，我们把这样的数据样本被称为定性数据，也称作分类数据。iris数据是R自带的内置数据集之一，其中的鸢尾花分类信息就是一个定性数据的实际例子。 str(iris)data.frame:150 obs. of 5 variables: $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 . $ Sepal.Width

2、 : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 . $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 . $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 . $ Species : Factor w/ 3 levels setosa,versicolor,.: 1 1 1 1 1 1 1 1 1 1 .定性数据该数据集的最后一列是属性Species，包含了对鸢尾花品种的分类。品种的名称分别用setosa、vers

3、icolor、virginica等不同级别的因子来表示，这就是一个定性数据的例子。 iris$Species 1 setosa setosa setosa setosa setosa setosa 7 setosa setosa setosa setosa setosa setosa 145 virginica virginica virginica virginica virginica virginica Levels: setosa versicolor virginica定性数据数据变量的频数分布是对数据在一组不重叠的类别中出现次数的概括。 species species.freq s

4、pecies.relfreq species.relfreqspecies setosa versicolor virginica 0.3333333 0.3333333 0.3333333定性数据使用柱状图可视化频数信息 barplot(species.freq)使用饼状图可视化频数信息 pie(species.freq)自定义饼图颜色 colors = c(red, yellow, blue) pie(species.freq, col=colors)定性数据 species s_species s_iris options(digits=3) mean (s_iris$Sepal.Len

5、gth)1 5.01 tapply (iris$Sepal.Length, iris$Species, mean) setosa versicolor virginica 5.01 5.94 6.59以setosa的花萼长度为例，可以分步骤得出其花萼长度的平均值。option()函数指定R与用户的交互方式。tapply()函数实现同样的效果。定量数据 head(faithful, 3) eruptions waiting1 3.600 792 1.800 543 3.333 74计算喷发持续时间的频数分布步骤：使用range()函数得到数据上下界将变量活动范围分成不重叠的区间对变量进行区间分类

6、得到统计信息 预览faithful数据集的前三行。定量数据使用range()函数得到数据上下界 duration range(duration) 1 1.6 5.1 将变量活动范围分成不重叠的区间 breaks breaks 1 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5查看数据活动范围 。以0.5为间隔的序列得到向量定量数据对变量进行区间分类 duration.cut str(duration.cut) Factor w/ 8 levels 1.5,2),2,2.5),.: 5 1 4 2 7 3 7 5 1 6 .得到统计信息 duration.freq dur

7、ation.freq colors hist (duration, right=FALSE, col=colors, main=“老忠实喷发次数”,xlab=“持续时间（分钟), ylab = 频数)定量数据累计频数分布 cumsum (1:10) 1 1 3 6 10 15 21 28 36 45 55 # breaks的长度比duration.freq多1，手动添加0。 cumfreq0 plot (breaks, cumfreq0, main=“老忠实喷发持续时间, xlab=持续时间（分钟）, ylab=累积喷发频数) lines (breaks, cumfreq0)定量数据茎叶图 d

8、uration stem (duration)The decimal point is 1 digit(s) to the left of the | 16 | 070355555588 18 | 000022233333335577777777888822335777888 20 | 00002223378800035778 48 | 00000022335800333 50 | 0370小数点位于“|”左侧一位，所以相邻两个树干的间隔为0.2。在树叶排列时按照升序排列，因此在最后一行“0370”中，第一个0表示的是5.0，而第二个0表示的则是5.1。第一行：1.60、1.67、1.70、1

9、.73、1.75、1.75、1.75、1.75、1.75、1.75、1.78和1.78分钟。定量数据 duration waiting head (cbind(duration, waiting), 3) duration waiting1, 3.600 792, 1.800 543, 3.333 74 plot (duration, waiting,xlab=喷发持续时间,ylab=等待时间)内容导航CONTENTS数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4均值观测样本的均值，等于数据值的和与数据个数的比值，是对数据值的中心位置的数值度量。如果数据的规模为

10、n，其样本均值定义为：函数mean ()返回的是一个数值向量的均值，可将na.rm设置为TRUE排除向量中的缺失值。 duration mean(duration) 1 3.487783中位值、四分位数和百分位数中位值 median (duration) #调用median ()函数1 4 四分位数 quantile(duration) #调用quantile ()函数 0% 25% 50% 75% 100% 1.60000 2.16275 4.00000 4.45425 5.10000百分位数 quantile(duration, c(.25, .5, .96) 25% 50% 96% 2.

11、16275 4.00000 4.83572 四分位距四分位距是一个观测变量上下四分位数之差 IQR (duration) #调用IQR ()函数1 2.2915 箱形图适用于直观地表示出四分位数和四分位距 par (mfrow=c(1,2) #绘图布局，用1行2列排列方式 boxplot (faithful$eruptions) boxplot (faithful$waiting) par (mfrow=c(1,1)summary()函数 # R语言中提供的summary ()函数可以直接对数据生成一些统计信息 summary (faithful) eruptions waiting Min.

12、 :1.600 Min. :43.0 1st Qu.:2.163 1st Qu.:58.0 Median :4.000 Median :76.0 Mean :3.488 Mean :70.9 3rd Qu.:4.454 3rd Qu.:82.0 Max. :5.100 Max. :96.0 方差与标准差方差是对数据相对于其均值的分散程度的一个数值度量。函数var ()计算给定数值向量的方差。观测变量的标准差是其方差的平方根，用函数sd ()得到。 var (duration) #调用var ()函数1 1.3027 sd (duration) #调用sd ()函数1 1.1414协方差数据集中

13、两个变量x和y的协方差用于度量两者之间的线性相关度。如果协方差为正数，表明变量之间存在着正相关的关系，负的协方差表明两者之间是负相关关系。cov ()函数可以计算出两个向量的协方差 cov(duration, waiting) #调用cov ()函数1 13.978相关系数两个变量的相关系数等于它们的协方差除以各自标准差的乘积。当相关系数接近于1时，意味着变量之间线性正相关，在散点图上表现为散点几乎沿着一条斜率为正的直线分布。如果相关系数接近-1，则表示变量存在线性负相关关系，散点几乎落在一条斜率为负的直线附近。如果相关 系数为0，表明变量之间线性相关性很弱。cor ()函数就可以计算出它们的

14、相关系数。 cor (duration, waiting) #调用cor ()函数计算相关系数1 0.901内容导航CONTENTS数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4相关系数R语言中提供了一组函数，分别以d、p、q和r开头，后面跟着概率分布的名称，用于返回一个给定参数的随机分布的概率密度、累积概率密度、分位数和按给定分布生成的（伪）随机数。以正态分布为例，其名称为norm，所对应的函数及参数默认值具有下列形式： dnorm (x, mean = 0, sd = 1, log = FALSE) pnorm (q, mean = 0, sd = 1, l

15、ower.tail = TRUE, log.p = FALSE) qnorm (p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) rnorm (n, mean = 0, sd = 1)分布 R语言名称参数Beta分布betashape1, shape2二项式分布binomsize, prob柯西分布cauchylocation, scale分布chisqdf指数分布exprateF分布fdf1, df2Gamma分布gammashape, scale几何分布geomprob逻辑分布logislocation, scale对数正态分布l

16、normmeanlog, sdlog正态分布normmean, sd泊松分布poislambda学生t分布tdf均匀分布unifmin, max威布尔分布weibshape, scale二项分布二项式分布是一种离散概率分布，描述的是在n次独立试验的最终结果。假定每次试验可以有两种结果，要么成功，要么失败。如果一次试验成功的概率为p，在n次独立试验中取得x次成功结果的概率如下，其中假设在一次考试中有10道多元选择题，每道题有4种可能的答案，其中只有一个答案是正确的。如果某个学生以随机方式回答了所有的问题，不及格（答对5道题及以下）的概率为： pbinom (5, size=10, prob=0.

17、25) #返回P(x ppois (11, lambda=5) #默认计算下尾概率(x ppois (11, lambda=5, lower.tail=FALSE) #上尾(x11)1 0.0055连续均匀分布连续均匀分布是在从a到b的连续区间中随机选择数值的概率分布。其概率密度函数的定义如下：例：例如，如果用户想在区间1,5中随机选取10个数，可以调用runif ()函数返回10个随机值： runif (10, min=1, max=5) 1 2.907 3.241 3.983 3.288 2.213 1.992 2.032 3.735 3.002 3.076指数分布指数分布表示了一系列随机

18、重复发生的独立事件的到达时间的分布。假设到下一个事件发生的平均等待时间为，指数分布的概率密度函数形式如下： x dmu colors labels plot (x, dexp(x,1), type=l, lty=1, col=colors1, xlab=x, ylab=概率密度, main=分布比较) for (i in 2:4) lines (x, dexp(x,1/dmui), lwd=2, lty=i, col=colorsi) legend (topright, inset=.05, title=均值, labels, lwd=1, lty=c(1, 2, 3, 4), col=col

19、ors)正态分布正态分布用下列的概率密度函数定义，其中为均值而是方差：不同方差下均值为0的正态分布概率密度函数：假设在一次考试中全班的成绩符合正态分布，平均值是71分，标准差是12.8。计算一下，考试成绩在85分以上的概率是多少？前面已经使用过一些分布的累积概率函数，对于正态分布，这样的函数是pnorm ().可以通过设置参数lower.tail=FALSE，也就是计算上尾概率求出所需结果。 pnorm (85, mean=71, sd=12.8, lower.tail=FALSE) #85分以上，上尾累积概率1 0.137学生t分布统计假设检验统计假设检验举例来看，假设有一个制造商宣称所生产

20、的灯泡平均寿命高达10000小时。检验时在30个灯泡样本中，发现其平均寿命只有9900小时。假定知道样本总体的标准差是120小时。给定5%的显著性水平，那么是否应该接受还是拒绝制造商的这一说法？。 xbar mu0 sigma n z alpha z.alpha z.alpha #临界值 1 -1.644854根据计算结果，检验统计量-4.5644小于临界值-1.6449，因此，在5%的显著性水平上应该拒绝平均寿命超过10000小时的宣传。统计假设检验统计假设检验假设有一种食品商在曲奇包装上标签称每一块曲奇最多含有2克的饱和脂肪酸。在共计35块的样本中，发现平均的饱和脂肪酸含量为2.1克。假定

21、已知总体标准差为0.25克，那么在5%的显著性水平上，是否应该拒绝或接受食品的标签？。 xbar mu0 sigma n z alpha z.alpha z.alpha #临界值 1 1.644854检验统计量2.367大于临界值1.645，因此在5%的显著性水平，可以拒绝假设每一块曲奇只含有2克的饱和脂肪酸。统计假设检验统计假设检验在南极洲发现了一群国王企鹅，去年它们的平均体重时15.4公斤。假设在今年的35只样本中，测量到的平均体重只有14.6公斤。如果已知总体标准差时2.5公斤。那么在5%的显著性水平下，能否拒绝企鹅平均体重与去年相同的假设？。 xbar mu0 sigma n z al

22、pha z.half.alpha c (z.half.alpha, z.half.alpha) 1 1.9600 1.9600从计算结果可知，检验统计量-1.8931位于临界值 -1.9600到1.9600之间。因此，在5%的显著性水平，无法拒绝这群企鹅平均体重与去年一样的假设。统计假设检验统计假设检验统计假设检验假设一个制造商宣称一种灯泡的平均使用寿命超过10000小时。如果真实的平均使用寿命只有9900小时，总体方差为120小时。假设现在掌握了30个灯泡的样本，如何计算在5%的显著性水平，犯下第二类错误的概率是多少呢？。 n - 30; sigma sem alpha - .05; mu0

23、 q mu pnorm (q, mean=mu, sd=sem, lower.tail=FALSE) #用SEM作为标准差，求上尾1 0.26196 内容导航CONTENTS数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4简单线性回归简单线性回归 eruption.lm coeffs coeffs #一元线性回归的参数：截距、斜率(Intercept) waiting -1.87401599 0.07562795 plot (eruptions waiting, faithful, #绘图变量col = blue, #绘图参数main = 老忠实线性回归结果, #

24、标题xlab = 等待时间, #x轴标签ylab = 持续喷发时间)#y轴标签 fit abline (fit, col=red) #画出回归模型lm参数：因变量eruptions ，自变量waiting，数据集faithful。使用coefficients ()来显示所得到的回归方程中的系数 。简单线性回归 waiting duration duration (Intercept) 4.1762 newdata predict(eruption.lm, newdata) 1 2 4.176220 1.907381 使用模型进行预测。简单线性回归简单线性回归简单线性回归简单线性回归qq图可用于直观验证一组数据是否来自于某个给定的分布，或者验证两组数据是否来自同一分布。根据所讨论的分布计算出每个数据点的理论预期值，如果数据确实遵循假定的分布，那么在qq图上的点将大致散落在一条直线上。正态概率图就是一种把数据集与正态分布进行比较的图形化工具。例如，可以比较线性回归模型的标准化残差来检验残差是否真正地符合正态分布规律。 qqnorm (eruption.stdres, ylab=标准化残差, xlab=正态得分