こんにちは。やまもとです。
統計分析をしていて、検定で有意になった相関係数だけを計算したい場合があります。というか、自分が必要でした。
Excelでも相関係数は計算できますが、検定は意外と面倒です。しかも、偏相関係数は計算できませんし、ましてや検定はできません。
そこで、pythonでサクッとプログラムしてみたのでコードを残しておこうと思います。
有意な相関だけを取り出す方法
有意な相関関係だけを取り出す方法は、下記のブログで説明されているので、これを参考にしました。
参考：【Python入門】pandasとscipyで相関係数とp値を算出し統計分析する方法【実践】
import itertools as itr
import pandas as pd
import scipy as sp
import numpy as np
# 有意な相関のリスト
# df : pandas dataFrame format
# alpha : max of p-value
def corlist(df,alpha=0.05):
clist=[]
for i,j in itr.combinations(df,2):
x = df.loc[:,[i]].values
y = df.loc[:,[j]].values
r,p = sp.stats.pearsonr(np.ravel(x),np.ravel(y))
if p < alpha :
cpair = { 'first':i, 'second':j, 'r':r, 'p':p }
clist.append(cpair)
return clist
戻り値は、相関のある要素名がclist[i].firstとclist[I].secondに格納され、相関係数はclist[i].rに、p値はclist[i].pに格納されて戻されます。
入力値のdfはpandas.DataFrame型を想定しています。また、alphaは有意水準の閾値で通常は0.05にしておけば問題ありません。
有意な偏相関だけを取り出す方法
相関係数は、本当は関係がないのに関係あると計算されてしまう（疑似相関の）場合があります。疑似相関の疑いがある場合には、疑似相関をある程度排除した偏相関係数を計算します。しかし、偏相関はライブラリに含まれていないので、自分でプログラミングする必要があります。
幸運なことに、すでに多くの人が偏相関係数を計算するプログラムを公開してくれているので、係数計算はあまり心配しなくても大丈夫です。
他にも検索すれば、多くの記事が見つかると思います。
でも、有意性を判定できるp値を同時に計算するpythonコードは、自分が探した限りは見つかりませんでした。
偏相関係数のp値付きで計算する方法
R言語で、p値付きの偏相関係数を計算する方法は、下記の記事で紹介されています。
参考：偏相関係数
上記のページで紹介されている方法から、p値は次のように計算していることが分かります。
- 偏相関の検定統計量tを計算する
- 自由度をn-3にして、正規分布を使って、t値をp値に変換する
偏相関の検定統計量tの計算式は、偏相関係数の検定で下記のように説明されています。
ここで、添字qは固定する変数の数ですが、今回の偏相関では1個の変数を固定する場合（１次偏相関の場合）を考えるので、q=1にします。そのため、自由度をn-3にしています。
また、計算したt値からp値を求めるには、scipy.stats.t.cdf(-t,df)を使って算出することができます。ここで、dfは自由度のことです。
これらを、コードにまとめると次のようになりました。
import itertools as itr
import pandas as pd
import scipy as sp
import numpy as np
#p値つき編相関係数行列の計算
# df : pandas.DataFrame
def pcorr(df):
cor = df.corr()
n = df.shape[0] # number of data
l = df.shape[1] # number of labels
# 偏相関係数の計算
cor_inv = sp.linalg.pinv(cor) # 従属変数が存在する場合に備えて、擬逆行列pinvを使用
denom = 1 / np.sqrt(np.diag(cor_inv))
numer = np.repeat(denom,l).reshape(l,l)
pcor = (-cor_inv) * denom * numer
np.fill_diagonal(pcor,0) #ゼロ割防止のため対角成分を０にする
# 統計検定量(t)の計算
t = np.abs(pcor) * np.sqrt(n-3) / np.sqrt(1-pcor*pcor)#３変数間の相関なのでn-3
np.fill_diagonal(pcor,1) #対角成分を自己相関値=1にする
# p値の計算
p = 2*sp.stats.t.cdf(-t,n-3) #両側、Array型
# pandas.DataFrameへ変換
pcor_df = pd.DataFrame(data=pcor,index=cor.index,columns=cor.columns)
p_df=pd.DataFrame(data=p,index=cor.index,columns=cor.columns)
return pcor_df,p_df
偏相関係数の計算で、相関行列の逆行列を計算する必要がありますが、一般には逆行列が存在しない場合があるため、逆行列を計算するscipy.linalg.inv()ではなく擬逆行列を計算するscipy.linalg.pinv()を使用しています。
cor_invの計算によって、データ型がpandas.DataFrame型からArray型に変わってしまうので、最後にpandas.DataFrame型に変換しています。
戻り値は、pandas.DataFrame型に変換した偏相関係数行列pcor_dfとp値行列p_dfを返します。
有意な偏相関を取り出す
上記の関数を使えば、次のようなコードで有意な偏相関関係だけを取り出すことができます。
import pandas as pd
import scipy as sp
import numpy as np
#有意な偏相関係数リスト
# df : pandas dataFrame format
# alpha : max of p-value
def pcorlist(df,alpha):
r,p = pcorr(df)
clist = []
n = df.columns.size
j = 0
while j < n:
i = 0
while i < j:
if( p.iloc[i,j] < alpha ):
first = r.columns[i]
second= r.columns[j]
rvalue= r.iloc[i,j]
pvalue= p.iloc[i,j]
dat = { 'first':first, 'second':second, 'r':rvalue, 'p':pvalue }
clist.append(dat)
i+=1
j+=1
return clist
戻り値は、上記のcorlist関数と全く同じ形式で返されます。
おまけ
上記のcorlist関数やpcorlist関数の戻り値clistを出力する関数のコード例を書いておきます。
# 有意な相関係数リストを表示する
def print_clist(clist):
for item in clist:
print("r=%.2f, p=%.4f : %-16s <-> %-16s"%(item['r'],item['p'],item['first'],item['second']))
まあ、確認用ですね。
