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

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

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」となるので、無事日本だけ抜き出せたことがわかる。
面倒でも狙った通りのデータが切り出せたかは、その都度確認しないといけない。

続きはまた後日。

前回までで、PISA2006のデータを落とすまではやった。 今日は実際にデータをRで読み込むところまで。

#########サンプルコード#############
setwd("C:/01 ブログ/01 PISA2006/INT_Stu06_Dec07")	#データのあるディレクトリに異動
moto<-readLines("INT_Stu06_Dec07.txt")	#データを読み込む
#########ここまで#############

最初のsetwdで、作業ディレクトリに移動させる。
そこにあるテキストファイルを「readLines」で読み込み 「moto」というオブジェクトに代入する。慣れると簡単。

Rはテキストデータをはじめ、様々なフォーマットのデータを読み込めるようになっている。
foreignパッケージで、SASSPSSなどのデータを読める。 xlsxパッケージで、エクセルのブックを読める。 (噂によると、xlsxファイルはxmlファイルを圧縮したものらしいが、噂の真偽はまだ試していない)
さて、PISA2006のデータは700MBと結構重たいので メモリの上に全部展開するのに、数分かかる。 バカでかいファイルを扱う時には、メモリを増設したり SQLなどでデータを切り出してからRに読ませるなどの工夫が必要になる。

データを読ませたら、本当に読み込めたか試してみる必要がある。
PISA2006のコードブックを手元に用意して、次のコマンドを実行する。

#########先頭から100行だけ切り出してくる############# atama100<-head(moto,100) #先頭から100行切り出す 
fix(atama100) #先頭から100行を表示する 
#########切り出す作業終わり########### 

これを実行するとこんな画面が出てくる。
f:id:namazu1945:20170523211708p:plain

・・・色々苦しいが耐えてくれ

さて、コードブックを見ながら、ちゃんと読み込めたかチェックする。
コードブックのP2に、16列目から18列目に国コードが3つのアルファベットで載っているのがわかる。
本当に16列目から18列が国コードになっているか、確かめる。
よかった、ちゃんと読めているようだ。
面倒くさいと思われるかもしれないが、 データをある程度の大きさで切り出して、
コードブックと合わせるという作業は 必ずやらなくてはならない。
コードブックと実際のデータが異なることはよくある話で
もし実際のデータがコードブックと異なっていたら
その部分についてどういう処置を行うか、考えてから分析しなくてはならない。
次はそんな時の対処の方法を予定(カンタンな例ですが)。

ブログのネタ用データ

さて、ブログのネタになるデータをどうするか思案しているのだが PISA2006のデータが良さそうだ。 まず個人が扱える規模の大きな個別データがある、というのが大きい。 何よりも子供たちに将来の夢なんか聞いていて、データとしてなかなか楽しそうである点が良い。 (コードブックを見ると良い) PISAについての詳細は、こちらをご覧いただきたい。 データのダウンロードの方法は、上記のリンクから Data sets in TXT format (compressed)の下にある Student questionnaire data file (151 MB) こいつを、右クリック→名前を付けて保存 あとは適当なフォルダにこれを解凍するだけだ。 解凍すると700MBくらいのデータになるので要注意。 対応する質問用紙はこれだ。 Questionnairesの下にある Student questionnaireのリンク先である。 ちなみにP23にある、Q30が将来就きたい職を聞いている箇所である。 必要なコードブックはこれだ。 100ページ以上もあるブ厚いものだが、可能ならば紙に出すことをお勧めする。 もしくはマルチ画面で分析ソフトを回す画面とは別の画面に映すと良い。 作業効率が格段に違う。(やはりワシはアナログ人間なのだ) 質問用紙もコードブックも英語だらけ(というか英語だけ)なのだが 読んでみると簡単な英語なので、意外となんとかなる。 TOEIC500点もないワシでもどうにかなる程度なので、是非見てみてほしい。 暫くはこのデータを弄って遊ぶことにしよう。