はじめに
目的： R を使ったデータ分析を通じ、データドリブンな意思決定力や論理的思考力を身につける
使わせていただいた書籍： Rによるやさしい統計学（著者：山田剛史氏、杉澤武俊氏、村井潤一郎氏）。こちらに記載の検定手法をベースに、Kaggleからデータセットを収集し分析してみる
GitHub： 今回作成したソースコードはGitHubに公開している
※注意１： データセットを DL して得られた .csv も GitHub リポジトリに入れていますが、７．主要都市の気温のデータセットのみサイズが大きすぎるので縮小版を GitHub にあげています。縮小前の本物の .csv はリンク先から入手できます
※注意２： 以下のコードブロックは必要な事前処理が省略されており、検定部分しか書いていません。事前処理の内容は GitHub をご参照ください
１．アイリスのデータセット
アイリスのデータセット（150行×6列）を使用
わかったこと
- アイリスの花びらの長さ（cm）と花びらの幅（cm）の標本統計量には強い相関がある
- アイリスのがく片の長さ（cm）の母平均は、95%の確率で5.71～5.98（cm）の間にある
- アイリスのがく片の長さ（cm）とがく片の幅（cm）の母相関は95%の確率で０である
- セトサのがく片の長さ（cm）とバージカラーのがく片の長さ（cm）の母平均には、95%の確率で1.11～0.75（cm）の差がある
iris.R（アイリスのがく片の長さの母平均）
> t.test(iris$SepalLengthCm, mu=6)
One Sample t-test
data: iris$SepalLengthCm
t = -2.3172, df = 149, p-value = 0.02186
alternative hypothesis: true mean is not equal to 6
95 percent confidence interval:
5.709732 5.976934
sample estimates:
mean of x
5.843333
- 検定方法
- t分布を用いた検定（１つの平均値の検定、母分散σ^2が未知）
- 分布
- 統計量は、帰無仮説のもとで自由度 df=n-1 の t分布 にしたがう
- 帰無仮説
- μ=6（アイリスのがく片の長さ（cm）の母平均は6である）
- 対立仮説
- μ≠6（母平均は6でない）
- 結果
- p値が0.02186なので帰無仮説は棄却される
- また、本当の母平均は95%の確率で5.709732～5.976934の間にある
- すなわち、母平均が６の可能性は低いと結論づけられる
iris.R（アイリスのがく片の長さとがく片の幅の母相関）
> cor.test(iris$SepalLengthCm, iris$SepalWidthCm)
Pearson's product-moment correlation
data: iris$SepalLengthCm and iris$SepalWidthCm
t = -1.3386, df = 148, p-value = 0.1828
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.26498618 0.05180021
sample estimates:
cor
-0.1093692
- 検定方法
- 相関係数の検定（無相関検定）
- 分布
- 統計量は、帰無仮説のもとで自由度 df=n-2 の t分布 にしたがう
- 帰無仮説
- ρ=0（アイリスのがく片の長さとがく片の幅の母相関は0である）
- 対立仮説
- ρ≠0（母相関は0ではない）
- 結果
- p値が0.1828なので帰無仮説は棄却されない
- また、本当の母相関は95%の確率で-0.26498618～0.05180021の間にある
- すなわち、母相関が０の可能性があると結論づけられる
iris.R（セトサのがく片の長さとバージカラーのがく片の長さの母平均の差）
> t.test(setosa$SepalLengthCm, versicolor$SepalLengthCm, var.equal = F)
Welch Two Sample t-test
data: setosa$SepalLengthCm and versicolor$SepalLengthCm
t = -10.521, df = 86.538, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.1057074 -0.7542926
sample estimates:
mean of x mean of y
5.006 5.936
- 検定方法
- 独立な２群のt検定
- 分布
- 統計量は、帰無仮説のもとで自由度 df = n1 + n2 -2 の t分布 にしたがう
- 帰無仮説
- セトサのがく片の長さとバージカラーのがく片の長さの母平均は等しい
- 対立仮説
- ２つの母平均は等しくない
- 結果
- p値が2.2e-16以下なので帰無仮説は棄却される
- また、本当の母平均の差は95%の確率で-1.1057074～-0.7542926の間にある
- すなわち、２つの母平均の差が０（母平均が等しい）の可能性は低いと結論づけられる
２．航空機利用客のデータセット
航空機利用客のデータセットを使用（103904×25）を使用
※注意１： このデータセットでは航空機のクラスが Business、Eco、EcoPlus と３種類ありますが、今回は Business、Eco の２種類に絞りました
※注意２： 今回は １．ロイヤル顧客 or 非ロイヤル顧客、２．ビジネスクラス or エコノミークラス、３．サービスに満足 or ニュートラルまたは不満足 の２値データとして分析しました
わかったこと
- 「ロイヤル顧客ほどビジネスクラスを選ぶ傾向がある」とはいえない
- 「ロイヤル顧客ほど満足する傾向がある」とはあまりいえない
- 「ビジネスクラスを選んだ人ほど満足する傾向がある」といえる
- サービスに対する顧客の５段階満足度表（９～２２列目）は、顧客の全体的な満足度を表す結果として一貫性があり信頼できる
サンプルサイズの大きさ
このデータセットはサンプルサイズが大きいため、全体的に検定結果が有意になりやすい
airplane.R（顧客タイプ-クラス）
> table(customerType, class)
class
customerType Business Eco
disloyal Customer 7356 10910
Loyal Customer 42309 35835
> cor(customerType01, class01)
[1] 0.1087623
- 検定方法
- クロス集計表とファイ係数
- 結果
- ファイ係数は約 0.109 なので、「ロイヤル顧客ほどビジネスクラスを選ぶ傾向がある」とはいえない
airplane.R（顧客タイプ-満足度）
> table(customerType, satisfaction)
satisfaction
customerType neutral or dissatisfied satisfied
disloyal Customer 13830 4436
Loyal Customer 39399 38745
> cor(customerType01, satisfaction01)
[1] 0.1993428
- 検定方法
- クロス集計表とファイ係数
- 結果
- ファイ係数は約 0.199 なので、「ロイヤル顧客ほど満足する傾向がある」とはあまりいえない
airplane.R（クラス-満足度）
> table(class, satisfaction)
satisfaction
class neutral or dissatisfied satisfied
Business 15185 34480
Eco 38044 8701
> cor(class01, satisfaction01)
[1] 0.5106618
- 検定方法
- クロス集計表とファイ係数
- 結果
- ファイ係数は約 0.511 なので、「ビジネスクラスを選んだ人ほど満足する傾向がある」といえる
airplane.R（クラス-満足度のカイ二乗検定）
> chisq.test(table(class, satisfaction), correct = F)
Pearson's Chi-squared test
data: table(class, satisfaction)
X-squared = 25141, df = 1, p-value < 2.2e-16
- 検定方法
- カイ二乗検定
- 分布
- 統計量は、帰無仮説のもとで自由度 (行の数-1)×(列の数-1) の カイ二乗分布 にしたがう（今回は自由度１）
- 帰無仮説
- ２つの変数（クラス - 満足度）は独立である
- 対立仮説
- ２つの変数には連関がある
- 結果
- p値が2.2e-16以下なので帰無仮説は棄却される
- すなわち、２つの変数が独立である可能性は低いと結論づけられる
※注意３： このデータはサンプルサイズが大きいので、検定の結果が有意になりやすくなっている
airplane.R（クラス-満足度のフィッシャーの正確確率検定）
> fisher.test(table(class01, satisfaction01), alternative="greater")
Fisher's Exact Test for Count Data
data: table(class01, satisfaction01)
p-value < 2.2e-16
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
9.681299 Inf
sample estimates:
odds ratio
9.927851
↑ 結果の見方はカイ二乗検定と同じだが、「フィッシャーの正確確率検定」はサンプルサイズが小さいときや 2×2 データのときに有効である
今回はサンプルサイズが大きいこともあり、カイ二乗検定と同様結果は有意（p-value < 2.2e-16）になっている
airplane.R（質問紙尺度）
> airplane <- read.csv("airplane.csv")
> library(psy)
> expairplane <- airplane[,9:22]
> cronbach(expairplane)
$sample.size
[1] 103904
$number.of.items
[1] 14
$alpha
[1] 0.7724205
- 検定方法
- 質問紙尺度
- 結果
- alpha が0.7以上あるため信頼性が高く、測定項目間に一定の一貫性がある可能性が高い。すなわち、項目が一貫して「その顧客から見た航空サービスの全体的な満足度」を測定できていると考えられる
airplane.R（因子分析）
> cor(airplane[,9:22])
> fa.parallel(airplane[,9:22])
> factanal(airplane[,9:22], factors = 5)
Call:
factanal(x = airplane[, 9:22], factors = 5)
Uniquenesses:
Wifi Departure.ArrivalTimeConvenient
0.326 0.553
EaseOfOnlineBooking GateLocation
0.190 0.541
Food.Drink OnlineBoarding
0.414 0.322
SeatComfort Entertainment
0.337 0.069
On.boardService LegRoomService
0.494 0.749
BaggageHandling CheckinService
0.420 0.696
InflightService Cleanliness
0.357 0.268
Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5
Wifi 0.144 0.460 0.648 -0.110
Departure.ArrivalTimeConvenient 0.650 0.127
EaseOfOnlineBooking 0.622 0.641
GateLocation 0.672
Food.Drink 0.763
OnlineBoarding 0.296 0.715 0.271
SeatComfort 0.764 0.157 0.230
Entertainment 0.807 0.500 -0.160
On.boardService 0.694 0.105
LegRoomService 0.481 0.124
BaggageHandling 0.750 0.121
CheckinService 0.102 0.234 0.485
InflightService 0.793 0.101
Cleanliness 0.841 0.124
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 2.647 2.244 1.482 1.411 0.480
Proportion Var 0.189 0.160 0.106 0.101 0.034
Cumulative Var 0.189 0.349 0.455 0.556 0.590
Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 5652.94 on 31 degrees of freedom.
The p-value is 0
- 検定方法
- 因子分析
- 結果
- さきほどの５段階満足度表に関して、顧客の選択の背景にどのような因子があるのか調べた。その結果の解釈として、以下のような５因子に分類されたのではと考えた
| 名前 | 説明 |
|---|---|
| １：席の快適さ | フードドリンク、席の落ち着き、娯楽、清潔さなど |
| ２：機内サービス | 娯楽、搭乗サービス、レッグルームサービス、手荷物の取り扱い、機内サービスなど |
| ３：空港での快適さ | 出発/到着時間の満足度、オンライン予約、ゲートの位置など |
| ４：システムの快適さ | Wi-fi、オンライン予約、オンライン搭乗など |
| ５：その他 | その他いろいろ |
３．ウォルマートのデータセット
ウォルマートのデータセット（6435×8）を使用
わかったこと
- 2010/02/05～2012/10/26までの約２年半で、売上は95%の確率で31830～155130（ドル）上昇した
- クリスマスやブラックフライデーに売上上昇する店舗が多い
- 大晦日に売上低下する店舗が多い（しかし1月や9月に売上低下する店舗が多いのはなぜだろう）
- ホリデーがある週/ない週で売上に有意差がある確率は、ランダムによるが10回中3回ぐらいである
walmart.R（約２年半での売上の上昇）
> feb2010 <- walmart[walmart$Date=="05-02-2010",]
> oct2012 <- walmart[walmart$Date=="26-10-2012",]
> t.test(feb2010$Weekly_Sales, oct2012$Weekly_Sales, paired = T)
Paired t-test
data: feb2010$Weekly_Sales and oct2012$Weekly_Sales
t = 3.0559, df = 44, p-value = 0.003805
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
31830.26 155130.81
sample estimates:
mean difference
93480.54
- 検定方法
- 対応のあるt検定
- 分布
- 統計量は、帰無仮説のもとで df=n-1 の t分布 にしたがう
- 帰無仮説
- μ=0（売上の変化の母平均は０である）
- 対立仮説
- μ≠0（売上の変化の母平均は０ではない）
- 結果
- p値が0.003805なので帰無仮説は棄却される
- また、本当の母平均の差は95%の確率で31830.26～155130.81の間にある
- すなわち、２つの母平均の差が０（母平均が等しい）の可能性は低い（＝売上は増えた）と結論づけられる
walmart.R（ホリデーがある週/ない週で売上に有意差があるか）
# データ整理を省略しているが、↓ の分析では店舗番号16、29、37の３店舗について、
# １店舗ごとにホリデー10週＆非ホリデー10週の売上数値を利用している
> sales_figure <- c(holidays_data$Weekly_Sales, not_holidays_data$Weekly_Sales)
> holi_not_holi <- factor(c(rep("Holi", 30), rep("notHoli", 30)))
> store123 <- factor(rep(c(rep("16",10),rep("29",10),rep("37",10)),2))
> summary(aov(sales_figure~holi_not_holi*store123))
Df Sum Sq Mean Sq F value Pr(>F)
holi_not_holi 1 4.703e+10 4.703e+10 6.226 0.0157 *
store123 2 1.583e+10 7.916e+09 1.048 0.3577
holi_not_holi:store123 2 3.634e+10 1.817e+10 2.405 0.0998 .
Residuals 54 4.079e+11 7.554e+09
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- 検定方法
- 二元配置分散分析（対応なし）
- 効果
- ・要因Aの主効果：休日週か非休日週かで売上に差がある
- ・要因Bの主効果：店舗によって売上に差がある
- ・要因Aと要因Bの交互作用効果
- 帰無仮説
- 主効果/交互作用効果がない
- 対立仮説
- 主効果/交互作用効果がある
- 結果
- 今回は要因Aの主効果のp値が0.0157、Bの主効果が0.3577、交互作用効果が0.0998なので、要因Aの主効果のみ帰無仮説は棄却される
- すなわち、有意差の有無でいうと、要因Aの主効果〇、要因Bの主効果×、要因Aと要因Bの交互作用効果×である
※注意１： しかしこれはどの非休日週をランダムに抽出するかに依り、今回は要因Aの主効果に有意差がある結果となったが、繰り返し行うと有意差がある確率はだいたい30%ぐらいである
４．バレンタインのデータセット
バレンタインのデータセット（6×9）を使用
わかったこと
- 対象商品カテゴリ（菓子、花、衣類など）や購入者の年代層によって、支出額に差がある
valentine.R（対象商品や年代による支出額の差）
# データ整理を省略しているが、↓ の分析では、
# spendingsに支出額、categoriesに対象商品、agesに年代を代入している
# 一元配置分散分析（対応なし）
> summary(aov(spendings~categories))
Df Sum Sq Mean Sq F value Pr(>F)
categories 6 5951 991.8 16.83 5.01e-09 ***
Residuals 35 2063 58.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 一元配置分散分析（対応あり）
> summary(aov(spendings~categories+ages))
Df Sum Sq Mean Sq F value Pr(>F)
categories 6 5951 991.8 29.358 2.82e-11 ***
ages 5 1049 209.9 6.212 0.000458 ***
Residuals 30 1013 33.8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> TukeyHSD(aov(spendings~categories))
# （結果省略）
- 検定方法
- 一元配置分散分析（対応なし/あり）
- 分布
- 統計量は、帰無仮説のもとで 分子の自由度（群の数-1）と分母の自由度（各群のデータ数-1）の F分布 にしたがう
- 帰無仮説
- 各商品カテゴリの平均支出額の母平均は等しい
- 対立仮説
- 各商品カテゴリの平均支出額の母平均は等しくない
- 結果
- 対応なし/ありに関わらず、p値が5.01e-09、2.82e-11、0.000458なので、商品カテゴリ間でも年代間でも、95%の確率で少なくとも１か所は群間で有意差がある。どの群の間で有意差があるかは TukeyHSD() で確認できる
疑問
今回は年代が６種類ある（6×8）一方、性別は２種類しかない（2×8）。そのため年代ごとの支出額のデータセットはある程度うまく分析できているようであるが、性別ごとの方は数値はほとんど同じなのに有意差なしとなっている。
そんなにサンプルサイズの影響があっていいものなのか？
５．車の価格予測のデータセット
車の価格予測のデータセット（22000×4）を使用
わかったこと
- 製造年と走行マイル数から価格予測をするのはあまり精度が良くない
car.R（回帰分析による価格予測）
# 製造年と走行マイルから予測
# Multiple R-squared: 0.2412なので、値段の24%ぐらいを説明できている
> data <- lm(price~year+miles)
> summary(data)
# （結果省略）
# 値段が25000以上のもの、製造年が2024以上のものを外れ値とみなし除いてみる
# 製造年と走行マイルから予測
# Multiple R-squared: 0.4035なので、値段の40%ぐらいを説明できている
> data2 <- lm(price2~miles2+year2)
> summary(data2)
Call:
lm(formula = price2 ~ miles2 + year2)
Residuals:
Min 1Q Median 3Q Max
-8607.5 -1799.6 -201.3 1771.4 9615.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.427e+05 1.525e+04 -55.26 <2e-16 ***
miles2 -4.701e-02 8.590e-04 -54.72 <2e-16 ***
year2 4.288e+02 7.553e+00 56.77 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2544 on 17434 degrees of freedom
Multiple R-squared: 0.4035, Adjusted R-squared: 0.4035
F-statistic: 5897 on 2 and 17434 DF, p-value: < 2.2e-16
# 回帰係数の等質性
summary(aov(price2~miles2*year2))
# （結果省略）
- 検定方法
- 回帰分析
- 結果
- 閾値を超えるデータを外れ値とみなし除いたところ、価格の40%ぐらいを製造年と走行マイル数で説明できるという結果になった。
- しかし回帰係数の等質性の検定により、miles2:year2のp値が2e-16以下という結果になったため、回帰線が平行でない（非線形）と考えられる。すなわち、そもそもこの回帰分析の信ぴょう性が低い疑いがある
６．A/Bテストのデータセット
A/Bテストのデータセット（90189×5）を使用
わかったこと
- 「バージョン30の人ほど1日後帰還率が高い」とはいえない
- 「バージョン30の人ほど7日後帰還率が高い」とはいえない
- 「1日後帰還率が高いほど7日後帰還率が高いがある」と少しいえる
ab.R（クロス集計表とファイ係数）
# クロス集計表とファイ係数（バージョン-1日後帰還率）
# ファイ係数は0.00594072なので、「バージョン30の人ほど1日後帰還率が高い」とはいえない
table(ver, ret1)
cor(ver01, ret1_01)
# クロス集計表とファイ係数（バージョン-7日後帰還率）
# ファイ係数は0.01053681なので、「バージョン30の人ほど7日後帰還率が高い」とはいえない
table(ver, ret7)
cor(ver01, ret7_01)
# クロス集計表とファイ係数（1日後帰還率-7日後帰還率）
# ファイ係数は0.3274012なので、「1日後帰還率が高いほど7日後帰還率が高いがある」と少しいえる
table(ret1, ret7)
cor(ret1_01, ret7_01)
# フィッシャーの正確確率検定
fisher.test(table(ret1_01, ret7_01), alternative="greater")
# p値が2.2e-16以下なので、1日後帰還率が高いほど7日後帰還率が高い傾向があると結論づけられる
↑ ２．航空機利用客のデータセットと同様に、クロス集計表とファイ係数による分析を行っている
７．主要都市の気温のデータセット
主要都市の気温のデータセット（2906327×8）を使用
わかったこと
- アメリカ/中国の主要５都市の、1995年8月と2019年8月の気温を比較することで、地域群によって気温変化の仕方に有意差があるとわかった
car.R（回帰分析による価格予測）
> # 共分散分析（lm）
> summary(lm(t2019~t1995+countryName))
# （略）
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.27498 3.94763 10.456 < 2e-16 ***
t1995 0.50202 0.04969 10.104 < 2e-16 ***
countryNameUS -3.43974 0.51393 -6.693 9.1e-11 ***
# （略）
# 回帰係数の等質性
> summary(aov(t2019~t1995*countryName))
# （結果省略）
- 検定方法
- 共分散分析
- 結果
- アメリカと中国の切片の差は-3.43974（アメリカの方が気温が低い）で、p値が9.1e-11のためその差は統計的に有意である
- また回帰係数の等質性の検定により、今回はt1995:countryNameのp値が0.0686なので、5%水準で回帰線が平行であると結論づけられる（＝共分散分析の結果に信ぴょう性があるといえる）
おわりに
R の使い方やデータの見方を学ぶことができて、とても勉強になった。
しかし基礎が足りないと感じた場面も多く、「これで合っているのか？そもそもどういう知見を得られるのが有益なのか？何かが分かったと言っても、それは少し考えれば自明なのでは…。」など、正解のなさや活用方法の不明さから、トータルとして難しいという印象であった。
参考にさせていただいた記事