平たく言えば「変数(variable)の種類」に相当する尺度水準(level of scales, scale level)は,測定水準(測定のレヴェル,levels of measurement)や測定尺度(measurement scales, scales of measurement)などとも言われ,level, scale, measurementの3つの語がやや入り組んでいて必ずしも厳密に使い分けられているとは言えない。
よく見られる変数の分類法と測定尺度の関係を表に纏めてみる。表の上段ほど尺度水準が低く,言い換えれば含まれる情報量が少なく,表の下段ほど水準が高い,つまり含まれる情報量が多い。
4尺度水準 | 変数分類A | 変数分類A’ | 変数分類B | 本書での提案 | ||
名義尺度 nominal |
質的変数 qualitative |
離散変数 | カテゴリカル変数 |
カテゴリカル変数 categorical |
名義尺度 | |
順序尺度 ordinal |
順序尺度 | |||||
間隔尺度 interval |
離散尺度 discrete |
量的変数 quantitative |
連続変数 | 連続変数 |
数量変数 numerical |
|
比例尺度 ratio |
連続尺度 continuous |
「質的/量的」の二分法を使用しない理由については『入門・社会調査法〔第3版〕』第2章を参照して欲しい。
区別が難しい場合があるのは,「順序尺度/間隔尺度」と「間隔尺度/比例尺度」であるが,後者の「間隔尺度/比例尺度」については,本書のレヴェルでは無理して区別しなくても支障ない。
尺度水準 | 例 | 可能な演算 |
---|---|---|
名義(nominal) | 性別,居住地,職業カテゴリ | 最頻値(mode) |
順序(ordinal) | 順位,意識調査項目への5件法回答 | 中央値(median) |
間隔(interval) | 試験の点数,西暦年号 | 算術平均(arithmetic mean),加・減 |
比例(ratio) | 時間,距離,金額 | 幾何平均(geometric mean),乗・除 |
interval scaleは「間隔尺度」以外に「距離尺度」,ratio scaleは「比例尺度」「比率尺度」「比尺度」の訳し方がある。
「間隔尺度/比例尺度」の区別の例としてよく見るのは,温度の単位の「摂氏温度℃」と「華氏温度℉」である。日本で20℃と表現する気温は,アメリカ合衆国では68℉と表現される。
摂氏温度と華氏温度の間には単純な関係があり,摂氏温度を1.8倍して32を加えれば華氏温度に変換できる。以下のRスクリプトで色々と試したりグラフ化したりすると良い。
摂氏温度(Celsius)から華氏温度(Fahrenheit)に変換
テクスト本文では華氏温度の変数名に F を使用しているが,後に不都合が生じるので Fd としておく。
(論理値の FALSE の 省略形 F と重複してしまう。)
C <- 20.0 # 任意の摂氏温度を入力
Fd <- 1.8*C + 32; names(Fd) <- c("Fahrenheit"); names(C) <- c("Celsius")
C; Fd
Celsius
20
Fahrenheit
68
華氏温度(Fahrenheit)から摂氏温度(Celsius)に変換
Fd <- 100.0 # 任意の華氏温度を入力
C <- 1/1.8*(Fd - 32); names(Fd) <- c("Fahrenheit"); names(C) <- c("Celsius")
Fd; C
Fahrenheit
100
Celsius
37.77778
摂氏を横軸,華氏を縦軸にしてグラフを書く
C <- -10:40 # 日本で比較的ありそうな気温の範囲に設定
Fd <- 1.8*C + 32
plot(C, Fd, bty="n", type="l", family="serif",
xlab="摂氏温度(Celsius)", ylab="華氏温度(Fahrenheit)",
main="摂氏温度と華氏温度の関係を表す直線")
segments(c(0, 25, 35), 0, c(0, 25, 35), 1.8*c(0, 25, 35)+32,
col="gray", lty=1)
segments(-12, 1.8*c(0, 25, 35)+32, c(0, 25, 35), 1.8*c(0, 25, 35)+32,
col="gray", lty=2)
摂氏が10℃から20℃に2倍になっても,華氏は50℉から68℉であり,2倍にはならない。実際には気温や温度という同一のものを測定している筈なので,2倍であるか・2倍では無いかのいずれかの筈であるが,測定法によって2倍になったりならなかったりと云う事は,(少なくとも一方の測定法は)比例尺度では無いと云う事を意味する。
摂氏の0℃は華氏では32℉であり,0℉ではない。逆に華氏の0℉は摂氏の-17.78℃で0℃ではない。華氏も摂氏も0には実体的な意味がなく,或意味自由に定められている。こうした,0に実体的な・他の値とは質的に異なる意味が無い測定方法(単位)は,比例尺度では無い。
では温度を表す比例尺度は何かと云うと,「絶対温度または熱力学温度(K; Kelvin)」がそれである。熱力学温度は分子の平均運動エネルギーを表しており,熱力学温度が2倍になると云う事は含まれる熱エネルギーが2倍になる事に対応しており,正しく比例尺度と考える事が出来る。0Kを絶対零度とも呼ぶが,これは古典力学的には原子が運動を停止している状態を指し,実体的な意味を持つ。そして,通常の意味においては,負の絶対温度は存在しないとされる(現在の統計力学では,例外的に「負の絶対温度」と云う概念が存在するらしいが)。
間隔尺度か比率尺度かを簡単に見分ける一つのポイントは,自然に・素直にマイナスの値を考える事が出来るかどうか,と言えよう。長さのメートル,質量のグラムなどは,0というのは「そもそも長さ/重さが存在しない」と云う特別な意味を持ち,「或物の長さがマイナス1メートルである」と云った表現は,素直な・自然な理解が出来ないものである。こうした長さや重さの単位は比例尺度である。
もう一つ区別が難しい,と云うか特に社会調査において厳密には問題となるのが「順序尺度/間隔尺度」である。
社会調査・世論調査では,例えば,「『夫は外で働き,妻は家庭を守るべきである』という考え方について賛成ですか,反対ですか」との質問に対して,「1 賛成,2 どちらかといえば賛成,3 どちらともいえない,4 どちらかといえば反対,5 反対」の選択肢から一つを選んで貰う回答方式(5件法)がよくある。
この回答を,1~5(もしくは逆順に5~1)の数字に変換して,平均を計算したり,平均の男女差や世代差を計算したりする事がよくあるが,平均を計算すると云う事はこれを,単なる順序尺度ではなく間隔尺度と見做している事になる。間隔尺度だと見做すと云う事は,「1 賛成」と「2 どちらかといえば賛成」の差は,「2 どちらかといえば賛成」と「3 どちらともいえない」の差に等しく,「3 どちらともいえない」と「5 反対」の差の2倍であると考える事を意味する。
こうした計算に対して,「この回答・選択肢は,選択肢間の距離(賛否の強度の違い?)が等間隔であると前提する事は出来ず,順序尺度でしかない」とする批判もある。そうであれば,例えば,その1~5の値をそのまま重回帰分析(第9章)の従属変数に使うと云った事は誤りである事になる。
回答選択肢としてはしばしば「3 どちらともいえない」が存在しない4件法も用いられる事を考えると,5件法の場合には1~5の値で等間隔と見做し,4件法の場合には1~4の値で等間隔と見做す事は確かに問題かも知れない。どちらかといえば賛成とどちらかといえば反対の間隔は,5件法では2になるが,4件法では1になってしまう。
しかし計量社会学では,こうした5件法(4件法)回答を間隔尺度扱いにして分析する習慣が長く存在している。必ずしもそれが認められないとは考えないが,議論の余地がある事は認識しておくのが良いだろう。因みに現在は,順序ロジット(第12章2-2)などのカテゴリカルデータ分析法が発達しているので,順序尺度のままで分析する事も十分に可能である。
本文で「xを変数のヴェクトルだとするとき」と書いているが,実際に演習できるように,模擬データを発生させてみよう。
n <- 125
x <- round(rnorm(n, mean=50, sd=10), 0)
取り敢えず今はこのコマンドの意味が分からなくても問題ない。演習用データとして変数(ヴェクトル)xが用意できたとだけ思っておけばよい。
このxに対して,データの個数,中央値,算術平均,そして最頻値を求める。(最頻値については若干の工夫が必要になる。)
【注意!】 最後の一行は,テクスト本文で見ると長い等号が一つ(“=”)に見えますが,半角の等号が二つ(“==”)です。御注意下さい。
length(x) # データの個数
[1] 125
median(x) # 中央値
[1] 49
算術平均の定義式
\[\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}=\color{red}{\frac{1}{n}\sum_{i=1}^{n}}x_{i}\]
mean(x) # 算術平均
[1] 49.976
sum(x)/length(x) # 合計をデータの個数で割って算術平均を求める
[1] 49.976
table(x) # 最頻値を求める準備の為に,度数分布表の全体を表示
x
31 32 33 36 37 38 39 40 41 42 43 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
1 1 3 1 3 2 5 7 8 5 6 6 6 2 4 4 3 7 1 3 5 4 4 5 4 5
60 61 62 63 65 67 69 70 71 74 75
2 3 1 4 1 2 2 1 1 1 2
table(x)[table(x)==max(table(x))] # 度数分布表から最頻値だけ抜粋
41
8
本文同様,幾何平均,調和平均,トリム平均も求めてみよう。
幾何平均の定義式
\[G=\sqrt[n]{x_{1}\cdot x_{2}\cdot x_{3}\cdots x_{i}\cdots x_{n}}=\left(\mathstrut \prod_{i=1}^{n}x_{i}\right)^\frac{1}{n}\]
prod(x)^(1/length(x)) # 幾何平均(相乗平均)①
[1] 49.0141
exp(sum(log(x))/length(x)) # 幾何平均②
[1] 49.0141
1/(sum(1/x)/length(x)) # 調和平均
[1] 48.07177
mean(x, trim=0.05) # 10%トリム平均
[1] 49.69027
通常の初歩的実習なら以上でも良いが,実践的な調査データの演習としては,この例は欠損値(NA)が含まれていない点で単純であり,実際のデータの集計の際に躓くことになってしまう。よって上の変数xの中の値を幾つか NA に書き換えて同じことを試して欲しい。
x[seq(1, n, by=round(n/6, 0))] <- NA # xのうち6つくらいの値をNAに書き換える
x # 置換した結果を確認してみる
[1] NA 38 60 58 63 40 51 45 56 42 42 58 33 71 47 56 54 59 46 41 54 NA 51 42 41
[26] 41 43 57 54 41 39 45 33 45 59 40 43 50 40 40 57 38 NA 50 37 49 51 48 40 39
[51] 49 55 51 61 75 46 50 41 61 57 40 48 58 NA 41 55 59 54 65 32 56 53 67 60 42
[76] 57 46 57 41 69 46 37 39 62 NA 36 61 58 56 37 41 53 31 70 55 63 45 54 43 43
[101] 39 51 53 67 49 NA 46 33 48 46 52 63 74 47 51 45 69 45 43 49 43 55 63 42 59
これで上の length( ), median( ), mean( ), sum( ), table( ), prod( )などを実行すると,length( )とtable( )以外は結果がNAになる。
length(x) # データの個数
[1] 125
median(x) # 中央値
[1] NA
mean(x) # 算術平均
[1] NA
sum(x)/length(x) # 合計をデータの個数で割って算術平均を求める
[1] NA
table(x) # 最頻値を求める準備の為に,度数分布表の全体を表示
x
31 32 33 36 37 38 39 40 41 42 43 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
1 1 3 1 3 2 4 6 8 5 6 6 6 2 3 4 3 6 1 3 5 4 4 5 4 4
60 61 62 63 65 67 69 70 71 74 75
2 3 1 4 1 2 2 1 1 1 1
table(x)[table(x)==max(table(x))] # 度数分布表から最頻値だけ抜粋
41
8
面倒に思うだろうが,以下のように常にNAを意識して処理することになれるのが良い。 ポイントはデータの個数(ヴェクトルの長さ)を求める時にもNAを除外することを忘れないということであり,最初の length( )の中の“!”は否定を意味し,is.na(x) は「xのうちNAであるもの」を意味する。合わせて !is.na(x) で「xのうちNAではないもの」を意味するので,x[!is.na(x)] でxからNAではないものだけを抽出し,length(x[!is.na(x)]) でその個数(長さ)を求めている。
length(x[!is.na(x)]) # NAを除いたデータの個数
[1] 119
median(x, na.rm=T) # 中央値
[1] 49
mean(x, na.rm=T) # 算術平均
[1] 49.87395
sum(x, na.rm=T)/length(x[!is.na(x)]) # 合計をデータの個数で割って算術平均を求める
[1] 49.87395
table(x, useNA="ifany") # table( )では逆に,NAを省略させないことがしばしば重要
x
31 32 33 36 37 38 39 40 41 42 43 45 46 47 48 49
1 1 3 1 3 2 4 6 8 5 6 6 6 2 3 4
50 51 52 53 54 55 56 57 58 59 60 61 62 63 65 67
3 6 1 3 5 4 4 5 4 4 2 3 1 4 1 2
69 70 71 74 75 <NA>
2 1 1 1 1 6
table(x)[table(x)==max(table(x))] # 最頻値にはNAは関係ないので省略して良い
41
8
prod(x, na.rm=T)^(1/length(x[!is.na(x)])) # 幾何平均(相乗平均)①
[1] 48.93268
exp(sum(log(x), na.rm=T)/length(x[!is.na(x)])) # 幾何平均②
[1] 48.93268
1/(sum(1/x, na.rm=T)/length(x[!is.na(x)])) # 調和平均
[1] 48.00489
mean(x, trim=0.05, na.rm=T) # 10%トリム平均
[1] 49.66972
上ではあえていちいち length(x[!is.na(x)]) を書いたが,実際にはこれは n1 <- length(x[!is.na(x)]) などとしてオブジェクト(n1)に代入しておいて,このn1を必要なところで使用するのが良い。
分散を求めるには var(x) とすればよい。しかし分散が何であるのかをより深く理解する為に,本文同様,var( )関数を用いずに分散を求めてみよう。
偏差平方和の定義式
\[SS(Sum\;of\;Squares)=\sum_{i=1}^{n}(x_{i}-\bar{x})^2\] 分散の定義式
\[var(x)=s_{x}^2=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2=\color{red}{\frac{1}{n}\sum_{i=1}^{n}}(x_{i}-\bar{x})^2\]
x # まずは生成した変数(ヴェクトル)そのものを表示させてみよう。
[1] NA 38 60 58 63 40 51 45 56 42 42 58 33 71 47 56 54 59 46 41 54 NA 51 42 41
[26] 41 43 57 54 41 39 45 33 45 59 40 43 50 40 40 57 38 NA 50 37 49 51 48 40 39
[51] 49 55 51 61 75 46 50 41 61 57 40 48 58 NA 41 55 59 54 65 32 56 53 67 60 42
[76] 57 46 57 41 69 46 37 39 62 NA 36 61 58 56 37 41 53 31 70 55 63 45 54 43 43
[101] 39 51 53 67 49 NA 46 33 48 46 52 63 74 47 51 45 69 45 43 49 43 55 63 42 59
x1 <- x[!is.na(x)] # NAが含まれていると面倒なので,NAを除外したものを新たに x1 とする
mean(x1) # 算術平均
[1] 49.87395
x1-mean(x1) # (算術平均からの)偏差(deviation)
[1] -11.8739496 10.1260504 8.1260504 13.1260504 -9.8739496 1.1260504
[7] -4.8739496 6.1260504 -7.8739496 -7.8739496 8.1260504 -16.8739496
[13] 21.1260504 -2.8739496 6.1260504 4.1260504 9.1260504 -3.8739496
[19] -8.8739496 4.1260504 1.1260504 -7.8739496 -8.8739496 -8.8739496
[25] -6.8739496 7.1260504 4.1260504 -8.8739496 -10.8739496 -4.8739496
[31] -16.8739496 -4.8739496 9.1260504 -9.8739496 -6.8739496 0.1260504
[37] -9.8739496 -9.8739496 7.1260504 -11.8739496 0.1260504 -12.8739496
[43] -0.8739496 1.1260504 -1.8739496 -9.8739496 -10.8739496 -0.8739496
[49] 5.1260504 1.1260504 11.1260504 25.1260504 -3.8739496 0.1260504
[55] -8.8739496 11.1260504 7.1260504 -9.8739496 -1.8739496 8.1260504
[61] -8.8739496 5.1260504 9.1260504 4.1260504 15.1260504 -17.8739496
[67] 6.1260504 3.1260504 17.1260504 10.1260504 -7.8739496 7.1260504
[73] -3.8739496 7.1260504 -8.8739496 19.1260504 -3.8739496 -12.8739496
[79] -10.8739496 12.1260504 -13.8739496 11.1260504 8.1260504 6.1260504
[85] -12.8739496 -8.8739496 3.1260504 -18.8739496 20.1260504 5.1260504
[91] 13.1260504 -4.8739496 4.1260504 -6.8739496 -6.8739496 -10.8739496
[97] 1.1260504 3.1260504 17.1260504 -0.8739496 -3.8739496 -16.8739496
[103] -1.8739496 -3.8739496 2.1260504 13.1260504 24.1260504 -2.8739496
[109] 1.1260504 -4.8739496 19.1260504 -4.8739496 -6.8739496 -0.8739496
[115] -6.8739496 5.1260504 13.1260504 -7.8739496 9.1260504
mean(abs(x1-mean(x1))) # 平均偏差
[1] 8.156486
(x1-mean(x1))^2 # 偏差平方
[1] 140.99067862 102.53689711 66.03269543 172.29319963 97.49488031
[6] 1.26798955 23.75538451 37.52849375 61.99908199 61.99908199
[11] 66.03269543 284.73017442 446.31000636 8.25958619 37.52849375
[16] 17.02429207 83.28479627 15.00748535 78.74698115 17.02429207
[21] 1.26798955 61.99908199 78.74698115 78.74698115 47.25118283
[26] 50.78059459 17.02429207 78.74698115 118.24277946 23.75538451
[31] 284.73017442 23.75538451 83.28479627 97.49488031 47.25118283
[36] 0.01588871 97.49488031 97.49488031 50.78059459 140.99067862
[41] 0.01588871 165.73857778 0.76378787 1.26798955 3.51168703
[46] 97.49488031 118.24277946 0.76378787 26.27639291 1.26798955
[51] 123.78899795 631.31840972 15.00748535 0.01588871 78.74698115
[56] 123.78899795 50.78059459 97.49488031 3.51168703 66.03269543
[61] 78.74698115 26.27639291 83.28479627 17.02429207 228.79740131
[66] 319.47807358 37.52849375 9.77219123 293.30160299 102.53689711
[71] 61.99908199 50.78059459 15.00748535 50.78059459 78.74698115
[76] 365.80580467 15.00748535 165.73857778 118.24277946 147.04109879
[81] 192.48647694 123.78899795 66.03269543 37.52849375 165.73857778
[86] 78.74698115 9.77219123 356.22597274 405.05790552 26.27639291
[91] 172.29319963 23.75538451 17.02429207 47.25118283 47.25118283
[96] 118.24277946 1.26798955 9.77219123 293.30160299 0.76378787
[101] 15.00748535 284.73017442 3.51168703 15.00748535 4.52009039
[106] 172.29319963 582.06630888 8.25958619 1.26798955 23.75538451
[111] 365.80580467 23.75538451 47.25118283 0.76378787 47.25118283
[116] 26.27639291 172.29319963 61.99908199 83.28479627
sum((x1-mean(x1))^2) # 偏差平方和
[1] 11303.11
sum((x1-mean(x1))^2)/length(x1) # 分散①
[1] 94.98411
mean((x1-mean(x1))^2) # 分散②
[1] 94.98411
var(x1) # Rの分散の関数.上の①や②とは一致しない.
[1] 95.78906
sum((x1-mean(x1))^2)/(length(x1)-1) # 不偏分散
[1] 95.78906
標準偏差の定義式
\[sd(x)=s_{x}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2}\]
sqrt(mean((x1-mean(x1))^2)) # 標準偏差
[1] 9.745979
sd(x1) # Rの標準偏差の関数.上とは一致しない.
[1] 9.787189
sqrt(sum((x1-mean(x1))^2)/(length(x1)-1)) # 不偏分散の平方根
[1] 9.787189
今度は,カテゴリカル変数らしい模擬データを生成してみよう。本文2-1の表と,全く同じとはならないが,よく似た変数になるように工夫している。
var0 <- sample(c(1, 2, 3, 4, NA), size=384, replace=T, prob=c(.474, .240, .099, .177, .010))
シンプルな度数分布表やそれを元にした集計を以下に示す。
table(var0) # シンプルな度数分布表
var0
1 2 3 4
179 99 41 64
prop.table(table(var0)) # 比率に変換
var0
1 2 3 4
0.4673629 0.2584856 0.1070496 0.1671018
cumsum(table(var0)) # 累積度数
1 2 3 4
179 278 319 383
本文2-1のような度数分布表を作成するスクリプトを以下で説明しよう。
var1 <- var0 # 集計したい変数をそのままvar1に代入。こうすることで汎用的スクリプトになる
t0 <- table(var1, useNA="always") # NAを含んだ総度数
dimnames(t0)[[1]] <- c("1 正規雇用","2 非正規雇用","3 自営","4 無職・学生","NA 欠損値")
t1 <- table(var1) # NAを除外した有効度数
p0 <- prop.table(t0) # NAを含んだ相対度数
ct0 <- cumsum(t0) # NAを含んだ累積度数
cp0 <- cumsum(p0) # NAを含んだ累積相対度数
p1 <- prop.table(t1) # 有効相対度数
cp1 <- cumsum(p1) # 累積有効相対度数
p1 <- c(p1, 0) # NAを含含むか含まないかによって行数が異なるので調整
cp1 <- c(cp1, 0) # NAの分行数が異なるので調整
ft0 <- cbind(t0, ct0, round(p0*100,1), round(cp0*100,1), round(p1*100,1), round(cp1*100,1))
dimnames(ft0)[[2]] <- c("度数","累積度数","%","累積%", "有効%", "累積有効%")
ft0 # 最終的な度数分布表を表示
度数 累積度数 % 累積% 有効% 累積有効%
1 正規雇用 179 179 46.6 46.6 46.7 46.7
2 非正規雇用 99 278 25.8 72.4 25.8 72.6
3 自営 41 319 10.7 83.1 10.7 83.3
4 無職・学生 64 383 16.7 99.7 16.7 100.0
NA 欠損値 1 384 0.3 100.0 0.0 0.0
二つの dimnames( )関数は,値ラベルや表頭の列名を付けるものである。
t0からcp1までパーツを準備し,それを「列結合(cbind)」でひとまとめにしたのが ft0 である。
上記のスクリプトで,var1 <- var0 (本文で言えば var1 <- data$job)の代入式の右辺だけを変更すれば,任意の変数の度数分布表が同じスクリプトで作成できる。
模擬データを生成して,ヒストグラムや箱ひげ図で図示してみよう。hist( )もboxplot( )もNAを除外して作図してくれるので,NAを含まない単純な模擬データで演習する。
一つずつ試してみると,主なオプションの働きが分かるだろう。ここでは煩雑になるので,主な描画結果のみ掲載する。それ以外は自分でスクリプトを実行して確認して欲しい。
x <- rnorm(n=400, mean=50, sd=10) # 模擬データ
hist(x) # 最も単純なヒストグラム
hist(x, xlim=c(10, 90), ylim=c(0, 150)) # x軸の範囲,y軸の範囲を指定
hist(x, xlim=c(10, 90), ylim=c(0, 150), main="変数xのヒストグラム", xlab="変数の値", ylab="度数") # タイトルと軸ラベルを指定
hist(x, xlim=c(10, 90), ylim=c(0, 100),
main="変数xのヒストグラム", xlab="変数の値", ylab="度数",
breaks=seq(0, 100, by=5)) # x軸の区切り位置を指定
hist(x, xlim=c(10, 90), freq=FALSE, # y軸を,度数ではなく密度に
main="変数xのヒストグラム", xlab="変数の値", ylab="密度",
breaks=seq(0, 100, by=5))
hist(x, xlim=c(10, 90), freq=FALSE, # x軸の分割を不均等に,y軸は密度
main="変数xのヒストグラム", xlab="変数の値", ylab="密度",
breaks=c(0, 20, 30, 40, 45, 50, 55, 60, 70, 80, 100))
これも一つ一つ実行して,それぞれのオプションの意味を理解しよう(ここでは途中の幾つかの描画は省略する)。
x <- round(rnorm(n=100, mean=50, sd=10), 0)
y <- round(rnorm(n=200, mean=50, sd=15), 0) # ケース数と散布度の異なるヴェクトルを用意
boxplot(x) # 最も単純な箱ひげ図
boxplot(x, ylim=c(0, 100), # y軸の範囲を指定
main="xの箱ひげ図", ylab="変数の値") # タイトルとy軸ラベル
boxplot(x, y, ylim=c(0, 100), # 2つのヴェクトルを並べて図示
main="xとyの箱ひげ図", ylab="変数の値")
boxplot(x, y, ylim=c(0, 100),
main="xとyの箱ひげ図", ylab="変数の値",
names=c("ヴェクトルx", "ヴェクトルy")) # x軸にラベル
boxplot(x, y, ylim=c(0, 100), varwidth=T, # 幅がサンプルサイズの平方根に比例
main="xとyの箱ひげ図", ylab="変数の値",
names=c("ヴェクトルx", "ヴェクトルy"))
boxplot(x, y, ylim=c(0, 100), varwidth=T,
main="xとyの箱ひげ図", ylab="変数の値",
names=c("ヴェクトルx", "ヴェクトルy"),
col=c("cyan", "pink"), border=c("blue", "red")) # 色を付けてみた
もう少し他の作図法も紹介しておこう。
時系列変化などを表すには折れ線グラフを用いる。ここでは,総務省の「労働力調査」の長期時系列データにおける「完全失業率」の年次変化をグラフ化してみる(長期時系列表2 就業状態別15歳以上人口-全国,列M: 完全失業率(%))。2011年の数値だけ他とは算出方法が異なるので色を変える。
year <- c(1973:2015)
unemployment <- c(1.3, 1.4, 1.9, 2.0, 2.0, 2.2, 2.1, 2.0, 2.2, 2.4, 2.6, 2.7, 2.6, 2.8, 2.8, 2.5, 2.3, 2.1, 2.1, 2.2, 2.5, 2.9, 3.2, 3.4, 3.4, 4.1, 4.7, 4.7, 5.0, 5.4, 5.3, 4.7, 4.4, 4.1, 3.9, 4.0, 5.1, 5.1, 4.6, 4.3, 4.0, 3.6, 3.4)
marker_col <- rep("black", length(year))
marker_col[year == 2011] <- c("red")
# 比較的シンプルなグラフ(コメントアウトしておいてここでは描画しない)
# plot(year, unemployment, bty="n", type="b", pch=16, col=marker_col)
# オプションなどを活用したグラフ
plot(year, unemployment, bty="n", type="b", pch=16, col=marker_col,
axes=FALSE, family="serif", ylim=c(0, ceiling(max(unemployment))),
main="総務省「労働力調査」による完全失業率の長期的推移", xlab="西暦年",
ylab="完全失業率(%)", cex.main=0.8)
axis(side=1, at=c(1973, seq(1980, 2015, by=5)), family="serif", las=2)
axis(side=2, at=c(0:ceiling(max(unemployment))), family="serif", las=1)
もう一つは幹葉図と呼ばれるものであるが,これは,テクスト(文字)のみで表現するのでグラフ機能を必要としない割には詳しい情報を表示すると云う利点がある一方で,データの数が多いと有効でなくなると云う欠点がある。ここでは,第3章1-2で使用する模擬年収データ(4,894人分)から200人分を無作為抽出して幹葉図で表現してみよう。比較のために同じデータをヒストグラムで描画する。
download.file("http://sgn.sakura.ne.jp/text/mock_income.csv", "mock_income.csv")
data02 <- read.csv("mock_income.csv")
str(data02)
'data.frame': 4894 obs. of 2 variables:
$ case : int 1 2 3 4 5 6 7 8 9 10 ...
$ income: int 799 387 605 301 505 243 771 289 359 210 ...
# 個人年収200人分を無作為抽出(単位: 万円)
x <- sample(data02$income, size=200, replace=FALSE)
stem(x)
The decimal point is 2 digit(s) to the right of the |
0 | 000000334588899
1 | 1223334444557788888999
2 | 00112233445566677888999
3 | 00112233444555566666677888999
4 | 00000122233334444566667777777899
5 | 001112222334444455566667788
6 | 01112235568
7 | 000112344455678
8 | 00122689
9 | 124479999
10 | 78
11 |
12 | 4
13 | 6
14 | 07
15 |
16 | 46
17 |
18 |
19 | 4
hist(x, breaks=seq(0, ceiling(max(x)/100)*100, by=100))
summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 255.8 432.0 473.1 607.2 1943.0
summary( )の結果やヒストグラムと比較して幹葉図を見ると読み取れると思うが,この例の場合には,年収の百万円以上の桁の部分がセパレータ(|)の左側(表側)に示され,その下の位である十万円の桁の数字が右側に昇順で羅列されている。
例えば, 10|4677 の行は,1040万円台が1名,1060万円台が1名,1070万円台が2名いる事を意味している。 この様に一つ一つのケースを数字一つで表す為に,大量のデータの表示には向かない。
テクストにはないが,グラフ描画に定評の高い ggplot2 パッケイジによる描画例を示しておく。
上のカラーの「xとyの箱ひげ図」と同様のグラフを,ggplot2 で描画して見る。
ggplot2を使うには,使用する変数はデータフレイムになっている必要があるので,それをgdataとして作成しておく。
x <- round(rnorm(n=100, mean=50, sd=10), 0)
y <- round(rnorm(n=200, mean=50, sd=15), 0)
score <- c(x, y)
group <- c(rep("ヴェクトルx", length(x)), rep("ヴェクトルy", length(y)))
gdata <- data.frame(group, score)
g10 <- ggplot(gdata, aes(x=group, y=score))
g11 <- g10 + geom_boxplot(aes(color=group), varwidth=T)
g11
上のヒストグラムで使用した data02 のデータフレイムを再利用する。
g20 <- ggplot(data02, aes(x=income))
g21 <- g20 + geom_histogram(alpha=.2)
g21
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
もっと様々なオプションや装飾があるが,基本形のみ紹介しておくので,興味のある人は自分で調べると良い。
(オプションが多過ぎて途方に暮れるかもしれないが。
正規分布\(N(\color{blue}{\mu}, \color{blue}{\sigma^2})\)の定義式 \[y=\frac{1}{\sqrt{2\pi\color{blue}{\sigma^2}}}e^{-\frac{(\color{red}{x}-\color{blue}{\mu})^2}{2\color{blue}{\sigma^2}}}\]
理論的な正規分布のグラフ描画については〈このページ〉で紹介しているので参照して欲しい。以下にその一部を微修正して紹介する。是非解読して試してほしい。
標準正規分布において,面積(%)を指定すると,対応するx座標を求めてグラフで示す
p <- .90 # ここは90%なら.90,95%なら.95,99%なら.99など,自分で好きに指定
q <- qnorm((1+p)/2, mean=0, sd=1)
curve(dnorm(x, mean=0, sd=1), from=-3.5, to=3.5,
bty="n", xaxp=c(-3, 3, 2),
xlab="", ylab="確率密度", main="標準正規分布の面積", col="#339900")
abline(h=0, col="#339900")
segments(x0 <- seq(-1*q, q, by=.01), 0, x0, dnorm(x0, 0, 1), col="#33990080")
axis(side=1, at=round(c(-1*q, q), 3), col.axis="#006600", cex=0.8, las=2)
text(0, 0.2, paste("面積", p*100, "%"))
次は,x座標を指定すると,対応する面積を計算してグラフで示すスクリプト。
標準正規分布において,x座標を指定すると,対応する面積(%)を求めてグラフで示す
q <- 2.0 # x座標を,プラスの値で自分で好きに指定
p <- 1-(1-pnorm(q, mean=0, sd=1))*2
curve(dnorm(x, mean=0, sd=1), from=-4, to=4,
bty="n", xaxp=c(-4, 4, 2),
xlab="", ylab="確率密度", main="標準正規分布の面積", col="#339900")
abline(h=0, col="#339900")
segments(x0 <- seq(-1*q, q, by=.01), 0, x0, dnorm(x0, 0, 1), col="#33990080")
axis(side=1, at=round(c(-1*q, q), 3), col.axis="#006600", cex=0.8, las=2)
text(0, 0.2, paste("面積", round(p*100, 1), "%"))
標準化\(standardization\) \[\color{red}{z_{i}}=\frac{\color{red}{x_{i}}-\color{blue}{\bar{x}}}{\color{blue}{s_{x}}},\;\;(i=1,2,3,...,n)\;\;\color{green}{\Rightarrow}\;\bar{z}=0,\;s_{z}=1\]
ここでは,試験の点数を模した模擬データを発生させて,それを標準化して標準化スコアを求め,更に受験業界でよく知られている「偏差値」に換算して比較してみよう(素点を発生させる部分は試行錯誤で余計な工夫をしているので余り気にしなくて良い)。
n <- 200 # 取り敢えず200人分とする
score <- rpois(n, 60) # それっぽい素点の生成
score[score < 0] <- 0 # 念の為の処理
score[score > 100] <- 100 # 念の為の処理
z <- (score - mean(score))/sd(score) # 事情により手動で標準化
Z <- 10 * z + 50 # 標準化スコアから,「偏差値」を換算
summary(cbind("素点"=score, "標準化スコア"=z, "「偏差値」"=Z)) # 素点
素点 標準化スコア 「偏差値」
Min. :43.00 Min. :-2.2249 Min. :27.75
1st Qu.:53.75 1st Qu.:-0.7426 1st Qu.:42.57
Median :58.00 Median :-0.1565 Median :48.43
Mean :59.13 Mean : 0.0000 Mean :50.00
3rd Qu.:64.00 3rd Qu.: 0.6709 3rd Qu.:56.71
Max. :79.00 Max. : 2.7393 Max. :77.39
hist(score, breaks=seq(0, 100, by=2), main="素点分布", ylab="人数", xlab="")
hist(z, breaks=seq(-4, 4, by=.2), main="標準化スコア", ylab="", xlab="")
hist(Z, breaks=seq(0, 100, by=2), main="「偏差値」", ylab="", xlab="")
スクリプトに#でコメントを付けてある様に,0点から100点の間の点数をn=200人分発生させ,それを標準化したzを求め,更にそこからいわゆる「偏差値」に換算したZを求めている。
次に,それぞれの要約統計量をsummary( )関数で表示させている。
その次に,それぞれのヒストグラム3つを描いている。。
尚,元の変数が正規分布ではない限り,標準化スコアも正規分布にはならない。標準化すれば必ず「平均が0に,分散と標準偏差が1に」なるが,標準正規分布になるかどうかは,元の変数が正規分布かどうかに依存する。上のスクリプトの score <- の行だけを例えば score <- runif(n, min=30, max=80) と書き換えて試してみるとはっきりするだろう(以下にヒストグラムのみ例示)。
score <- runif(n, min=30, max=80) # それっぽい素点の生成を工夫
z <- (score - mean(score))/sd(score) # 事情により手動で標準化
Z <- 10 * z + 50 # 標準化スコアから,「偏差値」を換算
hist(score, breaks=seq(0, 100, by=2), main="素点分布", ylab="人数", xlab="")
hist(z, breaks=seq(-4, 4, by=.2), main="標準化スコア", ylab="", xlab="")
hist(Z, breaks=seq(0, 100, by=2), main="「偏差値」", ylab="", xlab="")
母分散
\[\sigma_{x}^2=\frac{1}{N}\sum_{i=1}^{N}(x_{i}-\bar{x})^2\]
標本分散
\[s_{x}^2=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2\]
不偏分散
\[\hat{\sigma}_{x}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2\]
人為的に発生させたデータについて,標本分散・標本標準偏差を定義式に従って求め,次にRの関数var( ),sd( )の結果と比べてみる。
n <- 50
x <- rnorm(n, mean=50, sd=10) # データの生成
# 平均からの偏差: x - mean(x)
# 偏差平方: (x - mean(x))^2
sum((x - mean(x))^2) # 偏差平方和
[1] 3801.441
sum((x - mean(x))^2)/n # 定義式に従った標本分散
[1] 76.02882
mean((x - mean(x))^2) # 偏差平方の平均=標本分散
[1] 76.02882
sqrt(mean((x - mean(x))^2)) # 標本標準偏差
[1] 8.719451
# Rの分散var( ),標準偏差sd( )
var(x); sd(x)
[1] 77.58043
[1] 8.807975
# 定義式に従った不偏分散
sum((x - mean(x))^2)/(n-1)
[1] 77.58043
定義式に従った標本分散はRの関数var( )の結果とは一致しない事,nではなくn-1で割った不偏分散がvar( )の結果と一致する事,当然ながら不偏分散の方が標本分散よりも少し大きな値になる事,を確認しよう。
n=50の場合は,50/49=1.020408で,不偏分散の方が2%程度大きくなるが,n=500になるとこれが0.2%の違いでしかなくなるので,実際の計算上はそれ程神経質にならなくても良い。ケース数が小さい時には注意しよう。
これも,模擬データを発生させて,代表値wからの偏差平方和が最小になるのが,wが算術平均である時である事を示してみよう。
n <- 100
x <- round(rnorm(n, mean=50, sd=15), 0) # 整数データの生成
# 最小値から最大値までの範囲で0.05刻みのヴェクトルを準備
w <- seq(min(x), max(x), by=.05)
L <- length(w) # ヴェクトルwの長さ(=数字の個数)をLとする
LS <- numeric(L) # Lだけ数字を格納出来る空のヴェクトルLSを準備する
# iからLまでiの値を一つずつ増やしつつ,あるwからの偏差の平方の合計(誤差平方和)を格納
for (i in 1:L) {
LS[i] <- sum((x - w[i])^2)
}
min(LS) # 誤差平方和の最小値(=最小二乗和)
[1] 19433.65
w[LS == min(LS)] # 誤差平方和が最小になる時のwの値
[1] 50.05
mean(x) # 実際のxの算術平均
[1] 50.06
var(x)*(n-1) # 実際のxの偏差平方和
[1] 19433.64
スクリプトを解読する事は未だ難しいかも知れないが,最終的に誤差平方和の最小値,その時のwの値,実際の算術平均と偏差平方和を見比べて,一致している事を確認して欲しい。
Rで実際の調査データを扱い始めた途端に困るのが,実際のデータには欠損値(missing)が含まれている事である。Rの入門的解説では,取り敢えず欠損値を度外視して説明している場合が多く,其の儘では使えなくなって困ってしまう事が多い。かなり面倒には感じるが,学習者は常に欠損値について処置しなければならない場合が多いと云う事を念頭においておこう。集計や分析をしていて「何か変だな?」と感じたら必ず「NAはきちんと扱われているかな?」と気にかけよう。
以下に,スクリプトに#コメントを付けたので,自分でも実際に試しながらよく理解しよう。結果は二つに分けて掲示する。
x <- c(1, 2, 3, 4, 5, 6, 7, NA, 8, 9, 10)
mean(x); median(x); sd(x); sum(x) # いずれもNAは不可
[1] NA
[1] NA
[1] NA
[1] NA
summary(x) # NAを含んでも可
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1.00 3.25 5.50 5.50 7.75 10.00 1
length(x) # NAも含んだ数字の個数を出力する
[1] 11
mean(x, na.rm=TRUE) # NA を除外して算術平均を計算。 NA を remove する。
[1] 5.5
median(x, na.rm=T) # TRUE は T と略する事が出来る。
[1] 5.5
var(x, na.rm=T); sd(x, na.rm=T)
[1] 9.166667
[1] 3.02765
table(x) # NAは自動的に除外
x
1 2 3 4 5 6 7 8 9 10
1 1 1 1 1 1 1 1 1 1
table(x, useNA="ifany") # NAがあるかどうかに注意
x
1 2 3 4 5 6 7 8 9 10 <NA>
1 1 1 1 1 1 1 1 1 1 1
y <- c(1, 1, 1, 2, 2, 3, 3, NA)
z <- c(7, NA, 5, 6, 7, 6, 5, 7)
table(y, z) # いずれか一方にでも NA があると除外
z
y 5 6 7
1 1 0 1
2 0 1 1
3 1 1 0
table(y, z, useNA="ifany") # 「もしNAがあれば含めよ」とのオプション
z
y 5 6 7 <NA>
1 1 0 1 1
2 0 1 1 0
3 1 1 0 0
<NA> 0 0 1 0
y1 <- y[complete.cases(y, z)] # いずれもNAでないものだけ
z1 <- z[complete.cases(y, z)]
y1; z1
[1] 1 1 2 2 3 3
[1] 7 5 6 7 6 5
table(y1, z1, useNA="always") # NAがあろうがなかろうが表示
z1
y1 5 6 7 <NA>
1 1 0 1 0
2 0 1 1 0
3 1 1 0 0
<NA> 0 0 0 0
特に table( )関数の様に,自動で削除する関数の場合は,ユーザはNAの事を忘れがちになるのでしばしば注意が必要である。
2変数の分割表(クロス表)では,いずれか一方にNAが含まれるケースはデフォルトでは全て除外される事,useNA=“always”を付けると実際にはNAが無くても分割表にNAの欄が明示される事を確認しよう。
最初に度数分布表や分割表を出力する場合には必ずuseNA=“always”またはuseNA=“ifany”として確認する習慣をつけておくのが良い。
pnorm(1.00, mean=0, sd=1) - pnorm(-1.00, mean=0, sd=1)
[1] 0.6826895
qnorm(.90, mean=0, sd=1)
[1] 1.281552
1-pnorm(70, mean=50, sd=10)
[1] 0.02275013
download.file("http://sgn.sakura.ne.jp/text/practice.csv", "practice.csv")
data01 <- read.csv("practice.csv")
summary(data01$income)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 75.0 300.0 345.1 500.0 2500.0 42
table(data01$income, useNA = "ifany")
0 25 75 125 200 300 400 500 600 700 800 925 1125 1375 1750 2500
45 24 28 29 40 35 31 29 27 18 12 13 6 1 3 1
<NA>
42
t01 <- table(data01$income)
t01[t01 == max(t01)]
0
45
x <- data01$income
x <- x[!is.na(x)]
n <- length(x)
sqrt(sum((x-mean(x))^2)/n)
[1] 338.2096
sd(x)*sqrt((n-1)/n)
[1] 338.2096
\[\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2=\frac{1}{n} \sum_{i=1}^n (x_i^2 - 2\bar{x}\cdot{}x_i+\bar{x}^2)\] \[=\frac{1}{n} \sum_{i=1}^n x_i^2 - 2\frac{1}{n} \sum_{i=1}^n\bar{x}\cdot{}x_i+\frac{1}{n} \sum_{i=1}^n\bar{x}^2\] \[=\frac{1}{n} \sum_{i=1}^n x_i^2 - 2\bar{x}\frac{1}{n}\sum_{i=1}^nx_i+\frac{1}{n}\cdot{}n\bar{x}^2\] \[=\frac{1}{n} \sum_{i=1}^n x_i^2 - 2\bar{x}\cdot{}\bar{x}+\bar{x}^2=\frac{1}{n} \sum_{i=1}^n x_i^2 - \bar{x}^2\]
\(\Sigma\)の中にある\(\bar{x}\)は添え字の\(i\)に関係なく一定の値なので定数扱いになり,\(\Sigma\)の外に出せる。定数を\(1\)から\(n\)まで足すと,元の定数の\(n\)倍になる。