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を絶対零度とも呼ぶが，これは古典力学的には原子が運動を停止している状態を指し，実体的な意味を持つ。そして，通常の意味においては，負の絶対温度は存在しないとされる3。
間隔尺度か比率尺度かを簡単に見分ける一つのポイントは，自然に・素直にマイナスの値を考える事が出来るかどうか，と言えよう。長さのメートル，質量のグラムなどは，0というのは「そもそも長さ／重さが存在しない」と云う特別な意味を持ち，「或物の長さがマイナス1メートルである」と云った表現は，素直な・自然な理解が出来ないものである。こうした長さや重さの単位は比例尺度である。
間隔尺度と順序尺度
もう一つ区別が難しい，と云うか特に社会調査において厳密には問題となるのが「順序尺度／間隔尺度」である。
社会調査・世論調査では，例えば，「『夫は外で働き，妻は家庭を守るべきである』という考え方について賛成ですか，反対ですか」との質問に対して，「1 賛成，2 どちらかといえば賛成，3 どちらともいえない，4 どちらかといえば反対，5 反対」の選択肢から一つを選んで貰う回答方式（5件法）がよくある。
この回答を，1～5（もしくは逆順に5～1）の数字に変換して，平均を計算したり，平均の男女差や世代差を計算したりする事がよくあるが，平均を計算すると云う事はこれを，単なる順序尺度ではなく間隔尺度と見做している事になる。間隔尺度だと見做すと云う事は，「1 賛成」と「2 どちらかといえば賛成」の差は，「2 どちらかといえば賛成」と「3 どちらともいえない」の差に等しく，「3 どちらともいえない」と「5 反対」の差の2倍であると考える事を意味する。
こうした計算に対して，「この回答・選択肢は，選択肢間の距離（賛否の強度の違い？）が等間隔であると前提する事は出来ず，順序尺度でしかない」とする批判もある。そうであれば，例えば，その1～5の値をそのまま重回帰分析（第9章）の従属変数に使うと云った事は誤りである事になる。
回答選択肢としてはしばしば「どちらともいえない」が存在しない4件法も用いられる事を考えると，5件法の場合には1～5の値で等間隔と見做し，4件法の場合には1～4の値で等間隔と見做す事は確かに問題かも知れない。「どちらかといえば賛成」と「どちらかといえば反対」の間隔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_{i\;\;\;\;\;\;\;\;\;\;i=1,2,3,...,n} \]
x [1] NA 49 55 51 59 46 56 36 50 44 44 72 51 33 52 41 46 42 60 45 67 NA 41 45 68
[26] 42 45 59 32 52 58 63 61 62 37 62 42 57 50 38 47 50 NA 46 44 50 61 40 50 44
[51] 50 47 78 45 48 45 46 48 42 34 56 43 39 NA 47 44 42 70 78 37 45 26 31 37 52
[76] 30 58 54 53 53 56 47 56 48 NA 50 64 53 62 46 50 56 50 65 40 55 45 49 57 49
[101] 29 32 48 38 46 NA 46 51 51 44 45 36 55 53 50 52 37 35 47 54 47 51 58 36 48NAが含まれていると面倒なので，NAを除外したものを新たに x1 としておく7。
x1 <- x[!is.na(x)]算術平均 \[ \bar{x}=\color{red}\frac{1}{n}\sum_{i=1}^{n}\color{black}x_{i} \]
mean(x1)[1] 48.82353
不偏分散(unbiased variance)の式
推測統計において母分散の推定値として用いるのは上記の分散ではなく，以下の不偏分散である事が多い。 \[ \hat{\sigma}_{x}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^2 \]
Rの分散の関数 var( ) の結果は，上記の①，②とは一致しない.
var(x1)[1] 95.46859不偏分散の定義式に従った以下の計算と一致する事を確認しよう。
sum((x1 - mean(x1))^2)/(length(x1) - 1)[1] 95.468591-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] 9.729663
統計ソフトの出力する標準偏差
sd(x1) # Rの標準偏差の関数．上とは一致しない．[1] 9.770803sqrt(sum((x1 - mean(x1))^2)/(length(x1) - 1)) # 不偏分散の平方根を定義に従って計算[1] 9.7708032 1変量の全体像の把握
2-1 カテゴリカル変数の度数分布表
今度は，カテゴリカル変数らしい模擬データを生成してみよう8。
var0 <- sample(c(1, 2, 3, 4, NA),
size = 384,
replace = T,
prob = c(.474, .240, .099, .177, .010))本文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 正規雇用 178 178 46.4 46.4 46.6 46.6
2 非正規雇用 95 273 24.7 71.1 24.9 71.5
3 自営 28 301 7.3 78.4 7.3 78.8
4 無職・学生 81 382 21.1 99.5 21.2 100.0
NA 欠損値 2 384 0.5 100.0 0.0 0.0 二つの dimnames( )関数は，値ラベルや表頭の列名を付けるものである。
t0からcp1までパーツを準備し，それを「列結合(cbind)」でひとまとめにしたのが ft0 である。
上記のスクリプトで，var1 <- var0 （本文で言えば var1 <- data$job）の代入式の右辺だけを変更すれば，任意の変数の度数分布表が同じスクリプトで作成できる9。
→のサイトでは，より整理して，新たな関数として作成する方法を紹介している。超初心者向けのRガイド d. 一つのカテゴリカル変数の度数分布表
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)
時系列変化などを表すには折れ線グラフを用いる。
ここでは，総務省の「労働力調査」の長期時系列データにおける「完全失業率」の年次変化をグラフ化してみる10。
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 | 0000113378999
1 | 00111223444555677779999999
2 | 01112233335555667777778888889999
3 | 0000000111222222333445555667777889
4 | 001222223334456667777899
5 | 000001122333556677777789
6 | 01112223456677899
7 | 112344567
8 | 0124557
9 | 012367
10 | 0078
11 |
12 |
13 | 5
14 | 9
15 |
16 |
17 |
18 | 8
19 | 9hist(x2,
breaks = seq(0, ceiling(max(x2)/100)*100,
by = 100))
summary(x2) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 247.8 369.0 436.5 572.2 1992.0 summary( )の結果やヒストグラムと比較して幹葉図を見ると読み取れると思うが，この例の場合には，年収の百万円以上の桁の部分がセパレータ(|)の左側（表側）に示され，その下の位である十万円の桁の数字が右側に昇順で羅列されている。ヒストグラムを時計回りに90°倒してみると幹葉図と似ているのがわかるだろう。
この様に一つ一つのケースを数字一つで表す為に，大量のデータの表示には向かない。
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
もっと様々なオプションや装飾があるが，基本形のみ紹介しておくので，興味のある人は自分で調べると良い11。
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')
g
確率分布を曲線で表現する
確率分布のグラフは曲線となっていて，それとx軸の間の面積を考える事がピンとこない人は，以下の図を見比べて理解してほしい。
グラフ作成のコード
x <- rnorm(n = (nx <- 100000), mean = 50, sd = 10) # 模擬データ
d01 <- data.frame("score" = x)
g0 <- ggplot(d01, aes(x = score)) +
labs(title = "正規分布") +
# theme_classic(base_family = "HiraginoSans-W3") + # Macだとこうした行が必要
theme(legend.position = 'none',
axis.title.y = element_text(size = 20))
g1 <- g0 + geom_histogram(position = "identity",
color = "black",
fill = "gray",
binwidth = 5,
alpha = 1) +
labs(subtitle = "縦軸が人数（度数）。合計すると総人数。")
g2 <- g0 + geom_histogram(position = "identity",
color = "black",
fill = "gray",
binwidth = 5,
alpha = 1,
aes(y = after_stat(density))) +
labs(subtitle = "縦軸が確率密度（比率や割合みたいなもの）。合計は 1.0（100%）")
g3 <- g0 + geom_histogram(position = "identity",
color = "black",
fill = "gray",
binwidth = 2,
alpha = 1,
aes(y = after_stat(density))) +
labs(subtitle = "棒の幅を狭くする（ケースが多いので可能）")
g4 <- g0 + geom_histogram(position = "identity",
color = "black",
fill = "gray",
binwidth = .5,
alpha = 1,
aes(y = after_stat(density))) +
labs(subtitle = "棒の幅を更に狭くする。棒の輪郭が密集してくる。")
g5 <- g0 + geom_histogram(position = "identity",
color = "gray",
fill = "gray",
binwidth = .5,
alpha = 1,
aes(y = after_stat(density))) +
labs(subtitle = "棒の輪郭線を省略する")
g6 <- g0 + geom_density() +
geom_hline( yintercept = 0) +
labs(subtitle = "外周の輪郭だけで描くことにする")g1
g2
g3
g4
g5
g6
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), "%"))
統計学の参考書の巻末によくあるような，標準正規分布表，t分布表，カイ二乗分布表へのリンク（このページもRで作成。）
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. :38.00 Min. :-2.70279 Min. :22.97
1st Qu.:54.00 1st Qu.:-0.75615 1st Qu.:42.44
Median :60.50 Median : 0.03467 Median :50.35
Mean :60.22 Mean : 0.00000 Mean :50.00
3rd Qu.:66.00 3rd Qu.: 0.70383 3rd Qu.:57.04
Max. :80.00 Max. : 2.40715 Max. :74.07 次に，それぞれの要約統計量を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] 4856.755sum((x - mean(x))^2)/n # 定義式に従った標本分散[1] 97.13511mean((x - mean(x))^2) # 偏差平方の平均＝標本分散[1] 97.13511sqrt(mean((x - mean(x))^2)) # 標本標準偏差[1] 9.855715# Rの分散var( )，標準偏差sd( )
var(x); sd(x)[1] 99.11746[1] 9.955775# 定義式に従った不偏分散
sum((x - mean(x))^2)/(n-1)[1] 99.11746 定義式に従った標本分散は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] 24537.6w[LS == min(LS)] # 誤差平方和が最小になる時のwの値[1] 51.4mean(x) # 実際のxの算術平均[1] 51.38var(x)*(n-1) # 実際のxの偏差平方和[1] 24537.56スクリプトを解読する事は未だ難しいかも知れないが，最終的に誤差平方和の最小値，その時の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\)倍になる。
Footnotes
「質的／量的」の二分法を使用しない理由については『入門・社会調査法〔第3版〕』第2章を参照して欲しい。↩︎
論理値の FALSE の 省略形 F と重複してしまう。↩︎
現在の統計力学では，例外的に「負の絶対温度」と云う概念が存在するらしいが↩︎
すなわち，性別役割分業に賛成する態度の強さの違い↩︎
NAを含む場合は以下の「 NAを含むデータの場合はこうせよ！」を参照。↩︎
【注意！】 以下ののコマンドの2行目は，テクスト本文で見ると長い等号が一つ（“＝”）に見えますが，半角の等号が二つ（“==”）です。御注意下さい。↩︎
定義式は一般的なxで表記する。↩︎
本文2-1の表と全く同じとはならないが，よく似た変数になるように工夫している。↩︎
上の例は「従業上の地位」で厳密には名義尺度なので，「累積」情報は意味がない事が多い点に注意。この例ではたまたま，「被雇用者」「就労者」と云った意味で累積の列も解釈できる。↩︎
長期時系列表２ 就業状態別15歳以上人口－全国，列M： 完全失業率（％）↩︎
オプションが多過ぎて途方に暮れるかもしれないが。↩︎