統計2
Tomoki Motegi
2024-09-26
統計の分類
要約統計量
データの持つ意味を理解するために情報を要約したもの。
Rでは非常に簡単に表現することができる。
以下のデータを用いて作成してみる
## Rに収載されているデータ"iris"を呼び出して"iris_data"に格納
iris_data = iris
## 最初の6行のみを表示
head(iris_data)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa## なお、以降このデータを使うため、一度Rstudioを閉じた場合は一番最初から順々に動かして、各データを取ると良い。“iris”のデータは３品種（Species）のアヤメ（setosa, versicolor,
virginica）のがくの長さ（Sepal.length）と幅（Sepal.Width）、花びらの長さ（Petal.Length）と幅（Petal.width）について50サンプル分調べたデータである。
無料で公開されているため、様々なところで使用される。
それでは要約統計量の作成である。
## 要約統計量を計算するsummayコマンドを実行
summary(iris_data)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## 左列からがくの長さ（Sepal.length）・幅（Sepal.Width）・花びらの長さ（Petal.Length）・幅（Petal.width）・品種（Species）の順に並んでおり、各要約統計量が下の行に書かれている。
各パラメータは以下のとおりである。
Min.：最小値
1st Qu.：第一四分位数
Median：中央値
Mean：平均値
3rd Qu.：第三四分位数
Max.：最大値
品種については品種ごとの総数が記載されている。
続いてデータを視覚化してみる。
## がくの長さ（Sepal.length）・幅（Sepal.Width）・花びらの長さ（Petal.Length）・幅（Petal.width）四つヒストグラムの作成
## このデータを25個の階級に分類している
hist(iris_data$Sepal.Length, breaks = 25)
hist(iris_data$Sepal.Width, breaks = 25)
hist(iris_data$Petal.Length, breaks = 25)
hist(iris_data$Petal.Width, breaks = 25)
横軸は1-25段階に分けた各階級、縦軸はデータの出現頻度を示している。
結果からわかるようにがくの長さと幅については比較的平均値がその数値を代表しているように見える。
一方で花びらの長さと幅については二つのグループに分かれており平均値・中央値共にその値が代表的な数字に見えない。
しかし、一般的なの感覚として、花びらの長さと幅の値がこんなにばらつくことは考えにくい。
そこでこの原因を探っていく。
今回のデータには3品種が混ざっていたため、各々に分けて調べてみる
## 花びらの長さ（Petal.Length）・幅（Petal.width）について各3品種でヒストグラムを作成してみる
## setosa, versicolor, virginicaの各グループごとに分けたデータを作成
iris_sub_set = subset(x = iris_data, subset = iris_data$Species == "setosa")
iris_sub_ver = subset(x = iris_data, subset = iris_data$Species == "versicolor")
iris_sub_vir = subset(x = iris_data, subset = iris_data$Species == "virginica")
## 花びらの長さ（Petal.Length)
## 元のデータ
hist(iris_data$Petal.Length, breaks = 25)
## setosaのヒストグラム
hist(iris_sub_set$Petal.Length, breaks = 25)
## versicolorのヒストグラム
hist(iris_sub_ver$Petal.Length, breaks = 25)
## virginicaのヒストグラム
hist(iris_sub_vir$Petal.Length, breaks = 25)
## 花びらの幅（Petal.width）
## 元のデータ
hist(iris_data$Petal.Width, breaks = 25)
## setosaのヒストグラム
hist(iris_sub_set$Petal.Width, breaks = 25)
## versicolorのヒストグラム
hist(iris_sub_ver$Petal.Width, breaks = 25)
## virginicaのヒストグラム
hist(iris_sub_vir$Petal.Width, breaks = 25)
結果を見ると花びらの長さについては、1-2のピークはsetosa、3-5のピークはversicolor、4.5-7までのピークはvirginicaが作っていることがわかる。
同様に花びらの幅は0.1-0.6のピークはsetosa、1.0-1.8のピークはversicolor、1.4-2.6までのピークはvirginicaが作っている。
花びらの長さは平均値が各品種の代表値に近いように見えるが、花びらの幅についてはvirginicaのみ平均値が代表値からやや外れているように見える。
これをもう少しわかりやすくなるように視覚化してみる。
## 花びらの長さ（Petal.Length)を箱ひげ図で表現
## ggplot2パッケージの呼び出し
library(ggplot2)## Warning: パッケージ 'ggplot2' はバージョン 4.2.3 の R の下で造られました## ggplot2の機能を使って箱ひげ図を作成
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Length, fill = Species)) +
geom_boxplot() +
labs(x = "", y = "Petal.Length")+
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
## 花びらの幅（Petal.Width)を箱ひげ図で表現
## ggplot2の機能を使って箱ひげ図を作成
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Width, fill = Species)) +
geom_boxplot() +
labs(x = "", y = "Petal.Width") +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
これを見ると花びらの幅は品種間で箱ひげ図の幅が異なるためデータの幅があることがわかる。
しかし分布が分かりづらいため、ここに各個体ごとのデータを重ねてみる。
まずはきれいに分かれている花びらの長さから。
## 花びらの長さを箱ひげ図で表現
## 今回は箱ひげ図の外れ値を表示しないようにしている
## この上に個体のデータをドットプロットで記載
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Length, fill = Species)) +
geom_boxplot(outlier.shape = NA) +
labs(x = "", y = "Petal.Length") +
geom_jitter(size = 1.0, width = 0.1) +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
続いて花びらの幅の視覚化。
## 花びらの幅（Petal.Width)を箱ひげ図で表現
## 今回は箱ひげ図の外れ値を表示しないようにしている
## この上に個体のデータをドットプロットで記載
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Width, fill = Species)) +
geom_boxplot(outlier.shape = NA) +
labs(x = "", y = "Petal.Width") +
geom_jitter(size = 1.0, width = 0.1) +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
重ね合わせてみると箱ひげ図が示す中央値に集まらず、測定値にばらつきがあり、その差が分かる。
しかしこれでもまだ視覚的に理解することが難しい。
このため分布情報を盛り込んだバイオリンプロットを作成してみる。
## 花びらの長さをバイオリンプロットで表現
## この上に個体のデータをドットプロットで記載
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Length, fill = Species)) +
geom_violin(trim = FALSE) +
labs(x = "", y = "Petal.Length") +
geom_jitter(size = 1.0, width = 0.1, alpha = 0.3) +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
## 花びらの幅バイオリンプロットで表現
## この上に個体のデータをドットプロットで記載
bp_gg = ggplot(data = iris_data,
aes(x = Species, y = Petal.Width, fill = Species)) +
geom_violin(trim = FALSE) +
labs(x = "", y = "Petal.Width") +
geom_jitter(size = 1.0, width = 0.1, alpha = 0.3) +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
バイオリンプロットを見ると花びらの長さではsetosaと異なってvirginicaの中央の大きなピークが存在しておらず、幅があることがわかる。また花びらの幅では、各品種で複数個所にピークが見えており、安易に中央値や平均値で各品種の代表的な値が取れないことがわかる。
標準偏差と信頼区間
四分位範囲は不定な分布を理解するためには便利だが、正規分布のようなきれいな分布になるものであれば平均と標準偏差で求めたほうが比較しやすい。そこでがくの長さ（Sepal.Length）と幅（Sepal.Width）について標準偏差と信頼区間を比べていく。まずは各データのヒストグラムを作成する。
## がくの長さ（Sepal.length）と幅（Sepal.Width）について各3品種でヒストグラムを作成してみる
## setosa, versicolor, virginicaの各グループごとに分けたデータを作成
iris_sub_set = subset(x = iris_data, subset = iris_data$Species == "setosa")
iris_sub_ver = subset(x = iris_data, subset = iris_data$Species == "versicolor")
iris_sub_vir = subset(x = iris_data, subset = iris_data$Species == "virginica")
## 花びらの長さ（Petal.Length)
## 元のデータ
hist(iris_data$Sepal.Length, breaks = 25)
## setosaのヒストグラム
hist(iris_sub_set$Sepal.Length, breaks = 25)
## versicolorのヒストグラム
hist(iris_sub_ver$Sepal.Length, breaks = 25)
## virginicaのヒストグラム
hist(iris_sub_vir$Sepal.Length, breaks = 25)
## 花びらの幅（Petal.width）
## 元のデータ
hist(iris_data$Sepal.Width, breaks = 25)
## setosaのヒストグラム
hist(iris_sub_set$Sepal.Width, breaks = 25)
## versicolorのヒストグラム
hist(iris_sub_ver$Sepal.Width, breaks = 25)
## virginicaのヒストグラム
hist(iris_sub_vir$Sepal.Width, breaks = 25)
いくつかのデータは正規分布していないように見えるが、データ数は十分にあるため正規分布すると仮定して進めていく。ちなみに正規分布するデータ生成すると以下のようなヒストグラムになる。乱数を固定していないためndを何回か実行すると見た目上あまり正規分布していないようなデータをが出現することがわかる。
## 正規分布する50個の数字を生成する
nd = rnorm(50)
hist(nd, breaks = 25)
しかしデータ量を増やすときれいなデータになっていく。
## 正規分布する500個の数字を生成する
nd = rnorm(500)
hist(nd, breaks = 250)
## 正規分布する500個の数字を生成する
nd = rnorm(50000)
hist(nd, breaks = 25000)
では本題のがくの長さと幅について調べていく。まず各データの平均と標準偏差を求めていく。
## がくの長さについて各グループの平均と標準偏差を表示
# setosaの平均
mean(iris_sub_set$Sepal.Length)## [1] 5.006# setosaの標準偏差
sd(iris_sub_set$Sepal.Length)## [1] 0.3524897# versicolorの平均
mean(iris_sub_ver$Sepal.Length)## [1] 5.936# versicolorの標準偏差
sd(iris_sub_ver$Sepal.Length)## [1] 0.5161711# virginicaの平均
mean(iris_sub_vir$Sepal.Length)## [1] 6.588# virginicaの標準偏差
sd(iris_sub_vir$Sepal.Length)## [1] 0.6358796## がくの幅について各グループの平均と標準偏差を表示
# setosaの平均
mean(iris_sub_set$Sepal.Width)## [1] 3.428# setosaの標準偏差
sd(iris_sub_set$Sepal.Width)## [1] 0.3790644# versicolorの平均
mean(iris_sub_ver$Sepal.Width)## [1] 2.77# versicolorの標準偏差
sd(iris_sub_ver$Sepal.Width)## [1] 0.3137983# virginicaの平均
mean(iris_sub_vir$Sepal.Width)## [1] 2.974# virginicaの標準偏差
sd(iris_sub_vir$Sepal.Width)## [1] 0.3224966この何気なく出したデータの凄い点はデータが1000点でも10000点でも平均±標準偏差で表せてしまうことである。そして、このデータの約68%はこの範囲に存在することがわかる。ではこの2群を棒グラフで比較してみる。
## まずSepal.Lengthの各データを格納する
## 平均のデータ作成
SL_mean = as.data.frame(c(mean(iris_sub_set$Sepal.Length),
mean(iris_sub_ver$Sepal.Length),
mean(iris_sub_vir$Sepal.Length)))
## sdのデータ作成
SL_sd = as.data.frame(c(sd(iris_sub_set$Sepal.Length),
sd(iris_sub_ver$Sepal.Length),
sd(iris_sub_vir$Sepal.Length)))
## お互いのデータを結合
SL_sum = cbind(SL_mean, SL_sd)
# 行名を指定
rownames(SL_sum) = c("setosa", "versicolor", "virginica")
# 列名を指定
colnames(SL_sum) = c("mean", "sd")
## 棒グラフを描く
library(ggplot2)
box_SL = ggplot(data = SL_sum,
aes(x = rownames(SL_sum), y = mean, fill = rownames(SL_sum))) +
labs(x="", y = "Sepal.Length") +
theme_classic() +
scale_fill_manual(values = c("red", "green", "blue"))+
geom_bar(stat = "identity", colour = "black")+
geom_errorbar(aes(ymax = mean + sd, ymin = mean - sd), width = 0.1)
box_SL
## 続いてSepal.Widthの各データを格納する
## 平均のデータ作成
SW_mean = as.data.frame(c(mean(iris_sub_set$Sepal.Width),
mean(iris_sub_ver$Sepal.Width),
mean(iris_sub_vir$Sepal.Width)))
## sdのデータ作成
SW_sd = as.data.frame(c(sd(iris_sub_set$Sepal.Width),
sd(iris_sub_ver$Sepal.Width),
sd(iris_sub_vir$Sepal.Width)))
## お互いのデータを結合
SW_sum = cbind(SW_mean, SW_sd)
# 行名を指定
rownames(SW_sum) = c("setosa", "versicolor", "virginica")
# 列名を指定
colnames(SW_sum) = c("mean", "sd")
## 棒グラフを描く
library(ggplot2)
box_SW = ggplot(data = SW_sum,
aes(x = rownames(SW_sum), y = mean, fill = rownames(SW_sum))) +
labs(x="", y = "Sepal.Width") +
theme_classic() +
scale_fill_manual(values = c("red", "green", "blue"))+
geom_bar(stat = "identity", colour = "black")+
geom_errorbar(aes(ymax = mean + sd, ymin = mean - sd), width = 0.1)
box_SW
このデータを見るとSepal.Lengthのsetosaとvirginicaは平均±標準偏差の範囲では重ならないように見える。ではこれを95％信頼区間の範囲で見てみる。95%信頼区間自体は簡単な計算で求めることができるが、より簡便にRの1標本のt検定を行うと（t.test）その内部で計算されるため、それを用いる。
## まず1標本t検定を行って結果を見てみる
## setosaのSepal.Lengthを見てみる
iris_sub_set = subset(x = iris_data, subset = iris_data$Species == "setosa")
iris_sub_ver = subset(x = iris_data, subset = iris_data$Species == "versicolor")
iris_sub_vir = subset(x = iris_data, subset = iris_data$Species == "virginica")
t.test(iris_sub_set$Sepal.Length)##
## One Sample t-test
##
## data: iris_sub_set$Sepal.Length
## t = 100.42, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 4.905824 5.106176
## sample estimates:
## mean of x
## 5.006## 各データの1標本t検定の結果から平均と95 percent confidence intervalを抽出する
set_95ci = c(t.test(iris_sub_set$Sepal.Length)$estimate,
t.test(iris_sub_set$Sepal.Length)$conf.int[1:2])
ver_95ci = c(t.test(iris_sub_ver$Sepal.Length)$estimate,
t.test(iris_sub_ver$Sepal.Length)$conf.int[1:2])
vir_95ci = c(t.test(iris_sub_vir$Sepal.Length)$estimate,
t.test(iris_sub_vir$Sepal.Length)$conf.int[1:2])
## 95%信頼区間のデータを結合する
all_95ci = as.data.frame(rbind(set_95ci,ver_95ci,vir_95ci))
# 行名を指定
rownames(all_95ci) = c("setosa", "versicolor", "virginica")
# 列名を指定
colnames(all_95ci) = c("mean", "Lower", "Upper")
## 95%信頼区間をグラフ化
library(ggplot2)
ci_gg = ggplot(all_95ci, aes(x = rownames(all_95ci), y = mean)) +
geom_point(size = 3) +
xlab("Species") +
ylab("Confidence interval") +
theme_bw() +
geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, size=1) +
coord_flip() #図を横に## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ci_gg
続いてがくの幅を見てみる
## 1標本t検定を行って結果を見てみる
## setosaのSepal.Widthを見てみる
t.test(iris_sub_set$Sepal.Width)##
## One Sample t-test
##
## data: iris_sub_set$Sepal.Width
## t = 63.946, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 3.320271 3.535729
## sample estimates:
## mean of x
## 3.428## 各データの1標本t検定の結果から平均と95 percent confidence intervalを抽出する
set_95ci = c(t.test(iris_sub_set$Sepal.Width)$estimate,
t.test(iris_sub_set$Sepal.Width)$conf.int[1:2])
ver_95ci = c(t.test(iris_sub_ver$Sepal.Width)$estimate,
t.test(iris_sub_ver$Sepal.Width)$conf.int[1:2])
vir_95ci = c(t.test(iris_sub_vir$Sepal.Width)$estimate,
t.test(iris_sub_vir$Sepal.Width)$conf.int[1:2])
## 95%信頼区間のデータを結合する
all_95ci = as.data.frame(rbind(set_95ci,ver_95ci,vir_95ci))
# 行名を指定
rownames(all_95ci) = c("setosa", "versicolor", "virginica")
# 列名を指定
colnames(all_95ci) = c("mean", "Lower", "Upper")
## 95%信頼区間をグラフ化
library(ggplot2)
ci_gg = ggplot(all_95ci, aes(x = rownames(all_95ci), y = mean)) +
geom_point(size = 3) +
xlab("Species") +
ylab("Confidence interval") +
theme_bw() +
geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, size=1) +
coord_flip() #図を横に
ci_gg
結果を見るとがくの長さも幅もすべての群で信頼区間に重なりが無い。このためこれらは有意に異なっていると言えそうだが、この次はこれを有意差検定で証明していく。