- 在线时间
- 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 多元线性回归 一、背景
4 @% I; r+ l5 j" q8 a* J& w数据集展示了X市外来人口的相关数据情况,包括出生年月、收入、初次来到X市的日期、迁离X市的日期和现在的朋友数量。现假设外来人口的年龄、在X市的居住时间和朋友数量影响他们的收入。试加以证明。 二、要求和代码 一、分析收入的影响因素 I8 O% B" e' H- l5 ]
#15 l! a. a/ C2 Q2 Q4 ?: O" o. B5 t
#展示数据集的结构) f$ D; f2 `, S2 K0 s9 d
data2 <- read.csv(file="F:/hxpRlanguage/homework2.csv",header=TRUE,sep=","), _; }9 D! o1 x: r
str(data2) #显示的结果有一列是多余的,需要删除
& ^4 c& G/ s0 y' Ldata2 <- data2[,1:9]0 Z: U8 o3 \- G8 V
str(data2) #删完之后的显示效果是正常的没有多余列1 m" `* r* I# [& t# a4 x- T8 D# k
5 K* x) u7 J" b4 M) }2 x- ]#2
" M# t, M) I! ~: C1 M. t: N; q#显示前10条数据记录! t* E& w7 U' A/ m
data2[1:10,]
6 [: W! P1 s6 L" q0 B; [
7 @0 h& m# Q# k# G. X t( Z0 w#3
0 T0 D; A. U( W. s4 t#将变量名重新命名为英文变量名
5 P* w8 Z$ g7 _( Z( ~cnames <- c("number","birthyear","birthmonth","salary","inyear","inmonth","outyear","outmonth","friends")8 V* G. Q# u# D Z
colnames(data2) <- cnames% ?6 T/ ?6 } k# }
View(data2)0 v% ]( R8 S8 z8 ?$ _
% d7 B& i- x9 h+ B* z! H
#4
1 j- @0 [. }" D" ?; \) |#查找数据集中居住时间小于等于0的异常记录,若存在,从数据集中删除这些异常记录3 K9 U4 {8 d& S1 j% c0 f% P" ?- u
x2 <- ((data2$outyear-data2$inyear)*12+(data2$outmonth-data2$inmonth))- V% ]- x; {* U2 Z3 I
#View(x2) #①先算出居住时间! j7 f- q8 [) s$ J
data3 <- cbind(data2,x2)+ O+ A+ @+ {. X& y
#View(data3) #②使用cbind函数把x2和原数据拼成新的矩阵,方便之后删除异常数据列,并且是127条. B `( }- Q5 _" a" e3 p
list <- which(x2<=0)
7 f- @0 V" z2 M3 Cdata3 <- data3[-list,]
; c6 Z0 B/ R( u6 b, c/ P1 kView(data3) #删除异常数据后是125条数据
. b2 @6 T+ j( z( A* U6 [. v" }2 n" z. E" k+ b
#5
! o! w! x3 W/ F/ m% u+ h! z; x( E#展示数据集中因变量与自变量的均值、最小值、中位数、最大值和标准差,要求保留2位小数。) r' F" H V- ?3 X7 T D
library(lubridate)( `0 u0 q6 }( c" \( I
date<-Sys.Date() #返回系统当前的时间
, Q" L* p, X6 \! N1 Z- |$ Ynowyear<-year(date) #提取年份
- q; L" Z, A7 a+ Z! Onowmonth<-month(date) #提取月份
6 C! [& f/ B$ U E- r#View(date) #查看现在的日期+ S" g6 V4 K( |+ Y \
#View(month(date)) #查看现在日期中的月份/ m: H: I& _, Q" M* V5 g
x1 <- array(1:nrow(data3),dim=c(nrow(data3),1))5 W5 x. d; T l$ @. V- w. g
for(i in c(1:nrow(data3)) ){
1 W3 R5 H3 ^9 P if(nowmonth-data3[i,"birthmonth"]<0){3 W" g, P8 _" |6 Y
x1[i,1] <- nowyear-data3[i,"birthyear"]-1
7 l* t; j$ a. }3 k( v. |7 l9 d+ q }else{
6 c8 y6 \/ o9 p7 ?9 h; T1 S x1[i,1] <- nowyear-data3[i,"birthyear"]
$ ?7 d. W# o3 U2 z2 K5 I6 i }
! E6 i$ h. ^( C7 e% ]: n}# b" P5 L' Q& p* j: h
#View(x1) #算出年龄x1,并加入到数据表中
, M4 g/ M; M4 X2 odata4 <- cbind(data3,x1) . b6 z/ G7 }4 ]* }. x0 b
View(data4) #加入x1年龄变量的新表展示1 m0 I9 z/ X* q% X9 ^2 c
x2 <- data4$x2/ `2 m1 y ~: x3 ?' d- f* s5 e& G
Mean.x2 <- round(mean(x2),2)
* z7 |% T6 h+ N5 AMin.x2 <- round(min(x2),2)$ i! u$ J6 |2 A
Max.x2 <- round(max(x2),2)% e5 c! n- \+ x: U: V
Median.x2 <- round(median(x2),2)
* s. w. D' A0 H7 a' l) t( v& ~/ nSd.x2 <- round(sd(x2),2)
9 M8 T2 k: j! ~cbind(Mean.x2,Min.x2,Max.x2,Median.x2,Sd.x2) #x2居住时间的相关结果! T8 @$ e! r' |& A6 _
Mean.x1 <- round(mean(x1),2)4 }: m4 K9 D4 Y' G
Min.x1 <- round(min(x1),2)
8 |9 k6 \) b) y4 EMax.x1 <- round(max(x1),2)( ?4 K! W0 C0 g" Y* P
Median.x1 <- round(median(x1),2); G/ k; D1 f( N* `
Sd.x1 <- round(sd(x1),2)
, W _# w2 ^8 a3 l" a' v, f) Hcbind(Mean.x1,Min.x1,Max.x1,Median.x1,Sd.x1) #x1年龄的相关结果
7 | H' Q/ @. N! Ux3 <- data4$friends
$ ^" z. [, B8 m7 |' R1 N6 I+ KMean.x3 <- round(mean(x3),2)
! K$ i/ Y+ P. k6 s. J1 v9 mMin.x3 <- round(min(x3),2)
, E5 u. ~* f# n$ XMax.x3 <- round(max(x3),2)0 l$ E5 E' s. g: @7 i! |6 R
Median.x3 <- round(median(x3),2)1 q# i3 B4 `9 b l5 H! l
Sd.x3 <- round(sd(x3),2)' |: t" B' ^. q, ^% c- w
cbind(Mean.x3,Min.x3,Max.x3,Median.x3,Sd.x3) #x3朋友数量的相关结果
) |6 I* c, O: U& g& m& x# ey <- data4$salary
3 T* h/ F3 ]9 s+ {: k, pMean.y <- round(mean(y),2)0 S- m7 X! i# {
Min.y <- round(min(y),2)
0 u s% I5 O5 {; P0 GMax.y <- round(max(y),2)" F& M- p/ ^4 D7 q R1 D
Median.y <- round(median(y),2). I" H$ |2 p) N* o0 {5 S
Sd.y <- round(sd(y),2)- S# l# ?! r9 A- P- o7 T
cbind(Mean.y,Min.y,Max.y,Median.y,Sd.y) #因变量y的相关结果( o# D- R, h% x9 f3 \! R0 u$ U
0 J2 a) H2 q7 \. h& I; ~. i9 R#63 S8 } d9 S4 [6 F
#计算数据集中因变量和自变量的相关系数,要求保留2位小数。- z4 r) _. w4 a7 J
round(cor(y,x1),2) #y和x1年龄 A* B# h6 J4 l2 _+ |& L6 m
round(cor(y,x2),2) #y和x2居住时间7 I9 P: f" b% [+ O5 \7 y5 v" x a
round(cor(y,x3),2) #y和x3朋友数量
$ E' D& t" N' K9 o5 p# x# T1 W( v& H$ D" [3 R* o2 f) b! Z, t
#7+ F1 X- J/ ^8 l; b
#分别绘制数据集中因变量与各个自变量的散点图
9 A+ K0 c) m3 r' spar(mfrow=c(1,3)) #布局,一行画3个图
% C- V! c h3 r4 b: i" Fplot(x1,y,xlab="年龄x1",ylab="工资y")5 A# s+ s' x. k) j2 y
plot(x2,y,xlab="居住时间x2",ylab="工资y")
. B) y3 Y) l* c& Z6 z& B$ ?plot(x3,y,xlab="朋友数量x3",ylab="工资y")
- s# T; _- w- g8 Q3 A6 n) S* P% T+ c& o( f. E! t
#8
3 t( v1 T7 B3 ], F9 k0 V#利用多元线性回归模型对数据集中因变量与自变量的关系进行拟合。( s0 y( }* F4 L7 p y
lm.xy <- lm(y~x1+x2+x3)) i$ J9 c7 T2 {6 H; Y; t! N7 s, }
lm.xy
9 b: U- B( s( n9 @, b9 vsummary(lm.xy) #得到的结果是方程是显著的具有线性关系,但是每一个系数不都是显著的: y: K+ r" ^+ s% q9 h
" U5 K T) b: K1 e
#9# J5 x7 f& ?4 k+ J
#对#8中的多元线性回归模型进行诊断,确定异常值记录。
9 T/ H9 R& W& M0 H, jpar(mfrow = c(2,2)) #生成四种模型诊断的图形,2行2列
; I$ c) O6 W9 V f#生成四种模型诊断的图形:①残差与真实值的关系图 ②qq图用来检测其残差是否是正态分布
, M9 x M/ Y; s6 f' O#③用来检查等方差假设的。在一开始我们的五大假设第二条便是,我们假设预测的模型里方差是一个定值。
' C) J, E, H/ i* o* u6 _' j! Y/ S#如果方差不是一个定值那么这个模型的可靠性也是大打折扣的。
' q, X% ?' L7 C1 D# e, |) n#④Leverage就是杠杆的意思。这种图的意义在于检查数据分析项目中是否有特别极端的点。5 _! @( E: }) _
plot(lm). ^6 K+ J# i' m0 U5 C% p* X
library(carData)8 _' Y3 x3 b4 N+ b9 w
library(car)
$ F7 e. X- s+ d: o" d6 @outlierTest(lm.xy) #显示离群点,Bonferroni校正,残差最大的点是136号点
B; N" H. U' E. f& K) J) E) ~# U' d( Z+ d& U1 {
#10
! k8 c9 i6 }3 O: ~% |6 C1 P+ h#删除异常值记录后重新利用多元线性回归模型拟合数据。9 g" V- h7 B v# b' H! P4 z! Y
data4 <- data4[-136,] #删除该点) M. y2 |) R/ a& a: F" @: M
x1 <- data4$x1
1 {* G9 M3 ^! D$ Q$ G& C9 g9 b/ Ex2 <- data4$x2$ A% t7 h1 R. x
x3 <- data4$friends
6 Y1 q4 T4 u8 @y <- data4$salary
4 _( Z3 i# A- p5 [1 Xlm.xy2 <- lm(y~x1+x2+x3) #重新拟合回归模型! ~1 Z. ~0 ?6 y" p2 `: B
lm.xy29 Z/ Z! L) K. D) E( C- b
: S- m% y' N: [( V6 P! T0 U- h7 Y#11
. `& _7 \1 i1 y! @2 T#对#10中的多元线性回归模型进行多重共线性检验,若存在多重共线性,删除相关变量后重新进行拟合。
0 Z( z u6 Y) ?/ i L9 [. nvif(lm.xy2) #p判断多重共线性0<VIF<10(不存在)% ~. x' B8 ?- ~! K1 T
7 h9 j) L: Q# w9 [( D# A/ \$ r#12
9 v( y) [9 u$ a: B3 P1 h* O5 k9 U#对#11中的结果进行解释,重点分析年龄、在北京的居住时间和朋友数量如何影响收入。
|9 ~3 W) A3 D# ~4 Zsummary(lm.xy2) #可知年龄和朋友数量对收入有影响,显著性*一颗星
# Q0 t: x' {+ U1 E3 M) X3 M \! c; \" ]; c" H
**********************************************************************3 X5 W+ j# I& `/ z. V( x
. B/ l' c) m# h" M4 t$ X! d二、利用多元线性回归模型预测收入 T! x# g/ q$ c) y. D4 s/ X+ F0 u
View(data4) #124条数据
0 W7 W0 k" o. _1 h4 Z+ T, T4 }#1
( L! a4 G3 I, e7 \ U#从数据集中随机抽取50条记录作为预测集,剩下的数据作为训练集。" A1 T2 `. `7 r, p3 e
train0 <- sample(nrow(data4),nrow(data4)-50) #训练集和测试集
& C& I# m4 l# T7 otrainData <- data4[train0,] #训练数据
! A* h8 r8 c" t) z0 h9 @testData <- data4[-train0,] #测试数据" s0 N$ q9 _4 _/ K- K% Q
& f A3 r' ]1 z3 d3 S1 ?#2
8 @& }+ w/ `2 X4 ~* S#针对训练集,利用多元线性回归模型拟合数据。* i3 u3 f5 ~& {& \5 U: n
lm.xy3 <- lm(trainData[,"salary"]~trainData[,"x1"]+trainData[,"x2"]+trainData[,"friends"])3 z6 _+ A3 }, ]$ u/ T/ J
0 g- Z9 e" z. o2 \6 A
#3" Z" K# Q; k9 \* k7 U `
#对(2)中的多元线性回归模型进行诊断,处理异常值。
* \* |. R# d P9 j7 M! h7 Q4 Osummary(lm.xy3)
' {0 i% `- i# S) o6 e- Npar(mfrow=c(2,2))
% K/ {/ f6 s: \; ?' i5 I* i3 ~plot(lm.xy3)
1 F0 a9 N7 @ E# N2 Z* JoutlierTest(lm.xy3)
# r8 d1 r* l! J" M% a2 v$ f( _trainData<-trainData[-c(150,32,82),] #删除异常值,随机的
# q/ a" a" e3 b$ C
: v! P+ P7 b( O) @; H" o ]#4: @& U5 d+ x/ x& |
#对(3)中的多元线性模型进行多重共线性检验并加以处理。
+ y* d4 t2 d+ z1 x, U% R4 nvif(lm.xy3) #p判断多重共线性0<VIF<10(不存在)
$ A4 j/ D$ v) ~/ z. Bsalary<-trainData[,"salary"] #引入的数据是训练集的数据$ Q6 ~$ n' f$ C, ]
x2<-trainData[,"x2"]3 u p& O/ x: w' W, i
x1<-trainData[,"x1"]
) h' @2 j% Q: [friends<-trainData[,"friends"]
8 Y* f3 |$ [6 c4 b" c" Wlm.xy3 <-lm(salary~x2+x1+friends)$ b$ w. N& V. o
% \: }: }' \$ i! A g. `#5
5 r" y/ [* r# L7 ?' c3 k) H6 T#针对(4)中的模型,分别利用AIC和BIC选择最优模型。( Q- i7 ]2 S9 g% w
#AIC检验,赤池信息准则,选择最小的
8 j ]8 } f" S) T$ E- dAIClm<-step(lm.xy3,direction="both")
7 X1 t1 t+ U: J. `) B7 _#BIC检验,贝叶斯信息准则,选择最小的5 U: q2 ]- ~3 W8 T- M
BIClm<-step(lm.xy3,k=log(nrow(trainData)),direction="both")
0 J- X- A9 O" g, A+ j1 W7 V" S, k, R8 G- _+ V- S3 B
#6
/ z! S; e% J# @ |/ j#利用预测集进行预测,比较全模型(包含所有自变量)、AIC选择的最优模型、BIC选择的最优模型
3 [+ D% u) ?" L8 U* `( I#这三个模型预测的准确性大小,并进行解释。& c) B8 f4 ]0 T) n$ I# ~
Allmodel<-predict(lm.xy3,testData): g; s; i; Z' E3 d. x; k
AICmodel<-predict(AIClm,testData)1 z) U, _3 M1 q3 j$ i
BICmodel<-predict(BIClm,testData)
7 j% P6 l+ d5 D1 O7 d1 c. x#均方误差检验,最小最优,分别计算全模型,AIC,BIC的均方误差/ a9 w0 w% f# u; U0 D9 q# ?
#均方根误差亦称标准误差,均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根# A9 L& \/ W6 x2 o! s
#标准误差能够很好地反映出测量的精密度, b# g0 ~ _: ~* B# m2 q
MSE <- function(x){
6 |# ^" U8 C- i( t. y7 E mse <- sum((testData[,"salary"]-x)^2)/500 P! h, r& q2 }( P& M# a- Q; V
return(mse)
' ^' |$ I& }! T; \: T# j$ [ I}
3 H# `, k; o4 K P3 c$ sMSE(AICmodel) #AIC/BIC/ALL是误差最小的
. Z! ~2 m; Y. a3 g d. }MSE(BICmodel)) V3 q) `& G% B0 j y
MSE(Allmodel)8 [2 K0 J) J9 M9 I& F' E
% [$ S, d! n X+ s2 [
8 d+ }% C. V) [, e
0 D1 h' \# i$ @5 ^8 u |
zan
|