次元の海で溺れる

Rとデータ解析と統計手法たちとわたし

ggplot2とクラスタリング:極めて個人的な疑問に基づいて都道府県データであそぶ。~k-means編~

こんばんは

巷で噂の「厚切りジェイソン」なる芸人さんが
「厚切りJSON」に脳内変換されたくらいにはJSONに汚染されてきてます。にちようび。

SappoRo.R

ところで

少し前に念願のSappoRo.Rに初参加してきました。
死ぬほど楽しかったです。ありがとうございました。
来年の開催時にLTやれるぐらいには成長してたいなあ・・・なんて

おかげさまで今まで全力で目を逸らしていたこの件と向き合っています。

http://rpubs.com/yamano357/85463rpubs.com

コーディングルール・・・

オブジェクトへの割り当てに「=」を使い続けていた私には耳の痛い話です。
「=」を使い続けようと誓った竹馬の友がいるんですが・・・まあ・・・背に腹は代えられないので今後従います

きょうの本題

コーディング規約準拠の練習もかねて、
前回記事と似たようなことをしようかなと。

・e-statの都道府県データで遊ぶ
クラスタリングしまくってみる
・ggplot2の練習を兼ねる

が基本です。

お題について考えてたら大学1年生の時の統計分析のレポートが出てきた。

その名も

「ネット普及率と大学進学率が睡眠時間に及ぼす影響」

???

当時18歳の私がレポートでやってたことを焼き直してみるとこんなん。

#データ
#sleep:都道府県別睡眠時間(分)・・・社会生活基本調査平成23年度
#net:都道府県別家庭ネット普及率・・・通信利用動向調査平成23年度
#school:都道府県別大学進学率(%)・・・日本の統計2013

> kako_data<-read.table("clipboard",header=TRUE)
> kako_data
       area sleep  net school
1    北海道   468 39.8   40.4
2    青森県   481 30.2   41.9
3    岩手県   472 30.6   41.2
4    宮城県   469 40.3   45.5
5    秋田県   482 33.5   44.5
6    山形県   477 40.4   46.3
7    福島県   475 38.2   42.3
8    茨城県   466 34.2   50.9
9    栃木県   465 38.8   54.3
10   群馬県   464 38.7   52.5
11   埼玉県   455 42.5   57.1
12   千葉県   454 39.4   54.8
13   東京都   457 42.1   65.5
14 神奈川県   451 49.8   60.8
15   新潟県   472 33.4   47.3
16   富山県   466 40.3   54.2
17   石川県   469 43.0   54.8
18   福井県   468 36.2   56.0
19   山梨県   465 39.3   57.9
20   長野県   471 36.3   49.6
21   岐阜県   463 41.5   57.1
22   静岡県   462 41.7   54.1
23   愛知県   462 42.1   58.9
24   三重県   464 43.7   52.2
25   滋賀県   462 45.2   58.0
26   京都府   464 42.5   66.1
27   大阪府   460 40.7   58.7
28   兵庫県   454 42.1   59.9
29   奈良県   453 42.8   58.1
30 和歌山県   470 35.0   49.3
31   鳥取県   468 35.4   43.9
32   島根県   472 29.6   47.2
33   岡山県   462 40.5   52.9
34   広島県   456 32.8   61.1
35   山口県   467 27.0   43.2
36   徳島県   466 38.2   53.0
37   香川県   465 34.7   51.1
38   愛媛県   465 25.2   52.7
39   高知県   480 29.3   45.4
40   福岡県   466 32.6   53.3
41   佐賀県   467 38.0   42.3
42   長崎県   464 24.0   42.2
43   熊本県   471 27.5   43.1
44   大分県   465 27.8   47.4
45   宮崎県   467 21.1   43.0
46 鹿児島県   471 26.4   42.0
47   沖縄県   462 28.6   36.7

> model<-lm(sleep~net+school,data=kako_data)
> summary(model)

Call:
lm(formula = sleep ~ net + school, data = kako_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.8256  -3.5658   0.0187   3.3553  12.1481 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 498.48390    5.78069  86.233  < 2e-16 ***
net          -0.04622    0.16341  -0.283 0.778594    
school       -0.60862    0.14572  -4.177 0.000138 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.389 on 44 degrees of freedom
Multiple R-squared:  0.4296,	Adjusted R-squared:  0.4037 
F-statistic: 16.57 on 2 and 44 DF,  p-value: 4.328e-06

シンプルでなんか愚かな感じの重回帰ですね。
マルチコの香りが・・・するよ・・・

とりあえず
「なんか睡眠時間少ないところって、ネット普及してて大学進学率が高いところと一致するんじゃない・・・?」
みたいな18歳の着想から始まったんでしょう、このレポートは。

ということで

お題

睡眠時間とネット普及率と大学進学率をそれぞれクラスタリングしてみて、
地域性についてなんとなく知見を得よう。

ってことで、過去の自分を全力でイジ・・・、18歳の私の無念を晴らそうかと。

早速やってみる

#データいろいろ見てみる

library(ggplot2)

kako_data2<-kako_data[1:47,-3]#睡眠×進学率
kako_data3<-kako_data[1:47,-2]#進学率×ネット
kako_data4<-kako_data[1:47,-4]#睡眠×ネット
#睡眠時間と進学率の関係
sp2<-ggplot(kako_data2,aes(x=school,y=sleep))+
  geom_point(colour="orange",size=4)
sp2+geom_text(aes(label=area),size=3,vjust=-1)

f:id:WAFkw:20150628194218p:plain

#進学率とネット普及率の関係
sp3<-ggplot(kako_data3,aes(x=net,y=school))+
  geom_point(colour="blue",size=4)
sp3+geom_text(aes(label=area),size=3,vjust=-1)

f:id:WAFkw:20150628194254p:plain


ほら、、、
やっぱマルチコっぽいじゃん、、

#睡眠時間×ネット普及率の関係
sp4<-ggplot(kako_data4,aes(x=net,y=sleep))+
  geom_point(colour="red",size=4)
sp4+geom_text(aes(label=area),size=3,vjust=-1)

f:id:WAFkw:20150628194330p:plain


ここは意外とバラけてますね。
ネット流行りまくり地域は夜寝ない、みたいな単純な仮説は
どうやら微妙な仮説だったようで、、、

#K-means法でクラスタリングしてみる

library(kernlab)
x2<-as.matrix(kako_data2[,2:3])
x3<-as.matrix(kako_data3[,2:3])
x4<-as.matrix(kako_data4[,2:3])


#クラスタ数3,初期値10で探索開始。
cl2<-kmeans(x2,3,nstart=10)
cl3<-kmeans(x3,3,nstart=10)
cl4<-kmeans(x4,3,nstart=10)

> #評価値の算出
> sum(cl2$withinss)
[1] 1184.379
> sum(cl3$withinss)
[1] 1085.985
> sum(cl4$withinss)
[1] 1466.956

> print(cl2)
K-means clustering with 3 clusters of sizes 9, 20, 18

Cluster means:
     sleep   school
1 456.0000 60.23333
2 471.3000 43.86500
3 464.7222 53.96111

Clustering vector:
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 
 2  2  2  2  2  2  2  3  3  3  1  1  1  1  2  3  3  3  3  2  3  3  3  3  3  1  1  1  1  2  2  2  3  1  2  3  3  3  2  3  2  2 
43 44 45 46 47 
 2  3  2  2  2 

Within cluster sum of squares by cluster:
[1] 233.5800 740.1455 210.6539
 (between_SS / total_SS =  74.3 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"        
[9] "ifault" 

####長くなりそうだったので割愛####
print(cl3)
print(cl4)
##################################

#plotします

#睡眠時間×進学率
cl22=data.frame(cl2$cluster)
cl23=cbind(kako_data$area,kako_data$sleep,kako_data$school,cl22)#ダサい

cl_plot2=ggplot(cl23,aes(x=kako_data$school,y=kako_data$sleep,colour=cl2.cluster))+
  geom_point(size=4)
cl_plot2+geom_text(aes(label=kako_data$area),colour="grey40",size=3,vjust=-1)

#進学率とネット普及率
cl32=data.frame(cl3$cluster)
cl33=cbind(kako_data$area,kako_data$school,kako_data$net,cl32)

cl_plot3=ggplot(cl33,aes(x=kako_data$net,y=kako_data$school,colour=cl3.cluster))+
  geom_point(size=4)
cl_plot3+geom_text(aes(label=kako_data$area),colour="grey40",size=3,vjust=-1)

#睡眠時間×ネット普及率
cl42=data.frame(cl4$cluster)
cl43=cbind(kako_data$area,kako_data$sleep,kako_data$net,cl42)

cl_plot4=ggplot(cl43,aes(x=kako_data$net,y=kako_data$sleep,colour=cl4.cluster))+
  geom_point(size=4)
cl_plot4+geom_text(aes(label=kako_data$area),colour="grey40",size=3,vjust=-1)

f:id:WAFkw:20150628194746p:plain

f:id:WAFkw:20150628194806p:plain

f:id:WAFkw:20150628194826p:plain


うーーーんなにかうまいこと言えそうで、、言えない

おまけとして
睡眠時間に関するクラスタの共通項とかないかなーと思って探ってみました

#subset

> kekka_sleep<-subset(all.cluster[,-3],cl2.cluster==cl4.cluster)
> order(kekka_sleep$cl2.cluster) %>%
+ kekka_sleep[.,]
   kako_data$area cl2.cluster cl4.cluster
2          青森県           2           2
3          岩手県           2           2
5          秋田県           2           2
15         新潟県           2           2
32         島根県           2           2
35         山口県           2           2
39         高知県           2           2
42         長崎県           2           2
43         熊本県           2           2
45         宮崎県           2           2
46       鹿児島県           2           2
47         沖縄県           2           2
21         岐阜県           3           3
22         静岡県           3           3
23         愛知県           3           3
24         三重県           3           3
25         滋賀県           3           3
33         岡山県           3           3

睡眠時間×進学率と
睡眠時間×ネットの
それぞれのクラスタをまとめてます。

みごとに首都圏が抜けてますね。。。

なんか言いたい・・・言えない・・・

ということでなんか遠回りしたくせにわりとふつーーのオチで
なんとも言い難い。

ごめんよ18歳のわたし・・・

追記

恐縮ながらお声掛け頂いたこともあり
今回使用したrawデータについて詳細を追記します。

(今回昔の加工データを流用するという横着をしてしまったため
前処理Script存在せずです・・・すみません・・・)


①睡眠時間
都道府県別睡眠時間(分)・・・社会生活基本調査平成23年度【表番号1-1 週全体】
統計表一覧 政府統計の総合窓口 GL08020103

※.xls

②ネット普及率(無線)
都道府県別家庭ネット普及率・・・通信利用動向調査平成23年度【世17 問1(2)】
http://www.e-stat.go.jp/SG1/estat/List.do?bid=000001049749&cycode=0

※.csv
※最新版が存在します
※今考えると【世22 問2(1)】の方がネットの普及率という点では適切なデータかもしれません・・・


③大学進学率
都道府県別大学進学率(%)・・・日本の統計2013【22-18】
統計局ホームページ/日本の統計 2013−第22章 教育

※.xls
※最新版が存在します

ggplot2実践編:都道府県別身長体重を使ってクラスタリング&結果の可視化

明日からお仕事~
来週はいよいよSapporo.R~

ggplot2とだんだん仲良くなってきた気がしなくもないので、
今週もしつこくggplot2をやります。

慣れたら超楽しいね、ggplot2!!!!
まだまだポンコツだけど!!!

今回のお題

【仮  説 】身長体重って、なんか地域によって特性ありそうじゃない?

【使用データ】総務省統計局/学校保健調査平成26年「都道府県別/身長・体重・座高」17歳男女

【分析の流れ】
 1.まずふつーーーに図示してみる。データ全体の把握。
 2.k-平均法でクラスタリングしてみる
 3.結果をうまいことggplot2に反映出来たらいいな

【動機】完全に思いつきです。

困ったときの教科書

いつもの2冊です。

パターン認識 (Rで学ぶデータサイエンス 5)

パターン認識 (Rで学ぶデータサイエンス 5)

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

でも今回全然開かなかった、成長。

早速やります

#統計局e-statより、学校保健調査平成26年度の調査表をお借り。
#17歳男女の身長・体重平均値です

> data
         NAME m_height m_weight w_height w_weight
1  北 海 道    170.9     63.3    157.9     53.4
2  青   森    171.1     63.8    158.5     54.6
3  岩   手    170.3     64.4    157.6     53.4
4  宮   城    170.9     64.9    158.2     53.7
5  秋   田    171.4     65.7    158.3     53.5
6  山   形    170.8     65.0    157.7     53.5
7  福   島    170.3     62.5    157.9     54.1
8  茨   城    170.2     63.0    157.8     53.0
9  栃   木    170.4     63.0    157.3     54.5
10 群   馬    170.8     63.0    158.0     53.5
11 埼   玉    171.0     62.0    158.2     52.9
12 千   葉    170.9     62.6    157.8     53.5
13 東   京    170.9     62.2    158.6     52.3
14 神 奈 川    171.1     62.4    158.0     52.3
15 新   潟    170.8     62.6    158.1     53.4
16 富   山    171.4     63.5    158.1     53.1
17 石   川    171.7     63.5    158.4     52.7
18 福   井    171.3     63.2    158.3     53.2
19 山   梨    170.6     63.3    158.0     52.6
20 長   野    170.5     62.5    157.5     52.5
21 岐   阜    171.2     62.4    157.7     52.2
22 静   岡    170.4     61.7    158.0     52.5
23 愛   知    170.7     62.3    157.9     52.7
24 三   重    170.7     63.1    157.6     52.4
25 滋   賀    171.1     62.9    158.1     52.6
26 京   都    171.0     62.3    157.9     52.7
27 大   阪    171.2     62.5    158.5     53.2
28 兵   庫    170.4     62.5    158.0     53.0
29 奈   良    170.9     62.5    157.7     51.9
30 和 歌 山    170.7     62.9    157.9     53.0
31 鳥   取    171.1     62.1    158.4     53.1
32 島   根    169.8     61.9    157.5     53.2
33 岡   山    170.3     62.5    157.1     52.7
34 広   島    169.8     62.6    157.6     52.7
35 山   口    170.9     62.2    157.1     52.4
36 徳   島    170.0     63.4    156.6     53.2
37 香   川    170.1     62.2    157.1     52.1
38 愛   媛    170.2     62.7    157.0     53.0
39 高   知    169.8     62.9    157.4     52.7
40 福   岡    169.3     61.3    157.7     52.8
41 佐   賀    170.1     62.5    157.5     52.5
42 長   崎    170.7     63.5    157.1     53.0
43 熊   本    170.7     63.3    157.5     53.3
44 大   分    170.0     62.7    157.2     53.4
45 宮   崎    170.0     63.2    157.1     52.7
46 鹿 児 島    170.0     62.2    157.5     53.3
47 沖   縄    168.8     61.1    155.8     51.3


#m_height,m_weight:男性平均身長・体重
#w_height,w_weight:女性平均身長・体重

#データはrawのままいじってません

library(ggplot2)

#いったん描いてみよう
ggplot(data,aes(x=m_weight,y=m_height))+
  geom_point()

#男性の身長をy,体重をxと置いてるだけ

f:id:WAFkw:20150607214420p:plain

うん、、
ラベル付けたい、、、

どうやらgeom_text()を使って自動的にラベル付けが出来そう。

> #プロットの色とかサイズとかを整えてオブジェクト化
> sp=ggplot(data,aes(x=m_weight,y=m_height))+
+   geom_point(colour="orange",size=4)

> sp+geom_text(aes(label=NAME),size=3,vjust=-1)

f:id:WAFkw:20150607214902p:plain

ちゃんと都道府県の名前が入りました。
ちなみに、このvjustという引数で若干ラベルの位置を調整出来ます。

色々試した結果、正の方向に動かすとラベルが下方向に
負の方向に増やすとラベルが上方向に動いてくれます。

今回はplotとラベルがかぶってたのでvjust=-1でちょい上にズラしてみました。


ここでplotを見てみると、

東北→体重が重め、ガタイよさげ
北陸・関東→身長高め、体重軽め~標準
九州・沖縄→小柄め

って感じに、なんとなく3グループに分けられそう。

ほんとかどうか確かめてみたい・・・

#なんか2、3個のグループに分けられそう、、、

#K-means法でクラスタリングしてみる

library(kernlab)
x=as.matrix(data[,2:3])
#(なんかNAMEの文字列処理に自信がなかったので切り出し。意味はないです)

#クラスタ数3,初期値10で探索開始。
cl=kmeans(x,3,nstart=10)


#評価値を算出。なかなか小さい値なのでいい感じ?
> sum(cl$withinss)
[1] 19.00378

#クラスタリング結果
> print(cl)
K-means clustering with 3 clusters of sizes 25, 18, 4

Cluster means:
  m_height m_weight
1  170.952 62.81600
2  170.000 62.41111
3  170.850 65.00000

Clustering vector:
 [1] 1 1 3 3 3 3 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 1 1 2 2 2 1 2 2 2 2 2 2 1 1 2 2 2 2

Within cluster sum of squares by cluster:
[1] 8.396000 9.137778 1.470000
 (between_SS / total_SS =  62.6 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"   
[7] "size"         "iter"         "ifault"      


> #ggplot2でのクラスタリング結果の反映を考える

> #クラスタを切り出してみた
> cl1=data.frame(cl$cluster)
> cl1
   cl.cluster
1           1
2           1
3           3
4           3
5           3
6           3
7           2
8           2
9           1
10          1
11          1
12          1
13          1
14          1
15          1
16          1
17          1
18          1
19          1
20          2
21          1
22          2
23          1
24          1
25          1
26          1
27          1
28          2
29          1
30          1
31          1
32          2
33          2
34          2
35          1
36          2
37          2
38          2
39          2
40          2
41          2
42          1
43          1
44          2
45          2
46          2
47          2

#うまいことplotに反映させたいよねー

#元データとくっつけてみた
cl2=cbind(data$NAME,data$m_height,data$m_weight,cl1)#ダサい

#みてみる
>cl2
    data$NAME data$m_height data$m_weight cl.cluster
1  北 海 道         170.9          63.3          1
2  青   森         171.1          63.8          1
3  岩   手         170.3          64.4          3
4  宮   城         170.9          64.9          3
5  秋   田         171.4          65.7          3
6  山   形         170.8          65.0          3
7  福   島         170.3          62.5          2
8  茨   城         170.2          63.0          2
9  栃   木         170.4          63.0          1
10 群   馬         170.8          63.0          1
11 埼   玉         171.0          62.0          1
12 千   葉         170.9          62.6          1
13 東   京         170.9          62.2          1
14 神 奈 川         171.1          62.4          1
15 新   潟         170.8          62.6          1
16 富   山         171.4          63.5          1
17 石   川         171.7          63.5          1
18 福   井         171.3          63.2          1
19 山   梨         170.6          63.3          1
20 長   野         170.5          62.5          2
21 岐   阜         171.2          62.4          1
22 静   岡         170.4          61.7          2
23 愛   知         170.7          62.3          1
24 三   重         170.7          63.1          1
25 滋   賀         171.1          62.9          1
26 京   都         171.0          62.3          1
27 大   阪         171.2          62.5          1
28 兵   庫         170.4          62.5          2
29 奈   良         170.9          62.5          1
30 和 歌 山         170.7          62.9          1
31 鳥   取         171.1          62.1          1
32 島   根         169.8          61.9          2
33 岡   山         170.3          62.5          2
34 広   島         169.8          62.6          2
35 山   口         170.9          62.2          1
36 徳   島         170.0          63.4          2
37 香   川         170.1          62.2          2
38 愛   媛         170.2          62.7          2
39 高   知         169.8          62.9          2
40 福   岡         169.3          61.3          2
41 佐   賀         170.1          62.5          2
42 長   崎         170.7          63.5          1
43 熊   本         170.7          63.3          1
44 大   分         170.0          62.7          2
45 宮   崎         170.0          63.2          2
46 鹿 児 島         170.0          62.2          2
47 沖   縄         168.8          61.1          2


#クラスタリング結果で先ほどのplotを色分けしてみる

#やり方はほぼさっきと一緒。カラーの引数をクラスタにしてみた

cl3=ggplot(cl2,aes(x=data$m_weight,y=data$m_height,colour=cl.cluster))+
  geom_point(size=4)
cl3+geom_text(aes(label=data$NAME),colour="grey40",size=3,vjust=-1)


f:id:WAFkw:20150607220037p:plain


できたーーーー

ちょっと色味がわかりにくいけど、
さっきなんとなく仮説立てた通り

東北→体重が重め、ガタイよさげ
北陸・関東→身長高め、体重軽め~標準
九州・沖縄→小柄め

って感じの3群になった・・・

k-means法の理論的背景とか、その辺も書きたかったし、
そもそもクラスタ結果は色分けじゃなくてshapeで△と〇とかに分けたかったけど

次回にしよう、、、(なぜかshape失敗したし・・・)

以上です。たのしかった。

ggplot2で信頼区間の描画をして遊びたい 練習編

面白いデータセットを探し続けているものの
なかなか見つからないので歌詞でも収集して形態素に切ってしまおうプロジェクトが
水面下で進行中です。

ネットワークでも描こうか、、とも思ったけど卒論で死ぬ思いした残像がまだ、、、


###

SappoRo.R前にggplot2で遊ぼうシリーズ第2弾、
信頼図を描こう。です。

テーマが地味。

テレビとかに出てくるお医者さんやら何かしらの専門家の人が、
「信頼区間95%」をうまーく表現丸めて喋ってるのを見ると尊敬します。

『○○にならない可能性が5%、つまりかなり高い確率で〇〇になるってことですね~』

絶妙なライン攻めてるかんじ。


ggplot2初心者なので、
相変わらず教科書はこれです。

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

いつもお世話になっております。

ではさっそく。

#---信頼区間を図示してみたい--------------------

#使用するデータセットについて
#gcookbookパッケージよりclimateデータセットをお借りします

install.packages("gcookbook")
library(gcookbook)
library(ggplot2)


> #データ構造

> head(climate)
    Source Year Anomaly1y Anomaly5y Anomaly10y Unc10y
1 Berkeley 1800        NA        NA     -0.435  0.505
2 Berkeley 1801        NA        NA     -0.453  0.493
3 Berkeley 1802        NA        NA     -0.460  0.486
4 Berkeley 1803        NA        NA     -0.493  0.489
5 Berkeley 1804        NA        NA     -0.536  0.483
6 Berkeley 1805        NA        NA     -0.541  0.475

> str(climate)
'data.frame':	499 obs. of  6 variables:
 $ Source    : chr  "Berkeley" "Berkeley" "Berkeley" "Berkeley" ...
 $ Year      : num  1800 1801 1802 1803 1804 ...
 $ Anomaly1y : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Anomaly5y : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Anomaly10y: num  -0.435 -0.453 -0.46 -0.493 -0.536 -0.541 -0.59 -0.695 -0.763 -0.818 ...
 $ Unc10y    : num  0.505 0.493 0.486 0.489 0.483 0.475 0.468 0.461 0.453 0.451 ...
> summary(climate)
    Source               Year        Anomaly1y          Anomaly5y         Anomaly10y      
 Length:499         Min.   :1800   Min.   :-0.60070   Min.   :-0.4995   Min.   :-1.01500  
 Class :character   1st Qu.:1884   1st Qu.:-0.21629   1st Qu.:-0.1053   1st Qu.:-0.28350  
 Mode  :character   Median :1926   Median :-0.02797   Median :-0.0042   Median :-0.07328  
                    Mean   :1923   Mean   : 0.01277   Mean   : 0.0555   Mean   :-0.07869  
                    3rd Qu.:1968   3rd Qu.: 0.15155   3rd Qu.: 0.1620   3rd Qu.: 0.05065  
                    Max.   :2011   Max.   : 0.96354   Max.   : 0.8953   Max.   : 0.88400  
                                   NA's   :207        NA's   :373       NA's   :20        
     Unc10y      
 Min.   :0.0110  
 1st Qu.:0.0430  
 Median :0.1040  
 Mean   :0.1452  
 3rd Qu.:0.2220  
 Max.   :0.5050  
 NA's   :294     

#変数は6つ。
#Anomaly10yは1950年~1980年までの平均気温からの偏差の10年移動平均
#Unc10yは95%信頼区間
#Sourceは出典というか、調査機関って感じですかね。NASAとか。
#Anomaly1yと5yは1年移動平均と5年移動平均ぽいですがSourceによって欠損値がバラバラなので今回は不使用

#使うデータ部分を切り出してきて格納
#Berkely調査の平均気温偏差10年移動平均、その信頼区間を使用

clim=subset(climate,climate$Source=="Berkeley",select=c("Year","Anomaly10y","Unc10y"))

> clim
    Year Anomaly10y Unc10y
1   1800     -0.435  0.505
2   1801     -0.453  0.493
3   1802     -0.460  0.486
4   1803     -0.493  0.489
5   1804     -0.536  0.483
6   1805     -0.541  0.475
7   1806     -0.590  0.468
8   1807     -0.695  0.461
9   1808     -0.763  0.453
10  1809     -0.818  0.451
11  1810     -0.842  0.437
12  1811     -0.912  0.406
13  1812     -0.963  0.395
14  1813     -1.015  0.389
15  1814     -0.949  0.369
16  1815     -0.935  0.368
17  1816     -0.918  0.361
18  1817     -0.802  0.330
19  1818     -0.743  0.349
20  1819     -0.703  0.382
21  1820     -0.625  0.395
22  1821     -0.527  0.363
23  1822     -0.392  0.339
24  1823     -0.332  0.351
25  1824     -0.321  0.343
26  1825     -0.278  0.344
27  1826     -0.270  0.295
28  1827     -0.378  0.284
29  1828     -0.398  0.254
30  1829     -0.388  0.243
31  1830     -0.429  0.247
32  1831     -0.506  0.247
33  1832     -0.625  0.257
34  1833     -0.689  0.246
35  1834     -0.667  0.224
36  1835     -0.689  0.220
37  1836     -0.727  0.221
38  1837     -0.662  0.218
39  1838     -0.612  0.222
40  1839     -0.613  0.256
41  1840     -0.638  0.274
42  1841     -0.542  0.267
43  1842     -0.472  0.282
44  1843     -0.407  0.269
45  1844     -0.402  0.285
46  1845     -0.390  0.285
47  1846     -0.351  0.284
48  1847     -0.346  0.280
49  1848     -0.362  0.273
50  1849     -0.357  0.260
51  1850     -0.290  0.244
52  1851     -0.321  0.239
53  1852     -0.372  0.228
54  1853     -0.363  0.226
55  1854     -0.351  0.228
56  1855     -0.348  0.216
57  1856     -0.374  0.203
58  1857     -0.418  0.187
59  1858     -0.420  0.186
60  1859     -0.427  0.178
61  1860     -0.473  0.172
62  1861     -0.458  0.154
63  1862     -0.412  0.149
64  1863     -0.393  0.145
65  1864     -0.373  0.141
66  1865     -0.344  0.148
67  1866     -0.317  0.147
68  1867     -0.282  0.145
69  1868     -0.269  0.145
70  1869     -0.241  0.149
71  1870     -0.233  0.149
72  1871     -0.255  0.146
73  1872     -0.271  0.142
74  1873     -0.247  0.140
75  1874     -0.248  0.132
76  1875     -0.271  0.128
77  1876     -0.252  0.124
78  1877     -0.235  0.120
79  1878     -0.257  0.115
80  1879     -0.297  0.119
81  1880     -0.321  0.124
82  1881     -0.314  0.123
83  1882     -0.325  0.117
84  1883     -0.376  0.108
85  1884     -0.373  0.111
86  1885     -0.368  0.110
87  1886     -0.386  0.109
88  1887     -0.397  0.107
89  1888     -0.394  0.108
90  1889     -0.364  0.104
91  1890     -0.342  0.102
92  1891     -0.323  0.099
93  1892     -0.298  0.100
94  1893     -0.273  0.107
95  1894     -0.286  0.107
96  1895     -0.251  0.106
97  1896     -0.201  0.104
98  1897     -0.176  0.104
99  1898     -0.153  0.104
100 1899     -0.174  0.106
101 1900     -0.171  0.108
102 1901     -0.162  0.109
103 1902     -0.177  0.108
104 1903     -0.199  0.104
105 1904     -0.223  0.105
106 1905     -0.241  0.107
107 1906     -0.294  0.106
108 1907     -0.312  0.105
109 1908     -0.328  0.103
110 1909     -0.281  0.101
111 1910     -0.247  0.099
112 1911     -0.243  0.097
113 1912     -0.257  0.100
114 1913     -0.268  0.100
115 1914     -0.257  0.097
116 1915     -0.249  0.095
117 1916     -0.214  0.096
118 1917     -0.201  0.096
119 1918     -0.176  0.096
120 1919     -0.182  0.097
121 1920     -0.193  0.097
122 1921     -0.167  0.098
123 1922     -0.128  0.096
124 1923     -0.075  0.097
125 1924     -0.064  0.098
126 1925     -0.065  0.100
127 1926     -0.050  0.100
128 1927     -0.020  0.099
129 1928     -0.018  0.099
130 1929     -0.026  0.100
131 1930     -0.014  0.101
132 1931     -0.047  0.098
133 1932     -0.035  0.096
134 1933     -0.017  0.093
135 1934      0.020  0.092
136 1935      0.053  0.089
137 1936      0.063  0.085
138 1937      0.048  0.081
139 1938      0.073  0.079
140 1939      0.113  0.076
141 1940      0.113  0.072
142 1941      0.134  0.071
143 1942      0.134  0.069
144 1943      0.127  0.070
145 1944      0.111  0.068
146 1945      0.072  0.066
147 1946      0.035  0.066
148 1947      0.042  0.064
149 1948      0.045  0.063
150 1949      0.013  0.062
151 1950      0.010  0.058
152 1951     -0.017  0.054
153 1952     -0.040  0.047
154 1953     -0.040  0.043
155 1954     -0.032  0.038
156 1955     -0.022  0.035
157 1956      0.012  0.031
158 1957      0.007  0.028
159 1958      0.002  0.027
160 1959      0.002  0.026
161 1960     -0.019  0.026
162 1961     -0.001  0.021
163 1962      0.017  0.018
164 1963      0.004  0.016
165 1964     -0.028  0.018
166 1965     -0.006  0.017
167 1966     -0.024  0.017
168 1967     -0.041  0.019
169 1968     -0.025  0.020
170 1969     -0.019  0.024
171 1970      0.010  0.026
172 1971      0.007  0.022
173 1972      0.015  0.015
174 1973      0.028  0.012
175 1974      0.049  0.014
176 1975      0.068  0.012
177 1976      0.128  0.011
178 1977      0.158  0.012
179 1978      0.167  0.013
180 1979      0.193  0.012
181 1980      0.186  0.016
182 1981      0.217  0.016
183 1982      0.235  0.014
184 1983      0.270  0.014
185 1984      0.318  0.014
186 1985      0.344  0.013
187 1986      0.352  0.012
188 1987      0.380  0.011
189 1988      0.370  0.013
190 1989      0.366  0.017
191 1990      0.433  0.019
192 1991      0.467  0.018
193 1992      0.496  0.017
194 1993      0.526  0.019
195 1994      0.554  0.020
196 1995      0.563  0.019
197 1996      0.565  0.022
198 1997      0.618  0.022
199 1998      0.680  0.023
200 1999      0.734  0.025
201 2000      0.748  0.026
202 2001      0.793  0.027
203 2002      0.856  0.028
204 2003      0.869  0.028
205 2004      0.884  0.029

これで準備はOKなはず。
コメントにあれこれ書いてますが、
Climateデータセットの一部を切り出してきて使用したってことです。

#いざ描画

#信頼区間は
#Anomaly10yにUnc10yを加算したものを上限(ymax)
#Anomaly10yからUnc10yを減算したものを下限(ymin)として設定
#網掛けの透過度は80%(alpha=0.2)

ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_ribbon(aes(ymin=Anomaly10y-Unc10y,ymax=Anomaly10y+Unc10y),alpha=0.2)+
  geom_line()

f:id:WAFkw:20150523161942p:plain

できたー
イメージどおりー

geom_ribbon(網掛け)の上に、geom_lineを重ねてるので
こういう順番でのコード記述になる。ふーん。

てことは

順番を逆にしたら網掛けが上に来てしまって見えなくなる・・・?

ためしてみた

#あえて網掛け部分と実線を逆に描画してみる
#しかも網掛けの透過度を0%にしてみる(alpha=1)
ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_line()+
 geom_ribbon(aes(ymin=Anomaly10y-Unc10y,ymax=Anomaly10y+Unc10y),alpha=1)

#なるほど見えない

f:id:WAFkw:20150523162208p:plain


コメントにも書いてますが、なるほど見えない。

一回透過80%でやってみたら実線が透けて見えちゃってたので、
最終的には黒く塗りつぶしてみた。


せっかくなので街で良く見る点線パターンでの描画もしてみる

#点線バージョン

ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_line(aes(y=Anomaly10y-Unc10y),colour="blue",linetype="dashed")+
  geom_line(aes(y=Anomaly10y+Unc10y),colour="red",linetype="dashed")+
  geom_line()

f:id:WAFkw:20150523162422p:plain

いろいろ遊んではみたけど、
linetypeはdottedよりdashedの方が見やすいと思う派でした。
colourも上限下限で色変えた方が自分のコードの反映状況が分かりやすいかなーとおもって。

またいろいろ遊んでみよう。


整体のじかんだー

どうしてもggplot2を使ってコロプレス図を描いてみたい。練習編

平日なので、というか仕事の日なので終始ねむい。
ねむいの。むかしよりもずっと。

最近、
日本の市町村のシェープファイルを取ってきてー
地図を描いてー
他のデータとマージさせてー
うまいこと色分けしたいー

みたいなことを考えて思うようにやったら失敗したので。笑

ちゃんと基礎から練習します。

シェープファイルから地図を描くことは出来たんだけどなあ)

私が思い描いているものは、
どうやらコロプレス図というらしい、、、

参考文献

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

友人に貸していたのが最近戻ってきた!
「やりたいこと」があった時に、辞書的に使ってます。

今回はこれの地図の章の例題をいじりながら使っていこうかなと。

#---まず地図を描いてみる----------------------

install.packages("mapproj")
library(maps)#地図データ読み込み
library(ggplot2)

#アメリカの地図データを取得
states_map=map_data("state")

#地図の描画(白く塗りつぶしたパターン)
ggplot(data=states_map,aes(x=long,y=lat,group=group))+
  geom_polygon(fill="white",colour="black")

##long=経度、lat=緯度。

f:id:WAFkw:20150519224923p:plain

まあまあまあ。
こないだ日本データでうまくいかなかったのでアメリカ地図なことは受け入れる、、、

#---いよいよ本番-------------------------------

#USArrestsデータセットの取得

#アメリカの犯罪発生件数というか、そんな感じのデータセットっぽいですね

crimes=data.frame(state=tolower(rownames(USArrests)),USArrests)

> crimes
                        state Murder Assault UrbanPop Rape
Alabama               alabama   13.2     236       58 21.2
Alaska                 alaska   10.0     263       48 44.5
Arizona               arizona    8.1     294       80 31.0
Arkansas             arkansas    8.8     190       50 19.5
California         california    9.0     276       91 40.6
Colorado             colorado    7.9     204       78 38.7
Connecticut       connecticut    3.3     110       77 11.1
Delaware             delaware    5.9     238       72 15.8
Florida               florida   15.4     335       80 31.9
Georgia               georgia   17.4     211       60 25.8
Hawaii                 hawaii    5.3      46       83 20.2
Idaho                   idaho    2.6     120       54 14.2
Illinois             illinois   10.4     249       83 24.0
Indiana               indiana    7.2     113       65 21.0
Iowa                     iowa    2.2      56       57 11.3
Kansas                 kansas    6.0     115       66 18.0
Kentucky             kentucky    9.7     109       52 16.3
Louisiana           louisiana   15.4     249       66 22.2
Maine                   maine    2.1      83       51  7.8
Maryland             maryland   11.3     300       67 27.8
Massachusetts   massachusetts    4.4     149       85 16.3
Michigan             michigan   12.1     255       74 35.1
Minnesota           minnesota    2.7      72       66 14.9
Mississippi       mississippi   16.1     259       44 17.1
Missouri             missouri    9.0     178       70 28.2
Montana               montana    6.0     109       53 16.4
Nebraska             nebraska    4.3     102       62 16.5
Nevada                 nevada   12.2     252       81 46.0
New Hampshire   new hampshire    2.1      57       56  9.5
New Jersey         new jersey    7.4     159       89 18.8
New Mexico         new mexico   11.4     285       70 32.1
New York             new york   11.1     254       86 26.1
North Carolina north carolina   13.0     337       45 16.1
North Dakota     north dakota    0.8      45       44  7.3
Ohio                     ohio    7.3     120       75 21.4
Oklahoma             oklahoma    6.6     151       68 20.0
Oregon                 oregon    4.9     159       67 29.3
Pennsylvania     pennsylvania    6.3     106       72 14.9
Rhode Island     rhode island    3.4     174       87  8.3
South Carolina south carolina   14.4     279       48 22.5
South Dakota     south dakota    3.8      86       45 12.8
Tennessee           tennessee   13.2     188       59 26.9
Texas                   texas   12.7     201       80 25.5
Utah                     utah    3.2     120       80 22.9
Vermont               vermont    2.2      48       32 11.2
Virginia             virginia    8.5     156       63 20.7
Washington         washington    4.0     145       73 26.2
West Virginia   west virginia    5.7      81       39  9.3
Wisconsin           wisconsin    2.6      53       66 10.8
Wyoming               wyoming    6.8     161       60 15.6

#データセットの統合
#地図の州とデータセットの州をうまいこと紐付ます

crime_map=merge(states_map,crimes,by.x="region",by.y="state")


#マージすると順番が入れ替わっちゃうみたいなので傍観

> head(crime_map)
   region      long      lat group order subregion Murder Assault UrbanPop Rape
1 alabama -87.46201 30.38968     1     1      <NA>   13.2     236       58 21.2
2 alabama -87.48493 30.37249     1     2      <NA>   13.2     236       58 21.2
3 alabama -87.95475 30.24644     1    13      <NA>   13.2     236       58 21.2
4 alabama -88.00632 30.24071     1    14      <NA>   13.2     236       58 21.2
5 alabama -88.01778 30.25217     1    15      <NA>   13.2     236       58 21.2
6 alabama -87.52503 30.37249     1     3      <NA>   13.2     236       58 21.2

#確かにregion(地域)順になってるっぽいことを確認。
#ここでのsubregionはregion以下の島のグループとかの細かい単位という認識です

#arrange関数でソートする
library(plyr)

#group,とorderで順にソート
crime_map=arrange(crime_map,order)

> #もっかい見てみる(orderのところが整った!)

> head(crime_map)
   region      long      lat group order subregion Murder Assault UrbanPop Rape
1 alabama -87.46201 30.38968     1     1      <NA>   13.2     236       58 21.2
2 alabama -87.48493 30.37249     1     2      <NA>   13.2     236       58 21.2
3 alabama -87.52503 30.37249     1     3      <NA>   13.2     236       58 21.2
4 alabama -87.53076 30.33239     1     4      <NA>   13.2     236       58 21.2
5 alabama -87.57087 30.32665     1     5      <NA>   13.2     236       58 21.2
6 alabama -87.58806 30.32665     1     6      <NA>   13.2     236       58 21.2


#変数の1つであるAssault(強盗事件)をfillにマッピング

ggplot(crime_map,aes(x=long,y=lat,group=group,fill=Assault))+
  geom_polygon()+
  coord_map("polyconic")

f:id:WAFkw:20150519225755p:plain

できたーーー
こないだはやっぱりデータのマージ部分orソートが良くなかったんだろうか・・・

別な変数も使って見てみます

#変数の1つであるAssault(強盗事件)をfillにマッピング
ggplot(crime_map,aes(x=long,y=lat,group=group,fill=Assault))+
  geom_polygon()+
  coord_map("polyconic")

f:id:WAFkw:20150519230045p:plain

おおお。

デフォルトの色を変えていないので
犯罪度数の多い方が「色が薄い」みたくなっていて
その点ちょっと直感的にはわかりにくくなってしまっておりますが。

どちらも南東側に州が犯罪度数高めって感じですかね、、

これをベースに色んなデータをマージさせて表示してみることからまず始めてみようかと。
アメリカ州別データならいろいろありそう~


今日は平日だということを言い訳にここまで。
ねむたい。

Rに関する主成分分析とSappoRo.Rのはなし

あと30分で月曜日がやってきてしまう。。

現実逃避がてら久々に最近ためになったSlideshareのはなし。

主成分分析


話の内容が、というより何よりも
スライドの作り方が非常に勉強になりました。

Rだったり機械学習だったり統計だったり
そこらへんのスライドは、たとえ

「データの可視化」の有用性をうたっていても
微妙に見栄えが悪かったり、

もちろんこのくらいはわかるよね・・・?
ってな具合の高尚な書かれ方をされていたり、

はたまた

あんまりにもやさしく書いて頂いていて
理論的根拠や数式や、実装イメージの無いものであったり

人に伝えるのって大変だよなあ。。って思わせてくることが多いんですが、

これは見せ方も理論的背景も実際の分析イメージもつきやすくて
わくわくしました。


主成分分析の話をすると、

次元をうまいこと圧縮したにも関わらず
他の分析手法に比べて解釈を主観に委ねる割合が高くない?

みたいな議論を半年前くらいにした記憶があるんですが、
このスライドでもそれに近いことがちらっと述べられていました。

主成分分析だけでデータに対するオチをつけるのは
かなりの技術力、というか語彙力を求められるよなあ、、という苦手意識があります。
いまのところ。


はてなブログだとこちらの記事も。

statsbeginner.hatenablog.com


私もこのくらい分かりやすい構成で話がしたい。
そして毎度のことながら金 明哲先生は本当にすごい。

SappoRo.Rのはなし

勇気を出して今年こそSappoRo.Rの参加を決めました。

ggplot2の話があるみたいなので

予習がてら一本テーマ決めてスライド作品作りでもしようかと思っていたものの
間に合いそうにないなあ、、、


身長体重データとか男女出生比率らへんでならなんか作れるかな。
やりかけたまま放置のコロプレスマップもあったりするけど、、、


聖書である高村本もしばらく読めてないし、
仕事を言い訳にしてグダグダしたらだめですね。


次は何をしようかなあ

【Rでのデータ操作】前処理のための基礎の基礎の基礎【マエショリストのススメ】

失踪を諦めて人混みと日常に戻ってきて数日が経ちました。

生まれてこのかた迷ってばかり。

必死になったってなんにもならない上にかえって損だって
わかっちゃいるんだけど。

                                • -

さて。
札幌に帰還したためR記事のリベンジをします。


定期的にマエショリスト(=狂ったように前処理をする人、たまにお金を稼ぐ)
を名乗っているので、前処理について。

データ分析における前処理とは


世間様の認識とはあまりかけ離れてないとは思うんですが、
私なりの前処理に対する考え方は、

表向きが
>>生のデータを、使用予定の分析の枠組みに適用できる形に整形し直すこと<<

心の中では
>>これが一番時間がかかるんだよもう疲れたよパトラッシュ<<

以上2点。

よくある具体的な前処理としては、

・生データの取得・・・WEBクローリングとか、e-statとか
・バラバラに取ってきた生データの結合・・・マージとか、キーを決めるとか
・外れ値の処理・・・外れ値の判断基準とか、外れ値を変換 or NA or そのままとか
・言語処理系・・・同一単語を漢字に統一するとか、置換とか、専門用語の辞書登録とか
・自動で結合したデータのうち、いらない列やら行やらを削除して投入可能状態にするとか

あとはケースによっていろいろです。


前処理をRでする意味

オープンデータだと、.xls形式とかでデータが取れたりするので、
そのままExcelで処理してしまったりした方が楽な場合も多いです。

が。

データ作業には手戻りが不可欠なので、
その場合のデータのバージョン管理だったり、長期に渡る作業だったりで、

「自分が以前どのような処理をどのような手順で行ったのか」

を残しておいて、「再現性を担保」することは超大事だと思ってます。

でもって、再現性を担保するにゃあ、コードが残ってた方がいいよねと。
でもって、自分はVBAが苦手なのでRでやるよ!と。
どうせ処理後のデータはRに投入しちゃうしね!と。

現場からは以上です。

RでExcelみたいなことがしたい

データ取り込みがしたい

data=read.csv("AAAAA.csv",header=TRUE) #csv形式の取り込み、変数名(ヘッダ)あり
data=read.table("clipboard",header=TRUE) #Excelとかでコピーしたものを取り込み

data=iris #Rにもともと入っているirisデータ読み込み
data=mtcars #Rにもともと入っているmtcarsデータ読み込み
<||

テキストデータはread.txtかな。
とりあえずデータをdataなりなんなり、任意のオブジェクトとして格納します。


**データの中身がみたい

ここからはRにもともと入っているデータである、
mtcarsデータセットを使います。irisでもいいです。

詳細は?mtcarsと打ち込めばヘルプが出ますが、
1973年あたりの車に関するデータが入ったデータセットになってます。

>|r|
mtcars
data=mtcars#mtcarsデータセットをdataとして格納

head(data,10) #格納したデータの上から10行を表示
length(data)  #データの行数を把握
dim(data)     #データの行列数を把握

summary(data) #各変数の分布状況(中央値、平均値やら)
str(data)     #データの構造や変数の型をチェック


#実行

> mtcars #mtcarsデータセット
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

> length(data)  #データの行数を把握
[1] 11

> dim(data)     #データの行列数を把握
[1] 32 11
> 
> summary(data) #各変数の分布状況(中央値、平均値やら)
      mpg             cyl             disp             hp             drat             wt             qsec      
Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0   Min.   :2.760   Min.   :1.513   Min.   :14.50  
1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5   1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89  
Median :19.20   Median :6.000   Median :196.3   Median :123.0   Median :3.695   Median :3.325   Median :17.71  
Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7   Mean   :3.597   Mean   :3.217   Mean   :17.85  
3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0   3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90  
Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0   Max.   :4.930   Max.   :5.424   Max.   :22.90  
       vs               am              gear            carb      
Min.   :0.0000   Min.   :0.0000   Min.   :3.000   Min.   :1.000  
1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
Median :0.0000   Median :0.0000   Median :4.000   Median :2.000  
Mean   :0.4375   Mean   :0.4062   Mean   :3.688   Mean   :2.812  
3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
Max.   :1.0000   Max.   :1.0000   Max.   :5.000   Max.   :8.000  

> str(data)     #データの構造や変数の型をチェック

'data.frame':     32 obs. of  11 variables:
$ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
$ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
$ disp: num  160 160 108 258 360 ...
$ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
$ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
$ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
$ qsec: num  16.5 17 18.6 19.4 17 ...
$ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
$ am  : num  1 1 1 0 0 0 0 0 0 0 ...
$ gear: num  4 4 4 3 3 3 3 4 4 4 ...
$ carb: num  4 4 1 1 2 1 4 2 2 4 ...

全データのうち、条件に合致するものを抽出したい

#抽出方法は2つ

#①subset関数を使用

> subset(data,data$mpg==23) #mpgが23のデータを全て抽出
[1] mpg  cyl  disp hp   drat wt   qsec vs   am   gear carb
<0> (または長さ 0 の row.names) #無かった・・・


#②dataオブジェクトをいじる。①と答えは同じになります。
data[data$mpg==23]  #mpgが23のデータを全て抽出

#①の応用

> subset(data,data$mpg>23)  #mpgが23より大きいもの

                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Fiat 128       32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic    30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla 33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Fiat X1-9      27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
Porsche 914-2  26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lotus Europa   30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2


> subset(data,data$mpg>=23) #mpgが23以上のもの

                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Fiat 128       32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic    30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla 33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Fiat X1-9      27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
Porsche 914-2  26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lotus Europa   30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
> 




*ほんとに基礎

「こんなの常識だろ」とか
「本読みゃ書いてあるだろ」とか

そんな感じではあるんだけど、個人的にググっても出てきにくくて
(出てきてもなんか七めんどくさい表現してあったり謎の英語だったり)

もうちょっとシンプルな記事あったら便利だよなって
むかし思ってたので。書いた。

〇〇がしたいけどRでどう書くっけ?ってなった時に、
せめて変数名が分かればその後ググりやすくなるし。

ぐらいの。

本当に詳しいことは関数ヘルプ見るなりもっと著名でスペシャルな方のサイト


とかを見て頂ければ。

次回は

**フィールドの追加や削除がしたい

**前処理における条件節の用途

**オブジェクトが増えてきて邪魔だから消したい

**各変数の型の変換


とかまとめる予定~

はてなブログでTex数式を書いてみる話(1)~自然言語処理シリーズ「言語処理のための機械学習入門」より~

WAF@失踪中です。

特に誰にも何も告げずに札幌を飛び出して数日が経ちました。
家出です。探さないでください。

「自分から逃げたい超逃げたいおなかいっぱい食べて寝たい!!!!」

こんなこと言ってる失踪中の身分にも関わらず
最新の論文が何かしら受賞したようでありがたかったです。
恐縮しきりです。
いいのかな。。貰っても。。


そんなこんなで家出中はいつものR環境も形態素解析環境も使えないので
バイブルの勉強をちまちまと。

聖書

言語処理のための機械学習入門 (自然言語処理シリーズ)

言語処理のための機械学習入門 (自然言語処理シリーズ)

My聖書こと、高村本です。

数式がとても丁寧でわかり易くて
頭が上がらない、ひれ伏してる。

おかげでめんどくさそうな数式が
怖くなくなりました!!(20代女性)

勉強したのをささっとまとめたい。しかし。

手書きの勉強ノート3冊を手に、はたと気が付いた。
私はLaTexに対してあまりにも無知であると。

だからいちいち手書きで数式を書き、
手書きで証明をし、
必要に迫られた時はofficeとかの数式エディタで
難を逃れてきたのだ。

論文製作でどれほど苦労したか思い出したくもない、、、

ということで、
バイブル読みながら今更ながらのLaTexデビューをします。

そもそもTexとかLaTexってなんぞや。

TeX入門 - TeX Wiki

元々文書を美しく効率的に記述するためのもの。
なので、数式みたい特殊な記述のオンパレードも効率良く記述出来るんですね。
数式書きまくる人には必須スキルというか出来て当たり前、と良く聞きます。

あたりまえ、、?(はて)

とりあえず自分に関しては、論文製作が佳境の頃

「うちらさ、、テフでやればよかったなんて言わないようにしようね」

と約束し合った友(計量経済専攻)が居るので過去は振り返りません。

はてなブログで出来るらしい

LaTex当時なんでやらなかったかって、
「だって環境構築めんどくさかったんだもん」
の一言に尽きます。

でもやってみたいなーということで。
こちらのサイトさんにちょいちょいコマンドを打ち込んで出力してみました。

TeX を使ってみよう

た、たのしい、、、

目的は統計関連の式をブログに適宜登場させる、というところなので
こちらのブログ様にお世話になりました。


はてなブログで数式を表示する方法。または、tex記法でLaTexを記述する方法。 - Pythonでも金融工学でもない。

聖書こと高村本のナイーブベイズ分類器周辺について、
LaTex方式で書いてみます。

実際にTex記法しながらナイーブベイズ分類器

・ナイーブベイズ分類器は確率に基づいた分類器
・事例dに対してP(c|d)がmaxとなるクラスcを出力する。

ベイズの定理
{ \displaystyle
P(c|d)=\frac{P(c)P(d|c)}{P(d)}
}

ベイズの定理における右辺が最大となるクラスcを出力するのが目的。
・しかし、P(d)はクラスcに依存しない。=文書dにのみ依存する

{ \displaystyle
c_{max}=arg max_c P(c)P(d|c)
}

上記により求めることが可能。
このP(d|c)の導出方法として、2種類いずれかのモデルが仮定される。

①多変量ベルヌーイモデル
②多項モデル

つまるところ、
P(d|c)に①②いずれかのモデルを仮定し、P(c)P(d|c)を最大化する
cの導出を行うのがナイーブベイズ分類器ってこって、、

はてなでのTex記法

唐突なナイーブベイズ。笑

[tex:{\displaystyle
数式
}]

みたいな感じで今回は書きました。
これでいつでも数式書けるようになるかなあ、、、