瓢箪鯰的な男がR初級者から中級者になるのを記すブログ

掴みどころのないオッサンのR勉強録。目標は中級者。 その他雑記なども

新・シルクロードを観ながら

国破山河在
城春草木深
杜甫『春望』)

**************************

 先日、家の中を片付けていたら
10年少し前に撮りためた「新・シルクロード」が出てきて
ちょっと片付けの手を止めて見入ってしまった。
その時妙な感じがしたので、ここに備忘的に書き留めておく。

**************************

 昭和55年に、NHKスペシャルで「シルクロードをやっていて
それを撮ったVHSテープが自宅にあり、好んでよく見ていた。
悠久とも言えるタリム盆地に、昔ながらのバザールをやっているさまをみて
自分もいつの日か西域に行きたいと思うようになった。
(その後治安が怪しいのと体調の加減で実現してはいないのだが)
歴史が好き、さらに「日本人はどこから来たのか?」を探るのが好きな自分にとって
シルクロードは歴史がそのまま保存されているような錯覚に陥らせてくれて
非常に心地がよかったのだ。

**************************

 タクラマカン砂漠、西端のカシュガル
昔からシルクロードの天山南路と西域南道が交わるオアシス都市だ。
40年前、この交差点は舗装もされておらず土埃の舞う
あまりにも牧歌的な道だった。

 それから30年。
道路は舗装され、近代的なビルが交差点を囲い
日干し煉瓦で作られた家が並んでいた通りは取り壊され
そこには瀟洒な商店が軒を連ねるようになっていた。
町は成長し、より豊かになっていた。

 それを見た時、本音を言うと非常に落胆してしまったのだ。

**************************

 長安の都が戦火に包まれても、ただ在ったのは山河だが
ここタクラマカン砂漠においては、河すらも不変ではない。
タリム川は雪解けの時を待たないと干上がったままであるし
楼蘭を支えたロプノール湖は、その位置を変えて旅人を悩ませた。

 敦煌の仏像、三蔵法師の物語だけを見れば
パミール高原以東のシルクロード仏教徒の道のように錯覚するが
ここは太古より戦乱が絶えない土地で、
目まぐるしく支配者が変わり、それごとに風俗も変わっていたのだ。

**************************

 悠久なるシルクロードとは幻想に過ぎず
絶えず変化するのがシルクロードである。
そこに僕は、勝手に「昨日的な今日」を期待してしまい
それが裏切られたからと言って落ち込むという、独り相撲を取っていたのだ。

**************************

 元来、成長とは変わることと不可分なはずである。
他者に対して豊かになってほしいとは、常に思っているのではあるが
それとは裏腹に変わらないことを望んでいたのだ…。

 この矛盾をどのようにうっちゃるか考えながら本棚を覗くと
学生の頃あれほど愛読したシュンペーターが埃をかぶって、
奥に押し込まれていたのだった。

ラーメン屋にて

 先日外出先で、中途半端な時間をうっちゃるがてら小腹がすいたので、
ラーメン屋の暖簾をくぐった。
食券を買ってそれを差し出してしばらくすると、
ラーメンが「へい!!お待ち!」という掛け声とともに出てきて
私はそれを湯気に隠れるようにすすりはじめた。

 不意にMy foolish heart が耳に入ってきた。
きっと店で有線契約しているであろう。洒落た調べであった。
その後も次々と品の良い曲がかかっていった。

 その時感じたのが、どうにも形容のし難い据わりの悪さである。
自分は現在、脂の浮いた色の濃いラーメンをすすっている腹の出た中年男であり
どう取り繕おうとも、Bill Evans の繊細な調べが似合う主体とは言えないのだ。

 その時はその違和感について不思議にも思わなかったのだが、
ラーメンの汁を飲み干して店を出た時、急に怪しげな気分になってきた。
さっき自分に「据わりの悪さ」を感じさせたものはなんだろう、と。

 よくよく考えてみると、ラーメンをすすりながらジャズを聴いてはならぬという法もなく
中年オッサンが洒落た場にいてはならぬという法もない。
さらに店内は店主と自分しかおらず、誰の目を気にするという場面でもなかった。
でも、やはり据わりが悪いのだ。

************************************

 こうしてみると、人が俗世間のいわゆる価値観から
完全に離れられる時間はあるのだろうかと、そう考えてしまう。
自分は自由気まま、勝手するままにふるまっているつもりでも
「俗世間のそれ」は、地球上のすべてが受ける重力加速度のように
まとわりついて離れないのではないか?
自分自身は勝手気ままにやっているつもりでも、
所詮は釈迦の掌で粋がっている孫悟空でしかないのか、と。

 「まっ、世の中そんなものだよな、、、」という乾いた笑いを自らに浴びせつつ
今宵もジャージ姿で、ショパンの調べに身を委ねるのである。

広告を非表示にする

PISA2006の分析 -親の属性と子の学力④-

前回からかなり時間が経ってしまったのだが、なんてことはない
あまり筆が乗らなかっただけの話である。

前回までの復習。

PISA2006のデータをダウンロードして、日本のデータまで切り出すところまで。

##########PISA2006の読み込み#############
setwd("C:/INT_Stu06_Dec07")	#データのあるディレクトリに移動(Cの直下にあると仮定)
moto<-readLines("INT_Stu06_Dec07.txt")	#データを読み込む
##########ここまで読み込み##########

jpn2006<-subset(moto,substr(moto,16,18)=="JPN") #日本だけ切り出す
kuni<-as.factor(substr(jpn2006,16,18))	#16列目から18列目までを切り出して、factor型でkuniへ代入
summary(kuni)	#日本だけ切り出せているか確認

さて、親の学歴のデータを切り出そう。
コードブックを見ると、379列目が母親の学力レベル、380列目が父親の学力レベルを表しているのが分かる。
381列目が「両親の学力レベル」とあるので、定義を確認する必要がある。
定義が分からないデータは、必ず定義を確認することが肝心だ。

さて、学力レベルがISCEDのレベルであらわされているので
適当なラベルを振ってやることにする。
ソースはウィキだが、レベル2が中卒、レベル3が高卒、レベル4とレベル5Bが難しいが、高専と短大が近い。レベル6は概ね大卒と見てよさそうだ。
レベル4は便宜的に高専、レベル5Bも便宜的に短大としておいて
あとから両方をくっつけても良い。

ちょっと両親の学力の定義を調べてみよう。
多分両親のうち、学力が高い方を選んでいると思われるのだが。。。
高専と短大は短大の方が処理の都合上、レベルが高い【ことにしている】だけである)

mg<-(substr(jpn2006,379,379))	#母親の最終学歴を切り出す
mg<-factor(mg, levels=c(2:6,9), labels=c("中卒", "高卒", "高専", "短大", "大卒", "不詳"))	#ラベル付け
fg<-substr(jpn2006, 380, 380)	#父親の最終学歴を切り出す
fg<-factor(fg, levels=c(2:6,9), labels=c("中卒", "高卒", "高専", "短大", "大卒", "不詳"))	#ラベル付け
rg<-substr(jpn2006, 381, 381)	#両親の最終学歴を切り出す
rg<-factor(rg, levels=c(2:6,9), labels=c("中卒", "高卒", "高専", "短大", "大卒", "不詳"))	#ラベル付け
gakureki<-data.frame(mg, fg, rg)	#学歴をデータフレーム化
head(gakureki, 10)	#gakurekiの頭10行を抜いてくる
   mg   fg   rg
1  高専 不詳 高専
2  大卒 大卒 大卒
3  大卒 大卒 大卒
4  大卒 大卒 大卒
5  大卒 大卒 大卒
6  短大 大卒 大卒
7  大卒 不詳 大卒
8  大卒 大卒 大卒
9  短大 大卒 大卒
10 大卒 大卒 大卒
11 大卒 大卒 大卒
12 短大 大卒 大卒
13 高専 大卒 大卒
14 高専 大卒 大卒
15 大卒 大卒 大卒
16 短大 高専 短大
17 短大 短大 短大
18 短大 大卒 大卒
19 短大 高専 短大
20 短大 大卒 大卒

やはりどちらか学力が高い方を取っているようである。

簡単に集計してみるとこんな感じ

summary(gakureki)
    mg          fg          rg      
 中卒: 247   中卒: 370   中卒: 142  
 高卒: 711   高卒: 782   高卒: 414  
 高専:2010   高専:1498   高専:1680  
 短大:1410   短大: 399   短大: 836  
 大卒:1309   大卒:2416   大卒:2702  
 不詳: 265   不詳: 487   不詳: 178  

親の学力のデータが抜けたので、次は点数のデータを抜いてみよう。

PISA2006の分析 -親の属性と子の学力③-

前回はコードブックや質問用紙を読みながら
使えそうな変数を取るところまでやった。
問題は点数が5つもあるというところの処置。
簡単に5つの平均値を出せばいいじゃん、と思われるかもしれないが。。。

平均は実は扱いにくい変数である

意外と知られていないのが、平均値の扱いづらさである。
平均値のマジックとしてよく知られているのが、「平均貯蓄」である。

【2017年家計調査】みんなの平均貯蓄額は?1820万円
当然「俺の家こんなに貯金ないよ」となるわけだ。

実はこの平均貯蓄は、一部の貯蓄が多い世帯が平均を引き上げてしまっているのだ。
(本文中にも言及があるが)実際に貯蓄額の分布を見てみると、
やはり貯蓄額の多い一部の世帯の影響が大きいことが分かる。
さらにこれは、負債は考慮していないので、貯蓄額は最低でもゼロ円だ。
平均値をある集団の特性を表す統計量だが
単純に平均しただけでは、その集団の特性を表さないことがある。
こちらから松坂大輔投手の
2005年の等級データを入手できる(提供元:(株)データスタジアム)。
平均球速を算出すると約138.1km/hとなって、全盛時の松坂投手からは考えられない程遅い。

ヒストグラムを書いてみよう。

setwd("C:/")	#ディレクトリ移動
moto<-read.csv("松坂球速.csv",header=TRUE)	#データの読み込み
hist(moto$kyusoku,main="松坂投手の球速分布",xlab="km/h",breaks = 18)	#球速のヒストグラムを表示。bureaksの値はいくつか実際に書いて決めればよい

f:id:namazu1945:20170604204726p:plain

明らかに2つ山があり、球速が何かしらの要因に左右されていることがわかる。
野球ファンならすぐにわかるだろうが、球種によって球速は変わってくる。
よって球速の平均値を出すには、全部を無暗に足して割るのではなく
球種で分類して、平均を出すという手間が生じてくる。

・・・なかなか本番の分析が進まないのぉ。。。

時系列分析「DECOMP法」の紹介

ちょっと今日は少し趣向を変えて。。。
たまには「これがデータ分析の花だ!!」見たいのを見せておかないと
やはり学習意欲が湧かないというのを、同僚を指導していて分かったんで。。

そんなわけで、今日はDECOMP法の紹介。


DECOMP法とは??

時系列データを扱うことは、仕事をしている方なら一度はあったかと思う。
日々の売り上げ、月の仕入れ額、月ごとの集客数etc。
時間で繋がっているデータのことを格好良く「時系列データ」と呼んでいて、
その分析を「時系列分析」と呼んでおります。

時系列データで興味があるのは、大きくわけて次の3つ。
①トレンドとして増えているOR減っている、上がっているor下がっている。またその水準。
②季節性(6月に傘が売れるとか、12月に餅が売れるとか、そういう季節特有の性質)
③イレギュラーな動き(ノイズ)

で、時系列データをこれら3つ(正確にはAR成分も入るので4つ)に分解してくれるのが
「DECOMP法」という便利なツール。
数学的な説明は、こちらの後ろの方をご覧いただきたい


さて、ちょっとやってみよう。

時系列データの例として、総務省の家計調査
2000年以降の時系列結果-二人以上の世帯の、消費支出を取ってくる。期間は2010年1月~2016年12月。
エクセルのブック形式で、時系列データが横に並んでいるというダメなフォーマットなので
必要なところを切り出して、手動でcsvに変換する。
ファイル名は「消費支出.csv」、データ名はsyouhiである。Cドライブの直下に置いておく。
Rで必要なコードを書く。
(本当は、1世帯当たり人員やCPIも使ってコントロールせにゃならんのだけど、それはまたの機会に)

setwd("C:/")	#ディレクトリ変更
moto<-read.csv("消費支出.csv",header=TRUE)	#消費支出のデータを取り込む
t.syouhi<-ts(moto$syouhi,start=c(2010,1),frequency=12)	#2010年1月からスタートし、12か月周期の時系列データに変換
ts.plot(t.syouhi,ylab="消費支出額",main="1世帯当たりの1か月の消費支出の推移")	#取りあえず折れ線グラフに

f:id:namazu1945:20170530201929p:plain

ジグザグしながら横ばい的な動き。

で、ここからがDECOMP法の実行だ。

library(timsac)	#DECOMPはTIMSACパッケージに含まれるので、それを呼び出す
kekka<-decomp(t.syouhi)	#DECOMPにかける

こんなグラフが表示される。
f:id:namazu1945:20170530202407p:plain

左上が原系列とトレンド成分、右上が季節成分
左下がノイズ成分、右下がAR成分である。
詳しい結果はkekkaというオブジェクトに格納されている。

このままだとグラフを読みづらいので、kekkaオブジェクトから
トレンド成分と、季節成分と、ノイズ成分を切り出して
同じように時系列データ化する。

t.trend<-ts(kekka$trend,start=c(2010,1),frequency=12)	#トレンド成分を時系列データ化
t.kisetsu<-ts(kekka$seasonal,start=c(2010,1),frequency=12)	#季節成分を時系列データ化
t.noise<-ts(kekka$noise,start=c(2010,1),frequency=12)	#ノイズ成分を時系列データ化

んで、グラフ化

ts.plot(t.trend,main="トレンド成分")	#トレンド成分をグラフ化
ts.plot(t.kisetsu,main="季節成分")	#季節成分をグラフ化
ts.plot(t.noise,main="ノイズ成分")	#ノイズ成分をグラフ化

トレンド成分
f:id:namazu1945:20170530203143p:plain
やっぱり消費税上げたあたりからトレンドとしても減ってきている感じですね。

季節成分
f:id:namazu1945:20170530203344p:plain
明確な季節性はありそう。

迷うのは、このノイズ成分。
f:id:namazu1945:20170530203642p:plain

2014年の3月から4月あたりで、大きく振れている。
これはどう考えても、消費税を上げた影響であって
ノイズ成分に分解できたからと言って、全部がノイズではないことを示している。
ノイズとは、ただの雑音ではなくて、中に何かしらの情報が潜んでいる場合があるので
注意深く扱ってやる必要がある。
こういうことが、データ分析では本当に多い。

DECOMP法は、アタリをつける上で非常に有用なツール。
ただ読み方は、他のツールと同じく注意して頂きたい。

手っ取り早くこのDECOMP法を試したい方は
こちらのWeb上でできるので試していただきたい。

PISA2006の分析 -親の属性と子の学力②-

 前回で、データをPCに読ますところまでできた。
さぁ張り切って○○分析を!!!とやりたいところだが・・・

分析の前にやることはたくさんある!!!*1

 まず用意するのは、
〇質問用紙(調査票)
〇コードブック
〇実際のデータ
この3つだ。いくつか本格的に分析にかける前にやるべきことを書いてみたい。


  • 3点セットから分析の内容をおおまかに決める

 ここで大事なのは、緩めに研究内容の枠を設定することである。
ガッチガチに「この変数とあの変数から、こういう分析をしたい」と考えても
実際のデータは、欲しい変数がない場合が多い(特に社会学系)。
従って緩めの枠を想定(ここでは何となく「親の状況と子供の学力」みたいな)しておいて
あとはどんなデータが実際に存在するのか、
存在するデータをみてから、研究の方向を絞り込んでいく。
やっぱりデータが何もないと、手も足も出ないということがあるので。



  • コードブックと質問用紙を見ながら、使えそうなデータを探す

 まず親に関する変数を挙げてみる。コードブック読み読み。
いくつか挙げてみよう。
〇母親の職業
〇母親の就業状態
〇母親の最終学歴
父親の職業
父親の就業状態
父親の最終学歴
〇両親の最終学歴(381行目)
これくらいは使えそうだが、質問用紙には「両親の最終学歴」なんて聞いていないので
この変数は後でどういう定義か後で確認しておく。
 
 ターゲットとなる変数は子供の学力なので、テストの点数を持ってくればよい、…のだが
数学、理科、国語とも、点数が5つある。
数学の点数は657~665列、666~674列、675~683列、684~692列、693-701列、この5つである。
理科も国語も同じ具合だ。さぁどうしたものか。。。

続きはまた次回。

*1:広い意味では分析だけどね

PISA2006の分析 -親の属性と子の学力①-

いつもお世話になっているフォロワーさんがこんなことを話していたので、
ちょっと簡単に分析してみましょうか、ってのがお題です。




PISAとは?

こちらのサイトが詳しいのだが、要するに世界的な学力調査試験である。
確か2006年は57の国と地域から、40万人くらい受けてた。
日本からは6000人が参加。

PISAのデータの何が素晴らしいか

 まず挙げられるのが、情報量の多さである。
ここに質問用紙があるが、本人への質問はともかく、
親の職業や学歴、家にある本の大まかな数、
携帯電話や車の数まで聞いている。
さらにこれらは、ひとりひとりのデータをOECDのサイトから
ダウンロードすることが可能で、各人が好きに分析できるという、
公的データとしては稀有なものなのだ。

分析の方法

 単純だが、最初は親の学歴などで子供をグループ分けして
そのグループ間で差があるか見ていくことにする。
で、分析をしていく過程で、何となく
「こういうことに気を付けてやってます」ということを書ければヨシかな、と。

まずはデータを読んで切り取る

 まずデータを読まないと仕方がないので、読ませてみる。
以下、Cドライブの直下にデータがあるものと仮定して、Rのコードを書いてみる。

##########PISA2006の読み込み#############
setwd(C:/)	#データのあるディレクトリに移動
moto<-readLines("INT_Stu06_Dec07.txt")	#データを読み込む
##########ここまで読み込み##########

motoというオブジェクトに全部のデータを読むことが出来たら、
次は日本のデータを切り出してくる。
こちらのコードブックを読むと、16列目から18列目に国のコードがあるので
ここが「JPN」のものだけ切り出してくる。

jpn2006<-subset(moto,substr(16,18)=="JPN") #日本だけ切り出す
kuni<-as.factor(substr(jpn2006,16,18))	#16列目から18列目までを切り出して、factor型でkuniへ代入
summary(kuni)	#日本だけ切り出せているか確認

結果は「JPN 5952」となるので、無事日本だけ抜き出せたことがわかる。
面倒でも狙った通りのデータが切り出せたかは、その都度確認しないといけない。

続きはまた後日。