はじめよう実験計画

実験を早く終わらせるための技術

Q-Qプロットの作成 - ステップ解説&R実装例

Q-Qプロットとは

Q-Qプロットとは、回帰分析の残差が正規分布に従うかをチェックする図であり、Rではqqnorm()で描くことができます(図1)。

図1  Rで描いたQ-Qプロットの例。qqnorm(rnorm(100))で作成。

図1はrnorm()で生成した100個の点をプロットしたものであり、おおむねy=xの直線上に乗っていることが、正規分布に従うことを示しています。

このQ-Qプロットを残差に対してプロットして、「だいたい一直線なら回帰分析はOK」、という非常にざっくりとしたイメージしか今まで持っていなかったので、このプロットがどのように作られているかを調べてみました。

 

Q-Qプロットの作り方

Q-Qプロットは、以下のフローで作ることができます。

■Q-Qプロットの作り方(図2)■

① パーセンタイルPi=(i-0.375)/(n-2×0.375+1), i=1,2,...n(=データ数)を計算※1

② パーセンタイルPiに対応する標準正規分布のパーセント点Z scoreを計算※2

③ データを昇順に並べる

④ ③のデータを縦軸に、②のZ scoreを横軸にプロットしてQ-Qプロットの完成

 

①のパーセンタイルの公式(※1)と、②のZ scoreの求め方(※2)については後述します。

図2はサイズn=10のデータ sample = (-0.463, -0.280, -0.414, 1.62, -0.72, -0.45,  0.01,  0.22, 0.19, -0.05) に対してQ-Qプロットを作図するフローを示したものです。

 

図2 sampleデータに対してQ-Qプロットを作成するフロー。
①公式を用いてパーセンタイルPiの計算する。②正規分布において累積確率がPiになるようなパーセント点(Z score)を計算する。③プロットしたいデータを昇順に並べる。④②のZ scoreと③のデータをプロットしてQ-Qプロットの完成。

 

パーセンタイルを求める公式

①パーセンタイルの公式(※1)について補足します。

データxiに対するパーセンタイルとは、x1~xiがx全体に占める割合のことで、このパーセンタイルPiの期待値は

 E[P_i]=i/(n+1)

と書くことができます。xが一様分布に従う場合、順序統計量x(i)がベータ分布に従うことから上式を導出できます*1。xが正規分布に従う場合は修正を加えて

 E[P_i]=\frac{i-a}{n+2a-1}

に対してa=0.375とすることで、パーセンタイルの期待値が計算されます。*2

 

Z scoreの求め方

次に、②Z scoreの求め方(※2)について説明します。

パーセンタイルPiを求めたら、正規分布において累積確率がPiになるような横軸の点ZがZ scoreとなります。図3左の累積分布では、Pi=0.84(縦軸)に対応するZ score(横軸)が1であることを示しています。図3右は全く同じことを確率密度関数を使って示したものですが、Z<1となる確率が0.84であることが分かります。

図3 パーセンタイルからZ scoreを求める様子。
左図は正規分布の累積分布において、縦軸がPiとなるような横軸の値としてZ scoreを求めている。右図は正規分布の確率密度において、下側確率(Z<Z score)がPiになるような横軸の値としてZ scoreを計算している。




Rコード

本記事で作成した自作のQ-Qプロットと、Rのデフォルトの関数であるqqnorm()を比較します。

#1.パーセンタイルの計算
percentlie <- (seq(10)-0.375)/(10-2*0.375+1)
#2.Z scoreの計算
zscore <- qnorm(percentlie, lower.tail = TRUE)

#今回プロットするデータの設定
set.seed(3)
#今回プロットするデータ10点
sample <- rnorm(10)
#3.データを昇順に並べる
plot(sort(sample), pch=16)

#Rの関数qqnorm()
qqnorm(sample, xlab = "", ylab="", pch=16, xlim=c(-2,2), ylim=c(-2,2))

#4. 2と3から自作のQ-Q plotを作成して重ね書き
par(new=T)
plot(zscore, sort(sample), col="red", pch=2, xlim=c(-2,2), ylim=c(-2,2))

#y=xの直線と凡例の追加
abline(a=0,b=1)
legend(x=-2,y=2, legend=c("qqnorm()", "self-made"), pch = c(16,2), col = c("black", "red"))

以下、上記のコードの出力結果です。自作のQ-Q plotとRが用意しているqqnorm()が一致していることを確かめられました!

図4 自作のQ-Q plotと、Rの関数qqnorm()の比較

 

おまけ - Half Normal Plot

実験計画法ではHalf Normal Plotという、データの絶対値に対するQ-Qプロットをよく使います。Half Normal Plotについてはこちらの記事をご覧ください。

www.doe-get-started.com

 

参考サイト