- 在线时间
- 514 小时
- 最后登录
- 2023-12-1
- 注册时间
- 2018-7-17
- 听众数
- 15
- 收听数
- 0
- 能力
- 0 分
- 体力
- 39388 点
- 威望
- 0 点
- 阅读权限
- 255
- 积分
- 12512
- 相册
- 0
- 日志
- 0
- 记录
- 0
- 帖子
- 1388
- 主题
- 1158
- 精华
- 0
- 分享
- 0
- 好友
- 15
TA的每日心情 | 开心 2023-7-31 10:17 |
---|
签到天数: 198 天 [LV.7]常住居民III
- 自我介绍
- 数学中国浅夏
|
【高级数理统计R语言学习】2 多元线性回归 一、背景0 Z9 n. \+ e" t" a
数据集展示了X市外来人口的相关数据情况,包括出生年月、收入、初次来到X市的日期、迁离X市的日期和现在的朋友数量。现假设外来人口的年龄、在X市的居住时间和朋友数量影响他们的收入。试加以证明。 二、要求和代码 一、分析收入的影响因素
* \4 D( d& Q! J2 v; g; M v#1
$ R+ P& ~. G$ H' u' r) i; i#展示数据集的结构1 G3 g( ~ P0 n% R! \) g
data2 <- read.csv(file="F:/hxpRlanguage/homework2.csv",header=TRUE,sep=",")
8 U- D) e+ i$ p: i+ u& Nstr(data2) #显示的结果有一列是多余的,需要删除' M4 S/ j: r% _' d1 T; n
data2 <- data2[,1:9]6 g5 ^ q. c" B2 r4 Y" C0 I* R( X
str(data2) #删完之后的显示效果是正常的没有多余列3 @- v$ @) J; b, [' i
3 P/ c5 |" @' y; f* a; k4 e& Y#2
3 S/ d5 M% w1 m5 r6 R#显示前10条数据记录
5 m8 v& F) Y/ H: L- [data2[1:10,]1 X) O, f: T' M6 C1 }
: d$ n+ m: c: W#3; J3 ~6 w( I" m1 Y' G: ?4 O
#将变量名重新命名为英文变量名* y# B: B0 {+ Q0 J- r( w
cnames <- c("number","birthyear","birthmonth","salary","inyear","inmonth","outyear","outmonth","friends")3 I4 J0 V5 R( K& S, n* g
colnames(data2) <- cnames
# P+ y7 r8 t8 N6 g9 sView(data2)6 D- L. y% F$ o V
9 R/ v( X% O- A$ e2 j6 F#4
6 z3 k2 i9 ^# w6 K, S9 @#查找数据集中居住时间小于等于0的异常记录,若存在,从数据集中删除这些异常记录" @" Q' R5 Z; u
x2 <- ((data2$outyear-data2$inyear)*12+(data2$outmonth-data2$inmonth))
5 q& G7 D7 ^, y0 Z#View(x2) #①先算出居住时间
2 H& v" ]+ q. O+ Kdata3 <- cbind(data2,x2)
2 J9 C5 z- ]+ T' m& C#View(data3) #②使用cbind函数把x2和原数据拼成新的矩阵,方便之后删除异常数据列,并且是127条
7 w: a: Q* B% h' D* Z+ Q! h. {list <- which(x2<=0)
; b+ P/ x4 H2 u% Zdata3 <- data3[-list,]
1 i: x0 C! M, n( i, d) VView(data3) #删除异常数据后是125条数据( _$ h, r6 A" s( M7 M. |& J
: e) R: { L- ^, Q: v1 g8 [#5
& g* l$ [4 c* r6 J#展示数据集中因变量与自变量的均值、最小值、中位数、最大值和标准差,要求保留2位小数。
* @, w8 N# Y5 T# V! k: L2 C( {library(lubridate)
4 U0 j0 c5 W6 @3 D2 P j* ddate<-Sys.Date() #返回系统当前的时间' `- O3 v# U* R
nowyear<-year(date) #提取年份
: Y& x |( `; U0 _. @nowmonth<-month(date) #提取月份" s# V7 M: Y/ [7 k
#View(date) #查看现在的日期6 f" a7 u' m$ t% M$ ]
#View(month(date)) #查看现在日期中的月份
) h6 W) L# T& {1 f$ w. Fx1 <- array(1:nrow(data3),dim=c(nrow(data3),1))
- P0 }9 P+ Y; `# ?1 _& G* V cfor(i in c(1:nrow(data3)) ){, q8 T8 }! l3 G
if(nowmonth-data3[i,"birthmonth"]<0){
( I0 y; g$ W8 j' g3 s% R( u x1[i,1] <- nowyear-data3[i,"birthyear"]-11 V c2 `+ P2 V; d
}else{6 H) `# o4 J+ v! l2 s% o/ w7 l- i) N
x1[i,1] <- nowyear-data3[i,"birthyear"]
; }& E3 {$ ~ [* q& l }3 w" u* m, A8 a. _& ^- ]$ Y7 U+ B0 Z
}4 ~+ R9 {6 V( c* s! I( q, S4 ?
#View(x1) #算出年龄x1,并加入到数据表中
3 f. Q6 J, f& w4 s0 jdata4 <- cbind(data3,x1)
! Q8 M2 J0 O( R: |3 z8 G, I; jView(data4) #加入x1年龄变量的新表展示
' F: p% ]3 ]$ M$ v4 ux2 <- data4$x2
6 k2 E% E0 V' k$ r$ \Mean.x2 <- round(mean(x2),2)
8 @, O! j/ a" yMin.x2 <- round(min(x2),2)
3 k9 P# A8 i JMax.x2 <- round(max(x2),2)
* ^7 A0 P& f7 o- mMedian.x2 <- round(median(x2),2)( `; F( p- M, a/ Z) O% k3 x3 L
Sd.x2 <- round(sd(x2),2)
G5 B* m2 c0 v% u* r1 Q; acbind(Mean.x2,Min.x2,Max.x2,Median.x2,Sd.x2) #x2居住时间的相关结果$ @& T2 }7 A! H
Mean.x1 <- round(mean(x1),2)2 B3 H8 n' ~' {( E: _' C9 ~% Y; }
Min.x1 <- round(min(x1),2)
7 {3 b& `% A3 t" \4 p, N, `! A* N1 qMax.x1 <- round(max(x1),2)+ a6 P/ B$ Y8 o6 T( q: Y
Median.x1 <- round(median(x1),2)
& a d' N1 r0 p2 X2 oSd.x1 <- round(sd(x1),2)8 V7 S; H" j$ }
cbind(Mean.x1,Min.x1,Max.x1,Median.x1,Sd.x1) #x1年龄的相关结果6 k6 j2 A/ l/ g$ [' F# L
x3 <- data4$friends7 q/ `) T5 S( x0 P0 E0 n' X
Mean.x3 <- round(mean(x3),2)- v5 \4 L$ f$ K% y. `
Min.x3 <- round(min(x3),2)$ I2 M, H* X8 T7 _
Max.x3 <- round(max(x3),2)4 N% d& A. p& ]: q
Median.x3 <- round(median(x3),2)
# [3 ? T) C, S# M* L9 Z$ q1 @Sd.x3 <- round(sd(x3),2)
# q! p1 s+ d' Y8 H) D) D+ lcbind(Mean.x3,Min.x3,Max.x3,Median.x3,Sd.x3) #x3朋友数量的相关结果, z( _' ~) e" X B1 v2 f
y <- data4$salary1 `- K/ O( O3 ~1 n, }, Q7 D
Mean.y <- round(mean(y),2)
8 \. r3 E9 c: z5 K O& i, LMin.y <- round(min(y),2)& s5 G* [. d( x# @. y0 P3 ]
Max.y <- round(max(y),2)4 S9 W9 y( ?+ [/ [5 {
Median.y <- round(median(y),2)
8 N% h- i2 m, gSd.y <- round(sd(y),2)$ z- Y: |: O# [! ]
cbind(Mean.y,Min.y,Max.y,Median.y,Sd.y) #因变量y的相关结果
# Y0 N* V& q, `" G0 Y' j7 r0 E
- a# a1 v% c- j#6
* ?1 x5 N& d2 i$ P/ P4 J#计算数据集中因变量和自变量的相关系数,要求保留2位小数。
0 l! J" _8 r$ n# I3 g; z- dround(cor(y,x1),2) #y和x1年龄' E9 b+ V6 e) }
round(cor(y,x2),2) #y和x2居住时间% t( Q5 t5 u) J V
round(cor(y,x3),2) #y和x3朋友数量
* n( ] u9 R! x+ F) |( ?; ]8 p) E J* W. S D
#7- D: i/ l8 Y! V! p* h ]
#分别绘制数据集中因变量与各个自变量的散点图. u: F8 I3 k/ \/ j( i- r
par(mfrow=c(1,3)) #布局,一行画3个图
& _: ] R- F+ fplot(x1,y,xlab="年龄x1",ylab="工资y")0 F: ^/ M& I1 ^9 I- }; ^
plot(x2,y,xlab="居住时间x2",ylab="工资y")1 U9 d: \! X( o+ T. m" J1 M
plot(x3,y,xlab="朋友数量x3",ylab="工资y")" q6 Q% _' p% o
, v+ {# M N% Y! s2 d#86 X5 H% y7 G! D' }9 J- i" _! h
#利用多元线性回归模型对数据集中因变量与自变量的关系进行拟合。
# Y' \+ }$ b4 |% w! F( j; {" Wlm.xy <- lm(y~x1+x2+x3)
1 J6 m; \/ y+ @- u, ~5 Vlm.xy
- I( h* [6 }3 {summary(lm.xy) #得到的结果是方程是显著的具有线性关系,但是每一个系数不都是显著的
7 v' ` W* {- d' I
X9 r) ^( e# p! q' |+ w#92 }$ o& g5 y/ P8 W$ L
#对#8中的多元线性回归模型进行诊断,确定异常值记录。
* W" f- F9 ^5 |! Y( vpar(mfrow = c(2,2)) #生成四种模型诊断的图形,2行2列+ d# b$ d; j- }
#生成四种模型诊断的图形:①残差与真实值的关系图 ②qq图用来检测其残差是否是正态分布
& ~' n% ?# G9 R$ E* j/ C#③用来检查等方差假设的。在一开始我们的五大假设第二条便是,我们假设预测的模型里方差是一个定值。
6 q8 o6 N2 G+ Z. W$ d+ P* o5 E' X7 \#如果方差不是一个定值那么这个模型的可靠性也是大打折扣的。, J9 `9 p; Z1 g
#④Leverage就是杠杆的意思。这种图的意义在于检查数据分析项目中是否有特别极端的点。5 I& [0 |- Z6 @3 R. M' W
plot(lm), e1 D+ M; R8 ]" g
library(carData)
6 X2 _4 x+ S" vlibrary(car)5 t U* |, ~# Z% g* [7 o
outlierTest(lm.xy) #显示离群点,Bonferroni校正,残差最大的点是136号点
1 E/ T& Q3 Y8 F9 ^. q0 p
7 X, M8 l6 M2 d. B4 U#10
6 O# q& c2 V4 w5 K; R. a#删除异常值记录后重新利用多元线性回归模型拟合数据。
/ B1 A S P6 E9 x3 L) o- {! l0 [data4 <- data4[-136,] #删除该点. \( i" t% O0 v) y" B+ j# {! m
x1 <- data4$x11 a: P; S B! k( m! t
x2 <- data4$x26 [- n4 j8 S5 C
x3 <- data4$friends$ Q4 m3 r& ]. y N
y <- data4$salary
. Z( Y& A- d4 F; Mlm.xy2 <- lm(y~x1+x2+x3) #重新拟合回归模型 y9 p3 ~1 m% w; s' G0 w
lm.xy2. {4 D, P" S3 P( U- N4 E
0 k+ `' {6 y) H+ m#11
" L9 g5 W+ R9 G* @: J#对#10中的多元线性回归模型进行多重共线性检验,若存在多重共线性,删除相关变量后重新进行拟合。) l6 O, c. S0 K% ?5 ~9 S2 j8 b8 i
vif(lm.xy2) #p判断多重共线性0<VIF<10(不存在)7 U6 z/ F1 \/ v6 @0 w( i9 ?6 g7 }
, B3 j4 F' `: t#12
. ]( F8 u0 s7 i0 ?3 }% Z#对#11中的结果进行解释,重点分析年龄、在北京的居住时间和朋友数量如何影响收入。 k- G" X6 q3 @ z2 ^, z! F
summary(lm.xy2) #可知年龄和朋友数量对收入有影响,显著性*一颗星
% G; h8 d& X. G+ R$ ^8 d) A9 A( _
5 e+ m: h/ I; ?5 {**********************************************************************, D, F0 i" C% G' q8 s
% Y ~7 ^' U2 R# Q D
二、利用多元线性回归模型预测收入 J# g+ z* Z: C. B5 d3 ^/ `
View(data4) #124条数据8 [4 z1 _1 `3 _: P- W6 S& v
#11 V9 h# ]& Z' d5 z5 ^/ Y
#从数据集中随机抽取50条记录作为预测集,剩下的数据作为训练集。
- O u" `7 l1 k& Y! j, b7 M! i9 W0 M$ atrain0 <- sample(nrow(data4),nrow(data4)-50) #训练集和测试集
; V- L# i9 H1 @( A! B* vtrainData <- data4[train0,] #训练数据
/ h p" l3 L5 X0 @. AtestData <- data4[-train0,] #测试数据
2 F* [1 S9 `6 Q' R3 R! E9 w( f) O1 Q5 G8 j6 ]
#2
( H3 t- b% y, c5 b; J6 [#针对训练集,利用多元线性回归模型拟合数据。( r S# T, A/ T2 J; X( ?
lm.xy3 <- lm(trainData[,"salary"]~trainData[,"x1"]+trainData[,"x2"]+trainData[,"friends"])" L8 l0 w7 {& O: D/ Z, c
1 t" H5 s) s# B$ N3 o+ H
#39 p! `! e/ Z5 F7 V& `
#对(2)中的多元线性回归模型进行诊断,处理异常值。
( y w' G- |7 N" E; f8 ]$ Esummary(lm.xy3)7 e0 z% l {; f& _/ J
par(mfrow=c(2,2))9 u' ^1 \8 `' j( @$ h1 o
plot(lm.xy3)" e( u% _3 B' c- k. B, ~ p/ X7 c
outlierTest(lm.xy3)( Z: b3 A1 O( ` N. c, j& e" L
trainData<-trainData[-c(150,32,82),] #删除异常值,随机的& y+ ]9 P0 k3 g/ Q
$ Q8 W# i0 A5 P3 ~1 F) w5 C% i
#4
0 f H2 K% Q' [+ }- ` [+ ]1 l#对(3)中的多元线性模型进行多重共线性检验并加以处理。
) x$ C% G' Q6 B0 r! m, H3 d9 ovif(lm.xy3) #p判断多重共线性0<VIF<10(不存在)5 y4 O: v7 ]( L& j
salary<-trainData[,"salary"] #引入的数据是训练集的数据
5 h3 v6 S. M/ k- E3 f2 G; ^+ e- b1 px2<-trainData[,"x2"]& A M" g3 m6 E, l7 K
x1<-trainData[,"x1"]
$ Q- m' B# G+ E: L+ C5 d9 p8 xfriends<-trainData[,"friends"]
( _2 N1 ^, V- @* ilm.xy3 <-lm(salary~x2+x1+friends), E5 q! @2 [5 Q. S" u) C: B' [
/ E+ q% }9 z# ~; ]4 Y. s
#5+ H. A% F7 B3 D2 l4 x
#针对(4)中的模型,分别利用AIC和BIC选择最优模型。
1 n; ]: @% t& [#AIC检验,赤池信息准则,选择最小的
+ G7 a/ z' a2 [+ Z l8 J" yAIClm<-step(lm.xy3,direction="both")
: O9 a4 S7 T d: D#BIC检验,贝叶斯信息准则,选择最小的
/ v7 b) ]/ @5 B5 {9 [% x3 e# v+ CBIClm<-step(lm.xy3,k=log(nrow(trainData)),direction="both")
1 q; ?# M; T, Q8 Q( A# n+ f% {6 O- ~7 Y$ K/ I
#6
+ @, T& h1 B% }% @" M0 X#利用预测集进行预测,比较全模型(包含所有自变量)、AIC选择的最优模型、BIC选择的最优模型3 V) }( m+ {+ f. h" I/ ~+ H+ f
#这三个模型预测的准确性大小,并进行解释。
% c0 {# Y' ?! l! r bAllmodel<-predict(lm.xy3,testData)
: X! N: Z ~. b9 E0 i5 nAICmodel<-predict(AIClm,testData)
% ]! ?( z" T7 G* {+ R, \BICmodel<-predict(BIClm,testData)5 r+ ?/ ?" V( c( X
#均方误差检验,最小最优,分别计算全模型,AIC,BIC的均方误差
' {' s! F# Q7 N" M: C/ |( L2 i#均方根误差亦称标准误差,均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根
7 Y! l; H+ Z: y#标准误差能够很好地反映出测量的精密度- k7 S4 A4 D2 C
MSE <- function(x){
" @. S, o- I1 L0 e5 T' W% I0 c mse <- sum((testData[,"salary"]-x)^2)/50
( R2 _4 b$ @2 c/ K/ G return(mse)0 n, n7 k8 }2 f, |& r3 V
}
! i8 Z5 E1 M- i1 w* c+ MMSE(AICmodel) #AIC/BIC/ALL是误差最小的' T2 `2 @) ]1 m3 e: j
MSE(BICmodel)
; Q8 H! T& s" YMSE(Allmodel)
N* O& s7 P- C) t! [
; [, w7 o. ` p& \" w% r& `' Y
/ h4 i# ]* [* J5 u$ m
/ D! _' E6 n Y& w |
zan
|