【R语言】回归分析案例:北京市商品房价格影响因素分析
' T, @& n! M: G8 k* h这一案例是王汉生老师《应用商务统计分析》方差分析章节的案例,主要对离散型变量进行了处理。" L, P) y# h3 |, m: f3 N! \
这里将连续型变量也加进来,进行协方差分析,建立完整的模型。 首先对房价进行对数变换,解决异方差问题: & c0 w8 z# O+ J5 Z8 Y: |! L
行描述性统计分析,各连续型变量之间的相关关系如下:
* A s: i# o! {& R7 H3 L![]()
: D0 P7 Y' X. b名义变量的EDA一般做箱型图。 模型按照全模型-变量处理(分箱等)-变量选择-回归诊断等步骤建立。
( X& {( |$ g2 |" t![]()
5 F5 z5 b; ?# C! b7 h( \ 9 B" ]0 a) M' N1 ]8 Q: ?
最终模型残差图:8 H* E% h% |0 X* x& o: E
+ y( \9 ^8 A- c- g, ?* B
通过模型分析结果可知,影响北京市商品房平均销售价格的主要因素有:& {9 ~' a3 E& o4 y+ u
属性变量:所在辖区、所在环线、物业类别、装修状况、容积率大小(新引入);连续变量:绿化率、停车位住户比* u k, ~: N: l' a
属性变量的具体影响在此处分析略去。
5 j; Y& F* P }: n, x8 _* x! j连续型变量的影响主要为:
) c! w; r& T( y1 L 绿化率:绿化率的影响十分显著,由系数估计值为正,说明对房价有正向影响,绿化率越高的楼盘房价越高;
9 }, _) A, i8 U; g8 w+ o! P1 {, @ 停车位住户比:有较显著的影响,停车位住户比越高,价格越高;% n* ]% L5 g- e7 t: M, U+ X
同时,原本为连续型变量的容积率经过离散化变为属性变量后:; L0 c. x {9 Z$ {/ R
容积率大小:容积率分组有较显著的影响,高容积率的小区商品房价格更贵;
# x" s! E7 ]4 W- x' ]1 L5 Z 容积率与环线之间存在着交互效应。8 W* k( O& X1 c( q# k
rm(list=ls()) #清空当前工作空间( z6 a U* c& A8 L
setwd("D:/回归分析")
; W, Q( y4 v/ G1 c; g! A xa=read.csv("real.csv",header=T) #读入csv格式的数据,赋值为a
. B* k. L+ S3 LView(a)& m2 H. @ N# r& |8 n* ?
attach(a)
$ |- p% `- P9 Snames(a)
+ T3 f8 T/ c3 }; Y) A0 h- L* h1 ]
9 V! S( ?2 p8 Y! o# s3 M
+ U& _* G! p3 |+ u8 S% P0 n##描述性统计
. q1 P) j. H* v2 o7 P
8 v& a: b! p @3 K9 w& c
8 l% G" `4 f. `2 Q- F& d' W% R#未做处理的响应变量分布情况' `. q1 B3 i: s
par(mfrow=c(1,1))
* t6 k. ]( d! c: G* P1 H( Vhist(price)
+ ` J* S: w; p! i+ k9 K6 O, ^. `summary(price) #查看响应变量的描述统计量
8 n! w! k t+ E7 u/ s- I! a7 [! k: u1 ?#连续型变量描述性统计2 P8 Y5 p% L3 p
windows()
, z# Q8 ^7 a- V( }" C9 ~pairs(a[,c(6:10)]) #所有连续型变量间的散点图' ?% R, j0 P7 A
par(mfrow=c(2,2))
: D+ w! f2 B/ \plot(rong,price) #每个连续型因变量与响应变量间的散点图) G. w6 l1 h8 N* }
plot(lv,price)
& v, v) C. R5 W; |4 M+ rplot(area,price)( ~9 K* e2 s7 x9 E" I1 v2 T
plot(ratio,price)
7 M- y' O' W5 t* c6 fsummary(a[,c(6:10)]) #查看连续型变量的描述统计量$ [" n0 k9 W: ?& C, L: d
cor(a[,c(6:10)]) #查看连续型变量的相关系数) {" R; E$ k9 \# R! l1 R, N
#属性变量描述性统计
7 P2 V9 i; Q7 [. a* X6 C3 Z. lwindows()
3 k8 x v5 U% H$ t4 Gpar(mfrow=c(2,3))
) K9 q" C B9 ]1 Oboxplot(price~dis) #每个属性变量关于响应变量的箱型图& Z( z/ i4 M; h1 R* U0 s
boxplot(price~wuye)
! E' }& \! M! p' @5 g9 M; \2 w" dboxplot(price~fitment)
; N* q8 L9 a4 ?" z8 Jboxplot(price~ring)
$ W3 A q9 _- ~9 t% jboxplot(price~contype)
8 Q" p. O y. f X. |% T) F9 N
0 V2 v i3 h2 Z7 ^4 g2 w" c+ W: _ Q7 @6 r5 Y/ @1 J, @$ g' a
2 Q: n! E' A/ j7 a
5 r5 t7 X! y! J3 N* W8 c
##模型建立+ F. F, m1 t) H% |" O
( U" O L- V; i+ y4 ~ b8 F. U9 t
$ z4 f! |8 G/ ] U#在方差分析模型基础上加入连续型变量4 `9 c) I% V/ G4 U, |. H3 T6 F3 P
lm1=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+rong+lv+area+ratio)6 ~9 N, E. i3 S; }8 W
anova(lm1) #方差分析1 p3 l& c/ |- u+ Z0 p5 ~, w6 U
summary(lm1) #模型参数估计等详细结果% F- U9 Q( q$ y& V/ F
windows()8 B" E5 [& P7 r" F7 I/ k
par(mfrow=c(2,2))& V {% h; M D
plot(lm1,which=c(1:4)) #回归诊断做残差图7 O4 ^$ K s% J8 l. V V0 O
* y9 a* F$ I' H, \9 P+ D2 Y4 J: u6 A
/ \) l, h9 i$ P, D% [
. [& ]0 W+ G+ l##变量处理
9 a- c, ~) t) I4 l$ s' p7 r; r; \% ^/ {
5 \; f% J/ }( e h! W5 E4 \7 A
###对不显著的变量采用分组的方式希望能达到显著的效果
2 V6 f" r$ h$ Z) U) }##对容积率的处理
& f5 W" z9 y* c5 M$ \6 V! Z9 g% Hwindows()
1 K S# L% d K0 A* A; `n = 4& o0 d- C6 H9 @$ Q
boxplot(price~ceiling(rong/n)) #容积率多分组下的箱型图
+ V# [. F1 |' c( }2 o& dtable(ceiling(rong/n)) #容积率各分组下的样本数
4 r6 e# D, m7 z# tronggrp=1*(rong>n) #进行二分类
8 M% ?& T! ?: q5 a1 s#ronggrp=ceiling(rong/n) * I- ^) Y* i% T2 q5 a0 t
table(ceiling(ronggrp)) #容积率二分类下的样本数
0 [2 K, m& ~, N0 Xwindows()+ S! V8 n* n" D/ f3 n) R7 U
boxplot(price~ceiling(ronggrp)) #容积率二分类下的房价箱型图
4 Q- q6 E- H7 D9 mwindows()
" d0 `) E( o! _0 ]# n- Zpar(mfrow=c(1,2))
& t. c- a3 ]9 @# Fboxplot(rong~ring) #容积率与环线箱型图 [+ u+ i3 A, o0 y3 e
boxplot(price~ring) #房价与环线箱型图
9 v# D1 ?$ w; g! e#加入容积率分组和容积率分组*所在环线交互因子的模型7 O/ J% N: q7 }' i0 G# g, 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)
$ Y; w( u$ S8 V1 Y' Banova(lm2) #方差分析
$ I3 j+ c: Y% O8 @7 Q1 dsummary(lm2) #模型参数估计等详细结果8 B. [% A, U1 y9 o. L# b% n5 @/ t
windows()( x0 Q# B0 F# Q! W& h, C4 R8 I
par(mfrow=c(2,2))+ N5 ^# Z' X. |
plot(lm1,which=c(1:4)) #回归诊断# Y- T( N' _- F+ d4 A
' f4 ~2 P9 c: X. w. i! W
4 T7 j7 k; V+ e( F8 N
##对小区面积的处理2 {* z; S& P: j- _) ~
summary(area)
6 q$ K& N5 S! ~( H- a, J( r1 n3 W* a# Gplot(area,price)
/ }; ]3 P% x/ nwindows()& y. f& e, e# ?1 n" i
n = 150000
3 S) }$ e$ B4 i: T0 N3 K! b1 qboxplot(price~ceiling(area/n))
0 H! i# ~5 z8 P2 a- V X2 btable(ceiling(area/n)) 0 e8 Q# ?' ]+ E: K5 L- o0 z
areagrp=1*(area>n)5 Y7 |5 }" ^1 P5 n& E5 p- X5 t1 X
table(ceiling(areagrp))
- D: N* ?# X$ h4 S4 ^boxplot(price~ceiling(areagrp)). J) ] i) B7 \/ x% {9 i! v
#加入小区面积分组的模型* F2 K) u& }( H
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)! A6 B' r' N- N( U( n/ X4 a7 W7 G4 r; }
anova(lm3) #方差分析
; ~$ ^9 Z$ F T* Vsummary(lm3) #模型参数估计等详细结果. P; }8 c; @4 ]6 A- s
windows()$ G6 D5 N" w X& [5 @' M. a
par(mfrow=c(2,2)), p. D. D- z2 H
plot(lm3,which=c(1:4)) #回归诊断0 [* O0 U0 o- \# L
) t; j3 M- v) L. J4 v/ r
% {, g/ }0 m& _9 X
##变量选择* f, U( R& d3 ~7 `# z( e
- \( m3 M, z5 A2 t/ ^0 ^
% I; p5 l$ A& r7 H##AIC准则下的变量选择3 F& B1 |; H9 b# C4 N1 s6 A
lm4.aic=step(lm3,trace=F) #根据AIC准则选出最优模型,并赋值给lm.aic4 w! j! J, ?, Y V) R, L" q) |
summary(lm4.aic) #给出模型lm.aic中系数估计值、P值等细节
+ l: F+ J6 \8 g5 g##BIC准则下的变量选择; y+ W# h& ^# [5 u) `4 T& \: \* ]
lm5.bic=step(lm3,k=log(length(a[,1])),trace=F) #根据BIC准则选出最优模型,并赋值给lm.bic' t( i/ O8 J) B6 {! i
summary(lm5.bic) #给出模型lm.bic中系数估计值、P值等细节
, S* V! ?# B1 U3 J0 w% }
& H! i+ _, j% S% |5 D) ^
( c' O1 `4 J' F$ ]5 f6 u: U#选用AIC准则下的模型进行回归诊断& }% M! p; y7 j1 f; p! F! o0 h
windows()4 j; f& b% `0 c ?
par(mfrow=c(2,2))4 P, u; z4 G/ _% s
plot(lm4.aic,which=c(1:4))
' h T" U3 U5 A+ d5 ^8 W. S6 N0 s/ n& s
7 P6 E* t5 K2 W% s" c* j) u g( E9 t# d9 {
T0 r, D+ f( G" f; e# K
##数据变换" Z' V! \7 ~9 y& ?4 s
) ^/ g3 u3 @ s: f0 [* k! R5 y7 v" [/ W' h& L" Y
#box-cox变换
( n7 z) e. [ F* \( I) T5 Ylibrary(MASS)$ {( C9 k+ e, N4 X& w7 u& f( l
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))" ^4 S9 ^, K" ~) s
I=which(b$y==max(b$y)) #定位似然函数最大的位置
# k# Z/ U7 a: t( _9 mlambda = b$x[I] #精确的λ值
2 O5 f5 ]: s( o1 k2 j' I7 _3 S#λ接近于0,为模型简洁性,可以直接进行对数变换! R7 Q. _9 ?; J( ^9 t
logprice <- log(price)
& w3 ~' ]% Q$ B' o& m* Lhist(logprice)1 h, g+ ~# z% j+ V0 t
8 y, d9 Z3 N2 }( f
' ~3 Z% m0 f) o. g+ T/ W##最终模型与诊断% G# d/ @" g& ]
% v1 I6 H9 p" D' i- w& R
4 x0 O" Z% \6 o- h0 t# \/ Plm6=lm(logprice ~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(ring)*as.factor(ronggrp)+lv+ratio)
$ [: w1 S, L, twindows()
- W! Q: G2 }" [! y7 \: xpar(mfrow=c(2,2)), M. A# \: |2 X; n7 z- T! P
plot(lm6,which=c(1:4)). r* c) K( S. A! A* s+ v" P2 C" `
anova(lm6)
b: u3 p8 D2 x% C3 @4 Hsummary(lm6)
0 p, ?; x4 w/ S" T! a7 G, k9 R/ d( X7 t+ U2 c+ t
8 A* S) D) l9 w( q# S
请关注数学中国网微博和数学中国公众号,联系QQ 3243710560; Y. f6 ]: P$ j
) p) N7 F7 J; l8 z% g, L+ A
3 L+ z1 X( \; A3 I! j5 F5 Q! i
% j2 w# p' A# u4 J. M) X" A* q" r5 T: Y
7 {% \- U& Q% w |