C <- 20.0 # 任意の摂氏温度を入力
Fd <- 1.8*C + 32; names(Fd) <- c("Fahrenheit"); names(C) <- c("Celsius")
C; FdCelsius
20 Fahrenheit
68 『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第1章 1変量の記述統計の基礎
0 全体の章構成
- 第1章 1変量の記述統計の基礎
- 第2章 2変数の関連の記述統計
- 第3章 推測統計の基礎
- 第4章 統計的検定の一般型
- 第5章 2変数の関連の推定と検定
- 第6章 2群の母平均の差のt検定
- 第7章 平均値の差の分散分析(ANOVA)
- 第8章 単回帰分析(SRA)
- 第9章 重回帰分析（Ⅰ）(MRA)
- 第10章 重回帰分析（Ⅱ）(MRA)
- 第11章 主成分分析と因子分析
- 第12章 ロジスティック回帰分析
- 第13章 階層構造のあるデータの分析
統計学の見取り図
1 一変量についての要約統計量
データファイルのイメージ
- ケース case, unit, element, observation
- 統計的に集計する時の対象の単位。世論調査なら個人(respondent)，都道府県データなら都道府県。国際比較なら国。
- 変数 variable
- 単位によってその内容が異なり得る性質・属性・特性(property)。年齢，性別，人種・エスニシティ，教育程度，職業，婚姻上の地位，収入など。
- 値（変数値） value
- ある単位におけるある変数の性質・属性・特性の内容(property)。男性，22歳，四年制大学卒，未婚など。
1-1 尺度水準と変数の分類法
平たく言えば「変数(variable)の種類」に相当する尺度水準(level of scales, scale level)は，測定水準(測定のレヴェル，levels of measurement)や測定尺度(measurement scales, scales of measurement)などとも言われ，level, scale, measurementの3つの語がやや入り組んでいて必ずしも厳密に使い分けられているとは言えない。
よく見られる変数の分類法と測定尺度の関係を表に纏めてみる1。表の上段ほど尺度水準が低く，言い換えれば含まれる情報量が少なく，表の下段ほど水準が高い，つまり含まれる情報量が多い。
| 4尺度水準 | 変数分類A | 変数分類A’ | 変数分類B | 本書での提案 | ||
|
名義尺度 nominal |
質的変数 qualitative |
離散変数 | カテゴリカル変数 |
カテゴリカル変数 categorical |
名義尺度 | |
|
順序尺度 ordinal |
順序尺度 | |||||
|
間隔尺度 interval |
離散尺度 discrete |
量的変数 quantitative |
連続変数 | 連続変数 |
数量変数 numerical |
|
|
比例尺度 ratio |
連続尺度 continuous |
|||||
区別が難しい場合があるのは，「順序尺度／間隔尺度」と「間隔尺度／比例尺度」であるが，後者の「間隔尺度／比例尺度」については，本書のレヴェルでは無理して区別しなくても支障ない。
| 尺度水準 | 例 | 可能な演算 |
|---|---|---|
| 名義(nominal) | 性別，居住地，職業カテゴリ | 最頻値(mode) |
| 順序(ordinal) | 順位，意識調査項目への5件法回答 | 中央値(median) |
| 間隔(interval) | 試験の点数，西暦年号 | 算術平均(arithmetic mean)，加・減 |
| 比例(ratio) | 時間，距離，金額 | 幾何平均(geometric mean)，乗・除 |
interval scaleは「間隔尺度」以外に「距離尺度」，ratio scaleは「比例尺度」「比率尺度」「比尺度」の訳し方がある。
間隔尺度と比例尺度
「間隔尺度／比例尺度」の区別の例としてよく見るのは，温度の単位の「摂氏温度℃」と「華氏温度℉」である。日本で20℃と表現する気温は，アメリカ合衆国では68℉と表現される。
摂氏温度と華氏温度の間には単純な関係があり，摂氏温度を1.8倍して32を加えれば華氏温度に変換できる。以下のRスクリプトで色々と試したりグラフ化したりすると良い。
摂氏温度(Celsius)から華氏温度(Fahrenheit)に変換
テクスト本文では華氏温度の変数名に F を使用しているが，後に不都合が生じる2ので Fd としておく。
華氏温度(Fahrenheit)から摂氏温度(Celsius)に変換
Fd <- 100.0 # 任意の華氏温度を入力
C <- 1/1.8*(Fd - 32); names(Fd) <- c("Fahrenheit"); names(C) <- c("Celsius")
Fd; CFahrenheit
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）などのカテゴリカルデータ分析法が発達しているので，順序尺度のままで分析する事も十分に可能である。
1-2 代表値（中心傾向; central tendency）
本文で「xを変数のヴェクトルだとするとき」と書いているが，実際に演習できるように，模擬データを発生させてみよう。
このxに対して，データの個数，中央値，算術平均，そして最頻値を求める。（最頻値については若干の工夫が必要になる。）
【注意！】 最後の一行は，テクスト本文で見ると長い等号が一つ（“＝”）に見えますが，半角の等号が二つ（“==”）です。御注意下さい。
本文同様，幾何平均，調和平均，トリム平均も求めてみよう。
中央値と算術平均
1-3 偏差平方和と分散
分散を求めるには var(x) とすればよい。しかし分散が何であるのかをより深く理解する為に，本文同様，var( )関数を用いずに分散を求めてみよう。
偏差平方和(SS; Sum of Squares)の定義式
\[ SS(Sum\;of\;Squares)=\sum_{i=1}^{n}(x_{i}-\bar{x})^2 \]
分散(variance)の定義式
\[ 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}\color{black}(x_{i}-\bar{x})^2 \]
x # まずは生成した変数（ヴェクトル）そのものを表示させてみよう。 [1] NA 55 48 47 58 49 42 54 69 51 45 50 52 50 73 63 52 39 28 46 54 NA 61 53 38
[26] 58 45 53 61 29 50 42 50 52 41 29 40 41 58 59 39 57 NA 51 50 54 47 49 61 60
[51] 45 46 52 50 53 33 46 47 56 51 59 47 53 NA 51 41 47 60 52 64 48 52 65 53 53
[76] 44 42 49 53 54 42 59 63 43 NA 61 58 62 38 54 37 51 58 46 30 69 53 43 54 48
[101] 53 55 48 53 42 NA 45 44 61 59 52 32 64 60 40 53 31 49 41 56 44 41 36 58 55x1 <- x[!is.na(x)] # NAが含まれていると面倒なので，NAを除外したものを新たに x1 とする
mean(x1) # 算術平均[1] 50.08403mean(abs(x1 - mean(x1))) # 平均偏差[1] 6.988066sum((x1 - mean(x1))^2) # 偏差平方和 SS[1] 9381.16sum((x1 - mean(x1))^2)/length(x1) # 分散の求め方①[1] 78.83327mean((x1 - mean(x1))^2) # 分散の求め方②[1] 78.83327不偏分散(unbiased variance)の式
\[ \hat{\sigma}_{x}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2 \]
var(x1) # Rの分散の関数．上の①や②とは一致しない.[1] 79.50135sum((x1 - mean(x1))^2)/(length(x1) - 1) # 不偏分散[1] 79.501351-4 標準偏差
標準偏差(standard deviation)の定義式
\[sd(x)=s_{x}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})^2}\]
sqrt(mean((x1 - mean(x1))^2)) # 標準偏差[1] 8.87881sd(x1) # Rの標準偏差の関数．上とは一致しない．[1] 8.916353sqrt(sum((x1 - mean(x1))^2)/(length(x1) - 1)) # 不偏分散の平方根[1] 8.9163532 1変量の全体像の把握
2-1 カテゴリカル変数の度数分布表
今度は，カテゴリカル変数らしい模擬データを生成してみよう。本文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
185 102 38 55 prop.table(table(var0)) # 比率に変換var0
1 2 3 4
0.4868421 0.2684211 0.1000000 0.1447368 cumsum(table(var0)) # 累積度数 1 2 3 4
185 287 325 380 本文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 正規雇用 185 185 48.2 48.2 48.7 48.7
2 非正規雇用 102 287 26.6 74.7 26.8 75.5
3 自営 38 325 9.9 84.6 10.0 85.5
4 無職・学生 55 380 14.3 99.0 14.5 100.0
NA 欠損値 4 384 1.0 100.0 0.0 0.0 二つの dimnames( )関数は，値ラベルや表頭の列名を付けるものである。
t0からcp1までパーツを準備し，それを「列結合(cbind)」でひとまとめにしたのが ft0 である。
上記のスクリプトで，var1 <- var0 （本文で言えば var1 <- data$job）の代入式の右辺だけを変更すれば，任意の変数の度数分布表が同じスクリプトで作成できる。
2-2 1変量のグラフ表現
模擬データを生成して，ヒストグラムや箱ひげ図で図示してみよう。hist( )もboxplot( )もNAを除外して作図してくれるので，NAを含まない単純な模擬データで演習する。
2-2-1 ヒストグラム(histogram)
一つずつ試してみると，主なオプションの働きが分かるだろう。ここでは煩雑になるので，主な描画結果のみ掲載する。それ以外は自分でスクリプトを実行して確認して欲しい。
x <- rnorm(n = 400, mean = 50, sd = 10) # 模擬データhist(x)
hist(x,
xlim = c(10, 90),
freq = FALSE,
main = "変数xのヒストグラム",
xlab = "変数の値",
ylab = "密度",
breaks = seq(0, 100, by = 5))
hist(x,
xlim = c(10, 90),
freq = FALSE,
main = "変数xのヒストグラム",
xlab = "変数の値",
ylab = "密度",
breaks = c(0, 20, 30, 40, 45, 50, 55, 60, 70, 80, 100))
2-2-2 箱ひげ図(box-and-whisker plot)
これもそれぞれのオプションの意味を理解しよう。
x <- round(rnorm(n = 200, mean = 40, sd = 10), 0)
y <- round(rnorm(n = 50, mean = 50, sd = 15), 0) # ケース数と散布度の異なるヴェクトルを用意boxplot(x)
boxplot(x, y, # 2つのヴェクトルを並べて図示
ylim = c(0, 100), # y軸の範囲を指定
varwidth = T, # 幅がサンプルサイズの平方根に比例
main = "xとyの箱ひげ図",
ylab = "変数の値", # タイトルとy軸ラベル
names = c("ヴェクトルx", "ヴェクトルy"), # x軸にラベル
col = c("cyan", "pink"),
border = c("blue", "red")) # 色を付ける
2-2-3 折れ線グラフ(line graph)，幹葉図(stem and leaf plot)
もう少し他の作図法も紹介しておこう。
2-2-4 折れ線グラフ(line graph)
時系列変化などを表すには折れ線グラフを用いる。
ここでは，総務省の「労働力調査」の長期時系列データにおける「完全失業率」の年次変化をグラフ化してみる3。
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,
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)
plot(year, unemployment,
bty = "n",
type = "b",
pch = 16,
col = marker_col)
2-2-5 幹葉図(stem and leaf plot)
もう一つは幹葉図と呼ばれるものであるが，これは，テクスト（文字）のみで表現するのでグラフ機能を必要としない割には詳しい情報を表示すると云う利点がある。
他方で，データの数が多いと有効でなくなると云う欠点がある。
ここでは，第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人分を無作為抽出（単位： 万円）
x2 <- sample(data02$income,
size = 200,
replace = FALSE)stem(x2)
The decimal point is 2 digit(s) to the right of the |
0 | 000000000044445777999
1 | 00011114444455567778
2 | 00111222223445566666777777888899
3 | 0001122333444555566677799
4 | 000011111222334445578899
5 | 001111233444445566777779
6 | 0556777777899
7 | 012224445668
8 | 012468889
9 | 4455566
10 |
11 | 0125
12 | 5
13 |
14 | 1267
15 | 2
16 | 6
17 |
18 | 1
19 | 3hist(x2,
breaks = seq(0, ceiling(max(x2)/100)*100,
by = 100))
summary(x2) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 217.5 398.0 467.0 658.5 1931.0 summary( )の結果やヒストグラムと比較して幹葉図を見ると読み取れると思うが，この例の場合には，年収の百万円以上の桁の部分がセパレータ(|)の左側（表側）に示され，その下の位である十万円の桁の数字が右側に昇順で羅列されている。
この様に一つ一つのケースを数字一つで表す為に，大量のデータの表示には向かない。
2-2-6 ggplot2
テクストにはないが，グラフ描画に定評の高い ggplot2 パッケイジによる描画例を示しておく。
2-2-6-1 箱ひげ図
上のカラーの「xとyの箱ひげ図」と同様のグラフを，ggplot2 で描画して見る。
ggplot2を使うには，使用する変数はデータフレイムになっている必要があるので，それをgdataとして作成しておく。
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 + theme_classic()
2-2-6-2 ヒストグラム
上のヒストグラムで使用した data02 のデータフレイムを再利用する。
g20 <- ggplot(data02, aes(x = income))
g21 <- g20 + geom_histogram(alpha = .2)
g21`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
もっと様々なオプションや装飾があるが，基本形のみ紹介しておくので，興味のある人は自分で調べると良い。
（オプションが多過ぎて途方に暮れるかもしれないが。)
3 正規分布と標準化
正規分布\(N(\color{blue}\mu\color{black}, \color{blue}\sigma\color{black}^2)\)の定義式
\[ y=\frac{1}{\sqrt{2\pi\color{blue}\sigma\color{black}^2}}e^{-\frac{(\color{red}x\color{black}-\color{blue}\mu\color{black})^2}{2\color{blue}\sigma\color{black}^2}} \]
ggplot2で正規分布を描画
x <- rnorm(n = (nx <- 1000), mean = 50, sd = 10) # 模擬データ
y <- rnorm(n = (ny <- 1000), mean = 100, sd = 20) # 模擬データ
d01x <- data.frame("group" = rep("x", nx), "score" = x)
d01y <- data.frame("group" = rep("y", ny), "score" = y)
d01 <- rbind(d01x, d01y)
g <- ggplot(d01, aes(x = score, fill = group))
g <- g + geom_histogram(position = "identity", alpha = 0.7)
g <- g + labs(title = "色々な正規分布") +
theme(legend.position = 'none')
3-1 正規分布(normal distribution)と標準正規分布(SND)
理論的な正規分布のグラフ描画については〈このページ〉で紹介しているので参照して欲しい。
以下にその一部を微修正して紹介する。是非解読して試してほしい。
標準正規分布において，面積(%)を指定すると，対応する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),
xlim = c(-3.5, 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))*2curve(dnorm(x, mean = 0, sd = 1),
xlim = c(-4, 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), "%"))
3-2 変数の標準化(standardization)
標準化\(standardization\)
\[ \color{red}z_{i} \color{black}=\frac{\color{red}x_{i}\color{black}-\color{blue}\bar{x}\color{black}}{\color{blue}s_{x}\color{black}},\;\;(i=1,2,3,...,n)\;\;\color{green}\Rightarrow\color{black}\;\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 # 標準化スコアから，「偏差値」を換算スクリプトに#でコメントを付けてある様に，0点から100点の間の点数をn人分発生させ，それを標準化したzを求め，更にそこからいわゆる「偏差値」に換算したZを求めている。
summary(cbind("素点"=score, "標準化スコア"=z, "「偏差値」"=Z)) # 素点 素点 標準化スコア 「偏差値」
Min. :39.00 Min. :-2.7086 Min. :22.91
1st Qu.:55.00 1st Qu.:-0.6663 1st Qu.:43.34
Median :59.50 Median :-0.0919 Median :49.08
Mean :60.22 Mean : 0.0000 Mean :50.00
3rd Qu.:64.25 3rd Qu.: 0.5144 3rd Qu.:55.14
Max. :87.00 Max. : 3.4182 Max. :84.18 次に，それぞれの要約統計量をsummary( )関数で表示させている。
以下では，それぞれのヒストグラム3つを描いている。
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に，分散と標準偏差が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 = "")
発展1 平均と分散の理解を広げる
発展1-1 分散と不偏分散(unbiased variance)
母分散
\[
\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] 4850.759sum((x - mean(x))^2)/n # 定義式に従った標本分散[1] 97.01518mean((x - mean(x))^2) # 偏差平方の平均＝標本分散[1] 97.01518sqrt(mean((x - mean(x))^2)) # 標本標準偏差[1] 9.849628# Rの分散var( )，標準偏差sd( )
var(x); sd(x)[1] 98.99508[1] 9.949627# 定義式に従った不偏分散
sum((x - mean(x))^2)/(n-1)[1] 98.99508 定義式に従った標本分散はRの関数var( )の結果とは一致しない事，nではなくn-1で割った不偏分散がvar( )の結果と一致する事，当然ながら不偏分散の方が標本分散よりも少し大きな値になる事，を確認しよう。
n=50の場合は，50/49=1.020408で，不偏分散の方が2%程度大きくなるが，n=500になるとこれが0.2%の違いでしかなくなるので，実際の計算上はそれ程神経質にならなくても良い。ケース数が小さい時には注意しよう。
発展1-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] 22189w[LS == min(LS)] # 誤差平方和が最小になる時のwの値[1] 49mean(x) # 実際のxの算術平均[1] 49.01var(x)*(n-1) # 実際のxの偏差平方和[1] 22188.99スクリプトを解読する事は未だ難しいかも知れないが，最終的に誤差平方和の最小値，その時のwの値，実際の算術平均と偏差平方和を見比べて，一致している事を確認して欲しい。
発展1-3 Rにおける欠損値(missing value)
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] NAsummary(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] 11mean(x, na.rm=TRUE) # NA を除外して算術平均を計算。 NA を remove する。[1] 5.5median(x, na.rm=T) # TRUE は T と略する事が出来る。[1] 5.5var(x, na.rm=T); sd(x, na.rm=T)[1] 9.166667[1] 3.02765table(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 0table(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 0y1 <- 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 5table(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”として確認する習慣をつけておくのが良い。
第1章の【練習問題】の解答
1) の答え
pnorm(1.00, mean=0, sd=1) - pnorm(-1.00, mean=0, sd=1)[1] 0.68268952) の答え
qnorm(.90, mean=0, sd=1)[1] 1.2815523) の答え
1-pnorm(70, mean=50, sd=10)[1] 0.022750134) の答え
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.2096sd(x)*sqrt((n-1)/n)[1] 338.20965) の答え
\[\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\)倍になる。