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()

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

survdiff コマンドを使い、ログランク検定を行う。また、coxph コマンドを使い、Cox 比例ハザードモデル分析を行う。

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

library(survival)

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

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

data の内容を確認する。

data

Kaplan-Meier 法により、抗がん剤 A 群と抗がん剤 B 群を比較する。Surv コマンド、survfit コマンドを使う。

kaplan = survfit(Surv(生存期間, 打ち切り) ~ 抗がん剤, data)
kaplan #生存期間中央値と信頼区間
summary(kaplan) #生存確率、標準誤差、信頼区間
plot(kaplan) #Kaplan-Meier 曲線の図示

Surv コマンド、survdiff コマンドを使い、ログランク検定を行う。

survdiff(Surv(生存期間, 打ち切り) ~ 抗がん剤, data)

coxph コマンドを使うと、ログランク検定の結果に加え、Cox 比例ハザードモデル分析による尤度比検定、Wald 検定の結果も得られる。

cox = coxph(Surv(生存期間, 打ち切り) ~ 抗がん剤, data)
summary(cox)

Cox 比例ハザードモデルの説明変量を増やし、多変量解析とすることもできる。ここでは「全身状態」を説明変量に加えている。

cox.multi = coxph(Surv(生存期間, 打ち切り) ~ 抗がん剤 + 全身状態, data) 
summary(cox.multi) #結果の表示
AIC(cox.multi) #赤池情報量基準

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

prcomp() コマンドを使い、主成分分析を行う。また、factanal() コマンドを使い、因子分析を行う。

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

data = data.frame(手術時間 = c(60,50,100,90,30,40,80,90), 出血量 = c(15,100,90,20,18,20,120,90), 結紮失敗回数 = c(0,1,2,3,2,1,5,1))

data の内容を確認する。

data

prcomp() コマンドを使い、主成分分析を行う。

PCA = prcomp(data, scale = TRUE) #scale = TRUE でデータを標準化
PCA #固有ベクトルの表示
summary(PCA) #寄与率、累積寄与率の表示
PCA$x #主成分得点の表示
biplot(PCA) #結果の図示

princomp() コマンドを使っても、主成分分析を行えるが、prcomp()コマンドとは結果が微妙に異なる。

cor() コマンドとeigen() コマンドを使って固有値と固有ベクトルを表示することもできる。

eigen(cor(data))$values #固有値
eigen(cor(data))$vectors #固有ベクトル

factanal() コマンドを使い、因子分析を行う。ここでは因子数は 1としている。デフォルトで Varimax 回転が選択されている。

factanal(data,1) 

R によるガンマ回帰分析

glm() コマンドを使い、ガンマ回帰分析を行う。一般化線形モデルとして行う。

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

data = data.frame(出血量 = c(10,200,300,5,5,3,250,500), 手術方法 = c("内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"), 手術時間 = c(60,50,60,60,30,40,60,50))

data の内容を確認する。

data

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

result = glm(出血量 ~ 手術方法 + 手術時間, data=data, family=Gamma)
summary(result) #結果の表示

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

result = glm(出血量 ~ 手術方法 + 手術時間, data=data, family=Gamma(log))
summary(result) #結果の表示

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 #赤池情報量基準

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

glm() コマンドを使い、ポワソン回帰分析を行う。一般化線形モデルとして行う。

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

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

data の内容を確認する。

data

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

result = glm(重大インシデント数 ~ 残業80時間以上の医師数 + タイムアウトの導入, data=data, family=poisson)
summary(result) #結果の表示

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 #赤池情報量基準

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

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

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

data = data.frame(再発 = c("yes","yes","yes","no","no","no","yes","no","yes","yes"), 化学療法回数 = c(6,5,4,3,3,3,1,0,0,0), 手術時腫瘍残存 = c("yes","yes","yes","no","no","no","no","no","yes","yes"))

data の内容を確認する。

data

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

result = glm(再発 ~ 化学療法回数 + 手術時腫瘍残存, data=data, family=binomial)
summary(result) #結果の表示
exp(coef(result)) #オッズ比
exp(confint(result)) #オッズ比の95%信頼区間
exp(confint(result, level = 0.99)) #オッズ比の99%信頼区間

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%信頼区間の数値をひとつひとつ入れて計算する。

R による重回帰分析

lm() コマンドを使い、重回帰分析を行う。

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

data = data.frame(出血量 = c(10,200,300,5,5,3,250,500), 手術方法 = c("内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"), 手術時間 = c(60,50,60,60,30,40,60,50))

data の内容を確認する。

data

lm() コマンドを使い、重回帰分析を行う。単回帰分析を行うこともできる。

result = lm(出血量 ~ 手術方法 + 手術時間, data)
summary(result) #結果の表示
AIC(result) #赤池情報量基準

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 #赤池情報量基準