事例を通して、R 言語を勉強しよう

はじめに

はじめまして。私は大学院生です。

私は、Python を使って、研究をしています。
今学期、TA で R 言語を使った授業を担当したので、復習を兼ねて、ここでアウトプットさせていただきます。

適当な事例を使い、初心者にもわかりやすいように善処します。

この記事が、「R 言語」の勉強の手掛かりになれば幸いです。

事例が思いついたら、追加で載せていきます。

環境構築

以下のサイトから R をダウンロードします。

ダウンロードの流れ

サイトの左側にある CRAN をクリックしてください。

CRAN ミラーは、「0-Cloud」あるいは「Japan」をリクックしてください。

その後、「Download R for」の部分をクリックすれば、ダウンロードできます。

ライブラリ

この記事で使用するライブラリです。

library(ggplot2)

データ

ここで、使用するデータは、乱数を生成させることで作るので、事前の準備などは要りません。
また、アヤメのデータセットなど簡単に使うことができるデータセットは使います。

シード

実験では、乱数を使用するので、再現性に問題があります。
そこで、シードを固定したい方は、以下のコードを最初に実行してください。

set.seed(42)

事例

学校のテストの点の分析

データの生成

まず、平均を $50$ 点、標準偏差を $15$ 点で、$120$ 人分のデータを正規分布（ガウス分布）に従い生成させます。

まず、rnorm() 関数を使って、テストの点を生成します。 rnorm() 関数で生成したデータは、少数が含まれます。そこで、round() を使って、四捨五入します。

score_list = round(rnorm(120, 50, 15))

score_list にどのような点数が入ったかは、以下のコードを実行することで、確認できます。

score_list

可視化

ヒストグラム

ヒストグラムにして可視化します。
（正規分布に従って、生成していますので、ヒストグラムは単峰性になると思います。）

hist(score_list)

軸ラベルを変更したいときは、以下のように xlab や ylab を指定してください。

hist(score_list, xlab="score", ylab="Frequency")

ヒストグラム + 確率密度関数

上で生成したヒストグラムに確率密度関数のグラフを図示します。
今回のデータは、正規分布に従っていると仮定します。

第一軸をテストの点数、第二軸を確率となるようにします。

par() 関数を使って、追加でプロットできるようにします。
plot() 関数で空のグラフを図示します。
axis() 関数で第二軸を追加します。
mtext() 関数で第二軸のラベルを追加します。
curve() 関数を使って、確率密度関数をプロットします。

hist_data = hist(score_list, xlab="score", ylab="Frequency", col="blue")
par(new=TRUE)
plot(0, 0, type="n", xlim=range(hist_data$breaks), ylim=c(0, max(density(score_list)$y)), axes=FALSE, xlab="", ylab="")
axis(side=4)
mtext("Density", side=4)
curve(dnorm(x, mean=mean(score_list), sd=sd(score_list)), add=TRUE, col="red")

箱ひげ図

今度は、箱ひげ図にして可視化します。

boxplot(score_list)

図中に丸があった場合は、外れ値のデータです。
デフォルトだと、外れ値は、第一四分位数あるいは第三四分位数から四分位範囲の $1.5$ 倍以上離れた点です。

軸ラベルを変更したいときは、以下のように ylab を指定してください。

boxplot(score_list, ylab="socre")

記述統計量

summary 関数を使うことで、最大値や四分位数などの記述統計量をまとめて出すことができます。

summary(score_list)

検定

Shapiro–Wilk 検定

Shapiro–Wilk 検定¹ は、標本の正規性（データが正規分布に従っているか）を評価する統計的検定です。

帰無仮説 H0：「データは正規分布に従う」
対立仮説 H1：「データは正規分布に従わない」

R で Shapiro–Wilk 検定を行う場合、以下を実行してください。

shapiro.test(score_list)

$p$ 値が、有意水準 $0.05$ より小さい場合、帰無仮説を棄却する。つまり、データは正規分布に従わない。
$p$ 値が、有意水準 $0.05$ より大きい場合、帰無仮説を棄却できない。つまり、データは正規分布に従う。

身長と体重の関係性

データの生成

まず、身長と体重をランダムに生成します。今回は BMI (Body Mass Index) を活用します。
BMI は、身長 ($\mathrm{m}$) と体重 ($\mathrm{kg}$) を基に計算することができます。

\mathrm{BMI} = \mathrm{体重} \div \mathrm{身長}^2

日本肥満学会の肥満度分類によると、BMI $22$ を最も疾患の少ない標準体重（理想体重）としています。

まず、身長のデータを runif() という関数を使って、最小値 $1.40 \mathrm{m}$、最大値 $2.00 \mathrm{m}$ で $100$ 人分のデータを一様分布になるように生成します。

height_list = runif(100, min=1.40, max=2.00)

体重の生成は、以下の数式に従い、生成します。
誤差項を加える理由は、標準体重に誤差を足すことでデータを散らばらせることが目的で足しています。
誤差項は、標準正規分布に従い生成した行列式 $\det(\mathbf{\epsilon})$ に身長を掛けた値とします。理由としては、身長が高い人ほど誤差項の絶対値が大きくなるようにしたいので、生成した誤差に身長の値をかけます。

\mathrm{体重} = 22 \times \mathrm{身長}^2 + \mathrm{誤差項} = \mathrm{体重} = 22 \times \mathrm{身長}^2 + \mathrm{身長} \times \det(\mathbf{\epsilon})

今回は、行列 $\mathbf{\epsilon}$ を $2$ 行 $2$ 列とし、行列式 $\det(\mathbf{\epsilon})$ を生成します。生成した $2$ 行 $2$ 列の行列 $\mathbf{\epsilon}$ の行列式 $\det(\mathbf{\epsilon})$ を det() 関数を使って，計算します。replicate() 関数で length(height_list) 回同じ処理をします。

n = 2
det_epsilon_list = replicate(length(height_list), det(matrix(rnorm(n^2), nrow=n, ncol=n)))

体重を計算します。

weight_list = 22 * height_list^2 + height_list * det_epsilon_list

ピアソンの積率相関係数

ピアソンの積率相関係数（相関係数）を調べることで、身長と体重がどの程度線形に関係があるのかを求めることができます。

cor(height_list, weight_list, method="pearson")

相関係数が $1$ に近いほど正の相関関係があり、$-1$ に近いほど負の相関関係があります。また、相関係数が $0$ に近いほど相関が無いと言えます。

可視化

散布図

散布図にして図示します。
身長と体重の関係性は、体重は身長の $2$ 乗に比例しています。

plot(height_list, weight_list, xlab="height", ylab="weight")

線形回帰

相関係数を調べたことで、身長と体重の線形関係があるかどうかを調べることができました。
そこで、線形回帰を用いて、フィッティングをしようと思います。今回は説明変数が $1$ つ（身長）だけなので、単回帰で予測します。

\mathrm{体重} = \beta + \alpha \times \mathrm{身長}

lm() 関数を使って、切片 $\beta$ と傾き $\alpha$ を求めます。

simple_regression_model = lm(weight_list ~ height_list)

モデルの詳細な要約が知りたい場合は、summary() 関数を使ってください。

summary(simple_regression_model)

この単回帰モデルで予測した直線を先ほどのグラフに重ねてみます。
まず、等差の身長のリストを新たに作成します。その後、predict() 関数で体重を予測します。

fitted_height_list = data.frame(height_list=seq(1.40, 2.00, length=100))
fitted_weight_list = predict(simple_regression_model, fitted_height_list)

元のデータの散布図と単回帰モデルによる予測を重ねて図示します。

plot(height_list, weight_list, xlab="height", ylab="weight", col="blue")
lines(fitted_height_list$height_list, fitted_weight_list, col="red")

非線形回帰

次の式にフィッティングさせます。

\mathrm{体重} = \alpha \times \mathrm{身長} ^ 2 + \beta \times \mathrm{身長} + \gamma

nls() 関数を使って、$\alpha$、$\beta$、$\gamma$ を求めます。$\alpha$、$\beta$、$\gamma$ の初期値は BMI の定義に従い、それぞれ $22$、$0$、$0$ に設定しました。

nonlinear_model = nls(weight_list ~ alpha * height_list^2 + beta * height_list + gamma, start=list(alpha=22, beta=0, gamma=0))

モデルの詳細な要約が知りたい場合は、summary() 関数を使ってください。

summary(nonlinear_model)

散布図に予測したモデルを重ねて図示します。
今回は、曲線なので、curve() 関数を使います。

plot(height_list, weight_list, xlab="height", ylab="weight", col="blue")
curve(predict(nonlinear_model, newdata=data.frame(height_list=x)), add=TRUE, col="red")

二つの会社の株価の関係性

データの生成

Kotler の競争戦略は、「リーダ」、「チャレンジャ」、「フォロワ」、「ニッチャ」の四つにマーケットシェアを分類しています。
今回は、「リーダ」と「チャレンジャ」の株価を想定し、データを生成します。（データを生成するために使う数式は現実のモデルとはかけ離れています。）

リーダ：トップシェアを維持する全方位戦略
チャレンジャ：リーダのトップシェアを奪取するための差別化戦略

誤差項を $\epsilon$ とし、実数 $k$ を $0 < k < 1$ とする。

リーダの式を以下のように定義します。

y_{l_1} = 0, \quad y_{l_2} = \epsilon_{l_2}, \quad y_{l_t} = k y_{l_{t-1}} + (1 - k) y_{l_{t-2}} + \epsilon_{l_t}

チャレンジャの式を以下のように定義します。

y_{c_1} = \epsilon_{c_1}, \quad y_{c_t} = k y_{c_{t-1}} + (1 - k) y_{l_{t-1}} + \epsilon_{c_t}

リーダとチャレンジャの式を R で実装していきます。数式に再帰が含まれているので、独自に関数を定義します。Recall() 関数を使うことで、再帰関数を実装できます。

leader = function(t, mean=0, sd=1, k=0.5) {
    if (t==1) {
        y_l = 0
    } else if (t==2){
        y_l = rnorm(1, mean=mean, sd=sd)
    } else {
        y_l = k * Recall(t-1, mean, sd, k) + (1 - k) * Recall(t-2, mean, sd, k) + rnorm(1, mean=mean, sd=sd)
    }
    return (y_l)
}

challenger = function(t, y_l, mean=0, sd=1, k=0.5) {
    if (t==1) {
        y_c = rnorm(1, mean=mean, sd=sd)
    } else {
        y_c = k * Recall(t - 1, y_l, mean, sd, k) + (1 - k) * y_l[t-1] + rnorm(1, mean=mean, sd=sd)
    }
    return (y_c)
}

作成した関数を使って、データを生成します。

t = 20
y_l_list = sapply(1:t, leader, mean=0, sd=1, k=0.9)
y_c_list = sapply(1:t, challenger, y_l=y_l_list, mean=0, sd=1, k=0.7)

可視化

plot(1:t, y_l_list, type = "o", col = "blue", xlab = "Time", ylab = "Value", ylim = range(c(y_l_list, y_c_list)))
lines(1:t, y_c_list, type = "o", col = "red")
legend("topright", legend = c("Leader", "Challenger"), col = c("blue", "red"), lty = 1, pch = 1)

検定

等分散性の検定

「リーダ」と「フォロワ」の株価の母集団の分散が等しいかどうかを検定します。

帰無仮説 H0：「株価の分散に等しい」
対立仮説 H1：「株価の分散は等しくない」

var.test(y_l_list, y_c_list)

$p$ 値が、有意水準 $0.05$ より小さい場合、帰無仮説を棄却する。つまり、等分散でない。
$p$ 値が、有意水準 $0.05$ より大きい場合、帰無仮説を棄却できない。つまり、等分散である。

t 検定

今回のデータは、「リーダ」と「フォロワ」という別の企業から得た株価なので、対応のない t 検定を使います。

帰無仮説 H0：「株価の平均に差がない」
対立仮説 H1：「株価の平均に差がある」

R で t 検定を行う場合、以下を実行してください。ここで、二つの株価（「リーダ」と「フォロワ」の株価）は異なる組織から得たデータなので、対応のない t 検定を行います。対応のない t 検定を行うためには、paired=FALSE と指定します。

等分散性を仮定した場合は、var.equal=TRUE としてください。

t.test(y_l_list, y_c_list, paired=FALSE, var.equal=TRUE)

等分散を仮定しない場合は、var.equal=FALSE としてください。等分散を仮定しない場合は、Welch の t 検定 ² を行います。

t.test(y_l_list, y_c_list, paired=FALSE, var.equal=FALSE)

$p$ 値が、有意水準 $0.05$ より小さい場合、帰無仮説を棄却する。つまり、株価の平均に差がある。
$p$ 値が、有意水準 $0.05$ より大きい場合、帰無仮説を棄却できない。つまり、株価の平均に差がある。

アヤメのクラスタリング

データの用意

「アヤメ」という機械学習でよく使われるデータセットを使います。
このデータセットには、「アヤメ」の品種の情報もありますが、この情報を使わずに学習します。つまり、教師なし学習を行います。

data = iris[, -5]

k-means 法

今回は、クラスタリングの代表的な手法の一つである k-means 法を使います。
kmeans() 関数を使って、クラスタの数は centers で指定します。このデータセットは $3$ つの品種が含まれているので、クラスタの数は $3$ にします。

result = kmeans(data, centers=3)

可視化

扱いやすい型にするためデータフレームに変換します。

iris$Cluster = as.factor(result$cluster)

各クラスタは色を変えて図示します。また、正解となる品種を点の形を変えて図示します。

ggplot(iris, aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(color=Cluster, shape=Species), size=3) + labs(x="Petal Length", y="Petal Width") + theme_minimal()

おわりに

この記事では、事例を使って、R 言語を勉強しました。

この記事が、R 言語の勉強の支えになれば幸いです。

参考文献

Samuel Sanford Shapiro and Martin B Wilk. An analysis of variance test for normality (complete samples). Biometrika, 1965. ↩
B. L. Welch. The Generalization of `Student's' Problem when Several Different Population Variances are Involved. Biometrika, 1947. ↩