- 在线时间
- 514 小时
- 最后登录
- 2023-12-1
- 注册时间
- 2018-7-17
- 听众数
- 15
- 收听数
- 0
- 能力
- 0 分
- 体力
- 39407 点
- 威望
- 0 点
- 阅读权限
- 255
- 积分
- 12518
- 相册
- 0
- 日志
- 0
- 记录
- 0
- 帖子
- 1388
- 主题
- 1158
- 精华
- 0
- 分享
- 0
- 好友
- 15
TA的每日心情 | 开心 2023-7-31 10:17 |
---|
签到天数: 198 天 [LV.7]常住居民III
- 自我介绍
- 数学中国浅夏
|
【高级数理统计R语言学习】2 多元线性回归 一、背景
& G/ W9 X# ^6 C+ s! b" Y" }数据集展示了X市外来人口的相关数据情况,包括出生年月、收入、初次来到X市的日期、迁离X市的日期和现在的朋友数量。现假设外来人口的年龄、在X市的居住时间和朋友数量影响他们的收入。试加以证明。 二、要求和代码 一、分析收入的影响因素
9 P, A1 d$ L6 r$ o# H#10 ]/ s" K; B- s( v- Q
#展示数据集的结构
& Y! i2 M; o! d" b* _) \- ?data2 <- read.csv(file="F:/hxpRlanguage/homework2.csv",header=TRUE,sep=",")5 Z% U. F4 R$ i- Y
str(data2) #显示的结果有一列是多余的,需要删除+ ?! x* h# Y( A8 P' q( D0 I
data2 <- data2[,1:9]/ ^. Q6 p6 M& P! m+ r
str(data2) #删完之后的显示效果是正常的没有多余列
- e$ [# m( u8 V$ e2 Q3 q
; V. D! V( ?4 N9 t- [#2
4 {3 C6 u' A+ y" [* D( N! }+ O#显示前10条数据记录
9 |( Z: q `. s# s$ ^0 I, gdata2[1:10,]
^! N$ G7 J' o5 }! s) L9 `( m8 v8 Y
#3
7 g3 @, v% Q8 z, X. S2 s#将变量名重新命名为英文变量名
1 ^& g" o( N `) scnames <- c("number","birthyear","birthmonth","salary","inyear","inmonth","outyear","outmonth","friends")
, J2 i4 z" B5 Bcolnames(data2) <- cnames
5 u' H3 S& B* u& ~$ k$ @View(data2)
@* K, p4 E) g' {$ s- f% o/ A- \% D6 ^ ~& @# _
#4
! q. ~7 A' X5 @9 Y" U& L8 b#查找数据集中居住时间小于等于0的异常记录,若存在,从数据集中删除这些异常记录
6 B* H- ^; P$ b+ d2 n% j6 `* }: ax2 <- ((data2$outyear-data2$inyear)*12+(data2$outmonth-data2$inmonth)), D6 b. M& A7 m6 v7 z
#View(x2) #①先算出居住时间
7 s; {5 w7 T( s+ ~4 Q" pdata3 <- cbind(data2,x2)7 [! I2 F% p. v% D. I
#View(data3) #②使用cbind函数把x2和原数据拼成新的矩阵,方便之后删除异常数据列,并且是127条
8 A+ g" L8 i2 r, R9 klist <- which(x2<=0)
; M: T4 l+ W c% U3 ?data3 <- data3[-list,]
; ?9 |( s* q7 o3 \/ M1 y4 g/ jView(data3) #删除异常数据后是125条数据
9 g) _- k( \- u' `. }! N
6 A" N6 A; a3 t- u4 o, U#5
, x# o4 ~2 K" }/ x' c#展示数据集中因变量与自变量的均值、最小值、中位数、最大值和标准差,要求保留2位小数。0 A% U4 q2 i8 y' K" A) V
library(lubridate)
, l: q7 f; F& }" z- y$ Mdate<-Sys.Date() #返回系统当前的时间$ J. D' ~* l K- K# W" r8 l
nowyear<-year(date) #提取年份$ ^4 @+ h4 \/ O: @6 e$ \# A. `
nowmonth<-month(date) #提取月份
/ n4 ^+ k1 Q8 X! a#View(date) #查看现在的日期6 ?/ `3 _& ]1 q6 ^ e: a2 I0 F2 ~
#View(month(date)) #查看现在日期中的月份) P b; R; \4 m
x1 <- array(1:nrow(data3),dim=c(nrow(data3),1))
4 _8 I W5 T( v$ Lfor(i in c(1:nrow(data3)) ){ q) ]+ B- H* s- J# s: `4 V) h
if(nowmonth-data3[i,"birthmonth"]<0){
- n: ^) E8 p' `/ n% O% B x1[i,1] <- nowyear-data3[i,"birthyear"]-1
) f0 O1 n# Y7 V3 K: t$ x) ]4 _/ R }else{. z3 ~: N2 _6 W7 n: r( e
x1[i,1] <- nowyear-data3[i,"birthyear"]. a$ n8 w, N+ Q0 L1 v8 H9 |; _
}
0 a" R$ n1 j% B/ P5 q3 A9 T: R- j6 ?}( }' s/ T/ j2 Q8 D
#View(x1) #算出年龄x1,并加入到数据表中
* B1 V- X3 U" u i% |data4 <- cbind(data3,x1)
# u r8 ]& w, J- Q0 z' E! i& zView(data4) #加入x1年龄变量的新表展示( r/ `8 |3 C7 V. f
x2 <- data4$x2" x V7 m1 v l+ l9 y5 V
Mean.x2 <- round(mean(x2),2)1 h. U5 R% m. _& v) w
Min.x2 <- round(min(x2),2)2 Q p; g, l# _. u. y7 z
Max.x2 <- round(max(x2),2)- U5 T( F2 c6 [9 |# o8 E) S5 g/ x' D
Median.x2 <- round(median(x2),2)* Y# L/ F* ^1 h! ?& f
Sd.x2 <- round(sd(x2),2)
9 H. P1 @5 H9 w8 J" M5 _+ Scbind(Mean.x2,Min.x2,Max.x2,Median.x2,Sd.x2) #x2居住时间的相关结果
4 \2 Z1 g) a4 v/ j) HMean.x1 <- round(mean(x1),2)+ E. ~/ C9 l, [- b/ J) ?) T
Min.x1 <- round(min(x1),2)8 i- ~" {, o5 }. |/ j
Max.x1 <- round(max(x1),2)
, a5 C0 ^8 |! H% YMedian.x1 <- round(median(x1),2)4 j% ~7 B& C0 ]. h; K- b4 j
Sd.x1 <- round(sd(x1),2)& _- ?% Q f6 J+ A! z
cbind(Mean.x1,Min.x1,Max.x1,Median.x1,Sd.x1) #x1年龄的相关结果1 q! o+ r; R6 ]3 N) J5 H
x3 <- data4$friends3 _, V: j) j( X; `7 M
Mean.x3 <- round(mean(x3),2)1 I& ]0 V$ w; Z1 \: n S
Min.x3 <- round(min(x3),2)
T) j2 |( B; `, ?0 JMax.x3 <- round(max(x3),2)/ P/ D! c) \4 K2 Q4 A1 i- @
Median.x3 <- round(median(x3),2)
9 n2 H* Y7 i. nSd.x3 <- round(sd(x3),2)! I3 k$ g, K5 ?/ ?
cbind(Mean.x3,Min.x3,Max.x3,Median.x3,Sd.x3) #x3朋友数量的相关结果
# p3 o& X: \7 _6 @- Ey <- data4$salary
3 M/ d8 k" T3 S% |Mean.y <- round(mean(y),2)) X9 _* V; Z8 \( A# @9 `4 j
Min.y <- round(min(y),2)2 W! k0 Z- ]8 ]7 L( K' s: q2 z4 n
Max.y <- round(max(y),2)5 L) D4 |' y9 [& S4 G/ f! q
Median.y <- round(median(y),2)8 N# f! R$ s6 R% x5 ]! p
Sd.y <- round(sd(y),2)
/ Y% n/ ]/ q" U7 scbind(Mean.y,Min.y,Max.y,Median.y,Sd.y) #因变量y的相关结果6 D) `- \0 K2 H) b7 v, _7 m
7 x' ?/ m& `6 R! y8 P: N2 M- w7 s
#6
9 }) d8 g/ u% x; X; f# \, |' L" N#计算数据集中因变量和自变量的相关系数,要求保留2位小数。
0 x) ~, ^' g" ^. g. V" _: D" {round(cor(y,x1),2) #y和x1年龄. q/ Y" |/ [, n1 ^' @
round(cor(y,x2),2) #y和x2居住时间
) \; s; J; _2 Xround(cor(y,x3),2) #y和x3朋友数量
/ F- f- Z2 u/ f, h$ O+ a
' _! g4 T) ~$ G# E- N3 p#71 C% Y4 Z, Y! _; U
#分别绘制数据集中因变量与各个自变量的散点图% J) W" j' o0 T0 ^
par(mfrow=c(1,3)) #布局,一行画3个图2 l0 K" }* {7 r8 o- E# m j
plot(x1,y,xlab="年龄x1",ylab="工资y")0 M7 j' Z# N5 {7 B1 \% M$ q1 l
plot(x2,y,xlab="居住时间x2",ylab="工资y")
* K/ g3 M0 [. Z9 c& M hplot(x3,y,xlab="朋友数量x3",ylab="工资y")1 z6 r& Z( U9 C4 h! k$ f
8 |+ ]0 k& {& h- z/ \3 E#8- s% l7 Y. n& R9 M% k
#利用多元线性回归模型对数据集中因变量与自变量的关系进行拟合。
0 b5 c! `3 g% G9 V8 ~9 v3 ^lm.xy <- lm(y~x1+x2+x3)
, k! i" |* n' m2 L5 Zlm.xy
+ @, O! ?% L1 ^# F) C- `* ysummary(lm.xy) #得到的结果是方程是显著的具有线性关系,但是每一个系数不都是显著的
( }3 J' N, F8 N! |: S, a0 L; p$ A! N1 x9 f. o; }3 n
#9
% v( T7 a! u0 `1 I6 M& G#对#8中的多元线性回归模型进行诊断,确定异常值记录。
) \2 x$ s$ W2 @) i t# l9 v( z' D6 ~7 \par(mfrow = c(2,2)) #生成四种模型诊断的图形,2行2列
}( Q1 @! {5 u. Z5 i' \#生成四种模型诊断的图形:①残差与真实值的关系图 ②qq图用来检测其残差是否是正态分布* Q$ f# Q- V3 D
#③用来检查等方差假设的。在一开始我们的五大假设第二条便是,我们假设预测的模型里方差是一个定值。4 Y; U! Z! H4 q6 J4 d, e
#如果方差不是一个定值那么这个模型的可靠性也是大打折扣的。# f% s5 }2 I& g1 f
#④Leverage就是杠杆的意思。这种图的意义在于检查数据分析项目中是否有特别极端的点。8 I+ ?. x5 } G- V+ @
plot(lm)$ L6 Y# V; i. h" z
library(carData)5 w8 g3 |% ?' I/ V
library(car): S2 f1 e) ~! L F
outlierTest(lm.xy) #显示离群点,Bonferroni校正,残差最大的点是136号点: c" H0 N8 j- |* p! H/ Q
# ~- e/ m; B( u#10
0 ~: o N2 C) ^#删除异常值记录后重新利用多元线性回归模型拟合数据。
" q) L8 I5 X# s, m1 fdata4 <- data4[-136,] #删除该点% I' n& V0 j1 }
x1 <- data4$x1
2 i) t& }4 J2 mx2 <- data4$x2
; q5 i5 j' R+ ^2 w5 n/ Ax3 <- data4$friends0 @1 Z% O( _1 o" K( Z
y <- data4$salary6 ]% ^% M2 Z1 `' q2 V7 \: K
lm.xy2 <- lm(y~x1+x2+x3) #重新拟合回归模型
5 y0 d+ @, h; L! Zlm.xy2
* V% R" k. D' s9 Y) P
8 M4 N5 F2 o+ S( J% a2 L#11
; }7 Y! J$ V$ e" r( b" N#对#10中的多元线性回归模型进行多重共线性检验,若存在多重共线性,删除相关变量后重新进行拟合。* u. K* K3 @7 v0 o; ~
vif(lm.xy2) #p判断多重共线性0<VIF<10(不存在)
, i9 \ i& O7 l5 L9 S( p8 o M6 |* J3 y' o' N
#123 d# r& a: U9 ^! i1 }% V
#对#11中的结果进行解释,重点分析年龄、在北京的居住时间和朋友数量如何影响收入。
~, t/ b. d" r0 a Ksummary(lm.xy2) #可知年龄和朋友数量对收入有影响,显著性*一颗星
9 Y$ l$ Y8 h1 O6 g6 T7 G S# c1 T7 b- @2 K' H8 K1 \# N; p( Q+ e, I
**********************************************************************6 C6 l- i9 Z5 ?8 s/ d; `1 J- ~4 m0 O/ T
( B1 e; Q3 P* Q, F
二、利用多元线性回归模型预测收入* C ]( ]+ H o9 S3 s- f
View(data4) #124条数据
1 \; }( l2 @9 p0 i#1
Z7 }2 L6 A8 S0 u#从数据集中随机抽取50条记录作为预测集,剩下的数据作为训练集。
/ A$ L2 X+ N+ S! Rtrain0 <- sample(nrow(data4),nrow(data4)-50) #训练集和测试集
) |" ~7 q- F$ d. ]trainData <- data4[train0,] #训练数据6 q* j1 h# W& e+ i% X" j2 Z
testData <- data4[-train0,] #测试数据- z& j3 ?$ U. S* T4 w: I; n
2 y3 K1 @! Q5 |6 e! s' r, Z0 l6 t' ]#2
' j `9 J: A- n! w6 X5 y#针对训练集,利用多元线性回归模型拟合数据。* |0 J. x. o! u1 K$ \
lm.xy3 <- lm(trainData[,"salary"]~trainData[,"x1"]+trainData[,"x2"]+trainData[,"friends"])
' K, o& _9 L8 P" H q8 B" p* w& z4 p# g6 R e/ r# K9 d
#3' D' q8 d- x( d; v) {: s
#对(2)中的多元线性回归模型进行诊断,处理异常值。% D4 y% p" ]8 I8 R" ^
summary(lm.xy3)! ^' }* a# T8 m& u' g+ v Y! W7 @
par(mfrow=c(2,2))
5 ~# K* }; M( m- E' G6 C; V9 xplot(lm.xy3)
7 C2 a* y* T8 [/ M- f' poutlierTest(lm.xy3)
7 t. \; u6 H, g: Q$ xtrainData<-trainData[-c(150,32,82),] #删除异常值,随机的
9 ^* q8 v7 G9 q% v% y! g+ H+ z
#43 x- O+ q, ^3 O! N+ S0 `, l
#对(3)中的多元线性模型进行多重共线性检验并加以处理。: o: ^) A1 F3 R6 y
vif(lm.xy3) #p判断多重共线性0<VIF<10(不存在)' b8 m- ~& m( r) ^+ ]$ U
salary<-trainData[,"salary"] #引入的数据是训练集的数据
`! D7 L" k% q( g( W1 _x2<-trainData[,"x2"]
% t: Q. b o b3 H: p; b( kx1<-trainData[,"x1"]
3 p8 m6 R. g3 J9 g6 ]friends<-trainData[,"friends"]; j6 X/ I2 t8 {, s8 K: D
lm.xy3 <-lm(salary~x2+x1+friends)
+ D) |" s: w! @% @2 ~" |- m; H" r" |
#5& ?. i# w. h4 H1 k( ^
#针对(4)中的模型,分别利用AIC和BIC选择最优模型。
, K$ g2 n7 Y$ m" g) \! _% p#AIC检验,赤池信息准则,选择最小的
. n" f i& R3 d8 CAIClm<-step(lm.xy3,direction="both")
6 ?, S2 u0 K6 P2 _: p7 C#BIC检验,贝叶斯信息准则,选择最小的% d) P8 t, e( n( x
BIClm<-step(lm.xy3,k=log(nrow(trainData)),direction="both")
5 N; i+ y2 b3 S
; Y8 K% C! p: Y" }#6
8 L* I! _+ T& f) w/ }6 H4 g#利用预测集进行预测,比较全模型(包含所有自变量)、AIC选择的最优模型、BIC选择的最优模型! q$ `) }9 z- w1 x( W K% `* f
#这三个模型预测的准确性大小,并进行解释。
# k1 ~6 E* f5 X( R6 F$ RAllmodel<-predict(lm.xy3,testData)
$ c+ I u- w6 J/ g6 QAICmodel<-predict(AIClm,testData)& H7 k. A% g. }4 [8 v
BICmodel<-predict(BIClm,testData)
- y: D. n. u/ h+ X5 _8 k# i! a#均方误差检验,最小最优,分别计算全模型,AIC,BIC的均方误差8 V& x$ f# W- p* p
#均方根误差亦称标准误差,均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根
$ U& D' d$ g( x4 F2 i) d#标准误差能够很好地反映出测量的精密度
1 G1 u6 r# U) i: l0 n I( C' aMSE <- function(x){
5 Y+ |* I" G- k3 B( r2 J mse <- sum((testData[,"salary"]-x)^2)/506 ^. r' L w9 K
return(mse)- Y% \( j& ^7 [0 A# o
}# g# D% d) e+ t! l `1 C$ c
MSE(AICmodel) #AIC/BIC/ALL是误差最小的 m" M9 y4 r: B- R2 S, q3 B/ g
MSE(BICmodel)
/ V; {$ b7 l$ f; f9 P4 [; YMSE(Allmodel)
- [, O2 y& b2 a5 r8 c% a' q3 y4 k% o- Z9 Y: e+ X0 \: T. |' R
+ `1 O. ~% l* T3 B8 L' z: ~' A
4 y; S% l9 H5 S: J2 H6 {
|
zan
|