今日はSQLとRをつないで解析を行うことを勉強しました。主に習ったのはRのパッケージのTableOneについてです。下のページのPDFを使いました。10ページ目にCreateTableOneがあります。Rのパッケージにはexampleがついているので、コピペして1行ずつやっていくのが練習には良さそうです。
https://cran.r-project.org/web/packages/tableone/tableone.pdf
TableOneは「論文などの患者背景を示すTable 1を作るためのパッケージ」です。PDFには以下のようなExampleが書いてありますので、それを実際にRで動かして1文ずつ見ていきました。##はTableOneを作られた方のコメントです。
intが整数値で、numは数値型(実数値)、Factorが因子型、ここには無いですがcharacterが文字列型です。本来Factorで表すべきstatusやtrt、accitesがintで入っていますね。
print(tableOne, nonnormal = c("bili","chol","copper","alk.phos","trig"),exact = c("status","stage"), cramVars = "hepato", smd = TRUE)
・printはただ出力するだけでなく、以下のような条件をつけて出力する。
・nonormal:非正規分布、平均・SDを出すのでなく中央値と四分位を出す。
・カテゴリー変数は何も設定しないとカイ二乗検定が行われます。Fisherの正確確率検定したいもの(標本数が少ないfactor)はexactにいれましょう。
ちなみにRのスクリプトはATOMで書いてからコピペすることを勧められました。Rで色々書くのは難しいですよね。RStudio使えうのも1つです。ATOMは便利なテキストエディタで、勉強するには以下のサイトのAtom入門などがいいかもしれません。
その後、ロジスティク回帰を試しにやってみました。
model1<-glm(death~sex+age+c_jcs+ambul+tPA,family=binomial,data=res)
・始めのdeathが従属変数。〜以下が説明変数
・familyの部分は引数選択:binomoal(二項分布)、正規分布(gaussian)、ポアソン分布(poisson)、逆正規分布(inverse.gaussian)、ガンマ分布(Gamma)など。
summary(model1)
・結果を見る
ModelName <- model1
alpha <- 0.05
x <- summary(ModelName)
y <- confint(ModelName, level=1-alpha)
OR <- exp(x$coefficients[,1])
LowerCL <- exp(y[,1])
UpperCL <- exp(y[,2])
pvalue <- summary(model1)$coefficients[,"Pr(>|z|)"]
z <- cbind(OR, LowerCL, UpperCL, pvalue)
round(z, 5)
・オッズ比に変える
<その他 R>
res2<-subset(res,res$sex==1)
・resというデータからsex==1のみres2という名前で抽出する。
・Rでは=は==になる。プログラム言語では=は代入。<-と一緒。
!=
・Rでの≠を意味する。
まだ大まかにしかわかっていないですが、今後勉強していこうと思います。
<自習>
・EFファイルとjoinする場合は事前にcodeで分けておく。
SELECT T1.qhosp
,T1.ptid
,T1.dtadm
,MAX(CASE WHEN T3.code = 'ICU' THEN 1 ELSE 0 END) AS ICU
,SUM(CASE WHEN T3.code = 'ICU' THEN 1 ELSE 0 END) AS ICU_stay
,MIN(CASE WHEN T3.code = 'ICU' THEN T2.ef24 END) AS first_ICU
,MAX(CASE WHEN T3.code = 'ER' THEN 1 ELSE 0 END) AS ER
,MIN(CASE WHEN T3.code = 'culture' THEN T2.ef24 END) AS first_culture
FROM #p2 AS T1
LEFT JOIN[~] AS T2
ON T1.qhosp = T2.qhosp AND T1.ptid = T2.ef2 AND T1.dtadm = T2.ef4
INNER JOIN [~]
AS T3 ON T2.ef9 = T3.ef9
GROUP BY T1.qhosp,T1.ptid,T1.dtadm
0,1でMAXにするとありなしになる。SUMにすると日数になる。日付でMINにすると最初の日にちになる。
・<>:等しくないことを示す。
・2つのグループの差を取る。
select * into #a4 from #a2
except
select * from #pp2
https://cran.r-project.org/web/packages/tableone/tableone.pdf
TableOneは「論文などの患者背景を示すTable 1を作るためのパッケージ」です。PDFには以下のようなExampleが書いてありますので、それを実際にRで動かして1文ずつ見ていきました。##はTableOneを作られた方のコメントです。
## Load
library(tableone)
・tableoneを起動する
library(tableone)
・tableoneを起動する
## Load Mayo Clinic Primary Biliary Cirrhosis Data
library(survival)
・survivalを起動する
data(pbc)
・pbcのデータ全てを見る
## Check variables
head(pbc)
・pbcのデータの先頭5つくらいを見る
str(pbc)
・pbcのデータの変数を見てみる。今回のデータ(pbc)では以下のようになります。i
library(survival)
・survivalを起動する
data(pbc)
・pbcのデータ全てを見る
## Check variables
head(pbc)
・pbcのデータの先頭5つくらいを見る
str(pbc)
・pbcのデータの変数を見てみる。今回のデータ(pbc)では以下のようになります。i
intが整数値で、numは数値型(実数値)、Factorが因子型、ここには無いですがcharacterが文字列型です。本来Factorで表すべきstatusやtrt、accitesがintで入っていますね。
## Make categorical variables factors
varsToFactor <- c("status","trt","ascites","hepato","spiders","edema","stage")
・連続変数でないカテゴリーのFactor(varsToFactor)にどの因子を入れるかを決める。
pbc[varsToFactor] <- lapply(pbc[varsToFactor], factor)
・pbcのデータの変数をfactor化する。str(pbc)で確認するとできてますね。
varsToFactor <- c("status","trt","ascites","hepato","spiders","edema","stage")
・連続変数でないカテゴリーのFactor(varsToFactor)にどの因子を入れるかを決める。
pbc[varsToFactor] <- lapply(pbc[varsToFactor], factor)
・pbcのデータの変数をfactor化する。str(pbc)で確認するとできてますね。
## Create a variable list
dput(names(pbc))
・変数リストを""付きで作る。下のlistvarsを作る際に使う。
listvars <- c("time","status","age","sex","ascites","hepato",
dput(names(pbc))
・変数リストを""付きで作る。下のlistvarsを作る際に使う。
listvars <- c("time","status","age","sex","ascites","hepato",
"spiders","edema","bili","chol","albumin",
"copper","alk.phos","ast","trig","platelet",
"protime","stage")
・listvarsに以下の変数を入れる
・listvarsに以下の変数を入れる
## Create Table 1 stratified by trt
tableOne <- CreateTableOne(vars = listvars, strata = c("trt"), data = pbc)
・tableOneをlistvarsを変数にtrtで層別化して作成する。
tableOne <- CreateTableOne(vars = listvars, strata = c("trt"), data = pbc)
・tableOneをlistvarsを変数にtrtで層別化して作成する。
## Just typing the object name will invoke the print.TableOne method
tableOne
・tableOneを見る。以下のように出力される。
全て平均とSDで表示されている。
tableOne
・tableOneを見る。以下のように出力される。
全て平均とSDで表示されている。
## Specifying nonnormal variables will show the variables appropriately,
## and show nonparametric test p-values. Specify variables in the exact
## argument to obtain the exact test p-values. cramVars can be used to
## show both levels for a 2-level categorical variables.
## and show nonparametric test p-values. Specify variables in the exact
## argument to obtain the exact test p-values. cramVars can be used to
## show both levels for a 2-level categorical variables.
print(tableOne, nonnormal = c("bili","chol","copper","alk.phos","trig"),exact = c("status","stage"), cramVars = "hepato", smd = TRUE)
・printはただ出力するだけでなく、以下のような条件をつけて出力する。
・nonormal:非正規分布、平均・SDを出すのでなく中央値と四分位を出す。
・カテゴリー変数は何も設定しないとカイ二乗検定が行われます。Fisherの正確確率検定したいもの(標本数が少ないfactor)はexactにいれましょう。
・基本設定では2つのカテゴリーの場合、1つの要素しか出してくれません。どちらも出したい場合はcramVarsにいれましょう。
・smd=trueにすると2つ以上のグループがある場合、すべてのペアワイズ比較の標準化された平均差が計算されます。
すると以下のような表が出力されます。
・smd=trueにすると2つ以上のグループがある場合、すべてのペアワイズ比較の標準化された平均差が計算されます。
すると以下のような表が出力されます。
## Use the summary.TableOne method for detailed summary
summary(tableOne)
・tableOneを作るのに裏で動いていた数値が全て出ます。
summary(tableOne)
・tableOneを作るのに裏で動いていた数値が全て出ます。
## See the categorical part only using $ operator tableOne$CatTable
summary(tableOne$CatTable)
・カテゴリー変数のみでます
summary(tableOne$CatTable)
・カテゴリー変数のみでます
## See the continuous part only using $ operator tableOne$ContTable
summary(tableOne$ContTable)
・連続変数が出ます、
summary(tableOne$ContTable)
・連続変数が出ます、
## If your work flow includes copying to Excel and Word when writing manuscripts,
## you may benefit from the quote argument. This will quote everything so that
## Excel does not mess up the cells.
## you may benefit from the quote argument. This will quote everything so that
## Excel does not mess up the cells.
print(tableOne, nonnormal = c("bili","chol","copper","alk.phos","trig"),exact = c("status","stage"), quote = TRUE)
・項目や数字全てに" "がつきます。これはエクセルに貼り付けてテキストファイルウィザードを使用する際にキレイに項目ごとに貼り付けるためのものです。
・項目や数字全てに" "がつきます。これはエクセルに貼り付けてテキストファイルウィザードを使用する際にキレイに項目ごとに貼り付けるためのものです。
## If you want to center-align values in Word, use noSpaces option.
print(tableOne, nonnormal = c("bili","chol","copper","alk.phos","trig"),exact = c("status","stage"), quote = TRUE, noSpaces = TRUE)
・中央揃えにする
・中央揃えにする
## If SMDs are needed as numericals, use ExtractSmd()
ExtractSmd(tableOne)
・SMDが数値として必要な場合に行う
ExtractSmd(tableOne)
・SMDが数値として必要な場合に行う
ちなみにRのスクリプトはATOMで書いてからコピペすることを勧められました。Rで色々書くのは難しいですよね。RStudio使えうのも1つです。ATOMは便利なテキストエディタで、勉強するには以下のサイトのAtom入門などがいいかもしれません。
その後、ロジスティク回帰を試しにやってみました。
model1<-glm(death~sex+age+c_jcs+ambul+tPA,family=binomial,data=res)
・始めのdeathが従属変数。〜以下が説明変数
・familyの部分は引数選択:binomoal(二項分布)、正規分布(gaussian)、ポアソン分布(poisson)、逆正規分布(inverse.gaussian)、ガンマ分布(Gamma)など。
summary(model1)
・結果を見る
ModelName <- model1
alpha <- 0.05
x <- summary(ModelName)
y <- confint(ModelName, level=1-alpha)
OR <- exp(x$coefficients[,1])
LowerCL <- exp(y[,1])
UpperCL <- exp(y[,2])
pvalue <- summary(model1)$coefficients[,"Pr(>|z|)"]
z <- cbind(OR, LowerCL, UpperCL, pvalue)
round(z, 5)
・オッズ比に変える
<その他 R>
res2<-subset(res,res$sex==1)
・resというデータからsex==1のみres2という名前で抽出する。
・Rでは=は==になる。プログラム言語では=は代入。<-と一緒。
!=
・Rでの≠を意味する。
まだ大まかにしかわかっていないですが、今後勉強していこうと思います。
<自習>
・EFファイルとjoinする場合は事前にcodeで分けておく。
SELECT T1.qhosp
,T1.ptid
,T1.dtadm
,MAX(CASE WHEN T3.code = 'ICU' THEN 1 ELSE 0 END) AS ICU
,SUM(CASE WHEN T3.code = 'ICU' THEN 1 ELSE 0 END) AS ICU_stay
,MIN(CASE WHEN T3.code = 'ICU' THEN T2.ef24 END) AS first_ICU
,MAX(CASE WHEN T3.code = 'ER' THEN 1 ELSE 0 END) AS ER
,MIN(CASE WHEN T3.code = 'culture' THEN T2.ef24 END) AS first_culture
FROM #p2 AS T1
LEFT JOIN[~] AS T2
ON T1.qhosp = T2.qhosp AND T1.ptid = T2.ef2 AND T1.dtadm = T2.ef4
INNER JOIN [~]
AS T3 ON T2.ef9 = T3.ef9
GROUP BY T1.qhosp,T1.ptid,T1.dtadm
0,1でMAXにするとありなしになる。SUMにすると日数になる。日付でMINにすると最初の日にちになる。
・<>:等しくないことを示す。
・2つのグループの差を取る。
select * into #a4 from #a2
except
select * from #pp2
コメント
コメント一覧 (1)
nagano1123
がしました