R による連続確率分布と累積分布の描画

正規分布、t 分布、カイ2乗分布、F 分布、ベータ分布、ワイブル分布、一様分布、ガンマ分布、コーシー分布、指数分布、対数正規分布、べき分布(パレート分布)を描画する(Distributions)。

正規分布

正規分布と標準正規分布を描く。

#正規分布。平均50、標準偏差5。
curve(dnorm(x,50,5),30,70) #確率分布
curve(qnorm(x,50,5)) #累積分布
#標準正規分布
curve(dnorm, -5,5) #確率分布
curve(qnorm) #累積分布

t 分布

ここでは自由度を1としている。

curve(dt(x,1),-1,1) #確率分布
curve(qt(x,1)) #累積分布

カイ2乗分布

ここでは自由度を3としている。

curve(dchisq(x, 3)) #確率分布
curve(qchisq(x, 3)) #累積分布

F 分布

ここでは自由度を5と10としている。

curve(df(x,5,10)) #確率分布
curve(qf(x,5,10)) #累積分布

ベータ分布

ここではパラメータを0.3と0.5としている。

curve(dbeta(x,0.3,0.5)) #確率分布
curve(qbeta(x,0.3,0.5)) #累積分布

ワイブル分布

ここでは shape = 0.8, scale = 2 としている。

curve(dweibull(x,0.8,2)) #確率分布
curve(qweibull(x,0.8,2)) #累積分布

一様分布

curve(dunif(x,-1,1),-1.1,1.1) #確率分布
curve(qunif(x,-1,1)) #累積分布

ガンマ分布

ここでは shape = 2, scale = 1.5 としている。

curve(dgamma(x,2,1.5),0,5) #確率分布
curve(qgamma(x,2,1.5)) #累積分布

コーシー分布

ここでは location = 0, scale = 1 の標準コーシー分布としている。

curve(dcauchy(x,0,1),-1,1) #確率分布
curve(qcauchy(x,0,1)) #累積分布

指数分布

ここでは rate = 2 としている。

curve(dexp(x,2)) #確率分布
curve(qexp(x,2)) #累積分布

対数正規分布

ここではパラメーターを 1, 2 としている。

curve(dlnorm (x,1,2),1,3, log = "xy") #確率分布
curve(qlnorm (x,1,2)) #累積分布

べき分布(パレート分布)

ここでは a を shape、b を scale として関数を定義し、a = 2, b = 1 としている。

#確率密度関数
Fd = function(x, a, b) a/b * (b/x)^(a+1) 
curve(Fd(x,2,1),1,10, log = "xy")
#累積分布関数
Fq = function(x, a, b) 1-(b/x)^a
curve(Fq(x,2,1))

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 による Tukey-Kramer test

TukeyHSD() コマンドで Tukey-Kramer test を行う。

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

data = data.frame(デバイス = c("デバイス A","デバイス B","デバイス A","デバイス A","デバイス C","デバイス C","デバイス B","デバイス C","デバイス A","デバイス A"), 出血量 = c(30,20,5,5,100,80,10,60,10,25))

data の内容を確認する。

data

デバイス間で出血量の平均に差があるかどうかを TukeyHSD() コマンドで検定する。

TukeyHSD(aov(data$出血量 ~ data$デバイス)) #95%信頼区間
TukeyHSD(aov(data$出血量 ~ data$デバイス), conf.level = 0.99) #99%信頼区間

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

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 での結果とは微妙に異なる

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 による基本統計量の計算

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

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

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

data の内容を確認する。

data

summary() コマンドで平均などの基本統計量を表示する。

summary(data)

tapply() コマンドで、手術時間の統計量を術者毎に出す。

tapply(data$手術時間, data$術者, summary)

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

nrow(data) #行数 = サンプルサイズ
table(data$術者) #それぞれの術者が何件の手術を行っているかカウント
sum(data$手術時間) #手術時間の合計
mean(data$手術時間) #手術時間の平均
var(data$手術時間) #手術時間の不偏分散
sd(data$手術時間) #手術時間の標準偏差(不偏分散の平方根)
max(data$手術時間) #手術時間の最大値
min(data$手術時間) #手術時間の最小値
median(data$手術時間) #手術時間の中央値

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

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

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

data = data.frame(術者 = c("太郎","太郎","花子","花子"), 手術方法 = c("内視鏡","開腹","内視鏡","開腹"), 回数 = c(60,5,20,40))

data の内容を確認する。

data

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

table = tapply(data$回数, list(data$術者, data$手術方法), sum)

table の内容を確認する。

table

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

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

R による Kruskal-Wallis test と ANOVA

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

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

data = data.frame(デバイス = c("デバイス A","デバイス B","デバイス A","デバイス A","デバイス C","デバイス C","デバイス B","デバイス C","デバイス A","デバイス A"), 出血量 = c(30,20,5,5,100,80,10,60,10,25))

data の内容を確認する。

data

デバイスによって出血量に偏りがあるかどうかを kruskal.test() コマンド、oneway.test() コマンドで検定する。

kruskal.test(出血量 ~ デバイス, data = data) #Kruskal-Wallis test
oneway.test(出血量 ~ デバイス, data = data) #one-way ANOVA (Welch 拡張)
oneway.test(出血量 ~ デバイス, data = data, var = T) #通常の one-way ANOVA

aov() コマンドや、anova() コマンドでも one-way ANOVA の計算ができる。

summary(aov(出血量 ~ デバイス, data = data))
anova(lm(出血量 ~ デバイス, data = data))

R によるクロス集計表

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

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

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

data の内容を確認する。

data

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

tapply(data$手術時間, list(data$術者, data$手術方法), mean) #手術時間の平均
tapply(data$出血量, list(data$術者, data$手術方法), sum) #出血量の合計

件数をみるだけの通常のクロス集計は table() コマンドを使って行う。

table(data$術者, data$手術方法)

R による相関係数の計算

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

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

data = data.frame(縫合結紮の練習時間 = c(100,150,200,250,300,1000), 手術時間 = c(210,220,150,100,90,45), 出血量 = c(50,70,60,30,20,10))

data の内容を確認する。

data

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

cor(data)

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

cor(data, method = "s") #スピアマン順位相関係数。
cor(data, method = "k") #ケンドール順位相関係数。

P 値も計算したい場合は cor.test() コマンドを使う。

x = 縫合結紮の練習時間
y = 手術時間
cor.test(x, y) #ピアソン積率相関係数。
cor.test(x, y, method = "s") #スピアマン順位相関係数。
cor.test(x, y, method = "k") #ケンドール順位相関係数。

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

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乗検定(イエーツ補正あり)

R による独立2群の比較

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

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

data = data.frame(手術方法 = c("内視鏡","開腹","開腹","内視鏡","内視鏡","内視鏡","開腹","開腹"), 手術時間 = c(60,50,100,90,30,40,80,90), 出血量 = c(15,100,90,20,18,20,120,90))

data の内容を確認する。

data

手術方法(内視鏡か開腹か)によって data を2群にグループ分けする。

A = data[data$手術方法 == "内視鏡",]
B = data[data$手術方法 == "開腹",]

A 群、B 群の手術時間と出血量を Student t test で比較する。P 値だけでなくA 群、B 群それぞれの平均や、その差の信頼区間等も計算されて表示される。

t.test(A$手術時間, B$手術時間, var.equal = T) #手術時間
t.test(A$出血量, B$出血量, var.equal = T) #出血量

古典的な Student t test 以外は、以下のようにする。

t.test(A$手術時間, B$手術時間) # Welch t test
wilcox.test(A$手術時間, B$手術時間) # Wilcoxon rank sum test

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

R による信頼区間の計算

t.test() コマンドで平均と信頼区間を計算する。

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

data = data.frame(手術時間 = c(60,50,100,90,30,40), 出血量 = c(15,25,35,55,45,25))

data の内容を確認する。

data

ここでは手術時間の平均の95%信頼区間を計算してみる。t 値や自由度や平均なども表示される。

t.test(data$手術時間)

以下のようにすれば99%信頼区間となる。

t.test(data$手術時間, conf.level = 0.99)

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() #中央値