在线时间 514 小时 最后登录 2023-12-1 注册时间 2018-7-17 听众数 15 收听数 0 能力 0 分 体力 40220 点 威望 0 点 阅读权限 255 积分 12777 相册 0 日志 0 记录 0 帖子 1419 主题 1178 精华 0 分享 0 好友 15
TA的每日心情 开心 2023-7-31 10:17
签到天数: 198 天
[LV.7]常住居民III
自我介绍 数学中国浅夏
【R语言】回归分析案例:北京市商品房价格影响因素分析
, L8 ~. v/ @+ t 这一案例是王汉生老师《应用商务统计分析》方差分析章节的案例,主要对离散型变量进行了处理。; d# h4 A `4 e B4 N d
这里将连续型变量也加进来,进行协方差分析,建立完整的模型。
首先对房价进行对数变换,解决异方差问题:
5 {: c9 V4 [& w4 Q( e
行描述性统计分析,各连续型变量之间的相关关系如下: , `, }" @2 b- e( @) G
& ~( |/ Y+ t/ k! \8 L- K 名义变量的EDA一般做箱型图。
模型按照全模型-变量处理(分箱等)-变量选择-回归诊断等步骤建立。
/ T8 f8 W9 @" g
, ^% u( ]& a. H/ O& y+ M7 N 0 ]* F- e! z9 M' c8 `8 a1 R
最终模型残差图: ( i, T4 C0 A% ^% b
% l7 }7 P. w5 X4 s4 f
通过模型分析结果可知,影响北京市商品房平均销售价格的主要因素有:
$ X, @( f8 ?/ S9 ~1 }" J) K 属性变量:所在辖区、所在环线、物业类别、装修状况、容积率大小(新引入);连续变量:绿化率、停车位住户比
, V* w0 u/ g+ y* x& A3 X; _ 属性变量的具体影响在此处分析略去。 . s3 i! |# I. Q% F0 C- i' M
连续型变量的影响主要为:
; M7 c/ [3 y$ P* A( j 绿化率:绿化率的影响十分显著,由系数估计值为正,说明对房价有正向影响,绿化率越高的楼盘房价越高; 4 o! _, V" T6 ^7 Z8 {5 V/ S
停车位住户比:有较显著的影响,停车位住户比越高,价格越高;
) ?! m4 Q( G, }8 `: W# k 同时,原本为连续型变量的容积率经过离散化变为属性变量后:
* n& i4 T$ `! | g" b: ? 容积率大小:容积率分组有较显著的影响,高容积率的小区商品房价格更贵; ) U8 z* q8 H% ]. X
容积率与环线之间存在着交互效应。
0 y$ j" a( r! e7 O7 @ rm(list=ls()) #清空当前工作空间
& b G0 V# r8 U" i( k6 p) b8 { setwd("D:/回归分析")
8 J+ j# ~+ R* T% N R a=read.csv("real.csv",header=T) #读入csv格式的数据,赋值为a ; z6 Y. u. S H8 b5 e9 d5 v
View(a)
8 f0 @/ ~7 P' D- \. O) y: ~1 p attach(a) ! N" p3 M( D! Z3 h
names(a) # ~' ~4 l9 u# H2 O) ?
; B' q' O- K5 b9 ?( ^/ {: X 9 V. C3 i+ ~. T7 K
##描述性统计 v3 N: s: ~% F
6 `2 {2 F4 j3 m2 ^) m' ~& b: l5 E3 u
1 B1 W h$ O- l! K7 Q #未做处理的响应变量分布情况 ) B3 U2 V: O0 | ]# j
par(mfrow=c(1,1)) , {( z" M- q+ H( l
hist(price)
4 E2 \+ g/ O+ F4 r+ j) B# q* X summary(price) #查看响应变量的描述统计量
$ \8 F5 J2 i; O/ L, _, U #连续型变量描述性统计 2 q2 R; h# { c$ t1 ]- B
windows()
- Z! P% f; i, A! \ pairs(a[,c(6:10)]) #所有连续型变量间的散点图 , B) N- _1 t* g: h( A# k/ s
par(mfrow=c(2,2)) ; k% t! Y, |- f I S/ B
plot(rong,price) #每个连续型因变量与响应变量间的散点图
Y7 ~# F; R: R3 n5 j6 V plot(lv,price)
; k! F- b8 C: ?! y plot(area,price) R' D% a: q' \. n% _- Y! B
plot(ratio,price)
7 K0 }$ Y7 E# q summary(a[,c(6:10)]) #查看连续型变量的描述统计量 - L- c2 b5 a8 [- s! Q0 t, h: e4 e
cor(a[,c(6:10)]) #查看连续型变量的相关系数
6 y# B. t* n% B g% l# y #属性变量描述性统计
: Y, k. Y1 ~7 I9 g windows()
: _0 D! X9 o! ?- o: { par(mfrow=c(2,3)) - Y0 `6 r+ ~( v ?! R& V. q, g; K) K5 I
boxplot(price~dis) #每个属性变量关于响应变量的箱型图 * A* d* O7 H; {1 [' W4 j
boxplot(price~wuye) . q% ^0 y Q- f5 }* B/ y2 ?* V
boxplot(price~fitment)
5 p, p: N8 p. S [$ w boxplot(price~ring) ( |, z1 D% D6 d* M W8 H
boxplot(price~contype)
& P5 }6 W( _' G- O1 V
0 }9 i* j- V9 r, {. ~$ m" J! o
& L3 Y0 m0 j) O3 x
1 c* c# L$ Q, O$ P. w! V" R
. X- {1 @6 a. A ##模型建立 % `+ M9 i5 y4 |* g( ]
6 p. c+ \# K% J; \* o5 v/ L3 T
, u- D( E1 \, V, F q; |8 J
#在方差分析模型基础上加入连续型变量
# l, |& e8 _. U- s) H" Y lm1=lm(price~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(contype)+rong+lv+area+ratio)
6 J5 c2 X7 V. _ anova(lm1) #方差分析
: D- a" s1 K+ }9 e8 G summary(lm1) #模型参数估计等详细结果
4 A, j; M3 h( j windows() & O( M) N; L$ j0 L- ?6 \- m
par(mfrow=c(2,2))
+ y" M% {# |, s% d" ~ plot(lm1,which=c(1:4)) #回归诊断做残差图 9 Q8 R# E4 d* l: E0 S
, F6 {/ f9 \4 Q
- ^1 R, U9 A( a
" t& U" m, t: A0 v. _4 y s9 i' U ; B/ u2 z, \* S+ p7 J8 A# }
##变量处理
% ?. ?$ K3 k) M% g7 r , d5 H1 R3 J/ b$ T
- k8 {0 |( U% P( Q7 o ###对不显著的变量采用分组的方式希望能达到显著的效果 0 j2 ?& y: z$ z/ X O3 t
##对容积率的处理
- t% H. j: {4 k windows()
@& N/ z/ J7 ~1 t n = 4
" q" U: {( _: P2 m4 o boxplot(price~ceiling(rong/n)) #容积率多分组下的箱型图 ' N3 L! r. c$ G5 b
table(ceiling(rong/n)) #容积率各分组下的样本数
( D2 y* O' R7 E0 n6 g$ D ronggrp=1*(rong>n) #进行二分类 - u6 }0 l! D2 q @% e' R
#ronggrp=ceiling(rong/n) * h9 k" B( }4 c p; P
table(ceiling(ronggrp)) #容积率二分类下的样本数 # k: n' j3 ?$ Q; J9 q
windows()
9 y& o' I- S% t boxplot(price~ceiling(ronggrp)) #容积率二分类下的房价箱型图
! W- ~/ k( ^2 @' [6 R- a8 x Q( B windows() $ s* F: j* I' R' u5 x
par(mfrow=c(1,2))
2 e, t! j6 B* f boxplot(rong~ring) #容积率与环线箱型图 1 u0 R: s5 w) Q, P5 a% j8 ~0 m9 Q
boxplot(price~ring) #房价与环线箱型图
$ a' _; d2 u2 ?: G) c5 | #加入容积率分组和容积率分组*所在环线交互因子的模型 ! q0 Q8 B, p" z
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) 7 q/ |, [* [- }6 c6 F. {
anova(lm2) #方差分析
, S1 L; P/ ~' N# a$ E% |7 @% Y summary(lm2) #模型参数估计等详细结果 # e$ F7 @, E: Q0 U7 d1 o
windows()
8 E! M+ W F. s' @ G7 Y" E par(mfrow=c(2,2)) X3 D- l9 I9 Y5 B% v( s( K7 G8 l
plot(lm1,which=c(1:4)) #回归诊断 / U/ s Z n" B( }$ C! W9 [
, C, }0 n7 I8 I" H
5 U% c! G* H' y- ~" U: X/ i' w N
##对小区面积的处理 ' J6 o& H" V; i
summary(area)
, x" I V) N- V& R( ?0 r plot(area,price) + k6 K& ^5 C" \' d! u
windows()
x0 k/ c5 j8 c. d n = 150000 + `7 q5 N x3 R9 r+ [7 ]1 ^& F
boxplot(price~ceiling(area/n)) 4 r2 a3 N: o4 [: c0 K7 D+ I
table(ceiling(area/n))
. l% B' D$ V S areagrp=1*(area>n) ' A: M4 b+ _4 q1 Y X0 v
table(ceiling(areagrp))
% @2 V) b p$ r1 `* c boxplot(price~ceiling(areagrp)) ) b0 v+ q+ E$ L- k; g
#加入小区面积分组的模型
2 J5 i8 D0 b! a, a* r 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) & \6 m1 j% m* w
anova(lm3) #方差分析 ( {( y0 }' c9 j$ F4 z6 D8 f
summary(lm3) #模型参数估计等详细结果
% [4 }" P$ s/ I7 d9 Z9 @ windows()
6 ~5 Z. @* \% e par(mfrow=c(2,2)) # |& i1 x' G, D9 F( s- n
plot(lm3,which=c(1:4)) #回归诊断
+ l8 O) d( `& p/ e* [0 U
Y' y6 v( v9 E: E( |/ e , ?2 j6 u# q7 [9 f9 k3 y' P% n
##变量选择
, ~2 `9 Y8 D/ ` 5 i' n9 `: x6 s" A D" V
7 G' |9 Q( P1 t0 w' [ ##AIC准则下的变量选择
/ b8 X; B& |( D' f8 i/ G lm4.aic=step(lm3,trace=F) #根据AIC准则选出最优模型,并赋值给lm.aic
0 _3 M h* `' D1 a summary(lm4.aic) #给出模型lm.aic中系数估计值、P值等细节
7 |, p8 _, k" w3 v. B6 s6 X6 U5 ^0 h ##BIC准则下的变量选择
5 e9 m0 O5 R4 k2 Q5 v8 D lm5.bic=step(lm3,k=log(length(a[,1])),trace=F) #根据BIC准则选出最优模型,并赋值给lm.bic ' m; U8 K0 L; ]+ f9 A
summary(lm5.bic) #给出模型lm.bic中系数估计值、P值等细节
; K& u2 M' {& p) X$ L$ |% s $ t2 G! ~2 ?2 Q/ j7 Q0 ]
7 W/ T6 D/ B/ ?* N5 D4 K2 j8 ]4 M #选用AIC准则下的模型进行回归诊断
: }$ }9 A' _3 h3 k+ P& Z% F windows() * M( ]$ s2 ?0 Q& @0 w
par(mfrow=c(2,2))
: p! k _. W& u# X. e& X, ^$ k plot(lm4.aic,which=c(1:4)) 8 A- c8 X Y& M L4 ` V
- J A1 ?! P; U7 W
9 f9 i9 e7 A. {6 F7 u
; | s& n/ b# |+ @5 u" l0 X2 T2 Y
9 c/ U& V8 T( ]1 i# ^ ##数据变换
( q( F$ i) e) c8 v% J2 T C" e
$ k! e8 c+ o+ K* z: F& l 8 U( U# G2 O3 S8 i: r
#box-cox变换 ; V% l% \# _5 e, d. L' T7 s
library(MASS) * a4 ]! s6 h7 i' O3 q) x5 i8 ?
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))
% L7 E% w5 N- h4 ]' t# v I=which(b$y==max(b$y)) #定位似然函数最大的位置 , R5 @( V X8 j- s1 c: _
lambda = b$x[I] #精确的λ值
) ~8 ?5 {, J; @" P3 E #λ接近于0,为模型简洁性,可以直接进行对数变换
3 e. L% ]+ O2 h- \ logprice <- log(price) 0 h( t- F4 Q/ U: Y1 A
hist(logprice) P' l! b" Y; J0 t. d/ P, v
7 i$ V, k9 W: A5 z4 j5 c6 e * Y$ Z, s/ C( K' T: E8 ~# Y
##最终模型与诊断
7 G5 h7 [( l4 t, `
0 }/ m. b4 w6 q2 D. a
( A8 S8 t6 k' j; g* }/ y. ~ } lm6=lm(logprice ~as.factor(dis)*as.factor(ring)+as.factor(wuye)+as.factor(fitment)+as.factor(ring)*as.factor(ronggrp)+lv+ratio) H3 }( |) c# ?3 {3 o$ R' ^
windows()
1 Y9 E$ E/ z: v Z( _ par(mfrow=c(2,2))
! w& \4 W0 `' u plot(lm6,which=c(1:4)) , A4 z7 i% s% ?7 J0 d. g" I# [) s5 _
anova(lm6) ' V/ l# U; ]9 ~. F, d* e$ E* W8 N
summary(lm6) , n1 U4 A0 c: Q( l$ F
; y+ p+ R4 Z. ]- m$ ]* u, z
J0 w, h. `5 u5 E5 @. ` 请关注数学中国网微博和数学中国公众号,联系QQ 3243710560$ _- f4 h5 \' h( h4 k
! v) c T9 A) f, w0 g
- Y) `6 A7 p8 b% ^$ `
0 j; I$ Z0 Y. W7 s- E9 r T
9 n, s5 \; m# {6 W5 e
zan