在线时间 514 小时 最后登录 2023-12-1 注册时间 2018-7-17 听众数 15 收听数 0 能力 0 分 体力 40339 点 威望 0 点 阅读权限 255 积分 12813 相册 0 日志 0 记录 0 帖子 1419 主题 1178 精华 0 分享 0 好友 15
TA的每日心情 开心 2023-7-31 10:17
签到天数: 198 天
[LV.7]常住居民III
自我介绍 数学中国浅夏
【R语言】回归分析案例:北京市商品房价格影响因素分析
/ t' A# j6 Z: K' U) S! f3 Z 这一案例是王汉生老师《应用商务统计分析》方差分析章节的案例,主要对离散型变量进行了处理。
2 r! g' c3 [! v1 d5 F 这里将连续型变量也加进来,进行协方差分析,建立完整的模型。
首先对房价进行对数变换,解决异方差问题:
S \% B; X/ T' G
行描述性统计分析,各连续型变量之间的相关关系如下: + h9 {; v- l# N
! x; X$ S/ w% k2 O& i) g 名义变量的EDA一般做箱型图。
模型按照全模型-变量处理(分箱等)-变量选择-回归诊断等步骤建立。
: y7 o8 a- m: U# Z9 R
' g2 r% s4 G( j
) b% B4 i" I' n1 [ L1 X 最终模型残差图: ?% c) b2 W* H$ N
3 D: G$ w; j! ~1 `' I6 U 通过模型分析结果可知,影响北京市商品房平均销售价格的主要因素有:
- Q* X4 u2 m# k5 \9 ~" f t 属性变量:所在辖区、所在环线、物业类别、装修状况、容积率大小(新引入);连续变量:绿化率、停车位住户比
. i) F5 b3 W/ A& m8 j' @ 属性变量的具体影响在此处分析略去。 6 X1 g# P( \- B0 x9 A
连续型变量的影响主要为: - c. w) g8 j! H7 j2 ]1 a
绿化率:绿化率的影响十分显著,由系数估计值为正,说明对房价有正向影响,绿化率越高的楼盘房价越高; 1 `" J3 ]% B0 G' I U( {% K
停车位住户比:有较显著的影响,停车位住户比越高,价格越高;
8 {# @# P/ ]2 Q9 Z 同时,原本为连续型变量的容积率经过离散化变为属性变量后:
8 v d- R4 K/ R# J$ G& ?" f% u 容积率大小:容积率分组有较显著的影响,高容积率的小区商品房价格更贵;
7 L. j$ H% y. r0 f 容积率与环线之间存在着交互效应。 + @: a, ~% o! {* e
rm(list=ls()) #清空当前工作空间
9 M( q8 f1 n1 t9 V3 L; g7 _( T setwd("D:/回归分析") + a5 D) _) O' d6 C3 e3 S+ y, v
a=read.csv("real.csv",header=T) #读入csv格式的数据,赋值为a + x2 G" k% Y- A, m" f2 n
View(a)
4 W, X% U. Q) Q1 r% s: F% t5 h7 H attach(a)
( q/ `# o0 \' o# { names(a) 2 E* w4 j, R. a
$ a; T. g- p' l1 {6 S
9 D* F' u! w& D0 v ##描述性统计
: l) a+ w- v1 M5 H
& ^- s- y+ _& ^, s1 s9 m
6 h0 q: M$ P, l5 S #未做处理的响应变量分布情况
$ i: p0 y* i# L# ?6 y5 Y- Y par(mfrow=c(1,1))
H8 q/ E/ q9 O; [ hist(price) * S( j+ R y% m. q' I
summary(price) #查看响应变量的描述统计量
& |! ^9 S5 k/ q( W4 L# S #连续型变量描述性统计
' h9 X: L2 J8 q* `8 d windows()
; P! W7 A ?7 ]7 h- N5 n pairs(a[,c(6:10)]) #所有连续型变量间的散点图 9 F+ E" P; L( P) R y& ^
par(mfrow=c(2,2)) ; l' j. c( c# F6 {
plot(rong,price) #每个连续型因变量与响应变量间的散点图 - T* s: g: `+ r4 f0 |
plot(lv,price)
- l; M* o3 Z# Y plot(area,price)
% }5 t% y7 W/ M4 a9 d plot(ratio,price)
; [3 y* K! \3 ^0 _) W | summary(a[,c(6:10)]) #查看连续型变量的描述统计量
4 S$ v2 t9 u0 i9 E cor(a[,c(6:10)]) #查看连续型变量的相关系数
& q9 \% I/ v! q4 Z9 X' V #属性变量描述性统计
" A+ l% T }5 Y F windows()
0 P; o: I% d$ V, }2 ] par(mfrow=c(2,3)) " e: M# r1 Y7 K2 F
boxplot(price~dis) #每个属性变量关于响应变量的箱型图
/ V: p7 ]+ R/ s% F2 l0 ~) J boxplot(price~wuye)
; {1 k' L/ ]" [# J, l boxplot(price~fitment) ! ]: V, k7 w0 x. k8 Q. x
boxplot(price~ring) 3 H; Y$ J6 H) H2 y
boxplot(price~contype)
0 f0 Z: C# H, C1 ~
# T( U$ u- C3 K. B0 X- u+ J5 R 8 `, S) M0 }5 q; {
+ O+ r W6 l( i0 r
' \: j% D9 s% b, Y$ Q ##模型建立
1 F9 h- l0 {( z7 G; G' v
3 x* J- B% F( ?" J+ \0 Y2 l9 g6 D4 @ - A0 W0 |9 ?1 E' z. k" ]
#在方差分析模型基础上加入连续型变量 / P2 |9 f7 M) |& \8 B6 n5 y* D
lm1=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+rong+lv+area+ratio)
! N3 P9 ^/ i/ e! A+ m anova(lm1) #方差分析 2 O; Z( }8 ?& f6 F3 k! w* c
summary(lm1) #模型参数估计等详细结果 * M- ^6 V& V( A) k( F
windows() 6 b" ~* s+ q7 [) t5 o. N, w; [1 P
par(mfrow=c(2,2)) 2 ? V( X( L( V% S4 _" ~' R
plot(lm1,which=c(1:4)) #回归诊断做残差图 ; ]$ T9 k& S: T+ |$ K4 U9 C: G, W
- ]# x9 a" S& `& q) W* T 3 R2 j3 M8 {, @. X" [, i& u8 {( E5 J
2 h3 F2 f, F9 j
4 x( y/ `3 Q# R2 l+ B ##变量处理 9 j: ?* r" Y/ a1 x8 u
2 f) b& K0 H6 f2 @
) m! c0 T) e5 Z* i9 {' [$ [ ###对不显著的变量采用分组的方式希望能达到显著的效果
( \( B7 @9 H. ^' X8 u ##对容积率的处理 0 c' z; M& S5 Q5 y7 S' y3 B- ]
windows()
6 g M6 j: _# R4 e2 h: |" `; G n = 4
) W0 y+ k9 \. ?" @ boxplot(price~ceiling(rong/n)) #容积率多分组下的箱型图
: {1 [# w* `" X/ v! o, ? table(ceiling(rong/n)) #容积率各分组下的样本数 * s5 A1 e( Q1 c: U8 S! p* N9 ]
ronggrp=1*(rong>n) #进行二分类
9 y1 q" b4 j" N% q; G, K0 ] #ronggrp=ceiling(rong/n)
) F: K# V8 }* y" l% f5 ~5 D& M2 B table(ceiling(ronggrp)) #容积率二分类下的样本数
2 |1 z+ k, q! W. f windows()
% W" T0 `- |1 t. R b/ e9 S0 e- b boxplot(price~ceiling(ronggrp)) #容积率二分类下的房价箱型图 $ y5 ^9 E! L' `/ Z2 ?6 B
windows() ; I0 u J: g/ L1 c1 a+ J
par(mfrow=c(1,2))
2 _8 v, l6 i4 |" H8 W* I$ C7 ~* M boxplot(rong~ring) #容积率与环线箱型图 . g* g7 T) h6 d0 G7 y' V
boxplot(price~ring) #房价与环线箱型图 7 Y1 K# O% R" u$ `
#加入容积率分组和容积率分组*所在环线交互因子的模型 + c" q; D# x- i Y! p! d5 h' T# s
lm2=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+as.factor(ring)*as.factor(ronggrp)+lv+area+ratio) # z7 O7 r7 ]$ {6 y8 d+ }! m
anova(lm2) #方差分析
, r( c5 @$ p! e summary(lm2) #模型参数估计等详细结果 ; ~8 B( ~& K0 |7 n" a$ C) M
windows()
8 J* K. S5 I. O& t! L) p8 } r% u, z par(mfrow=c(2,2))
0 \& s9 j; F' y3 \ F) { plot(lm1,which=c(1:4)) #回归诊断 9 @8 b# Z- i' R
' I# \8 A' y7 D( Z- Y5 K9 [2 d1 Q + H& N$ o9 x6 ]. }1 w
##对小区面积的处理
6 a6 n4 O2 l& u# c0 p summary(area) 1 w4 Y. N) O8 J
plot(area,price)
% @& ^+ w7 o I t& l( f! l7 j' n windows()
5 t+ I9 _% K1 I9 F1 I n = 150000 " r+ V( l- a4 q9 o+ s& q
boxplot(price~ceiling(area/n)) * V+ j. r( \% w4 T* m* G
table(ceiling(area/n))
0 d6 z' f4 }0 Q- t! t/ u2 p4 n- v( g areagrp=1*(area>n) ( @- G& G0 v1 q
table(ceiling(areagrp))
: N) S# n- I# D% n: a2 c boxplot(price~ceiling(areagrp))
& U% ?5 L- R! M7 d# ^+ q) b/ s #加入小区面积分组的模型 $ w" C3 d2 [4 h) P
lm3=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+as.factor(ring)*as.factor(ronggrp)+lv+as.factor(areagrp)+ratio)
2 m5 X2 ?: n1 g% w anova(lm3) #方差分析 # H x# u: @; N8 y' ?
summary(lm3) #模型参数估计等详细结果 % c1 t" a8 u9 j7 T6 t7 H& w
windows() " u# V1 t- u2 r
par(mfrow=c(2,2))
2 T7 u/ u8 \9 ^; M plot(lm3,which=c(1:4)) #回归诊断
: ] ~9 `# k' l2 @* b: X
* E5 s: G2 B$ v- ], r / d: N1 r) c' l9 [
##变量选择
1 I$ q3 @! K. K- ~! d2 D 6 ? t$ m; m9 {* p8 U4 H
' j* q4 Z" Y' J! D9 M- v
##AIC准则下的变量选择 ! `4 N4 E5 S9 @6 s5 F5 F6 w( a* f
lm4.aic=step(lm3,trace=F) #根据AIC准则选出最优模型,并赋值给lm.aic
- x+ d. x% e- c! Q summary(lm4.aic) #给出模型lm.aic中系数估计值、P值等细节
& \9 R1 K: P( M- O( `5 u$ x3 \% h ##BIC准则下的变量选择 $ e1 D' r# t# d( b1 @4 X
lm5.bic=step(lm3,k=log(length(a[,1])),trace=F) #根据BIC准则选出最优模型,并赋值给lm.bic 5 m) p5 r5 ]8 E) V1 w
summary(lm5.bic) #给出模型lm.bic中系数估计值、P值等细节
: p2 s0 ^. A6 i& K j, l6 U) _% X0 ? ( I( ?) q; Q! E+ H1 U
; z6 D, Z3 A7 n T" d( H( u #选用AIC准则下的模型进行回归诊断 : G1 s+ }1 E, w9 Z. f
windows() 4 u2 l' m7 e% V1 K j. v8 L
par(mfrow=c(2,2)) 0 j8 [1 S# b) W% S8 C5 H3 e
plot(lm4.aic,which=c(1:4)) & _8 X- A* f1 z; O$ [6 I% u
4 R" z+ r$ F2 y3 L
U& X6 W3 l7 ?
3 g3 O2 `$ a* ^( z @5 {: v3 |. x! E4 L; O( g
##数据变换
- O; y* d( e( S }4 V& ^9 O6 P7 { ! R0 y6 `8 u8 s! ]; H4 d9 V' b
& g# @8 R. O0 ~- ?& ^ #box-cox变换
' y* T% b& S8 z$ F/ w0 l8 j; s7 g3 h library(MASS) - y' Z" S5 M5 T. }0 y, `* a6 S; s/ w7 ]
b=boxcox(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(ring)*as.factor(ronggrp)+lv+ratio, data=a,lambda=seq(-3, 3, by=0.1)) * }7 V4 l5 i2 U" a! A
I=which(b$y==max(b$y)) #定位似然函数最大的位置 2 M8 m/ C/ G2 @$ v2 `: C" N
lambda = b$x[I] #精确的λ值 ' N5 e& [0 z3 N8 W' S
#λ接近于0,为模型简洁性,可以直接进行对数变换 % c! O5 k) E; C
logprice <- log(price)
8 f! ]9 ~9 [) Q- Q5 ^! p hist(logprice) - t0 V) O: P1 J4 V+ _6 W; @
# U0 a3 e# X. @1 q8 h
# f. E* Q2 F: {# K9 q& l* r
##最终模型与诊断
0 h4 j3 [- M, ^ . k1 |0 \7 O; G
$ J/ w* G- r! |4 U lm6=lm(logprice ~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(ring)*as.factor(ronggrp)+lv+ratio)
3 b2 Z) g4 d) X6 ?0 ~& f windows()
& l) M. K$ d; B! P5 W& @+ g1 P par(mfrow=c(2,2)) 8 W) u; ?; m- w O4 j6 n3 I
plot(lm6,which=c(1:4)) + h& n' t% b+ [6 e0 Q
anova(lm6) 4 ~% r6 x0 a5 ~7 p9 ?: K! ^
summary(lm6)
7 o4 A! x" r* u8 G8 f! N5 f d4 p " [& @" {# u# a, I$ T# F0 g/ s0 n
4 T2 m7 L: [* m. z) [! {1 I4 X 请关注数学中国网微博和数学中国公众号,联系QQ 3243710560& u8 V A2 S* R7 U
, M. L, ^( v$ J
7 ^( V; \( L0 t3 f) V
4 J W {# O1 k% J9 w
0 a3 L7 N- W. |; ?( a- |: q
zan