Rによる統計分析ガイド(Windows対応版)

概要

Rでは、記述統計量の計算、クロス集計、各種検定、グラフ描画を行える。R標準のパッケージ(base/stats/grDevices)による基本処理から、momentsパッケージによる分布解析まで対応する。データはデータフレーム(行と列からなる2次元の表)を基盤として扱う。

本資料の前提と準備

本資料はWindowsパソコンにインストールしたRで動作することを前提とする。利用前に次を準備する。

本資料のコード例は、各節を上から順に実行することを想定している。節単独で実行する場合は、その節で使うデータ作成コードも実行する。

統計的検定は次の流れで進める。本資料の各節はこの流れのいずれかの段階に対応する。

① データの確認・記述統計 → ② 正規性の確認 → ③ 手法の選択(正規性を棄却しない=パラメトリック/正規性を棄却する=ノンパラメトリック)→ ④ 検定の実行 → ⑤ 結果(p値)の解釈

目次

  1. 記述統計量の概要
  2. SPSS・Python・Rの処理対応表
  3. Rによる分析の実装例
  4. 用語説明

記述統計量の概要

記述統計量は、データの特徴を数値で要約する指標である。データ全体の特性を把握し、分析手法の選択に用いる。上記の流れの①に対応する。

主な記述統計量には、平均値、標準偏差、中央値、四分位数、最大値、最小値、分散、歪度、尖度がある。Rのsd関数とvar関数は、n − 1 で割る不偏推定に基づく標本標準偏差と標本分散を計算する。各語の定義は「用語説明」を参照する。

SPSS・Python・Rの処理対応表

次の表は、SPSS、Python、Rにおける主要な統計処理の対応を示す。

処理内容 SPSS Python (pandas/scipy) R
記述統計量 FREQUENCIES df.describe()stats.skew()stats.kurtosis() summary()sd()var()skewness()kurtosis()
頻度表 FREQUENCIES value_counts() table()
クロス集計表 CROSSTABS pd.crosstab() table()
集約 AGGREGATE groupby().agg() aggregate()
Welchのt検定 T-TEST stats.ttest_ind(equal_var=False) t.test(var.equal = FALSE)
一元配置分散分析 ONEWAY stats.f_oneway() oneway.test()
Mann-Whitneyの検定 NPAR TESTS /M-W= stats.mannwhitneyu() wilcox.test()
Kruskal-Wallis検定 NPAR TESTS /K-W= stats.kruskal() kruskal.test()

Rによる分析の実装例

次のデータセットを使用する。

科目 受講者 得点
国語A90
国語B80
算数A95
算数B90
理科A80

データフレームの作成

ベクトルでデータを定義し、data.frame関数でデータフレームに変換する。stringsAsFactors = FALSEを指定し、文字列を文字型として扱う。

# データの作成
科目 <- c("国語", "国語", "算数", "算数", "理科")
受講者 <- c("A", "B", "A", "B", "A")
得点 <- c(90, 80, 95, 90, 80)

# データフレームの作成
df <- data.frame(科目 = 科目, 受講者 = 受講者, 得点 = 得点,
                 stringsAsFactors = FALSE)
print(df)

個別の統計量の計算

R標準の関数で、基本統計量を個別に算出する。sd関数は標本標準偏差、var関数は標本分散を返す。

scores <- df$得点

print(paste("平均値:", mean(scores)))
print(paste("標準偏差:", sd(scores)))
print(paste("分散:", var(scores)))
print(paste("中央値:", median(scores)))
print(paste("最大値:", max(scores)))
print(paste("最小値:", min(scores)))
print(paste("第1四分位数:", unname(quantile(scores, 0.25))))
print(paste("第3四分位数:", unname(quantile(scores, 0.75))))

summary関数とカテゴリ別集計

summary関数は複数の統計量を一括で算出する。カテゴリ別の代表値を求める場合はaggregate関数に単一の関数(平均など)を指定すると、結果がデータフレームとして返る。複数の統計量をカテゴリ別に見る場合はtapply関数を使う。

# 全体の統計量
print(summary(df$得点))

# 科目別の平均得点(結果はデータフレーム)
print(aggregate(得点 ~ 科目, data = df, FUN = mean))

# 科目別に複数の統計量を見る場合
print(tapply(df$得点, df$科目, summary))

グラフによる可視化

R標準のグラフィック機能でデータの分布を可視化する。箱ひげ図は分布の概要を、ヒストグラムは頻度分布を示す。png関数で画像ファイルへの出力を開始し、dev.off関数で出力を終了する。

Windowsで日本語ラベルが文字化けしないようにするには、windowsFonts関数でWindows標準の日本語フォント(例:MS Gothic)を登録する。画像デバイスを開いた後、par(family = "JP")でそのフォントを使用する。

# Windows標準の日本語フォントを登録(名前 JP として登録)
windowsFonts(JP = windowsFont("MS Gothic"))

# 箱ひげ図
png("score_distribution.png", width = 800, height = 600, type = "windows")
par(family = "JP")
boxplot(得点 ~ 科目, data = df,
        main = "科目別得点分布",
        xlab = "科目", ylab = "得点")
dev.off()

# ヒストグラム
png("score_histogram.png", width = 800, height = 600, type = "windows")
par(family = "JP")
hist(df$得点, breaks = 5,
     main = "得点分布",
     xlab = "得点", ylab = "頻度")
dev.off()

歪度と尖度

momentsパッケージで分布の形状を計算する。歪度が0に近いほど左右対称である。尖度は、正規分布と比べた分布の裾の重さや尖り具合を表す指標である。パッケージ未導入の場合はinstall.packages("moments")でインストールする。

尖度の定義はパッケージにより異なる。momentsパッケージのkurtosis正規分布で値が3になる定義である。一方、e1071パッケージなど過剰尖度(kurtosis − 3)を返す実装では正規分布で0になり、値が3だけ小さくなる。どちらの定義かを確認して比較する。

# 初回のみ: install.packages("moments")
library(moments)

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布に従う乱数を10000個生成
normal_data <- rnorm(10000)

print(paste("歪度:", round(skewness(normal_data), 3)))
print(paste("尖度:", round(kurtosis(normal_data), 3)))

クロス集計表

table関数で2つのカテゴリ変数の関係を集計する。組み合わせごとの度数を表形式で示す。ここでは、前節までのdfを上書きしないよう、クロス集計用のデータフレーム名をcross_dfとする。

グループ1 <- c("a", "b", "c", "a", "b")
グループ2 <- c("d", "d", "e", "e", "e")

cross_df <- data.frame(グループ1 = グループ1, グループ2 = グループ2,
                       stringsAsFactors = FALSE)

print(table(cross_df$グループ1, cross_df$グループ2))

検定手法の選択

統計的検定では、データの特性に応じて手法を選択する(流れの③)。まず各群の正規性を確認し、その結果に応じて次のように選ぶ。

ここで扱う例は、互いに独立した群の比較である。対応のあるデータ(同じ人を前後で測る場合など)では、対応のある検定を選ぶ必要がある。

正規性の検定

shapiro.test関数でデータの正規性を検定する(流れの②)。p値が有意水準(通常0.05)以上の場合、正規性を棄却しない。この結果に基づき、続く節のパラメトリック検定(t検定・分散分析)またはノンパラメトリック検定(Mann-Whitney検定・Kruskal-Wallis検定)を選択する。

使用上の注意:shapiro.testは標本サイズが3〜5000の範囲でのみ実行できる。大標本では、実用上問題のない小さな逸脱でも正規性なしと判定されやすい。検定結果だけに頼らず、ヒストグラムやQ-Qプロット(理論分布と比べる散布図。qqnorm関数で描画)も併用する。

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布に従う乱数を100個生成
normal_sample <- rnorm(100)

result <- shapiro.test(normal_sample)

print(paste("検定統計量:", round(result$statistic, 3)))
print(paste("p値:", format.pval(result$p.value, digits = 3, eps = 0.001)))

# Q-Qプロット
windowsFonts(JP = windowsFont("MS Gothic"))

png("normal_qqplot.png", width = 800, height = 600, type = "windows")
par(family = "JP")
qqnorm(normal_sample, main = "Q-Qプロット")
qqline(normal_sample, col = "red")
dev.off()

# 解釈:
# p値が0.05以上なら、正規性を棄却しない。
# p値が0.05未満なら、正規性を棄却する。

Welchのt検定(2群・正規性を棄却しない場合)

t.test関数で2群間の平均値の差を検定する。Welchのt検定は等分散を仮定しないため、2群の分散が異なる場合にも適用できる。t.test関数では、var.equal = FALSEを指定するとWelchのt検定を実行する。

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布に従う2群のデータを生成
group1 <- rnorm(100, mean = 50, sd = 10)
group2 <- rnorm(100, mean = 55, sd = 15)

test_result <- t.test(group1, group2,
                      paired = FALSE,
                      var.equal = FALSE)

print(paste("t値:", round(test_result$statistic, 3)))
print(paste("p値:", format.pval(test_result$p.value, digits = 3, eps = 0.001)))

# 解釈:
# p値が0.05未満なら、2群の平均値に有意差があると判断する(帰無仮説を棄却)。

Mann-Whitneyの検定(2群・正規性を棄却する場合)

2群の正規性を棄却する場合は、wilcox.test関数でMann-Whitneyの検定(2標本のWilcoxon順位和検定)を行う。母集団の正規性を仮定せず、順位に基づいて2群の差を評価する。2群の分布形状が大きく異ならない場合は、中心位置の差として解釈できる。

本例では、正規分布ではない指数分布のデータを使う。また、exact = FALSEを指定し、正規近似によるp値を計算する。標本サイズが大きい場合や同順位を含むデータでも、同じ計算方式で実行できる。

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布ではない2群のデータを生成
group1 <- rexp(60, rate = 1.0)
group2 <- rexp(60, rate = 0.5)

test_result <- wilcox.test(group1, group2,
                           paired = FALSE,
                           exact = FALSE)

print(paste("W値:", round(test_result$statistic, 3)))
print(paste("p値:", format.pval(test_result$p.value, digits = 3, eps = 0.001)))

# 解釈:
# p値が0.05未満なら、2群の順位分布または中心位置に有意差があると判断する。

一元配置分散分析(3群以上・正規性を棄却しない場合)

oneway.test関数で3群以上の平均値の差を検定する。var.equal = FALSEを指定すると等分散を仮定しないWelch型の一元配置分散分析を実行し、var.equal = TRUEを指定すると等分散を仮定する通常の一元配置分散分析を実行する。有意差が認められた場合、事後検定(どの群間に差があるかを特定する検定。多重比較ともいう)を行う。

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布に従う3群のデータを生成
group_a <- rnorm(30, mean = 50, sd = 10)
group_b <- rnorm(30, mean = 55, sd = 12)
group_c <- rnorm(30, mean = 60, sd = 9)

# 値とグループ名を1つのデータフレームにまとめる
# factor() はグループ名をカテゴリ変数として扱うための関数
anova_df <- data.frame(
  値 = c(group_a, group_b, group_c),
  グループ = factor(c(rep("A", 30), rep("B", 30), rep("C", 30)))
)

result <- oneway.test(値 ~ グループ,
                      data = anova_df,
                      var.equal = FALSE)

print(paste("F値:", round(result$statistic, 3)))
print(paste("p値:", format.pval(result$p.value, digits = 3, eps = 0.001)))

# 解釈:
# p値が0.05未満なら、少なくとも1組の群間に平均値の有意差があると判断する。
# どの群間かは事後検定で特定する。

Kruskal-Wallis検定(3群以上・正規性を棄却する場合)

3群以上の正規性を棄却する場合は、kruskal.test関数でKruskal-Wallis検定を行う。一元配置分散分析のノンパラメトリック版であり、順位に基づいて3群以上の差を評価する。有意差が認められた場合、ノンパラメトリックの事後検定(例:pairwise.wilcox.test関数)でどの群間に差があるかを調べる。

本例では、正規分布ではない指数分布のデータを使う。

# 乱数シードの設定(同じ値を指定すると同じ乱数列を再現できる)
set.seed(42)

# 正規分布ではない3群のデータを生成
group_a <- rexp(30, rate = 1.0)
group_b <- rexp(30, rate = 0.7)
group_c <- rexp(30, rate = 0.4)

# 値とグループ名を1つのデータフレームにまとめる
# factor() はグループ名をカテゴリ変数として扱うための関数
kw_df <- data.frame(
  値 = c(group_a, group_b, group_c),
  グループ = factor(c(rep("A", 30), rep("B", 30), rep("C", 30)))
)

result <- kruskal.test(値 ~ グループ, data = kw_df)

print(paste("カイ二乗値:", round(result$statistic, 3)))
print(paste("p値:", format.pval(result$p.value, digits = 3, eps = 0.001)))

# 解釈:
# p値が0.05未満なら、少なくとも1組の群間に順位分布の有意差があると判断する。

用語説明