【R实验.7】回归分析
解法并不单一，下列方法带有璇子个人的偏好，因此仅供参考。如有错误，欢迎在评论区斧正！
7.1 为了估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X 与当年灌溉面积Y，测得连续10 年的数据如表7-6 所示：
（1）试画出相应的散点图，判断Y 与X 是否有线性关系；
> x <- data$X.米.
> y <- data$Y.公顷.
> plot(x,y)
由散点图，x与y具有一定的线性关系
（2）求出Y 关于X 的一元线性回归方程；
> fit1 <- lm(y~x)
> summary(fit1)
Call:
lm(formula = y ~ x)Residuals:Min 1Q Median 3Q Max
-128.591 -70.978 -3.727 49.263 167.228 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 140.95 125.11 1.127 0.293
x 364.18 19.26 18.908 6.33e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 96.42 on 8 degrees of freedom
Multiple R-squared: 0.9781, Adjusted R-squared: 0.9754
F-statistic: 357.5 on 1 and 8 DF, p-value: 6.33e-08
由fit1的描述统计结果，得y关于x的一元线性回归方程为：
Y=364.18∗x+140.95Y = 364.18*x + 140.95Y=364.18∗x+140.95
（3）对方程作显著性检验；
由结果：x系数对应t统计量P值<<0.05，F统计量对应P值小于0.05，即t检验与F检验均通过。同时修正后的拟合优度为0.9754，方程的解释较强，故该模型通过显著性检验。
（4） 现测得今年的数据是X=7 米，给出今年灌溉面积的预测值和相应的区间估计（alpha=0.05）
> predict(fit1,data.frame(x=7),interval='confidence')fit lwr upr
1 2690.227 2613.35 2767.105
得到今年灌溉面积的预测值为2690.227，相应的区间估计为[2613.35, 2767.105]
7.2 研究同一地区土壤所含可给态磷的情况，得到18 组数据如表7-7 所示，表中X1 为土壤内所含无机磷浓度，X2 为土壤内溶于K2CO3 溶液并受溴化物水解的有机磷，X3 为土壤内溶于K2CO3 溶液但不溶于溴化物水解的有机磷。
（1）求出Y 关于X 的多元线性回归方程
> data7.2 <- read.table('data7.2.txt',header = T)
> fit2 <- lm(Y~X1+X2+X3,data=data7.2)
> summary(fit2)
Call:
lm(formula = Y ~ X1 + X2 + X3, data = data7.2)Residuals:Min 1Q Median 3Q Max
-28.349 -11.383 -2.659 12.095 48.807 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.65007 18.05442 2.418 0.02984 *
X1 1.78534 0.53977 3.308 0.00518 **
X2 -0.08329 0.42037 -0.198 0.84579
X3 0.16102 0.11158 1.443 0.17098
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 19.97 on 14 degrees of freedom
Multiple R-squared: 0.5493, Adjusted R-squared: 0.4527
F-statistic: 5.688 on 3 and 14 DF, p-value: 0.009227
由回归结果,y关于x的多元线性回归方程为：
Y=1.78534x1–0.08329x2+0.16102x3+43.65007Y = 1.78534x1 – 0.08329x2 + 0.16102x3 + 43.65007Y=1.78534x1–0.08329x2+0.16102x3+43.65007
（2）对方程做显著性检验
易知x1,x2,x3系数的t统计量对应P值，有且仅有x1系数的对应P值<0.05显著，x2,x3的系数不显著。同时，虽然F统计量对应P值<0.05，但方程修正后的拟合优度只有0.4527，即方程的解释力度较弱。因此，我们还需进行逐步回归分析，来进一步筛选自变量，进而提升方程的解释力度。
（3）对变量做逐步回归分析
> library(MASS)
> stepAIC (fit2, direction = "backward")#向后逐步回归
Start: AIC=111.27
Y ~ X1 + X2 + X3Df Sum of Sq RSS AIC
- X2 1 15.7 5599.4 109.32
<none> 5583.7 111.27
- X3 1 830.6 6414.4 111.77
- X1 1 4363.4 9947.2 119.66Step: AIC=109.32
Y ~ X1 + X3Df Sum of Sq RSS AIC
<none> 5599.4 109.32
- X3 1 833.2 6432.6 109.82
- X1 1 5169.5 10768.9 119.09Call:
lm(formula = Y ~ X1 + X3, data = data7.2)Coefficients:
(Intercept) X1 X3
41.4794 1.7374 0.1548
开始时模型包含3个（全部）预测变量，然后每一步中，AIC 列提供了删除一个行中变量后模型的AIC 值，中的AIC 值表示没有变量被删除时模型的AIC。第一步，x2被删除，AIC 从111.27降低到109.32；然后再删除变量将会增加AIC，因此终止选择过程。
向后逐步回归的最终模型为：
Y=1.7374x1+0.1548x3+41.4794Y = 1.7374x1 + 0.1548x3 + 41.4794Y=1.7374x1+0.1548x3+41.4794
逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，因为不是每一个可能的模型都被评价了。为克服这个限制，故我们用全子集回归法来验证一下结果。
library(leaps)
> data1 <- as.data.frame(data7.2[,c("Y", "X1","X2", "X3")])
> leaps <-regsubsets(Y ~ X1 + X2 + X3, data=data1, nbest=4)
> plot(leaps, scale="adjr2")
图形表明，双预测变量模型（x1和x3）是最佳的模型，与向后逐步回归的结果一致。
7.3 已知如下数据，由表7-8 所示：
（1）画出数据的散点图，求回归直线y=β0^+β1^xy= \widehat{\beta_0}+\widehat{\beta_1}xy=β0+β1x，同时将回归直线也画在散点图。
> data7.3 <- read.table("data7.3.txt",header = T)
> plot(data7.3$X,data7.3$Y)
> fit3 <- lm(Y~X,data=data7.3)
> abline(fit3)
>
> summary(fit3)
Call:
lm(formula = Y ~ X, data = data7.3)Residuals:Min 1Q Median 3Q Max
-7.5173 -1.7928 -0.0147 1.1396 16.6652 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.7228 1.7066 -1.010 0.322
X 1.7175 0.2778 6.184 1.54e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 4.86 on 26 degrees of freedom
Multiple R-squared: 0.5952, Adjusted R-squared: 0.5797
F-statistic: 38.24 on 1 and 26 DF, p-value: 1.536e-06
回归模型为：y=1.7175x–1.7228y = 1.7175x – 1.7228y=1.7175x–1.7228
（2）分析T 检验和F 检验是否通过
Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.7228 1.7066 -1.010 0.322
X 1.7175 0.2778 6.184 1.54e-06 ***
F-statistic: 38.24 on 1 and 26 DF, p-value: 1.536e-06
由上面的描述统计，x的t统计量和方程的F统计量对应P值均<<0.05，故T检验和F检验通过。
（3）画出残差（普通残差和标准化残差）与预测值的残差图，分析误差是否是等方差
par(mfrow=c(2,2))
> plot(fit3)
在“残差图与拟合图”（Residuals vs Fitted）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0 的正态分布。正态Q-Q 图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率。显然，图像表明残差违反了正态性的假设。
若满足不变方差假设，那么在位置尺度图（Scale-Location Graph）中， 水平线周围的点应该随机分布。由生成图，这里似乎不满足此假设。（略，其和上面的图是一起生成的）
（4）修正模型，对响应变量Y 作开方，再完成(1)-(3)的工作
> plot(data7.3$X,sqrt(data7.3$Y))
> fit4 <- lm(sqrt(Y)~X,data=data7.3)
> abline(fit4)summary(fit4)
Call:
lm(formula = sqrt(Y) ~ X, data = data7.3)Residuals:Min 1Q Median 3Q Max
-1.13695 -0.42839 0.01512 0.29452 1.94385 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.76627 0.24422 3.138 0.0042 **
X 0.31150 0.03975 7.837 2.6e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.6954 on 26 degrees of freedom
Multiple R-squared: 0.7026, Adjusted R-squared: 0.6911
F-statistic: 61.41 on 1 and 26 DF, p-value: 2.598e-08
回归模型为：y(1/2)=0.3115x+0.76627y^(1/2) = 0.3115x + 0.76627y(1/2)=0.3115x+0.76627
与前一问原理相同，此时模型通过t检验与F检验，同时拟合优度有了一定的提升。
> par(mfrow=c(2,2))
> plot(fit4)
正态Q-Q 图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率。显然，图像表明残差大致满足了正态性的假设。
而在位置尺度图（Scale-Location Graph）中， 水平线周围的点应该随机分布。该图似乎勉强满足同方差假设。
7.4 一位饮食公司的分析人员想调查自助餐馆中的自动咖啡售货机数量与咖啡销售量之间的关系，她选择了14 家餐馆来进行实验，这14家餐馆在营业额、顾客类型和地理位置方面都是相近的，放在试验餐馆的自动售货机数量从0（这里由咖啡服务员端来）到6 不等，并且是随机分配到每个餐馆的，表7-9 所示的是关于试验结果的数据。
（1）作线性回归模型
data7.4 <- read.table('data7.4.txt',header = T)
> y <- data7.4$售货机数量
> x <- data7.4$咖啡销售量
> fit5 <- lm(y~x)
> summary(fit5)Call:
lm(formula = y ~ x)Residuals:Min 1Q Median 3Q Max
-0.44423 -0.23923 0.06288 0.22868 0.44527 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.2722554 0.5322309 -17.42 6.94e-10 ***
x 0.0178252 0.0007632 23.36 2.26e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.3169 on 12 degrees of freedom
Multiple R-squared: 0.9785, Adjusted R-squared: 0.9767
F-statistic: 545.5 on 1 and 12 DF, p-value: 2.265e-11
由统计结果，模型通过t检验和F检验，且具有较好的解释力度，回归方程式为：
Y=0.0178252∗x–9.2722554Y = 0.0178252*x – 9.2722554Y=0.0178252∗x–9.2722554
（2）作多项式回归模型
fit6 <- lm(y ~ x + I(x^2))
> summary(fit6)Call:
lm(formula = y ~ x + I(x^2))Residuals:Min 1Q Median 3Q Max
-0.17759 -0.12396 -0.02883 0.10082 0.20601 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.535e+00 1.622e+00 0.946 0.36425
x -1.542e-02 4.945e-03 -3.118 0.00978 **
I(x^2) 2.484e-05 3.685e-06 6.740 3.2e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.1461 on 11 degrees of freedom
Multiple R-squared: 0.9958, Adjusted R-squared: 0.995
F-statistic: 1305 on 2 and 11 DF, p-value: 8.422e-14
同理，模型均通过了t检验与F检验，且具有较好的解释力度。模型表达式为：
Y=−0.01542x+0.00002484(x2)+1.535Y = -0.01542x + 0.00002484(x^2) + 1.535Y=−0.01542x+0.00002484(x2)+1.535
（3）画出数据的散点图和拟合曲线
plot(x,y)
abline(fit5)
plot(x,y)
lines(x, fitted(fit6))
7.5 一位医院管理人员想建立一个回归模型，对重伤病人出院后的长期恢复情况进行预测，自变量是病人住院的天数（X），应变量是病人出院后长期恢复的预后指数（Y），指数的数值越大表示预后结局越好，为此，研究了15 个别病人的数据，这些数据列在表7-10 中，根据经验表明，病人住院的天数（X）和预后指数（Y）服从非线性模型
Yi=θ0*exp(θ1Xi)+εi, i=1,…, 15
（1）用内在线性模型方法计算其各种参的估计值；
> data <- read.table('data7.5.txt',header=T)
> x <- data$X; y <- data$Y; lny <- log(y)
> fit <- lm(lny ~ x)
> summary(fit)
Call:
lm(formula = lny ~ x)
Residuals:Min 1Q Median 3Q Max
-0.37241 -0.07073 0.02777 0.05982 0.33539
Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.037159 0.084103 48.00 5.08e-16 ***
x -0.037974 0.002284 -16.62 3.86e-10 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1794 on 13 degrees of freedom
Multiple R-squared: 0.9551, Adjusted R-squared: 0.9516
F-statistic: 276.4 on 1 and 13 DF, p-value: 3.858e-10
plot(x,y)
> y1 <- exp(4.037159-0.037974*x)
> lines(x,y1)
（2）用非线性方法（nls()函数和nlm()函数）计算其各种参的估计值
pop.mod1 <- nls(y ~ theta1*exp(theta2*x),start=list(theta1=56,theta2=-0.038),trace=T)
62.32926 : 56.000 -0.038
49.4633 : 58.559055 -0.039553
49.4593 : 58.60589607 -0.03958537
49.4593 : 58.60654476 -0.03958642> summary(pop.mod1)
Formula: y ~ theta1 * exp(theta2 * x)
Parameters:Estimate Std. Error t value Pr(>|t|)
theta1 58.606545 1.472160 39.81 5.70e-15 ***
theta2 -0.039586 0.001711 -23.13 6.01e-12 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.951 on 13 degrees of freedom
Number of iterations to convergence: 3
Achieved convergence tolerance: 5.558e-06
判断拟合效果
> library(ggplot2)
> d <- data.frame(x,y)
> p <- ggplot(d,aes(y, x))
> p+geom_point(size=3)+geom_line(aes(x,fitted(pop.mod1)),col='red')
由图，拟合效果还可以，但明显没有使用内在线性模型方法得到的拟合曲线与原数据更贴近。
残差诊断
为了检测这些假设是否成立我们用拟合模型的残差来代替误差进行判断
plot(fitted(pop.mod1) , resid(pop.mod1),type='b')
fitted是拟合值(predict是预测值) resid和residuals表示残差
【R实验.7】回归分析相关推荐
- SAS实验04 ——回归分析
实验04 回归分析 一.实验目的 通过实验进行对回归分析的学习,并有效掌握回归分析数据样本的解读和整理并从SAS输出结果中得到相关结论 二.实验内容 ①我近些日子复习英语单词的个数和每天的单词学习时间 ...
- 《数据分析实战》--用R做多元回归分析
<数据分析实战>--用R做多元回归分析 本文参考的是<数据分析实战>的第六章. 背景:针对某公司对产品的不同广告平台投放,基于过去的新增用户数据和投放数据,希望获得更好的广告投 ...
- R语言 线性回归分析实例
y,X1,X2,X3 分别表示第 t 年各项税收收入(亿元),某国生产总值GDP(亿元),财政支出(亿元)和商品零售价格指数(%). (1) 建立线性模型: ① 自己编写函数: > librar ...
- R语言与回归分析计算实例
6.1.7 计算实例 这里用Forbes数据为例,全面展示一元回归模型的计算过程. 例 6.5 Forbes数据 在十九世纪四.五十年代,苏格兰物理学家James D. Forbes,试图通过水的沸点 ...
- 实验4-EXCEL回归分析
EXCEL回归分析 通过数据间的相关性,我们可以进一步构建回归函数关系,即回归模型,预测数据未来的发展趋势. 相关分析与回归分析的联系是:均为研究及测量两个或两个以上变量之间关系的方法.在实际工作中, ...
- matlab 定义一个有自变量的方程_Eviews、Stata、Python、Matlab、R描述+相关+回归分析教程汇总...
1描述统计分析简介 基本统计分析,又叫描述性统计分析,描述性统计主要包括数据的集中趋勢分析.数据的离散程度分析.频数分布分析等. 通常对收集来的数据进行直接的频率.频数等描述,描述性统计分析一般对样本 ...
- r语言中残差与回归值的残差图_用R语言做回归分析_iris数据集/longley数据集
机器学习课程2 回归分析 [题目1] 使用R对内置鸢尾花数据集iris(在R提示符下输入iris回车可看到内容)进行回归分析,自行选择因变量和自变量,注意Species这个分类变量的处理方法. 解答 ...
- 菜鸟学R语言（回归分析1）
学会用R做回归分析 在统计学中,回归分析(regression analysis)指的是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法.回归分析按照涉及的变量的多少,分为一元回归和多元回 ...
- 太原理工大学linux与python编程r实验报告_太原理工大学算法设计与分析实验报告...
<太原理工大学算法设计与分析实验报告>由会员分享,可在线阅读,更多相关<太原理工大学算法设计与分析实验报告(12页珍藏版)>请在人人文库网上搜索. 1.本科实验报告课程名称: ...
最新文章
- 都说百度前端牛，来看看百度前端工程化之H5性能优化
- 服务器告警其一：硬盘raid问题
- MySQL中的information_schema
- jq mysql二级联动_jq+php+mysql 实现二级菜单联动
- Linux I2C核心、总线与设备驱动
- python mysql异常处理_python-处理PyMySql异常-最佳做法
- MQTT教程（二）：MQTT中的可变报头
- Axure授权码，2021年11月11日亲测有效
- 1.0 信息化与信息系统
- 三调与二调图斑叠加分析，筛选不同地类面积占比，筛选举证图斑
- bzoj 4082: [Wf2014]Surveillance 倍增
- 如何使用计算机计算平方面积,尺平方米换算计算器(面积单位换算器)
- 网站友情链接交换的方法
- docker 入门优质文章
- android触屏对焦_Android相机开发(五): 触摸对焦,触摸测光,二指手势缩放
- 20篇高质量程序人生文章分享，做开发不仅仅只有代码
- 易语言html5播放器问题,易语言媒体播放器 - 已处理问题存放区 - 中国红客联盟 - Powered by HUC...
- 五邑大学计算机学院白明,五邑大学政法学院来访我院
- 输入直角三角形的两个直角边，求三角形的周长和面积，以及两个锐角的度数。结果均保留一位小数。
- 电阻电位器位移电子尺信号隔离变送器