タグ: R

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

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) #結果の表示

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) #結果の表示

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%信頼区間

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

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

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)