Python による主成分分析と因子分析

sklearn.decomposition.PCA を使い、主成分分析を行う。また、sklearn.decomposition.FactorAnalysis を使い、因子分析を行う。

Python を起動し、必要なライブラリを import する。

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, FactorAnalysis

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"手術時間":[60,50,100,90,30,40,80,90], "出血量":[15,100,90,20,18,20,120,90], "結紮失敗回数":[0,1,2,3,2,1,5,1]})

data の内容を確認する。

data

data を標準化する。

sc = StandardScaler()
sc.fit(data)
data_s = sc.transform(data)

sklearn.decomposition.PCA を使い、主成分分析を行う。

pca = PCA()
pca.fit(data_s)
data_t = pca.transform(data_s)
 
data_t  #主成分得点

pca.explained_variance_ratio_  #寄与率
np.cumsum(pca.explained_variance_ratio_)  #累積寄与率

pca.explained_variance_ #固有値
pca.components_ #固有ベクトル

第1主成分(PC1)と第2主成分(PC2)で散布図を描きたい場合は matplotlib を使う。

from matplotlib import pyplot as plt
x = data_t[:, 0] #第1主成分を抽出
y = data_t[:, 1] #第2主成分を抽出
plt.figure()
plt.scatter(x, y)
plt.grid()
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()

sklearn.decomposition.FactorAnalysis を使い、因子分析を行う。ここでは因子数は 1、varimax 回転を選択している。

FA = FactorAnalysis(n_components = 1, rotation = "varimax") 
FA.fit_transform(data_s) #因子得点
FA.components_ #因子負荷量

Python による Kaplan-Meier 法、ログランク検定、Cox 比例ハザードモデル

lifelines を使い、ログランク検定および Cox 比例ハザードモデル分析を行う。

必要なライブラリを import する。

!pip install lifelines #必要な場合のみ
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test
import pandas as pd

通常は Excel 等で作成した CSV ファイルを R に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。打ち切り = 0 なら打ち切りありとしている。

data = pd.DataFrame({"患者番号":[1,2,3,4,5,6,7,8,9,10], "抗がん剤":["A","A","A","A","A","B","B","B","B","B"], "生存期間":[20,30,40,50,60,55,75,105,115,120], "打ち切り":[1,1,1,1,0,1,1,1,1,0], "全身状態":[2,3,2,3,2,0,1,0,0,1]})

data の内容を確認する。

data

data から A 群、B 群をそれぞれ抽出し、Kaplan-Meier モデルに適合させる。

#A 群
A = data.query('抗がん剤 == "A"')
kmf_A = KaplanMeierFitter()
kmf_A.fit(A.生存期間, A.打ち切り)

#B 群
B = data.query('抗がん剤 == "B"')
kmf_B = KaplanMeierFitter()
kmf_B.fit(B.生存期間, B.打ち切り)

Kaplan-Meier 法により、抗がん剤 A 群と抗がん剤 B 群を比較する。

#生存期間中央値
kmf_A.median_survival_time_
kmf_B.median_survival_time_

#生存確率、信頼区間
kmf_A.confidence_interval_survival_function_
kmf_B.confidence_interval_survival_function_

#5年 (60カ月) 生存率
kmf_A.predict(60)
kmf_B.predict(60)

#Kaplan-Meier 曲線の図示
kmf_A.plot() and kmf_B.plot()

ログランク検定を行う。

result = logrank_test(A.生存期間, B.生存期間, A.打ち切り, B.打ち切り)
result.print_summary()

p 値をもっと正確に知りたい場合は以下のようにする。

result.p_value

次に、Cox 比例ハザードモデル分析を行う。lifelines で Cox 比例ハザードモデル分析を行う場合、data に “A”、”B” などの文字が含まれているとうまく解析できない。また、変量名に “生存期間” などの日本語(2バイト文字)を使うとグラフのラベルに表示されない。

以上を踏まえ、あらためて以下の簡単なデータフレームを使う。

data_cox = pd.DataFrame({"Group":[0,0,0,0,0,1,1,1,1,1], "Duration":[20,30,40,50,60,55,75,105,115,120], "Censor":[1,1,1,1,0,1,1,1,1,0], "PS":[2,3,2,3,2,0,1,0,0,1]})

data_cox の内容を確認する。

data_cox

Cox 比例ハザードモデル分析を行う。

cph = CoxPHFitter()
cph.fit(data_cox, "Duration", "Censor")
cph.print_summary()

各変量の信頼区間をグラフに描出することもできる。

cph.plot()

Python によるガンマ回帰分析

statsmodels.formula.api.glm() コマンドでガンマ回帰分析を行う。一般化線形モデルとして行う。

Python を起動し、必要なライブラリを import する。

import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"出血量":[10,200,300,5,5,3,250,500], "手術方法":["内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"], "手術時間":[60,50,60,60,30,40,60,50]})

data の内容を確認する。

data

statsmodels.formula.api.glm() コマンドでガンマ回帰分析を行う。リンク関数は逆関数が自動的に選択される。

result = smf.glm(formula = "出血量 ~ 手術時間 + 手術方法", data = data, family = sm.families.Gamma()).fit()
result.summary() #結果の表示
result.aic #赤池情報量基準

リンク関数を対数関数にしたい場合は以下のようにする。

result = smf.glm(formula = "出血量 ~ 手術時間 + 手術方法", data = data, family = sm.families.Gamma(sm.families.links.log)).fit()
result.summary() #結果の表示
result.aic #赤池情報量基準

Python による Tukey-Kramer test

statsmodels.stats.multicomp.pairwise_tukeyhsd() コマンドで Tukey-Kramer test を行う。p 値も示す。

Python を起動し、必要なライブラリを import する。

import pandas as pd
from statsmodels.stats.multicomp import pairwise_tukeyhsd

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"デバイス":["デバイス A","デバイス B","デバイス A","デバイス A","デバイス C","デバイス C","デバイス B","デバイス C","デバイス A","デバイス A",], "出血量":[30,20,5,5,100,80,10,60,10,25]})

data の内容を確認する。

data

デバイス間で出血量の平均に差があるかどうかを statsmodels.stats.multicomp.pairwise_tukeyhsd() コマンドで検定する。結果を表示するのに print() コマンドだけでは p値が示されないが、summary() コマンドを使用することで p 値も示される。

pairwise_tukeyhsd(data["出血量"], data["デバイス"]).summary() #95%信頼区間
pairwise_tukeyhsd(data["出血量"], data["デバイス"], alpha = 0.01).summary() #99%信頼区間

計算結果は、R での結果とは微妙に異なる

Python によるポワソン回帰分析

statsmodels.formula.api.glm() コマンドでポワソン回帰分析を行う。一般化線形モデルとして行う。

Python を起動し、必要なライブラリを import する。

import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"年":["X","X+1","X+2","X+3","X+4","X+5","X+6","X+7","X+8","X+9"], "重大インシデント数":[5,3,4,4,5,0,1,0,0,0], "残業80時間以上の医師数":[50,40,41,30,31,7,8,7,6,5], "タイムアウトの導入":["no","no","yes","yes","yes","yes","yes","yes","yes","yes"]})

data の内容を確認する。

data

statsmodels.formula.api.glm() コマンドでポワソン回帰分析を行う。リンク関数は対数関数が自動的に選択される。

result = smf.glm(formula = "重大インシデント数 ~ 残業80時間以上の医師数 + タイムアウトの導入", data = data, family = sm.families.Poisson()).fit()
result.summary() #結果の表示
result.aic #赤池情報量基準

Python によるロジスティック回帰分析

statsmodels.formula.api.glm() コマンドでロジスティック回帰分析を行う。一般化線形モデルとして行う。

Python を起動し、必要なライブラリを import する。

import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"再発":["yes","yes","yes","no","no","no","yes","no","yes","yes"], "化学療法回数":[6,5,4,3,3,3,1,0,0,0], "手術時腫瘍残存":["yes","yes","yes","no","no","no","no","no","yes","yes"]})

data の内容を確認する。

data

statsmodels.formula.api.glm() コマンドでロジスティック回帰分析を行う。リンク関数はロジット関数が自動的に選択される。

result = smf.glm(formula = "再発 ~ 化学療法回数 + 手術時腫瘍残存", data = data, family = sm.families.Binomial()).fit()
result.summary() #結果の表示
result.aic #赤池情報量基準

オッズ比は NumPy をインポートし、 numpy.exp() コマンドの () 内に各パラメータとその95%信頼区間の数値をひとつひとつ入れて計算する。

Python による重回帰分析

statsmodels.formula.api.ols() コマンドで重回帰分析を行う。

Python を起動し、必要なライブラリを import する。

import pandas as pd
import statsmodels.formula.api as smf

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"出血量":[10,200,300,5,5,3,250,500], "手術方法":["内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"], "手術時間":[60,50,60,60,30,40,60,50]})

data の内容を確認する。

data

statsmodels.formula.api.ols() コマンドで重回帰分析を行う。単回帰分析を行うこともできる。手術方法のような2値データはダミーに変換される。

lm = smf.ols("出血量 ~ 手術時間 + 手術方法", data).fit()
lm.summary() #結果の表示。

statsmodels.formula.api.glm() コマンドを使い、一般化線形モデルとして重回帰分析を行うこともできる。この場合、検定は t 検定ではなく、Wald 検定が行われる。リンク関数は恒等関数が自動的に選択される。また、summary() コマンドで AIC(赤池情報量基準)は表示されない。

import statsmodels.api as sm #確率分布指定のために必要
glm = smf.glm(formula = "出血量 ~ 手術時間 + 手術方法", data = data, family = sm.families.Gaussian()).fit()
glm.summary() #結果の表示
glm.aic #赤池情報量基準

Python による Kruskal-Wallis test と ANOVA

scipy.stats.kruskal() コマンドで Kruskal-Wallis test を行う。また、scipy.stats.f_oneway() コマンドで one-way ANOVA (一元配置分散分析)を行う。

Python を起動し、pandas と scipy.stats を import する。

import pandas as pd
from scipy import stats

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"デバイス":["デバイス A","デバイス B","デバイス A","デバイス A","デバイス C","デバイス C","デバイス B","デバイス C","デバイス A","デバイス A",], "出血量":[30,20,5,5,100,80,10,60,10,25]})

data の内容を確認する。

data

query() コマンドを使い、デバイスによって data を3群にグループ分けする。

A = data.query('デバイス == "デバイス A"')
B = data.query('デバイス == "デバイス B"')
C = data.query('デバイス == "デバイス C"')

それぞれのグループの平均や中央値等の基本統計量は describe() コマンドによって得られる。

A.describe()
B.describe()
C.describe()

A, B, C の出血量を Series 型に変換する。

x = A["出血量"]
y = B["出血量"]
z = C["出血量"]

x, y, z に偏りがあるかどうかを scipy.stats.kruskal() コマンド、scipy.stats.f_oneway() コマンドで検定する。

stats.kruskal(x, y, z) #Kruskal-Wallis test
stats.f_oneway(x, y, z) #ANOVA

Python による Fisher exact test とカイ2乗検定

scipy.stats.fisher_exact() コマンドで Fisher exact test を行う。また、scipy.stats.chi2_contingency() コマンドでカイ2乗検定を行う。

Python を起動し、pandas と scipy.stats を import する。

import pandas as pd
from scipy import stats

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"術者":["太郎","太郎","花子","花子"], "手術方法":["内視鏡","開腹","内視鏡","開腹"], "回数":[60,5,20,40]})

data の内容を確認する。

data

pandas.pivot_table() コマンドで、data を2×2クロス集計表にし、table という変数に格納する(別記事)。

table = pd.pivot_table(data, values = "回数", index = "術者", columns = "手術方法", aggfunc = "sum")

table の内容を確認する。

table

scipy.stats.fisher_exact() コマンド、scipy.stats.chi2_contingency() コマンドで術者によって手術方法(内視鏡か開腹か)に偏りがあるかどうかを検定する。

stats.fisher_exact(table) #Fisher exact test
stats.chi2_contingency(table) #カイ2乗検定(イエーツ補正あり)

Python による独立2群の比較

scipy.stats.ttest_ind() コマンドを使い、独立2群を比較する。

Python を起動し、pandas と scipy.stats を import する。

import pandas as pd
from scipy import stats

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"手術方法":["内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"], "手術時間":[60,50,100,90,30,40,80,90], "出血量":[15,100,90,20,18,20,120,90]})

data の内容を確認する。

data

query() コマンドを使い、手術方法(内視鏡か開腹か)によって data を2群にグループ分けする。

A = data.query('手術方法 == "内視鏡"')
B = data.query('手術方法 == "開腹"')

それぞれのグループの平均等の基本統計量は describe() コマンドによって得られる。

A.describe()
B.describe()

A 群、B 群の手術時間と出血量を Student t test で比較する。

stats.ttest_ind(A["手術時間"], B["手術時間"]) #手術時間の比較 
stats.ttest_ind(A["出血量"], B["出血量"]) #出血量の比較

古典的な Student t test 以外は、以下のようにする(ここでは手術時間を比較)。

stats.ttest_ind(A["手術時間"], B["手術時間"], equal_var = False) # Welch t test
stats.mannwhitneyu(A["手術時間"], B["手術時間"]) #マン・ホイットニー u test

Python による信頼区間の計算

scipy.stats.t.interval() コマンドを使い、平均の信頼区間を計算する。

Python を起動し、pandas と scipy.stats を import する。

import pandas as pd
from scipy import stats

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"手術時間":[60,50,100,90,30,40], "出血量":[15,25,35,55,45,25]})

data の内容を確認する。

data

ここでは手術時間の平均の信頼区間を計算してみる。data から手術時間を抽出し、time という変数に格納しておく。

time = data["手術時間"]

平均の信頼区間の計算のためには、(1) 信頼水準 (alpha)、(2) 自由度 (df)、(3) 標本平均 (loc)、(4)標準誤差 (scale) の情報が必要なため、あらかじめこれらを計算する。

a = 0.95 #信頼水準。ここでは 95%。
d = len(time)-1 #自由度 (= サンプルサイズ-1)
m = time.mean() #標本平均
s = stats.sem(time) #標準誤差。不偏分散から計算。

分布を t 分布と仮定して stats.t.interval() コマンドで信頼区間を計算する。

stats.t.interval(alpha = a, df = d, loc = m, scale =s)

scipy.stats.norm.interval() コマンドを使えば正規分布を仮定することもできる。

stats.norm.interval(alpha = a, loc = m, scale = time.std()) #time.std() は標準偏差。

Python による相関係数の計算

pandas.DataFrame.corr() コマンドをつかい、相関係数を計算する。P 値も計算したい場合には scipy.stats.pearsonr() を使う。

Python を起動し、pandas を import する。

import pandas as pd

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"縫合結紮の練習時間 (時間)":[100,150,200,250,300,1000], "手術時間 (分)":[210,220,150,100,90,45], "出血量 (ml)":[50,70,60,30,20,10]})

data の内容を確認する。

data

corr()コマンドでピアソン積率相関係数(いわゆる相関係数)を表示する。

data.corr() #ピアソン積率相関係数。

以下のようにすると順位相関係数も計算できる。

data.corr(method = "spearman") #スピアマン順位相関係数。
data.corr(method = "kendall") #ケンドール順位相関係数。

Pandas だけでは P 値は自動計算されない。これも計算したい場合は以下のように scipy.stats を import して使う。ここでは、x = 縫合結紮の練習時間 (時間)、y = 手術時間 (分)を data から抽出し、x と y の相関係数と両側検定による P 値を計算している。

from scipy import stats
x, y = data["縫合結紮の練習時間 (時間)"], data["手術時間 (分)"]
stats.pearsonr(x, y) #ピアソン積率相関係数。
stats.spearmanr(x, y) #スピアマン順位相関係数。
stats.kendalltau(x, y) #ケンドール順位相関係数。

Python によるクロス集計表

pandas.pivot_table() コマンドを使い、クロス集計表を作成する。

Python を起動し、pandas を import する。

import pandas as pd

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"術者":["太郎","花子","太郎","花子","太郎","花子"], "手術方法":["内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡"], "手術時間":[60,50,100,90,30,40], "出血量":[15,25,35,55,45,25]})

data の内容を確認する。

data

術者と手術方法で、手術時間の平均がどう異なるかを集計する。

pd.pivot_table(data, values = "手術時間", index = "術者", columns = "手術方法", aggfunc = "mean")

手術時間だけでなく、出血量の平均も集計する。

pd.pivot_table(data, values = ["手術時間", "出血量"], index = "術者", columns = "手術方法", aggfunc = "mean")

aggfunc は以下のように変えられる。
mean: 平均
count: データの個数
sum: 合計
var: 分散
std: 標準偏差
max: 最大値
min: 最小値
median: 中央値

margins = True を付け加えると、小計を表示する。また、aggfunc = “count” とすることで件数を表示する通常のクロス集計になる。

pd.pivot_table(data,values = "手術時間", index = "術者", columns = "手術方法", aggfunc = "count", margins = True)

Python による基本統計量の計算

pandas.DataFrame.describe() コマンドを使い、平均や標準偏差などの基本統計量を計算する。

Python を起動し、pandas を import する。

import pandas as pd

通常は Excel 等で作成した CSV ファイルを pandas に読み込んで分析する。しかし、ここでは便宜的に以下の簡単なデータフレームを作成し、data という変数に格納しておく。

data = pd.DataFrame({"術者":["太郎","花子","太郎","花子","太郎","花子"], "手術方法":["内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡"], "手術時間":[60,50,100,90,30,40], "出血量":[15,25,35,55,45,25]})

data の内容を確認する。

data

pandas.DataFrame.describe() コマンドで平均や標準偏差などの基本統計量を表示する。

data.describe()

グループ毎に統計量を出す。

group = data.groupby("術者") #術者でグループ分けする。 
group.describe() #術者毎に手術時間や出血量の平均等を出す。

describe() コマンド以外の方法による基本統計量の計算。

len(data) #サンプルサイズ
data.sum() #合計
data.mean() #平均
data.var() #分散。()内を ddof=0 なら標本分散。省略または ddof=1 なら不偏分散。
data.std() #標準偏差。()内を ddof=0 なら標本分散の平方根。省略または ddof=1 なら不偏分散の平方根。
data.max() #最大値
data.min() #最小値
data.median() #中央値