記述統計量と推測統計量

1 記述統計量

あるデータセットの特徴を把握するときに用いられる統計量。元のデータから算出される。なお、ここではデータとは数値のみで構成されているものとする。

以下では実際に次のデータ【04_example.csv】を使って例を示していく。このデータを右クリックメニューからダウンロードしたうえで，自分自身のPosit CloudのProjectに前回の要領でアップロードしておくこと．

Shohoku <- read.csv("practice/04_example_Shohoku.csv", header = TRUE, fileEncoding="UTF8")
Shohoku

##       nam number score
## 1    赤木      4     7
## 2    桜木     10     4
## 3    三井     13    10
## 4    宮城      7     7
## 5    流川     11     9
## 6    木暮      5     6
## 7    潮崎      8     6
## 8    安田      6     7
## 9    角田      9     7
## 10   石井     12     5
## 11 佐々岡     13     5
## 12   桑田     15     4

CSVファイルの文字コード（fileEncording）について

なお，ダウンロードしたcsvファイルはUTF-8という形式で保存したファイルであり，そのまま読み込めば問題なくRで読み込めるが，このファイルをExcelで開いた場合には，文字化けが起こることがある． このファイルをExcelで開きたい場合には，一旦このファイルを右クリックメニューからメモ帳で開き，「名前を付けて保存」で「文字コード」を「ANSI」に変更して保存し直すことで，Excelで開いても文字化けが起こらないようになる． ただし，今度はRで読み込んだ場合に文字化けが起こる可能性がある．その場合には，Rで読み込む際にfileEncoding="SJIS"というオプションを追加して読み込むことで文字化けを防ぐことができる．

以降は、前々回紹介したRMarkdownを作成し、上記の内容を含めて自分で入力・実行して出力結果を確認しながら進めていくこと。

1.1 代表値（平均、中央値、最頻値）

あるデータセットの全体的な特徴を一番端的に表現する値であり、文字通りデータセットを「代表」する値。平均、中央値、最頻値の3種類がある。

1.1.1 平均

すべてのデータ（数値）を足し合わせて、データの数で割った値のこと。

\[ \bar{x} =\frac{1}{n}\sum^{n}_{i=1}x_{i} \] Rで求めるにはmean()関数を用いる。

example<-c(1,2,3,4,5) # 例として1,2,3,4,5という5つの数値からなるベクトルを使った
mean(example)

## [1] 3

mean(Shohoku$score) #先ほどの湘北高校バスケ部のフリースロー10本勝負のスコアの平均

## [1] 6.416667

なお、実際のデータを扱う際には、CSVファイルやExcelファイルからデータを読み込むことが多くなる。そうしたデータファイルにはたまに「空白」が含まれているケースがある（例えばアンケートで無回答だったり、実験で計測に失敗したりした場合）。その場合、Rにそれらのデータを読み込むと、空白だった箇所にはNA(Not Available：利用不可能という意味)という記号が自動的に入力される。

mean()関数を用いる際に、与えたベクトルにNAが含まれていた場合、Errorにはならないが出力はNAとなる。

このようなデータの場合には、以下のようにna.rmオプションをTrue（真）に設定する。するとNAを除いた上で平均が算出される（当然、データの総数はNAの数の分だけ減ることになる）。

example<- c(1,2,3,4,NA)
mean(example) #NAが含まれているのでNAが出力される

## [1] NA

example<- c(1,2,3,4,NA)
mean(example, na.rm = T) 

## [1] 2.5

1.1.2 中央値

全てのデータ（数値）を数の大小で順に並べ、ちょうど中央の順位に来る数値。データ点数が偶数だった場合には、中央の位置の前後の数の平均値を中央値とする。

Rではmedian()を用いる

median(Shohoku$score)

## [1] 6.5

先の湘北高校のデータは12人分なので、中央値は6位と7位の間をとることになるため、小数の値が出てきている。

1.1.3 最頻値

全てのデータ（数値）について、度数分布をとり、最も度数が大きい（出現頻度が最も多い）数値。

Rでは最頻値を直接求める関数は存在しない。そのため、まずはtable()コマンドを使って、対象となるデータの度数分布表を求める。 次いで、その度数分布表で度数が最も高いものを求める、という手続きで算出する。 ただし、度数分布表で最も度数が高いものを取り出したときのデータは名前付きベクトルで出力され、最頻値を表す値は度数値に対する「名前」（すなわち文字列）として出力れる。そこでその名前を数値として読みなおす必要がある。 以下は例である。行っている一連の処理をよく確認しておくこと。

table_data<-table(Shohoku$score) # 度数分布表を作成
print(table_data)

## 
##  4  5  6  7  9 10 
##  2  2  2  4  1  1

d<-which.max(table_data) ##度数分布表から最も度数が高い値を取り出す
print(d) # 7と4が表示されている。わかりにくいが、これは"7"という名前の付いた列に4という度数データが入っている、ということを表している。ここで"7"は数値ではなく文字列である。

## 7 
## 4

mode<-as.numeric(names(d)) #欲しいのは値としての7なので、names()を掛けた後に、as.numeric()を行っている。
print(mode)

## [1] 7

ちなみに、上記のスクリプトファイルは細かく一つ一つ分けて書いているが、以下のように関数の中に関数を入れることで、シンプルに書くことができる。このように関数が入れ子になっている場合、内側の関数から順次実行される。

as.numeric(names(which.max(table(Shohoku$score))))

## [1] 7

さらには、tidyversというライブラリを読み込んでおくことによって、以下のような書き方もできる。なお、初めてTidyverseを使う場合には、library(tidyverse)を読み込んだときにパッケージのインストールをするかどうか聞かれるので、installを選択しておく。tidyverseについては後ほど別章で詳しく説明する。

library(tidyverse)
Shohoku$score %>%
  table() %>%
  which.max() %>%
  names() %>%
  as.numeric() 

## [1] 7

1.1.4 平均・中央値・最頻値の関係

正規分布のように左右対象の一山の度数分布をしていた場合には、平均値・中央値・最頻値は全てほぼ一致する。

左右対称の一山の分布(正規分布)

一方、次の図のように左右非対称であった場合には、平均値・中央値・最頻値は一致しなくなる。

左右非対称の分布

1.2 ばらつき（分散・標準偏差、四分位数）

あるデータセットがどの程度ばらつき（広がり）を持っているのかを示す指標。平均系の平均偏差・分散・標準偏差と、中央値系の四分位数とがある。

1.2.1 平均偏差

平均偏差とは、データ全体の平均値を予め求めておいた上で、各データと平均値との距離（差の絶対値）を全てのデータに対して求め、その平均をとったもの。

平均からの距離

\[ M.A.D =\frac{1}{n}\sum^{n}_{i=1}|x_{i} - \bar{x}| \] Rではmad()関数を用いる。

MAD <- mad(Shohoku$score)
print(MAD)

## [1] 1.4826

1.2.2 分散

平均偏差は数式のなかに絶対値記号が入ってきて、数値計算が非常に面倒になる。絶対値を加えたのは単純に\(x_{i}-\bar{x}\)を求めると\(x_{i}\)によっては負の数字になるからである。負を正にするためだけなら、\(x_{i}-\bar{x}\)を2乗するということでもよい。そこで、「差の絶対値」ではなく「差の2乗」の平均値をばらつきの指標としたものが分散(Variance)である。一般に \(S^{2}\)と表現する。

「2乗の平均」にしたことによって数値自体の絶対的な意味の解釈は難しくなる。すなわち、例えば元々は平均値からの距離が1，2，3のデータはそれぞれ等間隔で平均から離れていることになるが、これを2乗すると1,4,9となり等間隔ではなくなる。このように、平均から離れれば離れるだけより大きな値になってしまって、数値自体の絶対的な意味が分かり難くなってしまう。ただ、少なくとも2つのデータセットがある時に、どちらのデータのほうがよりばらついているかは分散の値を比較することによって知ることができる。

\[ S^{2} = \frac{1}{n}\sum^{n}_{i=1}(x_{i}-\bar{x})^2 \]

Rではvar()関数とnrow()関数を用いて以下のように書く。nrow()はデータフレームの行数、すなわちデータフレームに含まれているサンプル数を返す関数である。複雑な式に見えるが、要するに、\(\frac{データ数-1}{データ数}\)をvar()に掛けているということである。これはvar()関数は不偏分散を算出する関数であり、データそのものの分散（標本分散）を算出するためには、補正をしてやる必要があるためである。不偏分散については後ほど改めて説明する。

s2<- var(Shohoku$score)
n <- nrow(Shohoku)
VAR <- (n-1)/n*s2
print(VAR) #こちらが標本分散

## [1] 3.076389

1.2.3 標準偏差

分散の場合、先に述べた通り、距離（差）の2乗にしてしまっているので、元のデータと次元が異なってしまい、出てきた数値が何を意味しているかがにわかには分からない。そこで分散の平方根を取ることによって、元のデータと同じ次元の数値として、元のデータと比較したり平均と組み合わせられるようにしたものが標準偏差(Standard Deviation)である。一般に\(S\)と表現する。

\[ S=\sqrt{S^{2}} = \sqrt{\frac{1}{n}\sum^{n}_{i=1}(x_{i}-\bar{x})^2} \]

Rではsd()関数を用いるが、データそのものの標準偏差を得たい場合には分散の時と同様にnrow()関数を用いた補正（標準偏差の場合にはさらに平方根を求める関数sqrt()も用いる）をしてやる必要がある。

s <- sd(Shohoku$score) # こちらのsは小文字
S <- sqrt((n-1)/n)*s # nは先に求めたnrow()関数の出力。
print(S)

## [1] 1.753964

標準偏差は「ばらつき」の単位（ものさし）として用いることが一般的である。データが完全に正規分布をしていた場合、次の図のように、平均値から±1Sだけ離れたところに境界を引くと、その境界内には全データの約68.27%のデータが含まれる。同様に±2Sの範囲内であれば約95.45%のデータが、±3Sになると約99.73%のデータが含まれることになる。

標準偏差とデータの分布

データの標準化（正規化）得点

平均\(\bar{x}\)を基準（すなわち0）とし、それぞれのデータの値\(x_{i}\)と平均との距離を標準偏差\(S\)を単位として表現した値。要するに、そのデータ全体の中で、あるデータ値が相対的にどのくらいの位置にあるのかを表した値のことを標準化得点、もしくは正規化得点と呼ぶ。一般に\(z_{i}\)と表現する。さらに、こうした処理を行うことを標準化、あるいは正規化と呼ぶ。 \[ z_{i} = \frac{x_{i}-\bar{x}}{S} \]

Rで標準化を行う場合、scale()関数を用いればよい。

z <- scale(Shohoku$score)
print(z)

##             [,1]
##  [1,]  0.3184211
##  [2,] -1.3191733
##  [3,]  1.9560156
##  [4,]  0.3184211
##  [5,]  1.4101508
##  [6,] -0.2274437
##  [7,] -0.2274437
##  [8,]  0.3184211
##  [9,]  0.3184211
## [10,] -0.7733085
## [11,] -0.7733085
## [12,] -1.3191733
## attr(,"scaled:center")
## [1] 6.416667
## attr(,"scaled:scale")
## [1] 1.831955

出力結果の上12個は各データ値を標準化した値である。続いて、attr()にあるのは、それぞれ標準化に用いた平均値と標準偏差である。先ほど求めた平均、標準偏差と比べてみると、平均は一致しているが標準偏差は異なる。これは、scale()関数が内部で計算に用いている標準偏差は、あくまで補正の掛けない不偏標準偏差、すなわちsd()関数の結果をそのまま計算に用いているからである。実際に、先ほどsd()関数の結果を格納したs（小文字）をプリントさせてみると、同じ値が出力される。

print(s)

## [1] 1.831955

偏差値

標準化得点が得られれば、受験の際にお世話になった「偏差値」というものを以下の式で算出することができる。 \[ 偏差値= z_{i}*10＋50 \] 要するに、偏差値とは標準化得点を10倍したものに50を足したものである。この式から得点が平均点と一致していた場合に偏差値が50になることが理解できる。また、データが正規分布をしている場合には、先に示した標準偏差とデータの分布の図より、偏差値が60から40の範囲内に約68%のデータが含まれていることがわかる。さらに、偏差値が70を超えていた場合には、そのデータ値は全体の上位2.5%以内にある、ということになる。

1.2.4 四分位数

全てのデータ（数値）を小さいものから順に並べ、最小値(0%)、下位25%、中央値(50%)、上位25%、最大値(100%)の順位に来る数値のことをまとめて四分位数と呼ぶ。特に下位25%を第1四分位数、上位25％を第3四分位数と呼ぶ。 分散と標準偏差が平均値を代表値とした場合のばらつきの指標であるのに対して、中央値を代表値とした場合のばらつきの指標が四分位数となる。

四分位数を求めるRの基本関数はquantile()である。quantile()関数は算出対象となるデータを与えた後に、どの四分位数を求めるかを小数で与えてやる。

Q_1 <- quantile(Shohoku$score, 0.25) # 第1四分位数　下位25%の値。
Q_2 <- quantile(Shohoku$score, 0.75) # 第3四分位数　上位25%、すなわち下から75%の値
print(Q_1)
print(Q_2)

## 25% 
##   5 
## 75% 
##   7

このように関数を用いるときに()に与えるデータやオプションのことを引数（ひきすう）と呼ぶ。

また、quantile()は以下のように2つ目の引数にc()を用いて複数のパラメータを指定してやることで、複数の四分位数を一度に求めることもできる。

Q <- quantile(Shohoku$score, c(0.0, 0.25, 0.50, 0.75, 1.0))
print(Q)

##   0%  25%  50%  75% 100% 
##  4.0  5.0  6.5  7.0 10.0

ちなみに、quantile()の出力は名前付きベクトルである。

箱ひげ図

四分位数はその数字そのものを示すより、一般に「箱ひげ図」と呼ばれるグラフを作成する際に用いられることが多い。箱ひげ図の具体例を以下に示す。箱の上端と下端で第1、第3四分位数を示し、箱の中の線で中央値を示す。さらに上下に伸びたひげで最大値、最小値を示す。場合によっては、箱の中に点線で平均値を示したり、一点鎖線で最頻値を示したりといったことも可能である。

箱ひげ図

このグラフの秀逸な点は、データが概ねどのような分布をしているかを示すことができる点である。例の箱ひげ図の場合、上にひげが長く伸び、さらに第3四分位数と中央値の間が、第1四分位数と中央値の間よりも広がっていることから、データの分布として、全体として点の低いほうに山が偏っており、点の高いほうに裾野が伸びた分布をしているということが読み取れる。

以前は棒グラフで平均値を示したグラフや、平均値の棒グラフに標準偏差の幅分だけ「ひげ」（エラーバーと呼ぶ）をつけたグラフを示すのがほとんどであったが、平均と標準偏差を示すだけでは、上記のような分布の形状までは分からない。このため、最近では、箱ひげ図を示すことも多くなってきている。

1.3 記述統計量をまとめて出力するコマンド

以上、個別に記述統計量を算出する方法について説明した。 算出したこれらの統計量を別の計算で用いたい場合には、個別にこれらのコマンドで算出する必要がある。 一方、レポートや論文などで記述統計量を報告したり、とりあえず全体的な傾向を把握したいといった場合には、summary()関数やdescribe()関数を利用することができる

1.3.1 summary()関数

summary()関数を使うと、ベクトルやデータフレームに含まれている各列について、四分位数と平均値を一度に得ることができる。また、ベクトルやデータフレームの列が文字列型だった場合にはデータの点数が表示される．

ベクトルやデータフレームの列が要因型（factor型）だった場合には要因ごとのデータ点数が表示される。

簡単にデータ全体の傾向を把握することができるので、データを読み込んだ場合にはsummary()をまずは実施する、というのがデータ分析の基本フローと言えるだろう。

summary(Shohoku) #データフレームを与えると、含まれているベクトルごとの平均と四分位数を出力される。

##      nam                number           score       
##  Length:12          Min.   : 4.000   Min.   : 4.000  
##  Class :character   1st Qu.: 6.750   1st Qu.: 5.000  
##  Mode  :character   Median : 9.500   Median : 6.500  
##                     Mean   : 9.417   Mean   : 6.417  
##                     3rd Qu.:12.250   3rd Qu.: 7.000  
##                     Max.   :15.000   Max.   :10.000

summary(Shohoku$score) #ベクトルを与えると、そのベクトルの平均と四分位数をが出力される。

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.000   5.000   6.500   6.417   7.000  10.000

1.3.2 describe()関数

psychパッケージに含まれるdescribe()関数を用いるとより詳しい記述統計量を得ることができる。psychパッケージを利用するためには、上部のメニューの中からTools->install.packagesを選択したうえで、以下のようにPackagesにpsychと入力してパッケージをインストールする。

パッケージのインストール

インストールが完了すれば、library()関数を使ってpsychパッケージを読み込み、describe()関数を実行すると、以下のような結果が表示される。

library(psych)

## 
## 次のパッケージを付け加えます: 'psych'

## 以下のオブジェクトは 'package:ggplot2' からマスクされています:
## 
##     %+%, alpha

describe(Shohoku)

##        vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
## nam*      1 12 6.50 3.61    6.5     6.5 4.45   1  12    11  0.00    -1.50 1.04
## number    2 12 9.42 3.50    9.5     9.4 4.45   4  15    11 -0.03    -1.47 1.01
## score     3 12 6.42 1.83    6.5     6.3 1.48   4  10     6  0.41    -0.88 0.53

それぞれ、以下のような指標である

名前	内容
mean	平均
sd	標準偏差（不偏標準偏差）
median	中央値
trimmed	トリム平均（データの上下5%ずつを取り除いて算出した平均）
mad	平均偏差
min	最小値
max	最大値
range	最大値－最小値
skew	歪度（データの分布がどの程度左右に偏っているか。負の場合は左に、正の場合には右に偏っている。また値が大きいほど偏りが大きい。）
kurtosis	尖度（データの分布がどの程度尖っているか。正規分布の尖り具合を0として、正の場合には正規分布より尖っており、負の場合には正規分布より平坦になっている）
se	標準誤差（詳しくは後述{#平均値の区間推定}。）

なお、上記のname列に対する結果の通り、describe()関数はたとえ列の型が文字列であっても数字を出力してくる。当然ながら意味は全くない。参考までに、この出力結果は各要素のインデックス番号をデータとして計算を行ったものであある。

この例のようにRでは標準的なコマンド以外にも様々なコマンドが用意されている。ただし、それらのコマンドを利用するためには、そのコマンドが定義された外部パッケージをTools->install.packagesによってインストールし、さらにスクリプト内でlibrary()関数を使ってそのパッケージを読み込む必要がある。 なおlibrary()関数での外部パッケージの読み込みはスクリプトの実行中1回だけで良い。1度読み込みば、その後はその外部パッケージのコマンドを入力するだけでコマンドを実行することができる。 またパッケージのインストールも一度行えばその後は行う必要はなくなる。

【課題ファイル】には大学共通テストの結果を基にシミュレーションした1000人分の各教科の得点を記載している．このデータについて，各教科の平均値，中央値，最頻値，分散，標準偏差，第1四分位数，第3四分位数を求めなさい．ただし，summary()関数とdescribe()関数は用いてはらならない．

2 推測統計量

2.1 全数調査とサンプリング調査

何らかの調査をするときに、その調査が想定している対象全体のことを母集団という。それに対して、母集団の中から一部だけを抜け出して調査を行った場合、実際に調査を行った対象をサンプル（標本）集団と呼び、サンプル集団に含まれる個々のデータのことをサンプル（標本）と呼ぶ。また、母集団に対して行う調査を母集団調査、あるいは全数調査と呼び、サンプル集団に対して行う調査をサンプル調査と呼ぶ。

母集団とサンプル

母集団調査（全数調査）を実施できればそれに越したことはないが、実際には費用的・時間的な制約であったり、現実に不可能であったりするために、母集団調査を実施できないことは多い。そういう場合には、現実的に実施可能な限られた数のサンプル調査を行い、その結果をもとに母集団の特徴を推測するということが行われる。

母集団調査（全数調査）の例

• 中学での校内の定期テスト：あくまでその中学の中での成績を把握するためのものであり、対象はその中学の生徒に限られ、さらに全員が受けることが前提のテストであることから、母集団にアクセスできていると言える。同様に大学入試共通テストもその年度の全受験生が受けている（もちろん受けない受験生もいないわけではないが、そのような受験生はそもそも対象とはなっていない）ので、これも母集団全体にアクセスできている例である。
• 国勢調査：国勢調査はその時の日本の状態を把握することを目的に、5年に一回、日本国に在住のすべての人を対象に行われる調査であり、母集団そのもの（その時に日本に在住している人全て）にアクセスしている例である。
• 小中学校で行われる健康診断：日本の場合、小学校と中学校は義務教育として対象者全員が通うことが義務付けられている。さらに健康診断は学校保健安全法により全ての学校で実施することが定められている。したがって、小中学校における健康診断はその年度の日本の6歳児から15歳児までの全員に対する診断であり、母集団にアクセスしている例である。一方、高校や大学でも全員に対する健康診断は行われているが、これらの学校は義務教育ではない。したがって、これらの学校での健康診断結果は、日本の16歳〜22歳の母集団全体に対する調査結果ではない。

母集団にアクセスできず、あくまでサンプル調査となる例

• 血液型による性格の違いの検証のための性格診断調査：実際に血液型によって性格が違っているかを全数調査によって検証するためには、全世界の全ての人に調査しなければならない。さらに現時点での全世界の人にアクセスできたとしても「単に今現在はそうなってるだけだろう」というツッコミが入り得るので、理論としての完全性を確保するためには、過去の全人類やこれから先の全人類を調査しなければならない。当然ながらそんなことは不可能である。したがって、血液型による性格の違いを母集団調査によって検証することは不可能であり、どれだけ多くのデータを集めたとしても、それはあくまでサンプル調査である。
• 新聞社等が実施する世論調査：世論調査は各自治体が持つ選挙人名簿から10人おきや20人おきにピックアップされた人に対してハガキや電話、あるいは訪問によって調査している。得たいのはその時点での日本全体の世論であり、その時の日本人への全数調査ができれば良いが、それをするにはコスト的にも時間的にも現実的ではないので、限られた数のサンプルだけに調査を行い、そこから得た結果をもとに、日本全体の世論を推測している。
• 視聴率調査：視聴率はビデオリサーチ社がモニター世帯に選ばれた世帯に対して「どの番組を視聴しているか」を記録する機械を設置し、その機械に記録された情報をもとに算出している。あくまでモニター世帯だけしか調査は行なっていないのでサンプル調査であり、その結果をもとに「日本全体でどの番組がどの程度視聴されているか」を推測している。

そのほか、マーケティング調査なども基本的に全てサンプル調査である。

2.2 母数と標本統計量

母集団全体での記述統計量を母数と呼ぶ。また、特に平均、分散については母平均、母分散と呼ぶこともある。 それに対してサンプル集団の記述統計量を標本統計量と呼び、平均を標本平均、分散を標本分散と呼ぶ。

母集団全体にアクセスできる場合には、得られたデータの記述統計量は母数そのものとなる。それに対してサンプル集団にしかアクセスできない場合には、そのデータの記述統計量はあくまで標本統計量であり母数ではない。

2.3 標本統計量からの母数の推測

先の例にも述べた通り、サンプル調査ではデータそのものの記述統計量（つまり標本統計量）が得たいのではなく、その結果を基に母数を推定することが目的である。従って、標本統計量から母数を推定する必要がある。

推定には大きく2種類あり、値を単刀直入に推定することを点推定と呼ぶ。一方、「〇%の確率で、この範囲内に母数がある可能性が高い」という範囲を推定することを区間推定と呼ぶ。点推定の値のことを不偏推定量、区間推定の範囲のことを〇%信頼区間と呼ぶ。なお、区間推定では「95%」が用いられることが多いが、場合によっては90%や99%が用いられることもああり、その場合には信頼区間の呼び名もそれぞれの確率値に変わる。

標本統計量と点推定値，信頼区間の関係を以下の図に示す。

標本統計量と推測統計量の関係

細かな導出過程や導出に当たっての考え方については省略するが、母数の推定値は以下の数式によって算出される。なお、以下では\(\mu\), \(\sigma^2\)は母平均と母分散（要するに真の値）、\(\hat\mu\)や\(\hat{\sigma^2}\)はデータから推定した母平均と母分散の推定量、\(\bar{x}\)と\(s^2\)は標本平均、標本分散を表す。

2.3.1 平均値の点推定

母平均の不偏推定量を式で表すと以下のようになる。 \[ \mu \approx \hat{\mu}=\bar{x}=\frac{1}{n}\sum^{n}_{i=1}x_{i} \] 要するに、母平均の推定量は標本平均そのものである。つまり、サンプル調査を行った時、得られたデータの平均を母集団全体の平均と見なそう、ということである。

したがってRで母平均の推定量を求めるときも、従来通りmean()関数を用いればよい。

mean(Shohoku$score)

## [1] 6.416667

2.3.2 分散の点推定

母分散の不偏推定量（不偏分散）は以下の式となる。

\[ \sigma^2 \approx \hat{\sigma^2}=\frac{n}{n-1}s^2=\frac{1}{n-1}\sum^{n}_{i=1}(x_i-\bar{x})^2 \]

要するに、標本分散を求める時には「データと平均との差の二乗の和」をnで割っていたのを、n-1で割る形にしているだけである（なぜ\(n-1\)になるかは結構面倒な話になるので省略。興味ある人は上の数式の導出過程をネットで調べてみてほしい。ちなみにこの\(n-1\)のことを自由度と呼ぶ）。

さらに母標準偏差の不偏推定量は、上記の式の平方根となる。

\[\begin{align} \sigma \approx \hat{\sigma}&=\sqrt{\frac{n}{n-1}s^2} = \sqrt{\frac{1}{n-1}\sum^{n}_{i=1}(x_i-\bar{x})^2} \\ &= \sqrt{\frac{n}{n-1}}s \end{align}\]

Rで不偏分散や標準偏差の不偏推定量を求める方法については、すでに上記で説明したとおり、var()やsd()をそのまま用いればよい。

すなわち、先ほどの湘北高校のデータを用いるならば、以下の通りとなる。

var(Shohoku$score)
sd(Shohoku$score)

## [1] 3.356061
## [1] 1.831955

ちなみに、この数値の解釈として、湘北高校のバスケ部員はフリースローで10本中、約6本（平均の不偏推定量）を入れることができるが、メンバー間でのばらつきとして1.8点ほどの標準偏差を持っている、ということである。これはあくまで母数の推定量をであり、赤木・桜木・三井・宮城・流川の5人の得点の平均や5人の中でのばらつきではない。5人をサンプルとした「湘北高校のレギュラーメンバー」という母集団のフリースローの成績の推定値ということである。つまり、赤木・桜木・三井・宮城・流川が入学する前の過去のレギュラーメンバーや、彼らが引退した後の未来のレギュラーメンバーも対象に入れた、湘北高校バスケ部のレギュラーメンバー一般の成績の推定値ということになる。

2.3.3 平均値の区間推定

続いて、区間推定であるが、平均の区間推定を行う際には標準誤差(Standard Error)と呼ばれる統計量を算出する必要がある。標準誤差は以下の数式で算出される。

\[ SE = \sqrt{\frac{\hat{\sigma^2}}{n}}=\frac{\hat{\sigma}}{\sqrt{n}} \] 要するに、標準偏差の不偏推定量を\(\sqrt{n}\)で割った値である。 この標準誤差（以下SE）を用いて、平均値の95%信頼区間（上限と下限）はそれぞれ以下の数式で得られる。

\[\begin{align} 上限 &= \hat\mu + t\times SE \\ 下限 &= \hat\mu - t\times SE \end{align}\]

ここで\(t\)とは自由度\(n-1\)（\(n\)はサンプル数）のT分布表の両側95%から得られる値である。T分布表は大抵の統計学の教科書の末尾に記載されている。具体的には以下のような表である。自由度n-1の行の\(\alpha＝0.025\)の値を読む。例えば、サンプル数が30の場合、自由度は29なので、t＝2.045となる。

T分布表

Rで平均値の信頼区間を求めるにあたっては、2つの方法がある。

まず一つ目はt.test()を使う方法である。本来t.test()関数は後日説明する平均の比較に用いる関数であるが、信頼区間を出力させるのに用いることができる。

res<-t.test(Shohoku$score)
res$conf.int

## [1] 5.252698 7.580636
## attr(,"conf.level")
## [1] 0.95

\(5.252698\)が下限、\(7.580636\)が上限となる。\(0.95\)は信頼区間の確率である。

ちなみに、res自体を出力させると以下のような出力になる。ごちゃごちゃと書かれているが、真ん中あたりに95 percent confidence interval:とあり、信頼区間が記載されている。

print(res)

## 
##  One Sample t-test
## 
## data:  Shohoku$score
## t = 12.133, df = 11, p-value = 1.038e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  5.252698 7.580636
## sample estimates:
## mean of x 
##  6.416667

もし信頼区間の確率設定を変えたい場合には以下のようにすればよい。

res<-t.test(Shohoku$score, conf.level = 0.99)
res$conf.int

## [1] 4.774192 8.059141
## attr(,"conf.level")
## [1] 0.99

二つ目の方法は、先に示した式を自分で打ち込む方法である。

#　平均と標準誤差の算出
n <-length(Shohoku$score)#length()はベクトルの長さ（要素数）を得る関数。もちろんnrow(Shohoku)としてもよい。
res_mean <- mean(Shohoku$score)
res_se <- sd(Shohoku$score)/sqrt(n)  

#ｔ値の取得
t_value <- qt(0.975, df=n-1) # 0.025でもよい。0.025とした場合、負の数字が出力される点に注意すること。

#信頼区間の算出
Band_Upper <- res_mean + t_value * res_se
Band_Lower <- res_mean - t_value * res_se

#出力
c(Band_Lower, Band_Upper)

## [1] 5.252698 7.580636

先ほどと同じ数値が算出されていることが確認できる。

2.3.4 分散の区間推定

分散についての95%信頼区間については不偏分散\(\hat{\sigma^2}\)を用いて次の式で与えられる。

\[\begin{align} 上限 &= \frac{(n-1)\hat{\sigma^2}}{\chi^2_{0.975}} \\ 下限 &= \frac{(n-1)\hat{\sigma^2}}{\chi^2_{0.025}} \end{align}\]

\(\chi^2\)(カイ2乗と呼ぶ)は、\(t\)と同様に自由度n-1（nはサンプル数）のカイ2乗分布表から得られる値である。分布表と同様に、カイ2乗分布表も統計学の教科書には記載されている。具体的には次々項のような表である。自由度n-1の\(\alpha＝0.025と0.975\)の値を読み取る。例えば、サンプル数が30の場合、自由度は29なので、上限算出に使う\(\chi^2\)値は16.047、下限算出に使う\(\chi^2\)値は45.722となる

\(\chi^2\)分布表

Rで分散の信頼区間を算出する方法も2つある。1つはEnvStatsパッケージに含まれているvarTest()関数を使う方法である。使うにはEnvStatsパッケージをインストールしたうえで、library()関数を実行する必要がある（こちら参照)。

library(EnvStats)

## 
## 次のパッケージを付け加えます: 'EnvStats'

## 以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     predict, predict.lm

res<-varTest(Shohoku$score)
res$conf.int

##      LCL      UCL 
## 1.684151 9.674817 
## attr(,"conf.level")
## [1] 0.95

res<-varTest(Shohoku$score, conf.level = 0.90) #90%信頼区間の場合
res$conf.int

##      LCL      UCL 
## 1.876310 8.069546 
## attr(,"conf.level")
## [1] 0.9

ちなみにvarTest()自体の出力結果は以下の通りである。一番下に信頼区間が記載されている。なお、上のコードで確率を90%に設定したものをresに格納しているので、90%信頼区間の結果が表示されている。

print(res)

## 
## Results of Hypothesis Test
## --------------------------
## 
## Null Hypothesis:                 variance = 1
## 
## Alternative Hypothesis:          True variance is not equal to 1
## 
## Test Name:                       Chi-Squared Test on Variance
## 
## Estimated Parameter(s):          variance = 3.356061
## 
## Data:                            Shohoku$score
## 
## Test Statistic:                  Chi-Squared = 36.91667
## 
## Test Statistic Parameter:        df = 11
## 
## P-value:                         0.0002379912
## 
## 90% Confidence Interval:         LCL = 1.876310
##                                  UCL = 8.069546

もう一つの方法は、平均値の区間推定と同様に先に示した式を自分で打ち込む方法である。

n <- nrow(Shohoku) # 今度はnrow()を用いた。当然length()を用いてもよい
res_var <- var(Shohoku$score)

# カイ二乗分布のクリティカル値を使用して信頼区間を計算
chi_upper <- qchisq(0.025, df = n - 1) #上で示した式の上限・下限の確率値が入れ替わっていることに注意
chi_lower <- qchisq(0.975, df = n - 1) 

# 分散の信頼区間
Band_Lower <- (n - 1) * res_var / chi_lower 
Band_Upper <- (n - 1) * res_var / chi_upper

#出力
c(Band_Lower, Band_Upper)

## [1] 1.684151 9.674817

なお、Rで計算するときの\(\chi^2\)値算出の関数qchisq()での確率設定値と、先ほど式で示したものの中での\(\chi^2\)値の確率設定値が入れ替わっている。これは確率を「残された方」で見るのか、「切り落とす方」で見るのかの違いである。入れ替わってしまったとしても数値自体が変わるわけではないので、実際に出力させてみて、上限と下限の値が逆順になっていれば、数値を入れ替えればよい。

平均値の区間推定にしろ、分散の区間推定にしろ、「95%」信頼区間なのに、\(\alpha＝0.05\)や\(0.95\)ではなく\(\alpha＝0.025\)や\(0.975\)の値を読むのはなぜなのか疑問に思う人もいると思う。これは「95％の確率でその区間の中に平均値がある」という範囲を取り出す、ということはその範囲よりも上下にある2.5%ずつを切り落とすということだからである。仮に5%と95%の値を読むことにしてしまうと、残された区間は90%となってしまう。

先ほどの【課題ファイル】にある大学共通テストのデータについて，各教科の平均と分散の点推定値および95%の推定区間を求めなさい．

データフレームの各列に対する平均・標準偏差

これまでに紹介してきたmean()やvar()、sd()はベクトルに対しての関数である。 データフレームの各列に対して平均や分散，標準偏差を算出しようとするとき，これらの関数では各列を一つ一つ$を用いてベクトルとして取り出して，これらの関数に与えていく必要がある． 取り出そうとしている列の数が少ないのであれば，それでも問題ないが，列の数が多くなってきたときには1つ1つベクトルで取り出していくのは骨が折れてしまう．

そういう場合には，先に紹介したdescribe()関数を用いるのが良いが，ほかに sapply()関数を用いるという手段もある．sapply()関数は，データフレームの各列に対して関数を適用する関数である．第1引数には対象となるデータフレーム，第2引数には適用したい関数を指定する。たとえば平均であればmean，標準偏差であればsd，分散であればvarといった具合である．注意点として，引数としてこれらの関数名を与える際には()はつける必要はなく，あくまで関数名だけを与えればよい．

さらに第3引数として，適用したい関数に引数がある場合にそれを指定することができる．たとえば，mean()で紹介した通り，na.rm=TRUEを第3引数として与えることで，NAが含まれる場合にそれを無視して計算することができる．

このデータは，ある大学で行われた主要5因子性格検査短縮版アンケート（外向性，協調性，勤勉性，神経症傾向，開放性）参考で得られた学生400人分のデータである．このデータから，各性格特性の平均と標準偏差を求めなさい．

出典：小塩 真司, 阿部 晋吾, Pino Cutrone, 日本語版Ten Item Personality Inventory（TIPI-J）作成の試み., パーソナリティ研究, Vol.21, No.1, pp.40-52 (2012)

data <- read.csv("practice/04_work_basicStatistics_sapply.csv")
sapply(data, mean)

##     外向性     協調性     勤勉性 神経症傾向     開放性 
##     7.8675     9.4125     6.4350     9.2200     8.1625

sapply(data, sd)

##     外向性     協調性     勤勉性 神経症傾向     開放性 
##   2.539896   2.087640   2.346175   2.276963   2.405319

3 グループごとの統計量の算出

これまでつかってきたデータ（湘北高校のデータ）には特にグループ（群）分けという概念は存在していないデータであった。 しかし実際のデータには大抵、属性項目（例えば、性別や、実験データにおいては実験条件、あるいは学年や年齢区分、国籍など）によるグループ分けというものが存在する。そういったグループ分けが存在するときに、グループごとに記述統計量（推測統計量）を算出するには、aggregate()関数を用いると良い。

以下では先の湘北高校のデータを、レギュラー組と控え組にグループ分けしたものにしたうえで、グループごとの平均と標準偏差を算出させてみる。

3.1 要因型変数の作成

まずは湘北高校のデータにグルプ分けのための変数Groupを追加する。

Shohoku["Group"] <- as.factor(c("レギュラー", "レギュラー", "レギュラー", "レギュラー", "レギュラー", "控え", "控え", "控え", "控え", "控え", "控え", "控え"))
Shohoku # 確認

##       nam number score      Group
## 1    赤木      4     7 レギュラー
## 2    桜木     10     4 レギュラー
## 3    三井     13    10 レギュラー
## 4    宮城      7     7 レギュラー
## 5    流川     11     9 レギュラー
## 6    木暮      5     6       控え
## 7    潮崎      8     6       控え
## 8    安田      6     7       控え
## 9    角田      9     7       控え
## 10   石井     12     5       控え
## 11 佐々岡     13     5       控え
## 12   桑田     15     4       控え

as.factor()関数に関しては先に軽く触れたが(参考)、ここで改めて詳しく説明する。 グループ分けに使う変数の型は「要因型」と呼ぶ。要因型に設定するためにはc()関数を使って作成したベクトルをas.factor()関数で括ってやる必要がある。内容はレギュラーや控えという文字列になっているが、列全体としては文字列型ではないことに注意してもらいたい。

文字列型と要因型とでは結果が異なる。summary()の実行結果として出力されるものが異なってくる。それぞれ以下で確認してみてほしい。

Shohoku$Group <- as.factor(c("レギュラー", "レギュラー", "レギュラー", "レギュラー", "レギュラー", "控え", "控え", "控え", "控え", "控え", "控え", "控え"))
summary(Shohoku)

##      nam                number           score               Group  
##  Length:12          Min.   : 4.000   Min.   : 4.000   レギュラー:5  
##  Class :character   1st Qu.: 6.750   1st Qu.: 5.000   控え      :7  
##  Mode  :character   Median : 9.500   Median : 6.500                 
##                     Mean   : 9.417   Mean   : 6.417                 
##                     3rd Qu.:12.250   3rd Qu.: 7.000                 
##                     Max.   :15.000   Max.   :10.000

Shohoku$Group <- c("レギュラー", "レギュラー", "レギュラー", "レギュラー", "レギュラー", "控え", "控え", "控え", "控え", "控え", "控え", "控え")
summary(Shohoku)

##      nam                number           score           Group          
##  Length:12          Min.   : 4.000   Min.   : 4.000   Length:12         
##  Class :character   1st Qu.: 6.750   1st Qu.: 5.000   Class :character  
##  Mode  :character   Median : 9.500   Median : 6.500   Mode  :character  
##                     Mean   : 9.417   Mean   : 6.417                     
##                     3rd Qu.:12.250   3rd Qu.: 7.000                     
##                     Max.   :15.000   Max.   :10.000

要因型にした場合には、それぞれのグループに含まれるデータの数が返されるのに対して、文字列型の場合には全体のデータ点数と文字列型であることを示すClass: characterやMode: charaterという表示が返されるだけである。

なお、要因型に設定するのであれば、その内容は数値でも構わない。ただ、数値をグループ分け変数に使用する場合には特にas.factor()を実行するのを忘れないようにしなければならない。文字列型だった場合には数値演算はできないので結果自体は要因型であってもなくても変わらないことはよくあるが、数値型のときにas.factor()を忘れてしまうと、グループ分け用の変数としていたつもりでも数値の変数として処理され、演算の中に組み込まれてしまうことがあり、誤った結果を出力するケースがある（例えば、回帰分析）。

具体的に例を示すと，まずは以下は正しく処理している例．

Shohoku$Group <- as.factor(c(1,1,1,1,1,2,2,2,2,2,2,2)) #1:レギュラー, 2:控え
summary(Shohoku)

##      nam                number           score        Group
##  Length:12          Min.   : 4.000   Min.   : 4.000   1:5  
##  Class :character   1st Qu.: 6.750   1st Qu.: 5.000   2:7  
##  Mode  :character   Median : 9.500   Median : 6.500        
##                     Mean   : 9.417   Mean   : 6.417        
##                     3rd Qu.:12.250   3rd Qu.: 7.000        
##                     Max.   :15.000   Max.   :10.000

一方，こちらはas.factor()を忘れてしまった例．

Shohoku$Group <- c(1,1,1,1,1,2,2,2,2,2,2,2)
summary(Shohoku)

##      nam                number           score            Group      
##  Length:12          Min.   : 4.000   Min.   : 4.000   Min.   :1.000  
##  Class :character   1st Qu.: 6.750   1st Qu.: 5.000   1st Qu.:1.000  
##  Mode  :character   Median : 9.500   Median : 6.500   Median :2.000  
##                     Mean   : 9.417   Mean   : 6.417   Mean   :1.583  
##                     3rd Qu.:12.250   3rd Qu.: 7.000   3rd Qu.:2.000  
##                     Max.   :15.000   Max.   :10.000   Max.   :2.000

Groupは数値型として処理されてしまい，summary()関数の結果も数値型として出力されてしまっている． しかし，元々は単なるラベルとして「1」や「2」としているだけであり，出力されている平均や中央値，四分位数は全く意味がないものである．

このように，あくまでラベルとして数値を用いるときには，Rがその列（ベクトル）を数値として扱わないように，as.factor()を忘れないようにすることが重要である．

3.2 グループ別での演算の実行

aggregate()関数の機能は、ベクトルデータに対して、グループ分けを行った上で、それぞれのグループのデータに対して設定した関数演算を実行するというものである。 第1引数には、演算対象データの列名~グループ分けデータの列名という書式で、それぞれのデータを指定する。第2引数には演算を行いたいデータが入っているデータフレームを与える。第3引数には、実行したい演算関数の名前を与える。入力するものはあくまで関数の名前だけでよく、()はつける必要はない。

結果はデータフレームで出力される。

#グループ分け変数の設定
Shohoku$Group <- as.factor(c("レギュラー", "レギュラー", "レギュラー", "レギュラー", "レギュラー", "控え", "控え", "控え", "控え", "控え", "控え", "控え"))

#平均を算出
res.mean<-aggregate(score~Group, Shohoku,mean)
print(res.mean)

##        Group    score
## 1 レギュラー 7.400000
## 2       控え 5.714286

#標準偏差を算出
res.sd<-aggregate(score~Group, Shohoku, sd)
print(res.sd)

##        Group    score
## 1 レギュラー 2.302173
## 2       控え 1.112697

次の【課題ファイル】は，Setosa, Versicolor, Verginica,という 3種類のアイリスの花の花びらとがくの長さと幅についてのデータである．このデータについて，花の種類ごとに花びらとがくの長さと幅の平均と標準偏差をそれぞれ求めなさい．

4 今回用いたR関数

• read.csv()
  • csvファイルを読み込む
  • 第1引数: 読み込むファイル名もしくはファイルパス
  • 第2引数: ヘッダー(header)がある(TRUE)かない(FALSE)か
  • 第3引数: 区切り文字．sep=“<<区切り文字>>”. 何も指定しなければカンマ．
  • 第4引数: 文字コード．fileEncoding= “UTF-8” もしくは”SJIS”など．
• mean()
  • 引数で指定されたベクトルが数値ベクトルの場合に，その平均値を返す
  • 第1引数: 数値データのベクトル
  • 第2引数(省略可): na.rm=TRUEを指定するとNAを無視して計算する．省略するとNAが含まれる場合はNAを返す．
• median()
  • 引数で指定されたベクトルが数値ベクトルの場合に，その中央値を返す
  • 引数: 数値データのベクトル
• table()
  • 与えたベクトルデータの度数分布表を算出する
  • 引数: ファクター型データのベクトル
• which.max()
  • ベクトルデータの中で最大値を持つ要素のインデックスを返す
  • 引数: ベクトルデータ
• names()
  • データフレームの列名を取得する
  • 引数: データフレーム
• as.numeric()
  • 文字列型を数値型に変換する
  • 引数: 文字列データ
• mad()
  • データの平均絶対偏差を算出する
  • 引数: ベクトルデータ
• var()
  • データの分散（不偏分散）を算出する
  • 引数: ベクトルデータ
• sd()
  • データの標準偏差（不偏標準偏差）を算出する
  • 引数: ベクトルデータ
• quantile()
  • データのパーセンタイル値を算出する
  • 引数: ベクトルデータ
• library()
  • インストールしたパッケージを現在実行しているRで読み込む
  • 引数: パッケージ名
• summary()
  • データの要約統計量を算出する
  • 引数: データフレーム
• describe()
  • データの記述統計量を算出する．実行するためにはpsychパッケージをインストールしたうえで，現在実行しているRで読み込む必要がある．
  • 引数: データフレーム
• t.test()
  • 本来は2群間の平均値の差の検定を行う関数だが，1つの引数だけを与えた場合には平均値の信頼区間を得ることができる
  • 引数: データフレーム
• varTest()
  • 分散の信頼区間を算出する．実行するためにはEnvStatsパッケージをインストールしたうえで，現在実行しているRで読み込む必要がある．
  • 引数: データフレーム
• aggregate()
  • データフレームの列ごとにグループ分けを行い，それぞれのグループに対して演算を行う
  • 引数: 演算対象データの列名~グループ分けデータの列名, データフレーム, 演算関数の名前