R Linear Regression (線形回帰) (2) 予測モデルの作成と評価

分析の題材はこのデータ：

#SPC:Sales per customer,顧客単価(JPY)
#OpenH:Hours open,営業時間
#NOHWT:Number of households within trade area,商圏内世帯数(世帯)
#NOPWT:Number of people with in trade area,商圏内人口(人)
#IPH:Income per households within trade area,商圏内世帯の平均世帯年収(JPY)
#PPT:Profit per Tubo,坪当たり営業利益(JPY)
#線形回帰で、顧客単価を予測するモデルを作成することが目標

SPC	OpenH	NOHWT	NOPWT	IPH	PPT
32512	8.5	6211	17360	3200000	10125
40122	9	6256	16891	5920000	13574
39213	9	6787	18865	3000000	18215
41211	8.5	6982	18852	5200000	18325
50211	10	4521	12207	5700000	15331
38122	8.5	6781	18309	4920000	15215
37543	11	7325	19778	5220000	17211
34511	8.5	6211	16700	3280000	9253
35223	11	6521	16878	4220000	18222
34322	8.5	5251	14178	3250000	16252
35871	8.5	5293	10292	3000000	10321
36921	8.5	5621	14177	3010000	13221
36211	10	6221	16797	3210000	17522
34841	8.5	5822	13719	3720000	13376
42355	10	6351	17148	4970000	14389
34211	8.5	5651	14296	4660000	13921
36662	9	5721	13447	4750000	15265
42121	8.5	5825	15726	5350000	10321
43212	9	5651	15258	6360000	14796
36672	8.5	5723	15452	3000000	15223
35213	10	5992	16178	5350000	16351
30102	8.5	5825	15728	3250000	13543
38218	8.5	6311	17067	4990000	15531
38225	9	6321	17040	5370000	14762
33421	8.5	6251	16888	3580000	13251
35611	9	5321	18367	4330000	15767
36213	8.5	5651	15367	4790000	12226
37322	10	7321	19767	4820000	17251
36276	10	5632	15206	4820000	14378
45321	8	5571	15042	4920000	12765
30221	8.5	5221	13997	4520000	11311
36337	11	6921	15987	5220000	15225
36452	8.5	6973	18730	4220000	15217
38628	9	5323	15372	4970000	15665

#########1.データの取り込み,整形#########
raw <- read.csv("4_Linear Regression.csv") #上記データ
head(raw)

#ヒストグラム
par(mfrow=c(6,1),mar=c(0,3,3,0),ps=10)　#レイアウト調整
for (i in 1:ncol(raw)) {
  #ヒストグラムの数が増えると余白の調整が必要。Zoomしてみると良い
  hist(raw[,i],main=names(raw[i]))
}

#欠損値の確認
sum(is.na(raw))

#全て量的変数の場合はplotで傾向を把握
plot(raw) 

#k-meansでplotも可能。クラスター分析の演習参照
par(mfrow=c(1,1),mar=c(3,3,3,3))　#レイアウト調整
cl <- kmeans(scale(raw), 3, 20) #時間があればWSSのグラフも描いてみて下さい。(各自コーディング)
plot(raw, col = cl$cluster)
points(cl$centers, col = 1:3, pch = 8)

#相関係数
cor(raw)

#########2.顧客単価を目的変数として線形回帰#########
lm_a <- lm(SPC~.,data=raw)
summary(lm_a)

#(a)説明変数の選択
lm_b <- step(lm_a) #AICに基づいた変数選択法(デフォルトはbackward)
#変数選択の観点は複数あり、吟味が必要。ここではとりあえずAICを選択。
summary(lm_b)

#(b)モデルの評価
#(b-1):残差分析
par(mfrow=c(2,2),mar=c(4,4,4,4),ps=15)
plot(lm_b) #残差分析。外れ値の行番号が表示される

#外れ値を外して再度回帰分析
raw_rv <- raw[-5,]
lm_c <- step(lm(SPC~.,data=raw_rv))

summary(lm_c)
plot(lm_c) #残差分析

#Residuals vs Fitted values:SPC:目的変数値に応じた残差の値を縦軸に表示。
#NormalQ-Q:観測値の標準化後残差が正規分布に従う場合の期待値をx軸、
#観測値の標準化後残差をy軸にとったプロット。
#Scale-Location:Residual vs Fitted valuesのy軸について平方根をとったもの。
#Residuals vs Leverage:Leverage(てこ比)とは回帰分析の観察点(サンプル)毎に
#説明変数のデータを変えずに目的変数yの値を1だけ変えたときの予測値の変化量
#cook's distanceは検証用に使われる距離

#stepの結果単回帰直線となったので、図示
par(mfrow=c(1,1),mar=c(5,5,5,5))
plot(raw[,1]~raw[,5],data=raw_rv,xlab="商圏内世帯の平均世帯年収(JPY)",
     ylab="顧客単価(JPY)",main="顧客単価の回帰分析",
     xlim=c(2800000,6500000),ylim=c(25000,55000))
abline(lm_c)
par(new=T)

#(b-2):95%信頼区間　(回帰直線の存在領域)
lm_con <- predict(lm_c,interval="confidence")

#図示出来るように昇順並び替え
con_lwr <- data.frame(raw_rv[,5],lm_con[,2])
con_lwr <- con_lwr[order(con_lwr[,1]),]
con_upr <- data.frame(raw_rv[,5],lm_con[,3])
con_upr <- con_upr[order(con_upr[,1]),]

par(ps=20) #フォント調整
lines(con_lwr,type="l",xlim=c(2800000,6500000),ylim=c(25000,55000),ylab="",xlab="",col="blue",lty=1)
par(new=T)
lines(con_upr,type="l",xlim=c(2800000,6500000),ylim=c(25000,55000),ylab="",xlab="",col="blue",lty=1)
par(new=T)

#(b-3):95%予測区間　(予測値の存在領域)
lm_prd <- predict(lm_c,interval="prediction")
lm_con <- predict(lm_c,interval="confidence")

prd_lwr <- data.frame(raw_rv[,5],lm_prd[,2])
prd_lwr <- prd_lwr[order(prd_lwr[,1]),]
prd_upr <- data.frame(raw_rv[,5],lm_prd[,3])
prd_upr <- prd_upr[order(prd_upr[,1]),]

lines(prd_lwr,type="l",xlim=c(2800000,6500000),ylim=c(25000,55000),ylab="",xlab="",col="red",lty=1)
par(new=T)
lines(prd_upr,type="l",xlim=c(2800000,6500000),ylim=c(25000,55000),ylab="",xlab="",col="red",lty=1)
par(new=T)

par(ps=10) #フォント調整
legend("topleft",legend=c("95%信頼区間","95%予測区間"),col=c("blue","red"),lty=c(1))

#########[Additional]3.ブートストラップ法の利用#########
intr <- numeric(10000)
coef <- numeric(10000)

for(i in 1:10000){
  trainNO <- sample(nrow(raw),nrow(raw)*0.67,replace=TRUE)　#replace=TRUE:非復元一様分布からの抽出
  train_dat <- raw[trainNO,]
  res <- lm(train_dat[,1]~train_dat[,5])  #train_datで回帰分析
  intr[i] <- res$coefficients[1]   # 定数項の取り出し
  coef[i] <- res$coefficients[2]   # 係数の取り出し
}

par(mar=c(3,3,3,3),mfrow=c(1, 2)) #余白幅の調整,グラフ画面の分割
hist(intr)                        #定数項のヒストグラム

#95%の信頼区間の境界値を書き加える
abline(v=quantile(intr, c(0.025, 0.975)), col="red")

hist(coef)                        #係数のヒストグラム
#95%の信頼区間の境界値を書く
abline(v=quantile(coef, c(0.025, 0.975)), col="red")

#########[Additional]4.主成分分析#########
raw_s <- scale(raw) #データの標準化
cor(raw_s)　#相関行列
PC <- princomp(raw_s)　#主成分分析
summary(PC)
attributes(PC)
PC$loadings

#第1,2主成分でプロット(69%を説明)
par(mfrow=c(1, 1))
par(ps=15,mar=c(5,3,3,3))
plot (data.matrix(raw_s) %*% unclass(loadings(PC))[,c(1,2)])
biplot(PC)