- 在线时间
- 514 小时
- 最后登录
- 2023-12-1
- 注册时间
- 2018-7-17
- 听众数
- 15
- 收听数
- 0
- 能力
- 0 分
- 体力
- 39380 点
- 威望
- 0 点
- 阅读权限
- 255
- 积分
- 12510
- 相册
- 0
- 日志
- 0
- 记录
- 0
- 帖子
- 1388
- 主题
- 1158
- 精华
- 0
- 分享
- 0
- 好友
- 15
TA的每日心情 | 开心 2023-7-31 10:17 |
---|
签到天数: 198 天 [LV.7]常住居民III
- 自我介绍
- 数学中国浅夏
|
【高级数理统计R语言学习】2 多元线性回归 一、背景
C8 m8 w. X. I数据集展示了X市外来人口的相关数据情况,包括出生年月、收入、初次来到X市的日期、迁离X市的日期和现在的朋友数量。现假设外来人口的年龄、在X市的居住时间和朋友数量影响他们的收入。试加以证明。 二、要求和代码 一、分析收入的影响因素
* R% V8 ?8 S4 X5 X( @$ b#1
3 J2 _$ F( [& \1 h- P9 \. B/ f#展示数据集的结构
/ X- }) ]( ^) p5 c. u- bdata2 <- read.csv(file="F:/hxpRlanguage/homework2.csv",header=TRUE,sep=","), T& U/ m# D: J& s9 h
str(data2) #显示的结果有一列是多余的,需要删除
# K# {6 o+ ?2 a p9 @( idata2 <- data2[,1:9]
! {6 ^5 A3 ~0 N% C+ Q/ |str(data2) #删完之后的显示效果是正常的没有多余列: c! e, d7 L! s. k
8 C$ v, C/ `% z/ O9 y2 u
#2$ f7 u0 _% u2 v( \
#显示前10条数据记录! x8 D2 V8 J- ^+ U
data2[1:10,]
! m( ~, |) ~( C3 {/ w9 |8 G) s/ t/ w6 n9 \# K/ ]; e" ~ E
#30 V2 @: h1 w. m8 Z7 @4 J
#将变量名重新命名为英文变量名
$ Q9 s( e/ h5 n' q# K/ ncnames <- c("number","birthyear","birthmonth","salary","inyear","inmonth","outyear","outmonth","friends")' h3 ^2 i4 v6 q/ m7 Y* o4 e& n: y
colnames(data2) <- cnames
& V' Q0 Y N; \5 k# }& tView(data2) i. K* x3 G9 m( N" ]
: H: t# W2 s% @$ ]4 W/ d! L
#49 T( ]3 c1 q( o' L
#查找数据集中居住时间小于等于0的异常记录,若存在,从数据集中删除这些异常记录+ _+ x2 ]" ^5 ?, W9 \3 M' K
x2 <- ((data2$outyear-data2$inyear)*12+(data2$outmonth-data2$inmonth))2 q, g W& ^$ g5 t5 R1 W4 S2 G' T
#View(x2) #①先算出居住时间
2 E7 E( q/ P2 W$ }. k2 `- mdata3 <- cbind(data2,x2)
6 T, d, A- u( K8 r: i; \# y#View(data3) #②使用cbind函数把x2和原数据拼成新的矩阵,方便之后删除异常数据列,并且是127条
2 V% O' w# K) D) \' H: o% r4 T. ulist <- which(x2<=0)5 i+ C4 m$ O! V7 L0 D
data3 <- data3[-list,]
0 t$ F5 G6 @# x, N$ EView(data3) #删除异常数据后是125条数据
! G8 q' f% D% s Q# m, N+ S* M- p0 ?$ P+ m5 H
#5$ S) c# n( _' x: o% e
#展示数据集中因变量与自变量的均值、最小值、中位数、最大值和标准差,要求保留2位小数。
5 K7 s0 O; o5 T- }/ ?- ]library(lubridate) f. V: |" H+ D; i
date<-Sys.Date() #返回系统当前的时间
2 ]0 d# o/ s+ u; ~! Dnowyear<-year(date) #提取年份: R6 y6 P. _$ v4 g* @
nowmonth<-month(date) #提取月份$ U- f* _& Q* c* x y
#View(date) #查看现在的日期
# q" ~4 z4 g; L* I" ]& c- W0 `#View(month(date)) #查看现在日期中的月份2 ?; v- r) |5 v
x1 <- array(1:nrow(data3),dim=c(nrow(data3),1))
- ~% {5 ]) X! j1 [1 f# Ffor(i in c(1:nrow(data3)) ){
@2 E; d# u% m3 V, d if(nowmonth-data3[i,"birthmonth"]<0){0 X) a% f: y3 f$ B p1 t9 z
x1[i,1] <- nowyear-data3[i,"birthyear"]-1
* D0 X0 V& Q w3 Q. \5 Y) T- C' y }else{2 ^$ ~$ Z$ u% l4 l, H! x
x1[i,1] <- nowyear-data3[i,"birthyear"]
" N5 J6 E" M( G }. m+ q2 w2 l. }6 I; e6 q# I
}" F& O9 C6 z# P5 p
#View(x1) #算出年龄x1,并加入到数据表中 ^: h1 {2 b8 a- N, w0 ^
data4 <- cbind(data3,x1) 7 D. ^. T" d9 x' k% m) r
View(data4) #加入x1年龄变量的新表展示4 |0 ~" R1 F! q3 ^3 X
x2 <- data4$x2: Z2 q8 C; F1 F" T D% i
Mean.x2 <- round(mean(x2),2)
5 B$ F2 F# _' w$ mMin.x2 <- round(min(x2),2)" b5 l- b, M/ J* ?
Max.x2 <- round(max(x2),2)
; v5 _8 a5 W( A; B: {3 UMedian.x2 <- round(median(x2),2)
7 |% T) V/ q: M+ _! oSd.x2 <- round(sd(x2),2)
f, |7 f; O' P Lcbind(Mean.x2,Min.x2,Max.x2,Median.x2,Sd.x2) #x2居住时间的相关结果
x0 V) Q& y0 [2 F! LMean.x1 <- round(mean(x1),2)
' R: R) M; i0 G( a7 BMin.x1 <- round(min(x1),2)1 ` c) x# \! x2 Z2 E( K
Max.x1 <- round(max(x1),2)! Y$ e2 S- j) _5 M/ L
Median.x1 <- round(median(x1),2)1 s& j1 \( M+ E6 F d
Sd.x1 <- round(sd(x1),2)
" |: x0 C6 ~. U, Y! _$ x( p/ rcbind(Mean.x1,Min.x1,Max.x1,Median.x1,Sd.x1) #x1年龄的相关结果
* ~) J. ?5 }5 [x3 <- data4$friends8 b' F8 m* ?0 o' C1 ?1 f
Mean.x3 <- round(mean(x3),2)
* |8 B' \0 i( H+ }2 k7 X+ O; W& R/ }Min.x3 <- round(min(x3),2)8 n: j( W1 S% V _# J
Max.x3 <- round(max(x3),2)( }- M1 ?# a+ ]7 {2 H b
Median.x3 <- round(median(x3),2)
' u( q( i/ J% t+ p9 ^% }8 O" Y$ iSd.x3 <- round(sd(x3),2)6 v! J% `" U- {/ r6 ^) A
cbind(Mean.x3,Min.x3,Max.x3,Median.x3,Sd.x3) #x3朋友数量的相关结果
4 i9 O7 n% }0 |/ ~y <- data4$salary
3 @- A k6 B1 S4 lMean.y <- round(mean(y),2), g3 C& l: x3 x- n7 T
Min.y <- round(min(y),2)
: h( w6 i7 F& lMax.y <- round(max(y),2)# t, A+ r/ T; j
Median.y <- round(median(y),2)6 C3 e( @9 }/ a" `' I; P0 P7 ?- y3 P+ R$ y
Sd.y <- round(sd(y),2): F( a7 B p9 r1 [1 B
cbind(Mean.y,Min.y,Max.y,Median.y,Sd.y) #因变量y的相关结果# X. F7 C& G o/ S- R5 y
! w$ S4 \$ }( }9 N4 B3 E#66 _! J9 W* w3 H. i8 R
#计算数据集中因变量和自变量的相关系数,要求保留2位小数。! E) d. W8 s( h& S: B9 C" k c
round(cor(y,x1),2) #y和x1年龄1 }5 V4 d! I8 \9 C( x) K- g
round(cor(y,x2),2) #y和x2居住时间: h B2 Q G9 f2 W( l
round(cor(y,x3),2) #y和x3朋友数量
8 a: a' s% h0 T- A; N2 i5 [& L5 w/ u6 J9 f: I
#7
- \# O5 T5 \9 m#分别绘制数据集中因变量与各个自变量的散点图' _9 f& Q* @7 U8 e: H
par(mfrow=c(1,3)) #布局,一行画3个图
5 ]7 v5 X) n/ z/ b/ T! ^plot(x1,y,xlab="年龄x1",ylab="工资y")
' |$ q ~/ f6 \# G0 pplot(x2,y,xlab="居住时间x2",ylab="工资y")) P1 N X; q- A- W; e, u, Z: X
plot(x3,y,xlab="朋友数量x3",ylab="工资y")) K$ S( m* Z T3 P) R
& c/ C: c9 m# q4 _#8& G/ z$ E" e4 {' Z( L' D/ M% \
#利用多元线性回归模型对数据集中因变量与自变量的关系进行拟合。
' M# e h; ^( Plm.xy <- lm(y~x1+x2+x3)8 g, h2 s7 q2 Z* F. N
lm.xy
* G8 Y1 u+ f6 p2 ]+ ~, ?- ~summary(lm.xy) #得到的结果是方程是显著的具有线性关系,但是每一个系数不都是显著的
: a, u) j" T0 R; `- @1 B! [0 e8 B1 R
#9
6 `, \5 K/ W. f* s#对#8中的多元线性回归模型进行诊断,确定异常值记录。
# g- ^$ B+ z# c8 k* _( V" [) xpar(mfrow = c(2,2)) #生成四种模型诊断的图形,2行2列
- z- ?- ~" ]: x8 S9 a; W p#生成四种模型诊断的图形:①残差与真实值的关系图 ②qq图用来检测其残差是否是正态分布9 K% T. P6 j( l R; [0 d3 R) X
#③用来检查等方差假设的。在一开始我们的五大假设第二条便是,我们假设预测的模型里方差是一个定值。
* h7 F4 w' v. i N#如果方差不是一个定值那么这个模型的可靠性也是大打折扣的。. {0 z4 v' _ m' i) v
#④Leverage就是杠杆的意思。这种图的意义在于检查数据分析项目中是否有特别极端的点。
& h& s) V% x+ |3 _plot(lm)
! v/ i! D6 E! @9 U' dlibrary(carData)
: j( j" }$ [% d4 P6 Xlibrary(car)
/ M) `" G! d1 k6 i" Y9 J3 T; {4 [outlierTest(lm.xy) #显示离群点,Bonferroni校正,残差最大的点是136号点% Z2 U& ]4 G( Q3 w3 D$ m( o
& }8 k/ o8 [+ r2 O# P2 I d
#10. F" ?# W) t/ p2 ?' G+ X' b) c1 H
#删除异常值记录后重新利用多元线性回归模型拟合数据。
9 N( ~/ a; \2 P% ^4 e- M% gdata4 <- data4[-136,] #删除该点
$ g& O" ^6 j# W3 K N* gx1 <- data4$x1
% R% n4 ^0 A5 ~5 tx2 <- data4$x2% ~' Y: R4 {( `, v3 N
x3 <- data4$friends0 n3 B+ N" |: R
y <- data4$salary
7 S; `9 D9 J& c+ Hlm.xy2 <- lm(y~x1+x2+x3) #重新拟合回归模型; b$ O. l( a4 a2 T
lm.xy2
9 Y" E# ^! f) N
3 E5 G& ]$ B6 T& q- Q#11
% Z7 |* f1 a* d1 ?/ f) h#对#10中的多元线性回归模型进行多重共线性检验,若存在多重共线性,删除相关变量后重新进行拟合。
& o! R) T8 F0 {2 X3 L* evif(lm.xy2) #p判断多重共线性0<VIF<10(不存在)/ g: k K1 ~' U( F
2 T M# x7 ]; o. q. y( r( T( F#120 o, i; T" N7 X9 T! D' b! x# ^
#对#11中的结果进行解释,重点分析年龄、在北京的居住时间和朋友数量如何影响收入。
/ E! L8 L* \( |& H+ ?% y8 Esummary(lm.xy2) #可知年龄和朋友数量对收入有影响,显著性*一颗星
& O5 V! F# A* } [1 O
; J% p) t U$ _* `8 Y, w, m**********************************************************************
( |' t T( T t! |3 |
1 {% M. k$ v7 i8 y% i. b- x! l二、利用多元线性回归模型预测收入. r& D' m/ R* @' y2 \6 v
View(data4) #124条数据4 L; b* ?4 D8 X( Q6 |; H
#1, S; t8 a! i3 }
#从数据集中随机抽取50条记录作为预测集,剩下的数据作为训练集。
0 }0 ^$ c5 g9 Itrain0 <- sample(nrow(data4),nrow(data4)-50) #训练集和测试集! h5 W4 g+ j' o0 e/ `
trainData <- data4[train0,] #训练数据5 G; Y& W7 }& l0 q8 M( W
testData <- data4[-train0,] #测试数据
, z2 b4 ]5 O! Z( f# ^" U
) ?, X+ }# A/ \. z0 Q" M9 R#2
$ P& a: P: v- R8 I5 i* E#针对训练集,利用多元线性回归模型拟合数据。
( a# Y t9 `" Ulm.xy3 <- lm(trainData[,"salary"]~trainData[,"x1"]+trainData[,"x2"]+trainData[,"friends"])
+ x3 D' W8 T4 o$ L5 Y; S# t4 G& T2 Q5 T' e# a4 Y5 i' I
#32 c$ G9 T/ m" }
#对(2)中的多元线性回归模型进行诊断,处理异常值。
- {" x( P' B$ J; Psummary(lm.xy3) S' @! R$ v$ u8 e m
par(mfrow=c(2,2)): q: N/ s9 m g& |
plot(lm.xy3). I. R U* X$ l
outlierTest(lm.xy3)
8 g0 L) G: o3 X, H1 h; KtrainData<-trainData[-c(150,32,82),] #删除异常值,随机的; y) D' s0 V* c, J1 y; U" o' `
& A2 L: ?9 W8 v8 H#4
4 D$ U( \" J$ W* E/ D( e3 u#对(3)中的多元线性模型进行多重共线性检验并加以处理。
1 Q0 j* i; p+ K( _vif(lm.xy3) #p判断多重共线性0<VIF<10(不存在)
X3 D# S* ]6 d: Rsalary<-trainData[,"salary"] #引入的数据是训练集的数据
8 A. q. U9 T- Tx2<-trainData[,"x2"]
& i% y) L2 O! S- R; ^# tx1<-trainData[,"x1"]) z" |9 e6 d, E! H& j6 Y; a7 {
friends<-trainData[,"friends"]
- ~3 g# {5 Z3 M# G- u, j% Tlm.xy3 <-lm(salary~x2+x1+friends)
* S' x5 @. g/ n) g' \9 E3 @- W- t0 g" c: e! g' M/ n
#5
' o7 |$ }3 `% M3 M- `#针对(4)中的模型,分别利用AIC和BIC选择最优模型。, \/ y0 y6 z) b9 Z
#AIC检验,赤池信息准则,选择最小的0 M, ~. ^5 y6 o* l
AIClm<-step(lm.xy3,direction="both")
+ U3 c' k3 F6 b. m% d#BIC检验,贝叶斯信息准则,选择最小的+ f' a0 Z: h+ @: f; Y
BIClm<-step(lm.xy3,k=log(nrow(trainData)),direction="both")4 _- [. C! A/ V' q" c6 @, A& m
7 }+ e3 \- `8 d: i. B8 A
#6 l3 t9 s- j( K/ w6 P* x
#利用预测集进行预测,比较全模型(包含所有自变量)、AIC选择的最优模型、BIC选择的最优模型
) c6 ?4 a5 v* k* e0 f) O4 G#这三个模型预测的准确性大小,并进行解释。' l' ~* A: S& p0 n. d
Allmodel<-predict(lm.xy3,testData)" c0 h5 U, t& I s$ p( z% v
AICmodel<-predict(AIClm,testData)
3 g& R; y4 V8 c8 C" zBICmodel<-predict(BIClm,testData): Q' G' a% a3 b. F: H
#均方误差检验,最小最优,分别计算全模型,AIC,BIC的均方误差( v1 L' J# o( }' g' n
#均方根误差亦称标准误差,均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根
& A* z; @1 r: ?1 o3 ?#标准误差能够很好地反映出测量的精密度7 d+ e, k# j7 }& q' ^- i
MSE <- function(x){" H3 @' z/ f2 }0 q9 E% n3 B/ c0 C' R
mse <- sum((testData[,"salary"]-x)^2)/50
' k- N, b ?: ? return(mse)$ _2 S: \& v7 l, V6 S1 T
}
, z$ Z; ?( _7 M1 AMSE(AICmodel) #AIC/BIC/ALL是误差最小的! U9 ]2 G% e% f7 o7 q1 _
MSE(BICmodel)
. x! x5 z4 S& v* T8 p" j8 fMSE(Allmodel)6 @. Q) v% }: H
/ j5 _) E' \& C$ A
9 U/ h+ a7 S+ G! t e0 q$ @$ d7 v/ @- U+ z Y3 J: Q/ n
|
zan
|