理系夫婦のうたたねブログ

理系夫婦が好きなことを書いていきます。たまに医学っぽいことを書いていますが、あくまで私見です。

Rを用いてt検定

はじめに

大学が変わったのと年度が変わったことで、統計ソフトとグラフ作成用のソフトを使えなくなった。
これを機に統計とグラフ作成を両方Rにしようと思う。
t testのやり方もよく分かっていないので、Helpを見ながら勉強してみた。

2標本のt検定

Rでのt testはt.test()を使えば良い。
例えば1~10の数列と7~20までの数列でt testしてみる。

t.test(c(1:10), y=c(7:20))
     Welch Two Sample t-test
 
 data:  c(1:10) and c(7:20)
 t = -5.4349, df = 21.982, p-value = 1.855e-05
 alternative hypothesis: true difference in means is not equal to 0
 95 percent confidence interval:
  -11.052802  -4.947198
 sample estimates:
 mean of x mean of y 
       5.5      13.5

出力にもあるとおり、デフォルトでwelch (等分散性を仮定しなくても使えるt test)になっているようだ。
あとはお馴染みの結果が返ってきている。
このようにt.test() に2群のベクタを渡せばt検定できるみたい。
簡単でよろしい。
df (自由度) が小数になっているのはwelchだから。

試しに、サンプルデータでもやってみる。
以下は前にどこかで使ったsleepデータ
勉強中のggplotで図示してみる。

testData <- sleep
ggplot(data = testData, aes(x=group, y = extra))+
  geom_boxplot()+
  theme_classic()

適当なプロット

前も

もう一つの書き方 fomula形式

t検定の時に、ベクタを渡す形式以外に、式で書くこともできる。 Fomula形式と呼ばれている。
従属変数 ~ 独立変数
というような形で記載することが可能であり、この際にはデータの指定が必要。

t.test(extra ~ group, data = testData)
    Welch Two Sample t-test

data:  extra by group 
t = -1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval: -3.3654832  0.2054832
sample estimates: 
mean in group 1 mean in group 2 
            0.75            2.33

このFomula形式で書くものといえば線形回帰が思いつく

res <- lm(extra ~ group, data = testData)
summary(res)
Call:
lm(formula = extra ~ group, data = testData)
 
Residuals:
   Min     1Q Median     3Q    Max 
-2.430 -1.305 -0.580  1.455  3.170 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.7500     0.6004   1.249   0.2276  
group2        1.5800     0.8491   1.861   0.0792 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
Residual standard error: 1.899 on 18 degrees of freedom
Multiple R-squared:  0.1613,    Adjusted R-squared:  0.1147 
F-statistic: 3.463 on 1 and 18 DF,  p-value: 0.07919

結果のGroup2のPrはStudentのt検定の結果と等しいはず。
t検定, ANOVA、線形回帰には本質的な違いはない (というか全て一般線形モデルの範疇内の) はずなのでこれで良い。

対応のあるt検定

次に対応のあるt検定をやってみる。
元のformatでは少しやりにくいが、練習と思ってやってみる。
dataの中からGourpが1, 2であるものをfilterで取り出して、格納し、対応のあるt検定

group1 <- testData %>% filter(group == "1")
group2 <- testData %>% filter(group == "2")
t.test(group1$extra, group2$extra, paired = T)
 Paired t-test

data:  group1$extra and group2$extra
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean difference 
          -1.58

対応のあるt検定は差の1標本のt検定なので以下のように書いても良いはずで、実際に一致している。

t.test(group1$extra - group2$extra)
    One Sample t-test

data:  group1$extra - group2$extra
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of x 
    -1.58

感想

これなら簡単に使えそう。
one way ANOVAとか他の検定も簡単に使えそうなので、使いそうな検定を勉強していこうと思う。

DeepLabCutのエラー「ValueError: too many values to unpack (expected 2)」

数年ぶりにDeepLabCutを使ってみたらいろいろ使い勝手が変わっていてびっくりしました。
そんな中でであったエラーの解決法を記しておきます。

経緯

数年ぶりにDeepLabCutを使ってみている途中、学習までは無事に終了して、トラッキングした動画を作ろうとしたら
ValueError: too many values to unpack (expected 2)
というエラーがでて、動画が作成されませんでした。

その前日に試したときはエラーなく動画が作成されたのですが、何が違うのだろう・・・?
となって1時間くらい悩みました。

原因と解決方法

さすがに使用している人口が多いためか、同じエラーが起こった人がいて、丁寧にgithubでissueになっていました。
https://github.com/DeepLabCut/DeepLabCut/issues/1081
こちらを読めば簡単にわかると思いますが、日本語解説はなかったので、一応まとめておきます。

原因はskeltonの指定方法の間違いでした。
DeepLabCutで検出するポイントはconfig.yamlを書き換えて作成しますが、その中には
bodyparts という部分と
skeleton という部分があります。 上はどの体の部位をトラックするかで、下はスケルトン化の際のパーツを示しています。

エラーが出た際私は、
skeltonのほうもBodypartsと全く同じにしていました。

どうやら、Skeltonのほうはskeltonになった場合の1本の線ごとに区切って書く必要があるらしく、
bodypartのように同定したいポイントをただ単に羅列するのではなく、
Nose-Center
Center-Tail
という風に分けて書く必要があったらしいです。今回の場合NoseからCenterに線をひいて、CenterからTailに線を引いているので、
結果として上のように書く必要があったということです。

結論

マニュアルをちゃんと読もう。
これにつきますね。

Endnoteの「描画の対象とならないオブジェクト」というエラーの対処法

調べてみると意外にもこのエラーに対する対処法が出てこないので記載しておきます。

経緯と症状

身近でEndnoteを使い始めた方が、このエラーに直面していました。
引用文献を入れようとすると上記のエラーがでて、うまく引用文献を挿入できていませんでした。

不思議なことに、新規ファイルを作成して挿入してみるとうまくいくのだそうです。

原因

OneDriveのせい (多分)
当該ファイルがOnedrive上に保存されていたので、ローカルに保存してみるとうまくいきました。
Endnote公式も、以下のようにオンラインストレージへの対応をしていないといっていますし、文章はローカルに保存するのが無難ですね。
google driveと併用していた時に、ほかにもいろいろ変な動きをしたことがあった気がします。
https://www.usaco.co.jp/faq/detail.html?pdid1=180

結論

Endnoteとオンラインストレージは相性が悪いので、注意
(Onedriveはいろいろ起こすので個人的にはあまりこのみではなかったり・・・)

Rでカプランマイヤー曲線をかく (ggsurvplot)

うだるような暑さですが、なんとか頑張ってます。
カプランマイヤーもggpubrで書けたら良かったのですが、ちょっと調べた感じは書けなそうだったので、ggsurvplotという違うライブラリを使ってみました。

はじめに

生存時間解析 (文字通りの生存時間以外にも、イベントが起こるまでの時間の推定などの時に使う) を行う際によく出てくるカプランマイヤー曲線を描いてみます。
カプランマイヤー曲線は生存率の推定値を表しており、曲線を比較することで議論をします。
Pythonで描いていた時はlifelines

lifelines.readthedocs.io

をよく使っていました。グラフ作成をRに統一するなら、これもRでかけたほうが良いなと思って探したところ、ggsurvplotで書けるようです。
少し使った感じは、やはり簡単に書けるなあと言ったところですが、カスタマイズが少しやりにくかったかな? あとドキュメントの充実具合はlifelinesの方が上な気もします。

データ

使用するデータはLungデータセットです。 肺がん患者さんの生存期間データですね。みた感じはこんな感じのデータです。

データの見た目

重要な変数は time:  Survival time in days status: censoring status 1=censored, 2=dead sex:   Male=1 Female=2 この辺りですね。 このデータを使ってカプランマイヤーを書いていきます。

処理と結果

まずは、ライブラリの読み込み

library(survival)
library(survminer)

続いて、カプランマイヤー推定量の計算のために、survfitを行います。 いくつも変数がありますが、簡単のために最初は性別のみで分けてみます。

fit<- survfit(Surv(time, status) ~ sex, data = lung)

最後に、グラフ作成を行います。

ggsurvplot(fit, data = lung)

シンプルなグラフ

これだけでも比較的綺麗なグラフになります。

ここからいくつかの要素を足したりして、整えていきます。

カプランマイヤー曲線に信頼区間を表示する。

ggsurvplot(fit, data = lung, conf.int = TRUE)

信頼区間つき

生存期間中央値のラインを引く

ggsurvplot(fit, data = lung,
           conf.int = TRUE,
           surv.median.line = "hv")

生存期間中央値を表示

リスク表を出す。 (各時点でのat riskの患者数)

ggsurvplot(fit, data = lung,
           conf.int = TRUE,
           surv.median.line = "hv",
           risk.table = TRUE)

risk table 付き

カプランマイヤー曲線とリスク表の大きさの比率を変えてみる。

ggsurvplot(fit, data = lung,
           conf.int = TRUE,
           surv.median.line = "hv",
           risk.table = TRUE,
           surv.plot.height = 0.6,
           tables.height = 0.4)

リスク表大き目

グラフの軸ラベルフォントを変更

ggsurvplotはggparの引数を受け取れるようなので、そのまま入れると変更可能

ggsurvplot(fit, data = lung,
           conf.int = TRUE,
           conf.int.style = "step",
           surv.median.line = "hv",
           surv.plot.height = 0.75,
           tables.height = 0.25,
           xlab = "Survival rate", 
           ylab = "Days",
           font.main = c(16, "bold"),
           font.x = c(16, "bold"),
           font.y = c(16, "bold"),
           font.tickslab = c(14,"bold"))

軸ラベル変更
本当はrisk tableのフォントも変えられるけど、面倒なのでやっていない。 →やってみた。

ggsurvplot(fit, data = lung,
           conf.int = TRUE,
           risk.table = TRUE,
           conf.int.style = "step",
           surv.median.line = "hv",
           surv.plot.height = 0.75,
           tables.height = 0.25,
           xlab = "Survival rate", 
           ylab = "Days",
           font.main = c(16, "bold"),
           font.x = c(16, "bold"),
           font.y = c(16, "bold"),
           font.tickslab = c(14,"bold"),
           tables.theme = theme_survminer(font.x = c(16, "bold"), 
                                          font.y = c(16, "bold")))

tableのフォントも変更

感想

他にも色々いじれそうだけど、短時間で分かったのはこの辺りです。 また分かったら追記しようと思います。

ggpubrで散布図

はじめに

何回か続いているggpubrのシリーズです。 今回は散布図について書きます。2つの連続変数の関係性を見る際の最も基本的なプロットですね。

データの紹介

今回はcarsというデータを用います。こちらは車の速度と制動距離 (止まるまでに必要な距離) のデータです。

処理と結果

data(cars)
ggscatter(cars, x = "speed", y = "dist",
          color = "black")

基本の散布図
こんな感じです。

例によって、ggparで見た目を変えることができて

p<-ggscatter(cars, x = "speed", y = "dist",
             color = "darkgray", size = 2)
ggpar(p,
      xlab = "Speed", 
      ylab = "Distance",
      font.main = c(18),
      font.x = c(18, "bold"),
      font.y = c(18, "bold"),
      font.tickslab = c(14,"bold"))

少し見た目を変えたグラフ

また、回帰直線を簡単に示すこともできます。

p<-ggscatter(cars, x = "speed", y = "dist",
             color = "black", size = 2,
             add = "reg.line", 
             add.params = list(color = "skyblue", fill = "lightgray"),
             conf.int = TRUE,
             cor.coef = TRUE)
ggpar(p,
      xlab = "Speed", 
      ylab = "Distance",
      font.main = c(18),
      font.x = c(18, "bold"),
      font.y = c(18, "bold"),
      font.tickslab = c(14,"bold"))

回帰直線と回帰係数

ここでは95%信頼区間 (conf.int)と回帰係数 (cor.coef) も出してみています。

感想

やはり簡単ですね。 もう少し複雑なデータもやってみたいです。

ggpubrで箱ヒゲ図、バイオリンプロット

はじめに

ggpubrで箱ヒゲ図とバイオリンプロットを書いてみます。 私としては箱ひげ図の方が好きです。バイオリンプロットはなんとなく、KDEでの推定があっているのかわからない時があって、誤解を与えそうな気がするからです。

データの紹介

今回使うのはToothGrowth というデータセットです。 モルモットにビタミンCを投与した際の歯の長さについてのデータセットだそうです。 投与方法が2つあって、オレンジジュースとして投与したデータと、精製したアスコルビン酸で投与したデータがありますが、今回はまとめて表示しています。

処理と結果

箱ヒゲ図

まずは何も考えずに箱ヒゲ図を作ってみます。

data("ToothGrowth")
df <- ToothGrowth
ggboxplot(df, x = "dose", y = "len", width = 0.8)

簡単に箱ヒゲ図

これだけだと、味気ないのと、データ点を全て出した方が誠実だと思われるので、その様にしてみます。 doseごとに色を変えるためにcolor = dose データ点を出すためにadd = jitter とすると

ggboxplot(df, "dose", "len",
          color = "dose", palette =c("blue", "red", "magenta"),
          add = "jitter", shape = "dose")

色を変えたプロット

となります。

塗りつぶしたい時はcolor = doseではなくfill = doseにすればよく

ggboxplot(df, "dose", "len",
          fill = "dose", palette ="Blues",
          add = "jitter", shape = "dose")

塗りつぶし

となります。こちらの方が良いかもしれません。

さらに、調整を行ってみます。

p<-ggboxplot(df, x = "dose", y = "len",
         add = "dotplot",
         add.params = list(size = 0.5, fill = "white"),
         fill = "dose",
         palette = "Blues") + theme(aspect.ratio = 1.3)
ggpar(p,
      xlab = "Dose", 
      ylab = "Tooth length",
      font.main = c(18),
      font.x = c(18, "bold"),
      font.y = c(18, "bold"),
      font.tickslab = c(14,"bold"),
      legend = "right",
      legend.title = "Condition")

ちょっと調整したグラフ

ここでは、Boxplotにaddでドットプロット (各データ点を白丸で付け加える) を追加しています。ドットプロットの調整はadd.paramsで行っており、丸の大きさと塗りつぶしの色を変更しています。

さらに+theme(aspect.ratio = 1.3) でいい感じの縦横比になるように調整を行っています。

ggparでは軸ラベルを調整しています。こないだの記事の通りです。

med-ruka.hatenablog.com

バイオリンプロット

上記のグラフは全て、ggboxplotをggviolinにすればバイオリンプロットに変更可能です。

p<-ggviolin(df, x = "dose", y = "len",
             add = "dotplot",
             add.params = list(size = 0.5, fill = "white"),
             fill = "dose",
             palette = "Blues") + theme(aspect.ratio = 1.3)
ggpar(p,
      xlab = "Dose", 
      ylab = "Tooth length",
      font.main = c(18),
      font.x = c(18, "bold"),
      font.y = c(18, "bold"),
      font.tickslab = c(14,"bold"),
      legend = "right",
      legend.title = "Condition")

バイオリンプロット
やっぱり私は箱ひげのほうが好みですが。

 感想

やっぱりggpubrは簡単でいいですね。 しばらくは使ってみようかと思います。ただ、ライブラリにないグラフは結局ggplot2を使う気がするので、ggplot2に習熟する方が楽??
あと、p-valueとかつける方法もあるみたいですが、個人的にはイラレでつける方が結局楽かなと思います。
なのでフローとしては
ミニマルな図をRで作成 → イラレで装飾をたす
みたいな感じですね今のところ。

ggpubrでのグラフの整え方

ggpubrで作成されるグラフの整え方を少し紹介します。

はじめに

ggpubrでデフォルトで出てくるグラフはある程度整ってはいますが、自分好みに整えるのが必要な部分もあります。 昨日のggpairedを題材に、少しいじってみます。

データの紹介

使うデータとコードは昨日の記事と同じです。

med-ruka.hatenablog.com

sleep dataを用いて少しやってみます。

処理と結果

# ライブラリのインポートとデータ取得
library(ggpubr)
data(sleep)
# 簡単なプロット
ggpaired(data = sleep, x = "group", y = "extra")

簡単なプロット

ggpairedにはいくつかの引数があり、指定することで見た目を変えることができます。 前回は色を変えたりしていましたが、それ以外にbox plotの幅を変えたり (width)、対応するサンプル間に引く線の太さ (line.size) を変えられます。

ggpaired(data = sleep,
         x = "group",
         y = "extra",
         width = 0.8,
         line.size = 0.1,
         )

幅と線を変更

一方で軸ラベルの字の大きさや、グラフタイトルの追加にはggparを用います。 使い方は、

p <- ggpaired(data = sleep,
              x = "group",
              y = "extra",
              width = 0.5,
              line.size = 0.5,
              linetype = "dashed")
ggpar(p,
      main = "Effects of drugs on sleep amount",
      xlab = "Group", 
      ylab = "Sleep amount",
      font.main = c(24,"bold.italic", "blue"),
      font.x = c(16, "bold"),
      font.y = c(16, "bold"))

のようにするみたいです。 これで以下のように、軸ラベルなどをいじったグラフが出せます。 font.mainなどのc() の中身は font size, font type, font colorになっています。

ラベルなどを変更したグラフ

p<- ggpaired(d,
             cond1 = "drug1",
             cond2 = "drug2",
             fill = "condition",
             line.color = "gray",
             palette = "Blues")
ggpar(p,
      xlab = "Group", 
      ylab = "Sleep amount",
      font.main = c(18),
      font.x = c(18, "bold"),
      font.y = c(18, "bold"),
      font.tickslab = c(14,"bold"),
      legend = "right",
      legend.title = "Condition")

Conditionの情報が冗長ですが、こんな感じで整えていくことが可能です。

ちょっと整えたグラフ

感想

ggpubrは簡単にグラフを作るのには便利だなーと思います。 同僚に説明したりする程度であればこのくらいのグラフでも全く問題ない気がしますし。 論文に となるとまたちょっと変わる気もしますが、もうちょっと修行したらもっとよくかけるかもしれません。 R はフリーなのが大変助かるところで、これからさらに慣れていきたいところです。