Rで正規QQ-Plotについてまとめてみた

はじめに

多くの統計モデルは正規分布（ガウス分布）を想定しており、繰り返しの観測などで得られたデータが正規分布であるかどうかを確認することが必要であると言えます。
そのような確認をするツールとして正規QQ-Plotというのがあります。
諸事情でこの正規QQ-Plotを説明する資料を作成したので、その資料の備忘録としてこの記事を投稿します。

正規QQ-Plotとは

正規QQ-Plotは以下の通り手元にあるデータが正規分布かどうかを確認するツールです。

• 観測値が正規分布に従う場合の期待値をX軸に、観測値そのものをY軸にプロットしたものを正規Q-Q plotという。
• X軸の期待値は以下の手順で求める。
  • 観測値を昇順に並べた順位から、それぞれのパーセンタイル（昇順に並べ変えてどこに位置するのかをパーセントで表示したもの）を求める。
  • そのパーセンタイルに対して、正規分布の確率密度関数の逆関数(qnorm)を用いて期待値となる統計値を求める。
• 正規Q-Qプロットが対角線上に並べば、観測値は正規分布に従っていると考えられる。

正規QQ-PlotをRで描いてみる

このような正規QQ-PlotをExcelで作成することも可能なようなのですが、Rで簡単に作成することができるので、紹介していきます。

手法１：Rに組み込まれている関数を作成する

Rには正規QQ-Plotを描画するのに便利な以下の関数が準備されています。
qqnorm(): 与えられたベクトル（観測値）から正規QQ-Plotを作成する
qqline(): qqnorm()で作成した正規QQ-Plotに対角線を引く

qqnorm()でもある程度観測値の当てはまり度合いは判断がつくのですが、qqline()で描画される補助線によって、観測値の分布が正規分布に当てはまっているかどうかの判断が容易になります。

それでは以下のコードで描画してみましょう。
今回は適当に準備した12点の測定値を用いて確認をしてみます。

x <- c(3.2, 3.2, 3.1, 2.3, 2.3, 2.3, 3.1, 3.0, 4.2, 3.0, 2.9, 3.0)

#-----Q-Q plotとQ-Q lineを描画（図１）-----
qqnorm(x,  main = "QQ Plot by qqnorm function")
qqline(x, col="red")

出力結果は以下の通りです。

ここで、X軸はこの12点のデータが正規分布であったときの統計量の期待値、Y軸はベクトルxに入れた実際の観測値になります。
観測値の分布が正規分布に近ければqqline()で引いた赤線上にプロットが配置され、分布の当てはまり度合いがわかります。
こちらの例では低い値の3点と高い値の1点のデータが赤線から外れており、正規分布から外れていることがわかります。

手法2：qqnorm()とqqline()をStep-by-stepで確認してみる

続いて、qqnorm()とqqline()で描画される結果を分解してStep-by stepで確認してみます。
qqnorm()では自動で観測値xをソートしてくれるのですが、手計算で確認するにはまず観測値を昇順に配置する必要があります。
その後、データの個数で分割し、昇順で低い値から高い値のパーセンタイル(fi)を以下のスクリプトで計算します。

x.sort <- sort(x)
i <- 1:length(x.sort)
fi <- (i - 0.5)/length(x.sort)
plot(fi)

fiの出力結果が以下の通りです。

今回の12点の観測点について12等分のパーセンタイルが作成できました。
正規分布でfiとなるような統計値は確率密度関数の逆関数qnorm()で計算できます。

x.norm <- qnorm(fi)
x.norm

> x.norm
 [1] -1.7316644 -1.1503494 -0.8122178 -0.5485223 -0.3186394 -0.1046335
 [7]  0.1046335  0.3186394  0.5485223  0.8122178  1.1503494  1.7316644

この値がQQ-PlotのX軸の値になります。
この値をX軸に、Y軸に観測値をプロットすることで正規QQ-Plotが描画できます。

QQ-Plotを描画する前にこのqnorm()で得られた値がどのようなイメージかを以下のスクリプトで説明します。

x <- seq(-3, 3, 0.05)
plot(x,dnorm(x), type="l",lty=2, xlab="x", ylab="density", ylim = c(0,1))
curve(pnorm(x), type="l", lty=1, col='blue', add=T)

まず１行目の列で-3から3まで0.05ステップのベクトルを作成します。
そのベクトルを用いてdnorm()関数で正規分布の関数（黒色、破線）を描画します。
続いて、curve()関数を使って、pnorm()で計算した正規分布の確率密度関数（青色、実線）を描画します。

黒色破線の正規分布の左側からの統計量の面積が青色実線の確率密度関数となっています。
正規分布では全ての面積値が1ですので、確率密度関数は統計量が-∞で0、0で0.5、＋∞で1の値を取ります。

今回の12点の観測値のデータが正規分布と仮定する場合、各点のパーセンタイルの値は確率密度関数を用いて求めることができます。
以下の線をプロットに上書きしてみましょう。

abline(h=fi, v=x.norm, lty=4, col=1:12)

それぞれ、Y軸のパーセンタイルの値と、先ほどqnorm()関数で求めた統計量が青色実線の確率密度関数上で交わっていることが確認できます。

このようにして以下の通り正規QQ-Plotが得られます。

plot(x.sort ~ x.norm, type = "p", xlab = "Normal quantiles", pch = 20)

続いて正規QQ-Plotに描画する対角線を求めます。
対角線はX軸とY軸の第一四分位（Q1)、第三四分位（Q3）を通る直線を計算します。

#Q1、Q3の観測値を求める。Rに搭載されているqqline()関数では、
#quantile()関数のアルゴリズムを指定するオプションとしてtype = 7を使用しているので
#type=7を指定する。
y <- quantile(x.sort, c(0.25, 0.75), type = 7)

# 正規分布のQ1およびQ3の統計量を求める
x <- qnorm( c(0.25, 0.75))

#Q1, Q3のxおよびyそれぞれの値から切片と傾きを求める
slope <- diff(y) / diff(x)
int <- y[1] - slope * x[1]

#QQlineを描画する（図４）
abline(a = int, b = slope)

Rに搭載されているqqnorm(), qqline()と同じ正規QQ-Plotが描画できました。

参考サイト

本記事は以下のサイトの情報を参考にさせていただきました。

https://bellcurve.jp/statistics/glossary/2071.html
https://mgimond.github.io/ES218/Week06a.html