こんにちは。やまもとです。
前回の記事の中でSEMCausalFinderクラスの説明を省略しましたが、実はこれは構造方程式モデリング（SEM）を使った因果探索を実装したクラスでした。とあるデータを使うとLiNGAMでうまく因果探索を計算できなかったため、代替手段として作成してみました。LiNGAMの代替なので、インターフェースはLiNGAMの各クラスに似せて作っています。
今回は、SEMを使った因果探索の実装について解説したいと思います。
考え方
考え方は、シンプルです。その方法は、まず２つのデータ（データAとデータB）について、因果関係A→BとB→AをそれぞれSEM計算を行います。そして、A→BとB→Aで適合度（例えばCFI）に差があれば、より適合度の高い方の因果関係（下図では右方向）を採用、もし適合度に差がなければ、相関関係だと判定するというものです。
もちろん、これは両方の因果が有意であることが前提になります。そのため、片方の因果だけが有意な場合は、適合度を比較せずに有意な因果関係だけを採用します。
この方法は、Chat-GPTに「SEMで因果の方向を調べるにはどうしたらよいか？」と聞いた場合、２番目の方法として提示される程度には一般的な方法だと推測されます。
SEM（構造方程式モデリング）は、複数の変数間の因果関係をモデル化することができます。SEMにおいて因果の方向を調べるには、以下のような方法があります。
- 理論的仮説を基に因果関係を決定する SEMは、ある変数が別の変数に影響を与えることを表現することができます。そのため、SEMを構築する前に、理論的仮説を明確にし、因果関係を決定することが重要です。理論的仮説は、過去の研究結果、専門家の意見、あるいは理論的枠組みに基づいて設定されます。
- 因果関係の方向を変えてモデルを構築し、モデルフィットを比較する SEMでは、因果関係の方向を逆転させたモデルも構築することができます。そのため、モデルフィット指標を比較することで、因果関係の方向がより適切なものを決定することができます。一般的には、モデルフィット指標の改善が見られた方が、より適切な因果関係の方向である可能性が高いと考えられます。
- イベントスタディ分析などの手法を用いる SEMは、様々な手法と組み合わせることで、因果関係の方向をより正確に決定することができます。例えば、イベントスタディ分析では、ある変数があるイベントによって引き起こされた影響を調べることができます。これにより、因果関係の方向をより正確に決定することができます。
以上のように、SEMで因果の方向を調べるには、理論的仮説を明確にし、モデルフィット指標を比較することや、他の手法と組み合わせることが重要です。
参照：Chat-GPT
なのですが、これを全ての要素について調べようとすると、組み合わせ回数の試行が必要になりとても大変です。
そこで、上記の方法を全要素に対して総当たりで実行するプログラムを作ってしまうことにしました。
使い方
テストデータ
テストデータは、CausalAnalysisクラスで使ったのと同じ下記の公開データをCSV形式で保存して使用します。
実行方法
SEMCausalFinderクラスのインターフェースは、LiNGAMライブラリに似せて作ってあるため、使い方は至ってシンプルです。
import numpy as np
import pandas as pd
import causal_analysis as ca
data = pd.read_csv("data08-01.csv",index_col=0,header=0) # データ読み込み
finder = ca.SEMCausalFinder() # インスタンス作成
finder.fit(data) # 探索
print(finder) # 分析結果の表示
多くの機械学習ライブラリで用いられるfit関数で探索を実行します。
パラメータ設定
因果係数の有意水準alphaと、適合度CFI, GFI, NFIについて差があると判定する閾値eps_cfi, eps_gfi, eps_nfiをインスタンス作成時に指定することができます。これらの規定値は、すべて0.05に設定しています。
finder = ca.SEMCausalFinder(alpha=0.05, eps_cfi=0.05, eps_gfi=0.05, eps_nfi=0.05) # インスタンス作成
例えば、eps_cfi=0.05は、CFIの値が左右差0.05以上であれば、相関関係ではなく因果関係があると判定します。
実行結果
実行結果は、因果行列を表すadjecency_matrix_メンバ変数に行列として格納しています。LiNGAMの因果行列との違いは、相関関係も含めている点で、必ずしも三角行列にはなりません。
また、結果の出力には、途中計算で使用した因果係数、p値、各適合度も行列形式で出力するようにしました。
テストデータの結果は次のとおりです。
Cutoff for p-value : 0.05
Cutoff for CFI diff. : 0.05
Cutoff for GFI diff. : 0.05
Cutoff for NFI diff. : 0.05
Est. Stds:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.000000 0.486598 0.597399 0.242352 0.276134 0.218465
社交性 0.486477 0.000000 0.557879 0.125003 0.316861 0.172234
積極性 0.597371 0.557998 0.000000 0.240791 0.037339 0.146476
知性 0.242667 0.125067 0.240701 0.000000 0.436245 0.627193
信頼性 0.276134 0.317011 0.037409 0.436322 0.000000 0.583388
素直さ 0.218892 0.172331 0.146600 0.627123 0.583625 0.000000
p-values:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 1.000000 0.012741 0.000864 0.263937 0.198839 0.316735
社交性 0.012770 1.000000 0.002645 0.573128 0.135168 0.434251
積極性 0.000865 0.002637 1.000000 0.267216 0.867291 0.507840
知性 0.263278 0.572929 0.267407 1.000000 0.030148 0.000317
信頼性 0.198839 0.134962 0.867045 0.030113 1.000000 0.001317
素直さ 0.315744 0.433983 0.507473 0.000317 0.001308 1.000000
CFIs:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.000000 1.293761 1.146420 -0.266988 -1.429839 0.017958
社交性 1.293760 0.000000 1.183162 0.406423 9.617953 0.283428
積極性 1.146419 1.183162 0.000000 -0.240504 0.492925 0.361111
知性 -0.266992 0.406423 -0.240503 0.000000 1.449982 1.125109
信頼性 -1.429839 9.618122 0.492925 1.449982 0.000000 1.157994
素直さ 0.017952 0.283429 0.361116 1.125109 1.157995 0.000000
GFIs:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.000000 1.000000 1.000000 0.999996 0.999999 0.999995
社交性 1.000000 0.000000 1.000000 0.999997 1.000000 0.999996
積極性 1.000000 1.000000 0.000000 1.000000 0.999999 0.999998
知性 0.999998 1.000000 1.000000 0.000000 1.000000 1.000000
信頼性 0.999999 0.999999 0.999995 1.000000 0.000000 0.999999
素直さ 1.000000 0.999995 0.999980 1.000000 1.000000 0.000000
NFIs:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.000000 1.000000 1.000000 0.999996 0.999999 0.999995
社交性 1.000000 0.000000 1.000000 0.999997 1.000000 0.999996
積極性 1.000000 1.000000 0.000000 1.000000 0.999999 0.999998
知性 0.999998 1.000000 1.000000 0.000000 1.000000 1.000000
信頼性 0.999999 0.999999 0.999995 1.000000 0.000000 0.999999
素直さ 1.000000 0.999995 0.999980 1.000000 1.000000 0.000000
AICs:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.0 4.0 4.000000 3.999999 4.0 4.0
社交性 4.0 0.0 4.000000 4.000000 4.0 4.0
積極性 4.0 4.0 0.000000 4.000000 4.0 4.0
知性 4.0 4.0 4.000000 0.000000 4.0 4.0
信頼性 4.0 4.0 4.000000 4.000000 0.0 4.0
素直さ 4.0 4.0 3.999999 4.000000 4.0 0.0
RMSEAs:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 1.0 0.0 0.0 0.0 0.0 0.0
社交性 0.0 1.0 0.0 0.0 0.0 0.0
積極性 0.0 0.0 1.0 0.0 0.0 0.0
知性 0.0 0.0 0.0 1.0 0.0 0.0
信頼性 0.0 0.0 0.0 0.0 1.0 0.0
素直さ 0.0 0.0 0.0 0.0 0.0 1.0
Adjecency Matrix:
外向性 社交性 積極性 知性 信頼性 素直さ
外向性 0.000000 0.486598 0.597399 0.000000 0.000000 0.000000
社交性 0.486477 0.000000 0.557879 0.000000 0.000000 0.000000
積極性 0.597371 0.557998 0.000000 0.000000 0.000000 0.000000
知性 0.000000 0.000000 0.000000 0.000000 0.436245 0.627193
信頼性 0.000000 0.000000 0.000000 0.436322 0.000000 0.583388
素直さ 0.000000 0.000000 0.000000 0.627123 0.583625 0.000000
Adjecency Matrixの結果を見ると、因果関係は存在せず、全て相関関係のようですね。
コード解説
探索のメインとなるコードは、左右の各要素についての単純な二重ループです。
def __make_sem_combination_matrix(self,data):
#パラメータの準備
n = len(data.columns)
alpha = self.alpha_
eps_cfi = self.eps_cfi_
eps_gfi = self.eps_gfi_
eps_nfi = self.eps_nfi_
#変数名リスト
varnames = data.columns.to_list()
# 行列の初期化
stdest_matrix = np.zeros((n,n))
pvalue_matrix = np.ones((n,n))
cfi_matrix = np.zeros((n,n))
gfi_matrix = np.zeros((n,n))
nfi_matrix = np.zeros((n,n))
aic_matrix = np.zeros((n,n))
rmsea_matrix = np.ones((n,n))
adj_matrix = np.zeros((n,n))
#判定ループ
for i in range(n):
lhs = varnames[i];
for j in range(i+1,n,1):
rhs = varnames[j];
#print(f"{lhs} ~ {rhs}")
# SEMモデルの作成
model_l = semopy.Model(f"{lhs} ~ {rhs}\n") # lhs <- rhs用SEMモデル
model_r = semopy.Model(f"{rhs} ~ {lhs}\n") # lhs -> rhs用SEMモデル
# SEMの計算
model_l.fit(data) # lhs <- rhsのモデリング
model_r.fit(data) # lhs -> rhsのモデリング
# SEMの計算結果の取り出し
try:
inspect_l = model_l.inspect(std_est=True)
stdest_l = inspect_l["Est. Std"][0]
pvalue_l = inspect_l["p-value"][0]
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : {lhs} ~ {rhs}")
stdest_l = 0.0
pvalue_l = 1.0
try:
inspect_r = model_r.inspect(std_est=True)
stdest_r = inspect_r["Est. Std"][0]
pvalue_r = inspect_r["p-value"][0]
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : {rhs} ~ {lhs}")
stdest_r = 0.0
pvalue_r = 1.0
# 適合度指標の取り出し
goodfit_l = semopy.calc_stats(model_l)
goodfit_r = semopy.calc_stats(model_r)
#
cfi_l = goodfit_l["CFI"][0]
cfi_r = goodfit_r["CFI"][0]
gfi_l = goodfit_l["GFI"][0]
gfi_r = goodfit_r["GFI"][0]
nfi_l = goodfit_l["NFI"][0]
nfi_r = goodfit_r["NFI"][0]
aic_l = goodfit_l["AIC"][0]
aic_r = goodfit_r["AIC"][0]
rmsea_l = goodfit_l["RMSEA"][0]
rmsea_r = goodfit_r["RMSEA"][0]
# 因果関係の判定、因果行列の作成
if pvalue_l < alpha and pvalue_r < alpha : # 両向き有意な場合
if abs(cfi_l-cfi_r) > eps_cfi or abs(gfi_l-gfi_r) > eps_gfi or abs(nfi_l-nfi_r) > eps_nfi : # 左右差がある場合
if cfi_l > cfi_r : # 左向き確定
adj_matrix[i,j] = stdest_l
else : # 右向き
adj_matrix[j,i] = stdest_r
else: # 左右差がない場合 --> 両向き
adj_matrix[i,j] = stdest_l
adj_matrix[j,i] = stdest_r
elif pvalue_l < alpha and pvalue_r >=alpha : # 左向きのみ有意な場合
adj_matrix[i,j] = stdest_l
elif pvalue_l >=alpha and pvalue_r < alpha : # 右向きのみ有意な場合
adj_matrix[j,i] = stdest_r
#else : # 左右どちらも有意ではない場合
#
# # do nothing
#
#print(f"i={i},j={j},n={n}")
# 使用した指標を行列に格納する
stdest_matrix[i,j] = stdest_l
stdest_matrix[j,i] = stdest_r
pvalue_matrix[i,j] = pvalue_l
pvalue_matrix[j,i] = pvalue_r
cfi_matrix[i,j] = cfi_l
cfi_matrix[j,i] = cfi_r
gfi_matrix[i,j] = gfi_l
gfi_matrix[j,i] = gfi_r
nfi_matrix[i,j] = nfi_l
nfi_matrix[j,i] = nfi_r
aic_matrix[i,j] = aic_l
aic_matrix[j,i] = aic_r
rmsea_matrix[i,j] = rmsea_l
rmsea_matrix[j,i] = rmsea_r
プログラム
GitHub
下記から、causal_analysis.pyをダウンロードしてください。SEMCausalFinderクラスは、causal_analysisパッケージの中に定義してあります。
https://github.com/jyam45/stat_analytics/tree/main/causal_analysis
コピペ用ソースコード
causal_analysis.py
import numpy as np
import pandas as pd
import lingam
import semopy
import graphviz
#import causal_discover as cd
####################################################################
###
### Scaing Functions for pandas.DataFrame
###
### 1. Use standardize() for normal-distribution data such as Likert N-item questionnaire.
### 2. Use normalize() for the data bounded between upper- and lower-limit.
### 3. Use scale_by() for dynamics data such as time-plot data.
###
####################################################################
def standardize(data,ddof=False):
stddat = ( data - data.mean() ) / data.std(ddof=ddof)
stddat.index = data.index
stddat.columns = data.columns
return stddat
def normalize(data):
stddat = ( data - data.min() ) / ( data.max() - data.min() )
stddat.index = data.index
stddat.columns = data.columns
return stddat
def scale_by(data,row=0):
stddat = data / data.iloc[row]
stddat.index = data.index
stddat.columns = data.columns
return stddat
####################################################################
###
### Causal Finder using SEM
###
### by comparing Fitness on Single Path Regression
###
####################################################################
class SEMCausalFinder:
def __init__(self,alpha=0.05,eps_cfi=0.05,eps_gfi=0.05,eps_nfi=0.05):
self.adjacency_matrix_ = None
self.p_values_matrix_ = None
self.stdests_ = None
self.pvalues_ = None
self.cfis_ = None
self.gfis_ = None
self.nfis_ = None
self.aics_ = None
self.rmseas_ = None
self.adjmax_ = None
self.alpha_ = alpha # significant level for p-value (default:0.05)
self.eps_cfi_ = eps_cfi # significant difference level for CFI (default:0.05)
self.eps_gfi_ = eps_gfi # significant difference level for GFI (default:0.05)
self.eps_nfi_ = eps_nfi # significant difference level for NFI (default:0.05)
######################
### STATIC METHODS ###
######################
@staticmethod
def print_full(df,title="No title"):
x = ""
if df is not None:
x = f"{title}:\n"
pd.set_option('display.max_columns',len(df.columns))
pd.set_option('display.max_rows',len(df.index))
x += str(df)
x += "\n\n"
return x
########################
### INTERNAL METHODS ###
########################
def __make_sem_combination_matrix(self,data):
n = len(data.columns)
alpha = self.alpha_
eps_cfi = self.eps_cfi_
eps_gfi = self.eps_gfi_
eps_nfi = self.eps_nfi_
varnames = data.columns.to_list()
stdest_matrix = np.zeros((n,n))
pvalue_matrix = np.ones((n,n))
cfi_matrix = np.zeros((n,n))
gfi_matrix = np.zeros((n,n))
nfi_matrix = np.zeros((n,n))
aic_matrix = np.zeros((n,n))
rmsea_matrix = np.ones((n,n))
adj_matrix = np.zeros((n,n))
for i in range(n):
lhs = varnames[i];
for j in range(i+1,n,1):
rhs = varnames[j];
#print(f"{lhs} ~ {rhs}")
model_l = semopy.Model(f"{lhs} ~ {rhs}\n") # lhs <- rhs
model_r = semopy.Model(f"{rhs} ~ {lhs}\n") # lhs -> rhs
model_l.fit(data)
model_r.fit(data)
try:
inspect_l = model_l.inspect(std_est=True)
stdest_l = inspect_l["Est. Std"][0]
pvalue_l = inspect_l["p-value"][0]
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : {lhs} ~ {rhs}")
stdest_l = 0.0
pvalue_l = 1.0
try:
inspect_r = model_r.inspect(std_est=True)
stdest_r = inspect_r["Est. Std"][0]
pvalue_r = inspect_r["p-value"][0]
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : {rhs} ~ {lhs}")
stdest_r = 0.0
pvalue_r = 1.0
goodfit_l = semopy.calc_stats(model_l)
goodfit_r = semopy.calc_stats(model_r)
#
cfi_l = goodfit_l["CFI"][0]
cfi_r = goodfit_r["CFI"][0]
gfi_l = goodfit_l["GFI"][0]
gfi_r = goodfit_r["GFI"][0]
nfi_l = goodfit_l["NFI"][0]
nfi_r = goodfit_r["NFI"][0]
aic_l = goodfit_l["AIC"][0]
aic_r = goodfit_r["AIC"][0]
rmsea_l = goodfit_l["RMSEA"][0]
rmsea_r = goodfit_r["RMSEA"][0]
if pvalue_l < alpha and pvalue_r < alpha : # 両向き有意な場合
if abs(cfi_l-cfi_r) > eps_cfi or abs(gfi_l-gfi_r) > eps_gfi or abs(nfi_l-nfi_r) > eps_nfi : # 左右差がある場合
if cfi_l > cfi_r : # 左向き確定
adj_matrix[i,j] = stdest_l
else : # 右向き
adj_matrix[j,i] = stdest_r
else: # 左右差がない場合 --> 両向き
adj_matrix[i,j] = stdest_l
adj_matrix[j,i] = stdest_r
elif pvalue_l < alpha and pvalue_r >=alpha : # 左向きのみ有意な場合
adj_matrix[i,j] = stdest_l
elif pvalue_l >=alpha and pvalue_r < alpha : # 右向きのみ有意な場合
adj_matrix[j,i] = stdest_r
#else : # 左右どちらも有意ではない場合
#
# # do nothing
#
#print(f"i={i},j={j},n={n}")
stdest_matrix[i,j] = stdest_l
stdest_matrix[j,i] = stdest_r
pvalue_matrix[i,j] = pvalue_l
pvalue_matrix[j,i] = pvalue_r
cfi_matrix[i,j] = cfi_l
cfi_matrix[j,i] = cfi_r
gfi_matrix[i,j] = gfi_l
gfi_matrix[j,i] = gfi_r
nfi_matrix[i,j] = nfi_l
nfi_matrix[j,i] = nfi_r
aic_matrix[i,j] = aic_l
aic_matrix[j,i] = aic_r
rmsea_matrix[i,j] = rmsea_l
rmsea_matrix[j,i] = rmsea_r
self.stdests_ = pd.DataFrame(stdest_matrix,columns=varnames,index=varnames)
self.pvalues_ = pd.DataFrame(pvalue_matrix,columns=varnames,index=varnames)
self.cfis_ = pd.DataFrame(cfi_matrix ,columns=varnames,index=varnames)
self.gfis_ = pd.DataFrame(gfi_matrix ,columns=varnames,index=varnames)
self.nfis_ = pd.DataFrame(nfi_matrix ,columns=varnames,index=varnames)
self.aics_ = pd.DataFrame(aic_matrix ,columns=varnames,index=varnames)
self.rmseas_ = pd.DataFrame(rmsea_matrix ,columns=varnames,index=varnames)
self.adjmat_ = pd.DataFrame(adj_matrix ,columns=varnames,index=varnames)
self.p_values_matrix_ = pvalue_matrix
self.adjacency_matrix_ = adj_matrix
#####################
### CLASS METHODS ###
#####################
def fit(self,data):
self.__make_sem_combination_matrix(data)
print(self)
return self
def get_error_independence_p_values(self,data):
return self.p_values_matrix_
def __str__(self):
x = "\n"
x += "Cutoff for p-value : " + str(self.alpha_) + "\n"
x += "Cutoff for CFI diff. : " + str(self.eps_cfi_) + "\n"
x += "Cutoff for GFI diff. : " + str(self.eps_gfi_) + "\n"
x += "Cutoff for NFI diff. : " + str(self.eps_nfi_) + "\n"
x += "\n"
x += self.print_full(self.stdests_,"Est. Stds")
x += self.print_full(self.pvalues_,"p-values" )
x += self.print_full(self.cfis_ ,"CFIs")
x += self.print_full(self.gfis_ ,"GFIs")
x += self.print_full(self.nfis_ ,"NFIs")
x += self.print_full(self.aics_ ,"AICs")
x += self.print_full(self.rmseas_ ,"RMSEAs")
x += self.print_full(self.adjmat_ ,"Adjercency Matrix")
return x
####################################################################
###
### Causal Analizer using LiNGAM and Semopy
###
### by connecting Causal Discovery (LiNGAM) & Causal Guess (Semopy)
###
####################################################################
class CausalAnalysis:
def __init__(self,style="lavaan",method="direct",tol=0.01,alpha=0.05,enable_cov=False):
self.model_ = None # SEM-model text
self.sem_ = None # SEM-model object
self.lingam_ = None # LiNGAM-model object
self.tol_ = tol # Threshold for Causality effectiveness (defualt=0.01)
self.alpha_ = alpha # Significant value level for p-value (defualt=0.05)
self.style_ = style # SEM-model format style ["lavaan"|"sem"] (default=lavaan)
self.method_ = method# Discovery medhod ["direct"|"ica"|"var"|"varma"|"rcd"|"sem"] (default=direct)
self.enable_cov_= enable_cov # Enable Covariance Path <--> (default=False)
######################
### STATIC METHODS ###
######################
@staticmethod
def print_full(df,title="No title"):
x = ""
if df is not None:
x = f"{title}:\n"
pd.set_option('display.max_columns',len(df.columns))
pd.set_option('display.max_rows',len(df.index))
x += str(df)
x += "\n\n"
return x
@staticmethod
def read_model(file,encoding='UTF-8'):
file = open(file,'r',encoding=encoding)
text = file.read()
file.close()
return text
@staticmethod
def write_model(file,data,encoding='UTF-8'):
file = open(file,'w',encoding=encoding)
file.write(data)
file.close()
@staticmethod
def lingam_to_sem(data,lingam_model,tol=0.01,alpha=0.05,sem_style="lavaan",enable_cov=False):
# 変数の数
n = len(data.columns)
# 変数名リスト
varnames = data.columns.to_list()
# 因果係数の行列
adjmat = lingam_model.adjacency_matrix_
# 変数の順序
if hasattr(lingam_model,"causal_order_"):
order = lingam_model.causal_order_ # 因果順序
else:
order = list(range(n)) # 0, 1, 2, ..., n
# p値の行列
if hasattr(lingam_model,"get_error_independence_p_values"):
pvalmat = lingam_model.get_error_independence_p_values(data)
else:
pvalmat = np.zeros((n,n)) # ゼロ行列
# モデルの書き方
if sem_style == "lavaan":
left = "~"
right= "=~"
both = "~~"
elif sem_style == "sem":
left = "<-"
right= "->"
both = "<>"
else:
print(f"lingam_to_sem: irregal option sem_style={sem_style}. ['lavaan','sem']")
sem_model=""
sem_model2=""
for i in range(n):
for j in range(i+1,n,1):
lhs = varnames[order[i]]
rhs = varnames[order[j]]
if adjmat[i,j] > tol and pvalmat[i,j] < alpha :
if adjmat[j,i] > tol and pvalmat[j,i] < alpha :
sem_model2 += f"{lhs} {both} {rhs}\n"
else:
sem_model += f"{lhs} {left} {rhs}\n"
else:
if adjmat[j,i] > tol and pvalmat[j,i] < alpha :
sem_model += f"{rhs} {left} {lhs}\n"
#print(sem_model)
if enable_cov : sem_model += sem_model2
return sem_model
########################
### INTERNAL METHODS ###
########################
def __lingam_to_sem(self,data,lingam_model):
return self.lingam_to_sem(data,lingam_model,tol=self.tol_,alpha=self.alpha_,sem_style=self.style_,enable_cov=self.enable_cov_)
#####################
### CLASS METHODS ###
#####################
# 因果分析を行う
def analyze(self,data):
model = self.discover(data)
#print(model)
#print(data.columns.to_list())
result = self.evaluate(data,model)
return result
# 因果探索を行う
def discover(self,data):
if self.method_ == "direct":
self.lingam_ = lingam.DirectLiNGAM()
self.lingam_.fit(data)
elif self.method_ == "ica":
self.lingam_ = lingam.ICALiNGAM()
self.lingam_.fit(data)
elif self.method_ == "var":
self.lingam_ = lingam.VARLiNGAM()
self.lingam_.fit(data)
elif self.method_ == "varma":
self.lingam_ = lingam.VARMALiNGAM()
self.lingam_.fit(data)
elif self.method_ == "rcd":
self.lingam_ = lingam.RCD()
self.lingam_.fit(data)
elif self.method_ == "sem":
self.lingam_ = SEMCausalFinder()
self.lingam_.fit(data)
#elif self.method_ == "lina":
# self.lingam_ = lingam.LiNA()
# self.lingam_.fit(data)
#elif self.method_ == "lim":
# self.lingam_ = lingam.LiM()
# dis_con = [0] * n
# self.lingam_.fit(data,dis_con)
#elif self.method_ == "resit":
# self.lingam_ = lingam.RESIT()
# self.lingam_.fit(data)
else :
print(f"discover : Invalid option self.method_ = {self.method_}")
self.model_ = self.__lingam_to_sem(data,self.lingam_)
return self.model_
# 因果推論を行う
def evaluate(self,data,model=None):
if model is None: model = self.model_
if model is None: raise ValueError("No model")
if str(model) == "":
if self.lingam_ is not None:
print("Adjacency Matrix:")
print(self.lingam_.adjacency_matrix_)
raise ValueError("Empty model")
self.sem_ = semopy.Model(model)
res = self.sem_.fit(data)
return res
# モデルや図をファイルに書き込む
def save(self,init_imgfile="found_causality.png",final_imgfile="evaled_causality.png",
init_dotfile="found_causality.gv", final_dotfile="evaled_causality.gv",
modfile="sem_model.txt",
basedir=".",show=False):
if self.lingam_ is not None:
dotdata = lingam.utils.make_dot(self.lingam_.adjacency_matrix_)
graph = graphviz.Source(dotdata,basedir+"/"+init_dotfile,format='png')
#graph.render(outfile=basedir+"/"+init_imgfile,view=show)
graph.render(view=show)
if self.model_ is not None:
self.write_model(basedir+"/"+modfile,self.model_)
if self.sem_ is not None:
self.write_model(basedir+"/"+final_dotfile,str(dotdata))
try:
dotdata = semopy.semplot(self.sem_,basedir+"/"+final_imgfile,plot_covs=True,std_ests=True,show=show)
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : Graph can not be drawn.")
# インスタンスをダンプする
def __str__(self):
x = "\n"
x += "SEM Model Style : " + str(self.style_) + "\n"
x += "LiNGAM Method : " + str(self.method_) + "\n"
x += "Cutoff Causality : " + str(self.tol_) + "\n"
x += "Cutoff p-value : " + str(self.alpha_) + "\n"
x += "\n"
if self.lingam_ is not None:
if hasattr(self.lingam_,"causal_order_"):
x += "Causal Order:\n"
x += str(self.lingam_.causal_order_)
x += "\n\n"
x += "Causal Matrix:\n"
x += str(self.lingam_.adjacency_matrix_)
x += "\n\n"
if self.model_ is not None:
x += "SEM Model:\n"
x += str(self.model_)
x += "\n\n"
if self.sem_ is not None:
try:
ins = self.sem_.inspect(std_est=True)
pd.set_option('display.max_columns',len(ins.columns))
pd.set_option('display.max_rows',len(ins.index))
x += "Infomation:\n"
x += str(ins)
x += "\n\n"
except np.linalg.LinAlgError as e:
print(f"NumPy Linear Algebra Error : Infomationcan not be shown.")
x += "Fitness Index:\n"
gfi = semopy.calc_stats(self.sem_)
pd.set_option('display.max_columns',len(gfi.index))
pd.set_option('display.max_rows',len(gfi.columns))
x += str(gfi.T)
x += "\n\n"
return x
