【R语言】回归分析案例：北京市商品房价格影响因素分析
" A5 Q. F& z0 u/ o这一案例是王汉生老师《应用商务统计分析》方差分析章节的案例，主要对离散型变量进行了处理。
* P- Q, J2 m) d, Q K这里将连续型变量也加进来，进行协方差分析，建立完整的模型。 首先对房价进行对数变换，解决异方差问题： 
: `' _. j1 x: I l) \ y0 Q行描述性统计分析，各连续型变量之间的相关关系如下：
* R( B) H: O/ ~* ~! b4 P$ q1 s* p ; A$ h! Y9 L: Q6 s
名义变量的EDA一般做箱型图。 模型按照全模型-变量处理（分箱等）-变量选择-回归诊断等步骤建立。
R2 Q4 h3 s3 {5 O* ]7 m+ j
$ d2 b" v! H# c2 J 4 t8 W! s* }1 B2 L, K: o
最终模型残差图：- C4 {* U1 v0 r4 t" `/ s
0 o( P3 B7 g6 U; Z2 Q/ k2 e& q通过模型分析结果可知，影响北京市商品房平均销售价格的主要因素有：- ~1 P) d. [9 N! W8 \1 W( A" K
属性变量：所在辖区、所在环线、物业类别、装修状况、容积率大小（新引入）；连续变量：绿化率、停车位住户比' B; S; b% [3 t8 q; t
属性变量的具体影响在此处分析略去。% B! t8 w+ X3 t, J; R
连续型变量的影响主要为：
' j* t, n( n2 c# `# { 绿化率：绿化率的影响十分显著，由系数估计值为正，说明对房价有正向影响，绿化率越高的楼盘房价越高；1 `. W( A! V5 d* m
停车位住户比：有较显著的影响，停车位住户比越高，价格越高；
0 }( U+ ~1 N7 a: q B4 M+ I同时，原本为连续型变量的容积率经过离散化变为属性变量后：3 N$ O% t; t3 j: ]# \7 V
容积率大小：容积率分组有较显著的影响，高容积率的小区商品房价格更贵；1 C, P* H& O; h' ?' G' W
容积率与环线之间存在着交互效应。
3 @9 a" T5 l" d/ _- B; G Z+ |0 \rm(list=ls()) #清空当前工作空间
/ R1 {2 P, N( @# n7 gsetwd("D:/回归分析")
" c% `1 u0 ~; B9 K z: l. ua=read.csv("real.csv",header=T) #读入csv格式的数据，赋值为a
" L$ o+ p8 _! u& Y% L% L" B5 l: R$ }View(a)
# U* L" W! q. ?' }$ b2 Nattach(a)
& O W) ^+ ?$ v: ~names(a)1 x. K3 A; _/ [1 G: @
8 I! a' O. r4 e! U
! ]+ M3 Q- L( c, C1 z##描述性统计" V) c8 r9 w) z/ H1 A, a2 Z
) `: D4 A$ l( I
Z! M' Q! U/ f1 b/ {$ r#未做处理的响应变量分布情况6 ~( g4 A4 ]: f) {" ?$ p
par(mfrow=c(1,1))
+ r7 E% k' h! c6 E; rhist(price)
2 n; ~ O% @& v4 x9 k7 J- Y3 o! O- qsummary(price) #查看响应变量的描述统计量1 u5 \( F3 Z, u% s6 m( W; y' f
#连续型变量描述性统计$ j! U) x* r, s( U; A* z' X
windows()
" s7 \7 _( O# Y9 K2 [/ o( Hpairs(a[,c(6:10)]) #所有连续型变量间的散点图4 H* F8 b4 ]' ]" |) c; u
par(mfrow=c(2,2)) " @( O7 S/ o* {% Y
plot(rong,price) #每个连续型因变量与响应变量间的散点图
2 d+ F% Y. g6 c! d: l E2 |plot(lv,price); l1 Q8 B; d( k( |+ U& A
plot(area,price)# {: L8 e% h0 X* f' X" }
plot(ratio,price)
( d0 j2 X) U4 z0 zsummary(a[,c(6:10)]) #查看连续型变量的描述统计量( K; X4 h1 V! p: T4 e+ A( ^
cor(a[,c(6:10)]) #查看连续型变量的相关系数* W, h+ F3 r. L
#属性变量描述性统计. [. f1 m& n- Z Z n# [
windows()
: Z# M1 }4 S X( h& i% Lpar(mfrow=c(2,3)) ( ^/ |$ f6 T3 h. U5 P' {
boxplot(price~dis) #每个属性变量关于响应变量的箱型图
" ]( b8 g7 G2 Gboxplot(price~wuye) ) m6 g5 z' q7 _3 g
boxplot(price~fitment) 8 p; t1 q# \3 O- W5 H# q
boxplot(price~ring) 6 ?! K3 l1 k0 }$ R9 q7 Z4 j
boxplot(price~contype)8 w! n3 ^/ c' b6 ?+ b
2 j( G6 Z6 J3 Q' G" B2 W
+ \, G2 j7 H2 J+ \1 @ P1 Y# s% y3 L' t3 D6 J
1 [2 J$ h+ }! z3 z" o& l
##模型建立
3 R) g/ U4 n& a# A: ]
9 f/ s# t, a. j* r
' m7 ^/ e! i) P& R- e#在方差分析模型基础上加入连续型变量- c# S; B/ c8 s- c
lm1=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+rong+lv+area+ratio)
0 i6 X1 Y( S1 ~# L: nanova(lm1) #方差分析; h$ m' b0 M l- w% r
summary(lm1) #模型参数估计等详细结果
8 V9 J4 s1 d) Y+ s8 S& B6 I" T' ewindows()
- Y, X! \) Y6 Q8 g+ E5 `par(mfrow=c(2,2))1 K2 }0 R6 i. Y3 X6 n% c/ ^9 d ~
plot(lm1,which=c(1:4)) #回归诊断做残差图8 L5 s: i: W1 H5 U* c# L
8 B' a5 P/ q" ~! V
+ j6 U e# X. h
$ L% `9 K( H: M; b6 T7 M5 ~
- f0 [: u* K. ~, f3 g1 P##变量处理1 H F5 b: l( x; X+ C) f) L( H
. f0 e" N# L* }& Z! A5 n: M
7 X/ _0 c7 A- L+ `8 b$ F
###对不显著的变量采用分组的方式希望能达到显著的效果$ y1 l( [, \" m& q, i* l
##对容积率的处理
6 ~: z4 N O9 s$ `6 T( ^windows()
" r) h( `! @) n- Q; M/ O- vn = 48 D; e! H8 @3 m3 k$ m* [* O A: v3 c
boxplot(price~ceiling(rong/n)) #容积率多分组下的箱型图 - Z% B* j3 o$ W/ c5 i
table(ceiling(rong/n)) #容积率各分组下的样本数 k R6 U5 I) |4 G' A# }. g( ]
ronggrp=1*(rong>n) #进行二分类" ^, x' T" x9 L* _: W
#ronggrp=ceiling(rong/n)
, x% S' L8 p) s( Stable(ceiling(ronggrp)) #容积率二分类下的样本数! K+ a0 p- K8 X- \/ \# H5 x. D
windows()9 h4 t3 d+ q- F; g2 |' J
boxplot(price~ceiling(ronggrp)) #容积率二分类下的房价箱型图4 ]3 V3 O0 L/ }4 E6 W6 U4 w% {
windows() O- Q; b; r4 l2 S7 \* g, z
par(mfrow=c(1,2))
* v! t: d# F( Jboxplot(rong~ring) #容积率与环线箱型图6 D: o8 h1 i2 e' X
boxplot(price~ring) #房价与环线箱型图
+ _7 I a8 J& X* k, c1 K# e$ `#加入容积率分组和容积率分组*所在环线交互因子的模型
& C- t0 [8 Q; slm2=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)
5 j" ~* t/ j7 k! L9 r4 lanova(lm2) #方差分析
0 V6 j: u' o" ~7 Isummary(lm2) #模型参数估计等详细结果
# z2 A! B C# Z t' x$ Cwindows()6 _5 i$ t3 q' L! V4 I; h8 p6 b
par(mfrow=c(2,2))* V7 |9 B8 R' G) | U- j
plot(lm1,which=c(1:4)) #回归诊断
( p2 A0 S, [0 r' `( p! M' u( N. c% ^% F
$ Y% x& p9 `3 @* u##对小区面积的处理
1 ~; L) J$ T1 q* j, `summary(area)3 u( ^& M. Q5 e [0 [& V
plot(area,price)8 q6 C& N* C8 N0 S, a: q
windows()1 L- t- J) z/ d+ l! ^# `
n = 150000! E8 a8 G E8 {3 p
boxplot(price~ceiling(area/n)) / U: E( ~; W# c/ P9 s- W C
table(ceiling(area/n)) ) S8 t1 C8 o# m* m5 ^/ M9 v$ K3 i" f( ~5 Y
areagrp=1*(area>n)
4 ^9 a: r$ b. G& Gtable(ceiling(areagrp))
8 B$ E8 O! H# e- } k. o, |. h1 \boxplot(price~ceiling(areagrp))+ I- l! ~ w% v
#加入小区面积分组的模型) w6 o1 I2 }3 U7 r, B3 B
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)
8 i" _! _9 H1 Sanova(lm3) #方差分析* r4 q6 U3 I5 s1 Y" {- z
summary(lm3) #模型参数估计等详细结果9 r$ z4 A, Q0 A! C- ]) e1 {8 A% o- D' }% k
windows() R$ l# g+ p. e2 C w8 z5 k
par(mfrow=c(2,2))
. d% ] [0 ~$ i6 e" h; p# V# mplot(lm3,which=c(1:4)) #回归诊断 u& B3 p0 [2 l9 U1 P4 k
4 l3 y3 |' M1 c
, c% P p/ c3 Y$ h
##变量选择8 I1 o. H" a) _- i8 F0 Y, s
" h; S: ?5 _ R: F
/ c' a. U/ P! i& Z##AIC准则下的变量选择 u6 F' I; w9 [- B: m. R) g9 C( G
lm4.aic=step(lm3,trace=F) #根据AIC准则选出最优模型，并赋值给lm.aic
$ W; _: C# r. x8 d" Jsummary(lm4.aic) #给出模型lm.aic中系数估计值、P值等细节
/ D8 W- ^. a3 j. ^, i##BIC准则下的变量选择4 J- O" e9 `% g% a
lm5.bic=step(lm3,k=log(length(a[,1])),trace=F) #根据BIC准则选出最优模型，并赋值给lm.bic+ P* [# x( H) i0 h) a) M4 B
summary(lm5.bic) #给出模型lm.bic中系数估计值、P值等细节
2 d$ b: k# ?9 v- D D
0 f5 }1 U+ D! _, ~2 F, i, l3 ^1 k; A7 c1 \- U. ?
#选用AIC准则下的模型进行回归诊断& J+ P0 ~/ H% H
windows()
5 o! R7 a1 Z+ ]0 q/ `5 ipar(mfrow=c(2,2))
8 _$ G* Q2 K5 F9 i6 G4 i8 ^2 \( x6 tplot(lm4.aic,which=c(1:4))
6 h9 k5 ?+ p) q* i$ P6 Q5 M% K; F% G
% J3 l6 F6 R2 ~) c# K1 [$ c9 l% s, p( o7 V; v$ E6 D3 k- h* j
$ c$ I5 X$ M7 y T9 b, v7 V, O7 l7 r! m, w4 D% g, k5 x1 k
##数据变换
5 s/ r" O2 t' t/ w @& w
' `& l" d( o) Q W; g: R8 w- H( e# \0 M; l+ p- a- x0 e
#box-cox变换+ ~) O8 q, k7 x* M
library(MASS)
" L8 C0 t- \8 ?6 c! Lb=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))- k- {+ G7 v# C! L4 q T
I=which(b$y==max(b$y)) #定位似然函数最大的位置
+ X- K! n& ]8 W9 _/ o) mlambda = b$x[I] #精确的λ值2 c3 Q/ _; K3 ~ H
#λ接近于0，为模型简洁性，可以直接进行对数变换
2 U4 }4 W/ ]* W$ mlogprice <- log(price)
" x" O+ I7 i! b6 ihist(logprice)
7 @/ Z/ h! t" g2 f1 i0 ]6 F; t( z9 h
- A$ t2 @" s& f##最终模型与诊断' @4 C1 a- V# |4 @+ a
p; p. \6 W6 E; s4 g
. N% t- x. g3 {3 N1 o( P& b- U0 t
lm6=lm(logprice ~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(ring)*as.factor(ronggrp)+lv+ratio)
+ p. h$ V% F7 X1 F8 A5 K( u% Ywindows()
+ U6 W' T$ O6 Zpar(mfrow=c(2,2))
( U- X0 Q* N6 n j0 hplot(lm6,which=c(1:4))/ \" Z& t" ^% I7 }3 I( A: i; L, V
anova(lm6) S* W$ C' p" k6 |' _
summary(lm6)) J0 ` H; T9 ^8 O
3 l0 j- H/ d1 Z& N6 o+ t2 A- |, Y v: l1 d# V& a3 t: q2 C& n
请关注数学中国网微博和数学中国公众号，联系QQ 3243710560' Q' f& v' E* V% Z, A i9 _
& M" A u" e# B6 g3 _. b* @# q
% R( v2 t. h9 u& ?3 W
9 D7 S( t0 T, ^8 \
4 T# s; n; _) R5 E |
IP卡
狗仔卡
发表于 2021-10-27 14:50
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
显身卡
发表于 2021-10-27 16:32
