こんにちは。やまもとです。
統計分析をしていて、検定で有意になった相関係数だけを計算したい場合があります。というか、自分が必要でした。
Excelでも相関係数は計算できますが、検定は意外と面倒です。しかも、偏相関係数は計算できませんし、ましてや検定はできません。
そこで、pythonでサクッとプログラムしてみたのでコードを残しておこうと思います。
（2024-01-01：自由度の記述を訂正）
有意な相関だけを取り出す方法
有意な相関関係だけを取り出す方法は、下記のブログで説明されているので、これを参考にしました。
参考：【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は固定する変数の数で、n-q-2は自由度dfです。上記の参考ページでは、１個の変数を固定した偏相関計数を考えていたので、自由度df=n-1-2=n-3になっていました。
今回は、l個のデータセットの偏相関行列を相関行列の逆行列を使って計算するので、偏相関を求める２データセット以外を全て固定することになります。そのため、固定数q=l-2、自由度df=n-q-2=n-lとなります。
また、計算したt値からp値を求めるには、scipy.stats.t.cdf(-t,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
dof = n - l #自由度
# 偏相関係数の計算
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(dof) / np.sqrt(1-pcor*pcor)#３変数間の相関なのでn-3
np.fill_diagonal(pcor,1) #対角成分を自己相関値=1にする
# p値の計算
p = 2*sp.stats.t.cdf(-t,dof) #両側、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']))
まあ、確認用ですね。
ｔ検定とｐ値の計算部分で自由度がn-３固定になっていますが正しいのでしょうか？
このコードは、除外する第3の変数が１つの場合の偏相関をまとめて計算するものです。そのため、各偏相関係数の計算に用いられた第3の変数は必ず１つです。結果、自由度n-q-2は、第3の変数の個数がq=1となるため、自由度はn-3に固定しています。間違っていたら、ごめんなさい。。。
お返事ありがとうございます
変数x1,x2,x3,x4があるとき、x1,x2の偏相関係数を求めるには他の変数を固定して残差の相関を求める・・プログラムでは逆行列でこれをやっている
で、この時に固定される変数は、２個なので、q=2になる
変数が８個の場合は、q=８－２なのかなと、思いましたが
その点どうでしょう？
公式の再導出や検算で時間がかかってしまいましたが、自由度の部分はこの記事が間違っていますね。
ご理解の通り、変数が４つの場合で逆行列を使って偏相関を求めた場合は、
相関を計算する２つの変数以外が全て固定され、残り２つの変数が制御変数になるため、q=2になります。
記事とコードを修正したいと思います。
教えていただき、ありがとうございました！