搁デモ03冲データ分析编2016
- 1. #-------------------------------------------------------
# R03_データ分析編_擬似ミクロ.r
#-------------------------------------------------------
rm(list=ls(all=TRUE)) # 作業領域のお掃除
#...............................................................................................
# i. 入出力編で作成した、擬似ミクロとコンスタントを収録したワークスペースファイルを読み込む
#setwd("Y:/PC123/giji/") #要修正%%%
setwd("D:/R 研修 2016/giji/") #要修正%%%
load("giji.RData")
ls() # ロードしたワークスペースファイルの内容表示
attach(dt1) # データフレーム dt1 を attach
#...............................................................................................
# ii. 主食分類符号の作成
# 主食に関するデータを確認。該当項目は、Youto039 から Youto043 で、データフレーム内では 53 から
# 57 番目にあることは、コンスタントデータで確認できる
cd1[53:57,]
# ren koumoku lvl type namae
#53 205 穀類 6 1 Youto039
#54 208 米 7 1 Youto040
#55 210 パン 7 1 Youto041
#56 212 めん類 7 1 Youto042
#57 214 他の穀類 7 1 Youto043
# 想定されない数字以外の文字が入っているか確認 => 全てなし
is.character(Youto039); is.character(Youto040); is.character(Youto041);
is.character(Youto042); is.character(Youto043);
# 総合計と内訳の和の差異は誤差の範囲
sum(Youto039) - sum(Youto040) - sum(Youto041) - sum(Youto042) - sum(Youto043)
# [1] -0.004349995
# NA は存在しない
length(which(is.na(dt1[,53:57])==TRUE)) #[1] 0
#...............................................................................................
# A. 主食分類符号の作成
# 主食インデックス（g1-g3）の作成
g1 <- which(Youto040 > Youto041 & Youto040 > Youto042 & Youto040 > Youto043)
g2 <- which(Youto041 > Youto040 & Youto041 > Youto042 & Youto041 > Youto043)
g3 <- which(Youto042 > Youto040 & Youto042 > Youto041 & Youto042 > Youto043)
g4 <- which(Youto043 > Youto040 & Youto043 > Youto041 & Youto043 > Youto042)
- 2. # g1～g4 の中身は、それぞれ主食が米、パン、麺、その他であるデータのインデックス番号
length(g1); length(g2); length(g3); length(g4);
#[1] 18535 [1] 12538 [1] 939 [1] 15
length(g1)+length(g2)+length(g3)+length(g4)
# [1] 32027 # データサイズと合致するので、同額（タイ）のデータは存在しない
# 主食分類符号（f.MD）作成
f.MD <- rep(NA, dim(dt1)[1])
f.MD[g1] <- 1
f.MD[g2] <- 2
f.MD[g3] <- 3
f.MD[g4] <- 4
table(f.MD, useNA="always")
# f.MD 1 2 3 4 <NA>
# 18535 12538 939 15 0
# ※ useNA オプションを指定しないと、NA があっても表示されない
# 比率を計算してみる
# 世帯数ベース
t1 <- matrix(NA, nc=4, nr=4)
colnames(t1) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与
rownames(t1) <- c("世帯数", " 世帯比率%", "金額(千円)", " 金額比率%")
# コンスタントファイルの漢字の項目名を付与
t1[1,] <- table(f.MD) # 世帯数ベース
t1[2,] <- t1[1,] / sum(t1[1,]) * 100 # 比率計算
t1[3,] <- apply(dt1[,54:57], 2, sum) /1000 # 金額ベース
t1[4,] <- t1[3,] / sum(t1[3,]) * 100 # 比率計算
round(t1, digits=2)
#...............................................................................................
# B. 円グラフ、帯グラフ、棒グラフ
#...............................................................................................
# 円グラフ作成(世帯数ベース)
main1 <- table(f.MD, useNA="ifany") # 比率計算は必要ない
names(main1) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与
#pie(main1, main="穀類支出の内訳(世帯数)", col=rainbow(4)) # 題目を付与 虹色指定
pie(main1, main="穀類支出の内訳(世帯数)", col=gray.colors(4)) # 題目を付与 白黒で
# 円グラフ作成(金額ベース)
#main2a <- c(sum(Youto040), sum(Youto041), sum(Youto042), sum(Youto043))
main2 <- apply(dt1[,54:57], 2, sum) # 上のコマンドでも同じ
names(main2) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与
- 3. #pie(main2, main="穀類支出の内訳(金額)", col=rainbow(4)) # 題目を付与 虹色指定
pie(main2, main="穀類支出の内訳(金額)", col=gray.colors(4)) # 題目を付与 白黒で
# 二つ横に並べてみよう
dev.new(width=16, height=10) # 新たなグラフィックウィンドウを作成
par(mfrow=c(1,2)) # 1 行 2 列にウィンドウを分割
pie(main1, main="穀類支出の内訳(世帯数)", col=rainbow(4)) # 題目を付与 虹色指定
pie(main2, main="穀類支出の内訳(金額)", col=rainbow(4)) # 題目を付与 虹色指定
#......................................................................
# 帯グラフ作成
#main1 と 2 の内容をパーセンテージにする
main1p <- main1 / sum(main1) * 100
main2p <- main2 / sum(main2) * 100
dev.new(width=16, height=7) # 新たなグラフィックウィンドウを作成
par(mfrow=c(1,2)) # 1 行 2 列にウィンドウを分割
barplot(as.matrix(main1p), horiz=T, main="穀類支出の内訳(世帯数)", col=gray.colors(4),
legend=names(main1),
args.legend = list(x = "topleft"), xlab="%")
barplot(as.matrix(main2p), horiz=T, main="穀類支出の内訳(金額)", col=gray.colors(4),
legend=names(main1),
args.legend = list(x = "topleft"), xlab="%")
#......................................................................
# 棒グラフ作成
dev.new(width=10, height=6) # 新たなグラフィックウィンドウを作成
main3 <- cbind(main1p, main2p)
colnames(main3) <- c("世帯数","金額")
barplot(main3, main="穀類支出の内訳", col=gray.colors(4), horiz=T, xlab="%")
#......................................................................
# C. プロット図の保存と加工（手作業の場合）
# ① グラフィックウインドウをクリックしてアクティブ化
# ② メニューバーから「ファイル」=> 「別名で保存」=> お好みのファイル形式を選択
# => ファイル名と置き場所を指定（ここでは png と emf を作成します）
# ③ 余白を削ったり、文字を加えたりなどの加工は、ペイントを使用する
# ※ 留意事項: メモリを食うのでグラフィックウィンドウは大量に開かない
#......................................................................
# D. 円グラフ二つを A4 横サイズの png 形式で保存（200 dpi）
png(filename="主食の割合_円グラフ.png", width=2339, height=1654, pointsize = 32)
- 4. par(mfrow=c(1,2)) # 1 行 2 列にウィンドウを分割
pie(main1, main="穀類支出の内訳(世帯数)", col=gray.colors(4))
pie(main2, main="穀類支出の内訳(金額)", col=gray.colors(4))
dev.off()
#...............................................................................................
# 詳細分析(1) エンゲル係数が大きいと米?
#...............................................................................................
# エンゲル係数は、支出総額(Youto037) に占める食料支出額（Youto038）の割合
#......................................................................
# E. 散布図作成 : 支出総額(Youto037) と食料支出額（Youto038）の関係
# E.1 基本の散布図
plot(Youto037[g1], Youto038[g1], main="主食別支出総額と食料費(実軸)", xlab="支出総額", ylab="食料費",
ylim=c(0, max(Youto038)))
# 確実に全データが収まるよう、y 軸は最大値までの幅を指定しているのがポイント!
points(Youto037[g2], Youto038[g2], col="green")
points(Youto037[g3], Youto038[g3], col="red")
points(Youto037[g4], Youto038[g4], col="yellow")
# データ量の多いものから順に描画すると良い
# E.2 常用対数軸で格子を追加
plot(log10(Youto037), log10(Youto038), main="主食別支出総額と食料費(常用対数軸)", xlab="支出総額",
ylab="食料費", ylim=c(log10(min(Youto038)), log10(max(Youto038))), col="white")
points(log10(Youto037[g1]), log10(Youto038[g1]))
points(log10(Youto037[g2]), log10(Youto038[g2]), col="green")
points(log10(Youto037[g3]), log10(Youto038[g3]), col="red")
points(log10(Youto037[g4]), log10(Youto038[g4]), col="yellow")
abline(h=seq(4, 5.5, by=0.5), v=seq(4.5, 7, by=0.5), lty=3, col="gray")
# E.3 透過プロットにする : 色は RGB 指定
# 黒:(0,0,0) 赤: (1,0,0) 緑:(0,1,0) 黄:(1,1,0) 第 4 パラメータは透過度
plot(log10(Youto037), log10(Youto038), main="主食別支出総額と食料費(常用対数軸)", xlab="支出総額",
ylab="食料費", ylim=c(log10(min(Youto038)), log10(max(Youto038))), col="white")
points(log10(Youto037[g1]), log10(Youto038[g1]), col=rgb(0,0,0,0.2), pch=20,)
points(log10(Youto037[g2]), log10(Youto038[g2]), col=rgb(0,1,0, 0.3), pch=20)
points(log10(Youto037[g3]), log10(Youto038[g3]), col=rgb(1,0,0, 0.2), pch=20)
points(log10(Youto037[g4]), log10(Youto038[g4]), col=rgb(1,1,0, 0.4), pch=20) # 黄は弱い色なので透過
度を下げる
abline(h=seq(4, 5.5, by=0.5), v=seq(4.5, 7, by=0.5), lty=3, col="gray")
#......................................................................
# F. 主食分類符号別エンゲル係数の基本統計量算出
# 支出総額(Youto037) に占める食料費（Youto038）の割合
- 5. eg1 <- Youto038 / Youto037 * 100
summary(eg1) # 基本統計量
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.9159 18.5500 23.9000 24.6000 29.8900 96.0800
# これを主食別に算出してくっつける
out1 <- cbind(summary(Youto038[g1]/Youto037[g1]), summary(Youto038[g2]/Youto037[g2]),
summary(Youto038[g3]/Youto037[g3]), summary(Youto038[g4]/Youto037[g4]),
summary(Youto038/Youto037)) * 100
colnames(out1) <- c("米", "パン", "麺", "その他", "合計")
round(out1, digits=3) # 小数点第 3 位で四捨五入
# 米 パン 麺 その他 合計
#Min. 1.413 0.916 3.571 4.196 0.916
#1st Qu. 18.810 18.250 17.840 15.620 18.550
#Median 24.370 23.350 22.830 17.770 23.900
#Mean 25.080 23.950 23.870 20.420 24.600
#3rd Qu. 30.640 28.930 29.330 25.900 29.890
#Max. 96.080 75.630 54.600 42.280 96.080
#......................................................................
# G. 複雑の棒プロット(主食分類符号別エンゲル係数の中央値と平均値
barplot(out1[2:3,], beside=TRUE) # 中央値と平均に着目
# y 軸改善, 凡例?目盛り追加
barplot(out1[2:3,], beside=TRUE, legend=c("中央値", "平均値"), ylim=c(0,25),
main="主食の違い別による中央値?平均値の差異") # 中央値と平均に着目
abline(h=0:25, col="gray", lty=3)
#.............................................................................................................
# 詳細分析(2) 世帯の人数が多いと米?
#.............................................................................................................
# 仮説: ご飯は一度に大量に炊けるので、世帯人員が多いと主食にする割合が多いのでは
#......................................................................
# 演習 H 世帯人員数を確認し、棒グラフを作成する
(tb1 <- table(SetaiJinin, useNA="ifany")) # useNA="ifany" がないと、NA が含まれていても表示されな
い。
# SetaiJinin
# 2 3 4 5 6 7 8 9 10
#7438 8537 9944 4405 1214 390 81 15 3
# 棒グラフ
barplot(tb1, ylim=c(0, 11000), main="世帯人員の数")
- 6. abline(h=seq(0, 11000, by=1000), col="gray", lty=3)
# 演習 I 世帯人員数と主食の別のクロス表を作成する
(tb2 <- table(SetaiJinin, f.MD))
# f.MD
#SetaiJinin 1 2 3
# 2 4605 2554 279
# 3 4814 3439 284
# 4 5152 4544 248
# 5 2703 1608 94
# 6 890 297 27
# 7 290 90 10
# 8 69 12 0
# 9 15 0 0
# 10 3 0 0
#......................................................................
# 演習 J. クロス表を棒グラフにして図示
barplot(t(tb2), beside=TRUE, legend=c("米", "パン", "麺", "その他"), # t()は転置
main="主食別世帯人員数")
abline(h=seq(0, 5000, by=500), col="gray", lty=3)
abline(h=seq(0, 5000, by=1000), col="gray", lty=1)
#l1 <- length(unique(SetaiJinin))
#leg1 <- apply(cbind(rep("世帯人員: ", l1), unique(SetaiJinin)) , 1, paste, collapse="")
#barplot(table(SetaiJinin, f.MD), beside=TRUE, legend=leg1)
#......................................................................
# 演習 K. クロス表を世帯人員の比率にする
# まず table(SetaiJinin, f.MD)に、合計欄を加える
addmargins(table(SetaiJinin, f.MD, exclude=NULL)) # exclude=NULL は NA もカウントする
# f.MD
#SetaiJinin 1 2 3 <NA> Sum
# 2 4605 2554 279 0 7438
# 3 4814 3439 284 0 8537
# 4 5152 4544 248 0 9944
# 5 2703 1608 94 0 4405
# 6 890 297 27 0 1214
# 7 290 90 10 0 390
# 8 69 12 0 0 81
# 9 15 0 0 0 15
# 10 3 0 0 0 3
# <NA> 0 0 0 0 0
- 7. # Sum 18541 12544 942 0 32027
tb2.p <- prop.table(table(SetaiJinin, f.MD), 2) # 縦（世帯人員別）の構成比
tb3.p <- prop.table(table(SetaiJinin, f.MD), 1) # 横（主食別）の構成比
round(tb2.p, 3) # 小数点 4 桁目で四捨五入
# f.MD 縦（世帯人員別）の構成比
#SetaiJinin 1 2 3
# 2 0.248 0.204 0.296
# 3 0.260 0.274 0.301
# 4 0.278 0.362 0.263
# 5 0.146 0.128 0.100
# 6 0.048 0.024 0.029
# 7 0.016 0.007 0.011
# 8 0.004 0.001 0.000
# 9 0.001 0.000 0.000
# 10 0.000 0.000 0.000
round(tb3.p, 3) # 横（主食別）の構成比
# f.MD
#SetaiJinin 1 2 3 4
# 2 0.619 0.343 0.038 0.000
# 3 0.564 0.402 0.033 0.001
# 4 0.518 0.457 0.025 0.000
# 5 0.614 0.365 0.021 0.000
# 6 0.733 0.245 0.022 0.000
# 7 0.744 0.231 0.026 0.000
# 8 0.852 0.148 0.000 0.000
# 9 1.000 0.000 0.000 0.000
# 10 1.000 0.000 0.000 0.000
# 凡例ラベル作成
l1 <- length(unique(SetaiJinin))
(leg1 <- apply(cbind(rep("世帯人員: ", l1), unique(SetaiJinin)) , 1, paste, collapse=""))
colnames(tb3.p) <- c("米", "パン", "麺", "その他") # x 軸のラベル用
barplot(tb3.p * 100, col=blues9, beside=T, legend=leg1, ylab="%",
main="主食毎の世帯人員数の構成割合")
#.............................................................................................................
# 詳細分析(3) 世帯収入が高いとパン?
#.............................................................................................................
# 収入総額の内訳
#Youto004 経常収入
#Youto018 特別収入
- 8. #Youto021 実収入以外の収入
# 演習 L. 散布図行列を眺めてみる
# まずプロットするデータ行列を作成
GSyunyu <- cbind(Youto004, Youto018, Youto021)
ix1 <- c(which(cd1$namae=="Youto004"), which(cd1$namae=="Youto018"),
which(cd1$namae=="Youto021"))
colnames(GSyunyu) <- cd1$koumoku[ix1]
pairs(GSyunyu, col=f.MD, pch=20, cex=0.2) # 実軸
# ファイル化したい場合は
png(filename="世帯収入と主食_散布図行列_実軸.png", width=2339, height=1654, pointsize = 32)
pairs(GSyunyu, col=f.MD, pch=20, cex=0.2) # 実軸
dev.off()
#------------------------------------------------------------------------------
# 演習 M. aov 関数を用いた一元配置分散分析
rs1 <- aov(f.MD ~ GSyunyu)
summary(rs1)
# Df Sum Sq Mean Sq F value Pr(>F)
#GSyunyu 3 29 9.636 31.54 <2e-16 *** <= 母平均は同じではない
#Residuals 32023 9783 0.306
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 有意だという結果は出ているが、この手法はデータが正規分布であることが前提
# データを対数化して、正規分布に近づける
#------------------------------------------------------------------------------
# 演習 N. 対数変換でデータを正規化する
png(filename="世帯収入と主食_散布図行列_常用対数.png", width=2339, height=1654, pointsize = 32)
pairs(log10(GSyunyu), col=f.MD, pch=20, cex=0.2) # 常用対数軸
dev.off()
rs2 <- aov(f.MD ~ log10(GSyunyu))
# データに欠測などがある、というエラーメッセージが出る。
# 0 値を対数化すると-Inf(マイナス無限大)、マイナス値を対数化すると NaN(非数)になり、
# 分析時に問題を起こすので、データを予めチェックし、元の数値が 0 以下のデータを除く
# チェック
length(which(is.na(Youto004)))
length(which(is.na(Youto018)))
- 10. plot(ct1) # 「根しかないよ」というエラーメッセージ。デフォルト設定では枝わかれしない
# 枝分かれを促進させるためには、cp の値を小さくする。デフォルト値は 0.01
ct2 <- rpart(f.MD ~., dat= as.data.frame(GSyunyu), method="class", control = rpart.control(cp = 0.001))
plot(ct2)
text(ct2)
# 分散分析の結果から、何らかの差異はあると思われるが、何か説明がつきそうにはない。
# 収入を一人あたり換算してみても、やはり説明は難しそう
ct3 <- rpart(f.MD ~., dat= as.data.frame(GSyunyu/SetaiJinin), method="class", control = rpart.control(cp
= 0.001))
plot(ct3)
text(ct3)
#.............................................................................................................
# III. 探索的な分析
#.............................................................................................................
# 演習 P. 全データで CART
ct4 <- rpart(f.MD ~., dat= dt1, method="class")
plot(ct4)
text(ct4)
# 説明変数に全データを放り込むと、米とパンの支出額での分類が提案される。
# ある意味当然だが、あまり意味がないので、説明変数から穀類の内訳を削る
#------------------------------------------------------------------------------
# 演習 Q. 説明変数を削って CART
cd1[53:57,] # dt1 の第 54～57 変数を除きたい
# ren koumoku lvl type namae
#53 205 穀類 6 1 Youto039 # 合計の穀類は残す
#54 208 米 7 1 Youto040
#55 210 パン 7 1 Youto041
#56 212 めん類 7 1 Youto042
#57 214 他の穀類 7 1 Youto043
dt2 <- dt1[,-(54:57)]
dim(dt1) #[1] 32027 197
dim(dt2) #[1] 32027 193 # 項目が 4 つ減っている
# dt2 で再び CART
ct5 <- rpart(f.MD ~., dat= dt2, method="class")
plot(ct5)
text(ct5)
# S1_Age は世帯主の年齢 5 歳階級
- 11. # Youto039 は穀類（主食）への支出
# Youto067 は菓子
#------------------------------------------------------------------------------
# 演習 R 世帯主の年齢階級
#--------------------------------------
# R.1 単純な表の作成
(t.age <- table(S1_Age))
#S1_Age
# 2 3 4 5 6 7 8 9 10 11 12 VV 1
# 887 2557 3558 3893 3938 4294 4188 1993 566 115 26 5944 68
# 年齢階級 1 が一番最後になっている!!!
#--------------------------------------
# R.2 因子水準の修正
# levels パラメータで因子水準の並び順を任意に指定する
table(factor(S1_Age, levels=c(" 1"," 2"," 3"," 4"," 5"," 6"," 7"," 8"," 9","10","11","12", "VV")) )
# 1 2 3 4 5 6 7 8 9 10 11 12 VV
# 68 887 2557 3558 3893 3938 4294 4188 1993 566 115 26 5944
# 因子の水準を直した S2_Age を作成 ※ 習慣として加工前のデータも残しておくと安全
S2_Age <-factor(S1_Age, levels=c(" 1"," 2"," 3"," 4"," 5"," 6"," 7"," 8"," 9","10","11","12", "VV"))
(t2.age <- table(S2_Age))
#--------------------------------------
# R.3 棒グラフ作成
barplot(t(t2.age), main="世帯主の年齢階級")
#--------------------------------------
# R.4 年齢階級と主食の別のクロス表作成
t3 <- table(S2_Age, f.MD)
t(t3)
# S2_Age
#f.MD 1 2 3 4 5 6 7 8 9 10 11 12 VV
# 1 23 275 796 1167 1850 2336 3052 3406 1656 480 105 21 3368
# 2 30 550 1641 2289 1946 1479 1159 730 313 80 7 5 2309
# 3 15 62 114 102 97 123 77 52 21 6 3 0 267
# 4 0 0 6 0 0 0 6 0 3 0 0 0 0
#--------------------------------------
# R.5 クロス表から複雑な棒グラフ
#barplot(table(S2_Age, f.MD)) # 棒は 3 つ。できれば積み上げでなく横に並べたい!
barplot(t(table(S2_Age, f.MD))) # 棒の数は世帯人員。できれば積み上げでなく横に並べたい!
barplot(t(table(S2_Age, f.MD)), beside=TRUE, legend=c("米", "パン", "麺", "その他"))
- 12. # 不詳 VV のバーが高く、凡例と重なるので locator を使って修正する
# カラーパレット作成
pcol <- c(blues9[1], blues9[5], blues9[7], blues9[9])
barplot(t(table(S2_Age, f.MD)), beside=TRUE, col=pcol, ylim=c(0, 3500))
abline(h=seq(0, 3500, by=500), col="gray", lty=3)
abline(h=seq(0, 3500, by=1000), col="gray", lty=1)
legend(locator(1), legend=c("米", "パン", "麺", "その他"), fill=pcol, bg="gray90")
# 主食が米の世帯とパンの世帯の年齢階級の分布が異なることが確認できた
q(save="no") # 作業スペースを保存せずに終了
#################################################
#################################################
#################################################
#################################################
![#-------------------------------------------------------
# R03_データ分析編_擬似ミクロ.r
#-------------------------------------------------------
rm(list=ls(all=TRUE)) # 作業領域のお掃除
#...............................................................................................
# i. 入出力編で作成した、擬似ミクロとコンスタントを収録したワークスペースファイルを読み込む
#setwd("Y:/PC123/giji/") #要修正%%%
setwd("D:/R 研修 2016/giji/") #要修正%%%
load("giji.RData")
ls() # ロードしたワークスペースファイルの内容表示
attach(dt1) # データフレーム dt1 を attach
#...............................................................................................
# ii. 主食分類符号の作成
# 主食に関するデータを確認。該当項目は、Youto039 から Youto043 で、データフレーム内では 53 から
# 57 番目にあることは、コンスタントデータで確認できる
cd1[53:57,]
# ren koumoku lvl type namae
#53 205 穀類 6 1 Youto039
#54 208 米 7 1 Youto040
#55 210 パン 7 1 Youto041
#56 212 めん類 7 1 Youto042
#57 214 他の穀類 7 1 Youto043
# 想定されない数字以外の文字が入っているか確認 => 全てなし
is.character(Youto039); is.character(Youto040); is.character(Youto041);
is.character(Youto042); is.character(Youto043);
# 総合計と内訳の和の差異は誤差の範囲
sum(Youto039) - sum(Youto040) - sum(Youto041) - sum(Youto042) - sum(Youto043)
# [1] -0.004349995
# NA は存在しない
length(which(is.na(dt1[,53:57])==TRUE)) #[1] 0
#...............................................................................................
# A. 主食分類符号の作成
# 主食インデックス（g1-g3）の作成
g1 <- which(Youto040 > Youto041 & Youto040 > Youto042 & Youto040 > Youto043)
g2 <- which(Youto041 > Youto040 & Youto041 > Youto042 & Youto041 > Youto043)
g3 <- which(Youto042 > Youto040 & Youto042 > Youto041 & Youto042 > Youto043)
g4 <- which(Youto043 > Youto040 & Youto043 > Youto041 & Youto043 > Youto042)](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-1-320.jpg)
![# g1～g4 の中身は、それぞれ主食が米、パン、麺、その他であるデータのインデックス番号
length(g1); length(g2); length(g3); length(g4);
#[1] 18535 [1] 12538 [1] 939 [1] 15
length(g1)+length(g2)+length(g3)+length(g4)
# [1] 32027 # データサイズと合致するので、同額（タイ）のデータは存在しない
# 主食分類符号（f.MD）作成
f.MD <- rep(NA, dim(dt1)[1])
f.MD[g1] <- 1
f.MD[g2] <- 2
f.MD[g3] <- 3
f.MD[g4] <- 4
table(f.MD, useNA="always")
# f.MD 1 2 3 4 <NA>
# 18535 12538 939 15 0
# ※ useNA オプションを指定しないと、NA があっても表示されない
# 比率を計算してみる
# 世帯数ベース
t1 <- matrix(NA, nc=4, nr=4)
colnames(t1) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与
rownames(t1) <- c("世帯数", " 世帯比率%", "金額(千円)", " 金額比率%")
# コンスタントファイルの漢字の項目名を付与
t1[1,] <- table(f.MD) # 世帯数ベース
t1[2,] <- t1[1,] / sum(t1[1,]) * 100 # 比率計算
t1[3,] <- apply(dt1[,54:57], 2, sum) /1000 # 金額ベース
t1[4,] <- t1[3,] / sum(t1[3,]) * 100 # 比率計算
round(t1, digits=2)
#...............................................................................................
# B. 円グラフ、帯グラフ、棒グラフ
#...............................................................................................
# 円グラフ作成(世帯数ベース)
main1 <- table(f.MD, useNA="ifany") # 比率計算は必要ない
names(main1) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与
#pie(main1, main="穀類支出の内訳(世帯数)", col=rainbow(4)) # 題目を付与 虹色指定
pie(main1, main="穀類支出の内訳(世帯数)", col=gray.colors(4)) # 題目を付与 白黒で
# 円グラフ作成(金額ベース)
#main2a <- c(sum(Youto040), sum(Youto041), sum(Youto042), sum(Youto043))
main2 <- apply(dt1[,54:57], 2, sum) # 上のコマンドでも同じ
names(main2) <- cd1$koumoku[54:57] # コンスタントファイルの漢字の項目名を付与](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-2-320.jpg)
![par(mfrow=c(1,2)) # 1 行 2 列にウィンドウを分割
pie(main1, main="穀類支出の内訳(世帯数)", col=gray.colors(4))
pie(main2, main="穀類支出の内訳(金額)", col=gray.colors(4))
dev.off()
#...............................................................................................
# 詳細分析(1) エンゲル係数が大きいと米?
#...............................................................................................
# エンゲル係数は、支出総額(Youto037) に占める食料支出額（Youto038）の割合
#......................................................................
# E. 散布図作成 : 支出総額(Youto037) と食料支出額（Youto038）の関係
# E.1 基本の散布図
plot(Youto037[g1], Youto038[g1], main="主食別支出総額と食料費(実軸)", xlab="支出総額", ylab="食料費",
ylim=c(0, max(Youto038)))
# 確実に全データが収まるよう、y 軸は最大値までの幅を指定しているのがポイント!
points(Youto037[g2], Youto038[g2], col="green")
points(Youto037[g3], Youto038[g3], col="red")
points(Youto037[g4], Youto038[g4], col="yellow")
# データ量の多いものから順に描画すると良い
# E.2 常用対数軸で格子を追加
plot(log10(Youto037), log10(Youto038), main="主食別支出総額と食料費(常用対数軸)", xlab="支出総額",
ylab="食料費", ylim=c(log10(min(Youto038)), log10(max(Youto038))), col="white")
points(log10(Youto037[g1]), log10(Youto038[g1]))
points(log10(Youto037[g2]), log10(Youto038[g2]), col="green")
points(log10(Youto037[g3]), log10(Youto038[g3]), col="red")
points(log10(Youto037[g4]), log10(Youto038[g4]), col="yellow")
abline(h=seq(4, 5.5, by=0.5), v=seq(4.5, 7, by=0.5), lty=3, col="gray")
# E.3 透過プロットにする : 色は RGB 指定
# 黒:(0,0,0) 赤: (1,0,0) 緑:(0,1,0) 黄:(1,1,0) 第 4 パラメータは透過度
plot(log10(Youto037), log10(Youto038), main="主食別支出総額と食料費(常用対数軸)", xlab="支出総額",
ylab="食料費", ylim=c(log10(min(Youto038)), log10(max(Youto038))), col="white")
points(log10(Youto037[g1]), log10(Youto038[g1]), col=rgb(0,0,0,0.2), pch=20,)
points(log10(Youto037[g2]), log10(Youto038[g2]), col=rgb(0,1,0, 0.3), pch=20)
points(log10(Youto037[g3]), log10(Youto038[g3]), col=rgb(1,0,0, 0.2), pch=20)
points(log10(Youto037[g4]), log10(Youto038[g4]), col=rgb(1,1,0, 0.4), pch=20) # 黄は弱い色なので透過
度を下げる
abline(h=seq(4, 5.5, by=0.5), v=seq(4.5, 7, by=0.5), lty=3, col="gray")
#......................................................................
# F. 主食分類符号別エンゲル係数の基本統計量算出
# 支出総額(Youto037) に占める食料費（Youto038）の割合](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-4-320.jpg)
![eg1 <- Youto038 / Youto037 * 100
summary(eg1) # 基本統計量
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.9159 18.5500 23.9000 24.6000 29.8900 96.0800
# これを主食別に算出してくっつける
out1 <- cbind(summary(Youto038[g1]/Youto037[g1]), summary(Youto038[g2]/Youto037[g2]),
summary(Youto038[g3]/Youto037[g3]), summary(Youto038[g4]/Youto037[g4]),
summary(Youto038/Youto037)) * 100
colnames(out1) <- c("米", "パン", "麺", "その他", "合計")
round(out1, digits=3) # 小数点第 3 位で四捨五入
# 米 パン 麺 その他 合計
#Min. 1.413 0.916 3.571 4.196 0.916
#1st Qu. 18.810 18.250 17.840 15.620 18.550
#Median 24.370 23.350 22.830 17.770 23.900
#Mean 25.080 23.950 23.870 20.420 24.600
#3rd Qu. 30.640 28.930 29.330 25.900 29.890
#Max. 96.080 75.630 54.600 42.280 96.080
#......................................................................
# G. 複雑の棒プロット(主食分類符号別エンゲル係数の中央値と平均値
barplot(out1[2:3,], beside=TRUE) # 中央値と平均に着目
# y 軸改善, 凡例?目盛り追加
barplot(out1[2:3,], beside=TRUE, legend=c("中央値", "平均値"), ylim=c(0,25),
main="主食の違い別による中央値?平均値の差異") # 中央値と平均に着目
abline(h=0:25, col="gray", lty=3)
#.............................................................................................................
# 詳細分析(2) 世帯の人数が多いと米?
#.............................................................................................................
# 仮説: ご飯は一度に大量に炊けるので、世帯人員が多いと主食にする割合が多いのでは
#......................................................................
# 演習 H 世帯人員数を確認し、棒グラフを作成する
(tb1 <- table(SetaiJinin, useNA="ifany")) # useNA="ifany" がないと、NA が含まれていても表示されな
い。
# SetaiJinin
# 2 3 4 5 6 7 8 9 10
#7438 8537 9944 4405 1214 390 81 15 3
# 棒グラフ
barplot(tb1, ylim=c(0, 11000), main="世帯人員の数")](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-5-320.jpg)
![#Youto021 実収入以外の収入
# 演習 L. 散布図行列を眺めてみる
# まずプロットするデータ行列を作成
GSyunyu <- cbind(Youto004, Youto018, Youto021)
ix1 <- c(which(cd1$namae=="Youto004"), which(cd1$namae=="Youto018"),
which(cd1$namae=="Youto021"))
colnames(GSyunyu) <- cd1$koumoku[ix1]
pairs(GSyunyu, col=f.MD, pch=20, cex=0.2) # 実軸
# ファイル化したい場合は
png(filename="世帯収入と主食_散布図行列_実軸.png", width=2339, height=1654, pointsize = 32)
pairs(GSyunyu, col=f.MD, pch=20, cex=0.2) # 実軸
dev.off()
#------------------------------------------------------------------------------
# 演習 M. aov 関数を用いた一元配置分散分析
rs1 <- aov(f.MD ~ GSyunyu)
summary(rs1)
# Df Sum Sq Mean Sq F value Pr(>F)
#GSyunyu 3 29 9.636 31.54 <2e-16 *** <= 母平均は同じではない
#Residuals 32023 9783 0.306
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 有意だという結果は出ているが、この手法はデータが正規分布であることが前提
# データを対数化して、正規分布に近づける
#------------------------------------------------------------------------------
# 演習 N. 対数変換でデータを正規化する
png(filename="世帯収入と主食_散布図行列_常用対数.png", width=2339, height=1654, pointsize = 32)
pairs(log10(GSyunyu), col=f.MD, pch=20, cex=0.2) # 常用対数軸
dev.off()
rs2 <- aov(f.MD ~ log10(GSyunyu))
# データに欠測などがある、というエラーメッセージが出る。
# 0 値を対数化すると-Inf(マイナス無限大)、マイナス値を対数化すると NaN(非数)になり、
# 分析時に問題を起こすので、データを予めチェックし、元の数値が 0 以下のデータを除く
# チェック
length(which(is.na(Youto004)))
length(which(is.na(Youto018)))](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-8-320.jpg)
![length(which(is.na(Youto021)))
length(which(is.infinite(Youto004)))
length(which(is.infinite(Youto018)))
length(which(is.infinite(Youto021)))
length(which(is.nan(Youto004)))
length(which(is.nan(Youto018)))
length(which(is.nan(Youto021)))
length(which(Youto004=="V"))
length(which(Youto018=="V"))
length(which(Youto021=="V"))
length(which(Youto004=="VV"))
length(which(Youto018=="VV"))
length(which(Youto021=="VV"))
length(which(Youto004<=0)) #[1] 17
length(which(Youto018<=0)) #[1] 951
length(which(Youto021<=0)) #[1] 168
# 0 以下ではないデータのインデックス f.G を作成
f.GS <- rep(1, dim(GSyunyu)[1])
f.GS[which(Youto004 <= 0)] <- 0; f.GS[which(Youto018 <= 0)] <- 0; f.GS[which(Youto021 <=
0)] <- 0
f.G <- which(f.GS == 1)
length(f.GS) #[1] 32027 総レコード数
length(which(f.GS == 1)) #[1] 30917 問題のないレコードの数
# aov 関数
rs2 <- aov(f.MD[f.G] ~ log10(GSyunyu[f.G, ]))
summary(rs2)
# Df Sum Sq Mean Sq F value Pr(>F)
#log10(GSyunyu[f.G, ]) 3 51 17.118 55.9 <2e-16 ***
#Residuals 30913 9467 0.306
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# やっぱり有意。ただしこの場合、対数化データを分析しているので、比較対象は平均値ではなく、
# 幾何平均であることに留意する。
#------------------------------------------------------------------------------
# 演習 O. CART
require(rpart) # rpart 関数 最初はインストールが必要
require(partykit) # 決定木の可視化に必要 最初はインストールが必要
ct1 <- rpart(f.MD ~., dat= as.data.frame(GSyunyu), method="class")
# method="class"で分類木を指定している](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-9-320.jpg)
![plot(ct1) # 「根しかないよ」というエラーメッセージ。デフォルト設定では枝わかれしない
# 枝分かれを促進させるためには、cp の値を小さくする。デフォルト値は 0.01
ct2 <- rpart(f.MD ~., dat= as.data.frame(GSyunyu), method="class", control = rpart.control(cp = 0.001))
plot(ct2)
text(ct2)
# 分散分析の結果から、何らかの差異はあると思われるが、何か説明がつきそうにはない。
# 収入を一人あたり換算してみても、やはり説明は難しそう
ct3 <- rpart(f.MD ~., dat= as.data.frame(GSyunyu/SetaiJinin), method="class", control = rpart.control(cp
= 0.001))
plot(ct3)
text(ct3)
#.............................................................................................................
# III. 探索的な分析
#.............................................................................................................
# 演習 P. 全データで CART
ct4 <- rpart(f.MD ~., dat= dt1, method="class")
plot(ct4)
text(ct4)
# 説明変数に全データを放り込むと、米とパンの支出額での分類が提案される。
# ある意味当然だが、あまり意味がないので、説明変数から穀類の内訳を削る
#------------------------------------------------------------------------------
# 演習 Q. 説明変数を削って CART
cd1[53:57,] # dt1 の第 54～57 変数を除きたい
# ren koumoku lvl type namae
#53 205 穀類 6 1 Youto039 # 合計の穀類は残す
#54 208 米 7 1 Youto040
#55 210 パン 7 1 Youto041
#56 212 めん類 7 1 Youto042
#57 214 他の穀類 7 1 Youto043
dt2 <- dt1[,-(54:57)]
dim(dt1) #[1] 32027 197
dim(dt2) #[1] 32027 193 # 項目が 4 つ減っている
# dt2 で再び CART
ct5 <- rpart(f.MD ~., dat= dt2, method="class")
plot(ct5)
text(ct5)
# S1_Age は世帯主の年齢 5 歳階級](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-10-320.jpg)
![# 不詳 VV のバーが高く、凡例と重なるので locator を使って修正する
# カラーパレット作成
pcol <- c(blues9[1], blues9[5], blues9[7], blues9[9])
barplot(t(table(S2_Age, f.MD)), beside=TRUE, col=pcol, ylim=c(0, 3500))
abline(h=seq(0, 3500, by=500), col="gray", lty=3)
abline(h=seq(0, 3500, by=1000), col="gray", lty=1)
legend(locator(1), legend=c("米", "パン", "麺", "その他"), fill=pcol, bg="gray90")
# 主食が米の世帯とパンの世帯の年齢階級の分布が異なることが確認できた
q(save="no") # 作業スペースを保存せずに終了
#################################################
#################################################
#################################################
#################################################](https://image.slidesharecdn.com/r03-161003143939/85/R-03_-2016-12-320.jpg)