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

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

Rグラフに注釈追加

はじめに

最近この本で勉強しています。

グラフを描画する際に、一緒に文字や直線、短径を描画したい時があると思います。
これらはデータそのものを描画しているわけではないのですが、グラフを解釈する際に有用になります。
今回はRでグラフないに文字や直線、短径を描画する方法を記載します。

目的

グラフ内に文字や直線、短径を描画する。
以下のようなグラフを作れるようになる。

使用するデータのチェック

お馴染みのiris (アヤメ) データを使用します。
データの中身を確認します。
irisデータセットはデフォルトで入っているので、一行で読み込めます。

library(tidyverse)
iris

irisデータの概観

各行が1サンプルに対応していて、Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Speciesが与えられています。1
Species (アヤメの種類) 以外は数値データです。
ここでは、Petal.Length とPetal.Widthの関係に着目してみます。 この二つに正の相関があるのではないかという仮説をもとに可視化してみることにします。

データの前処理

最初に Spealのデータは今回使用しないので、selectで最初に除去します。

dat <- iris %>% 
  select(-c(Sepal.Length, Sepal.Width))

シンプルなグラフの描画

まずは単純に描画してみます。

p <- ggplot(dat, 
            mapping = aes(x = Petal.Width,
                          y = Petal.Length,
                          color = Species)) 
p + 
  geom_point() +
  theme_classic()

シンプルにプロット

ThemeはClassicが好みです。 このグラフに、文字や直線、短径を追加してみます。

グラフ要素の追加

グラフ要素の追加はannotateを使って行います。
annotateの最初の引数には、追加する要素の種類を指定します。 - "test": 文字
- "segment": 線分
- "pointrange": 点と線分
- "rect": 四角形

文字の追加

最初は文字を追加してみます。

p + 
  geom_point() +
  theme_classic() +
  annotate("text", x = 0.5, y = 5, label = "Scatter plot", color = "black", size = 5)

文字追加

x, y は文字の位置を指定します。文字の位置は軸の値で指定しています。
その他色やサイズも指定可能です。

線分の追加

適当な線分を追加する場合にはsegmentを用います。 この線分の使い道としては、グラフ上に閾値や平均値を示すといったところでしょうか。

p + 
  geom_point() +
  annotate("segment", x = 0, xend = 3, y = 4, yend = 4, color = "gray", size = 1)+
  annotate("segment", x = 1, xend = 1, y = 0, yend = 8, color = "gray", size = 1)+
  theme_classic()

線の追加

そのまま描画すると軸が浮いてしまって不格好なので、scale_x_continuous, scale_y_continuousでexpand = c(0, 0)を指定してきちんとくっつけます。

p + 
  geom_point() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  annotate("segment", x = 0, xend = 3, y = 4, yend = 4, color = "gray", size = 1)+
  annotate("segment", x = 1, xend = 1, y = 0, yend = 8, color = "gray", size = 1)+
  theme_classic()

軸を調節したグラフ

もしくは今回のデータではあまり意味がありませんが、y=xの直線を引くことで縦軸と横軸の値がどちらが大きいのかということを示すのにも使えます。

p + 
  geom_point() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  annotate("segment", x = 0, xend = 3, y = 0, yend = 3, color = "gray", size = 1)+
  theme_classic()

y=xを追加

短径の追加

rectでは四角形を描画できます。
文字と合わせれば、グラフ上で注目して欲しいところを表すのに使えそうです。

p + 
  geom_point() +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 2.7)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 7.5)) +
  annotate("text", x = 0.8, y = 2.7, label = "Petal of setosa is small.", color = "black", size = 5)+ 
  annotate("rect", xmin = 0, xmax = 1, ymin = 0, ymax = 2.5, fill = "gray", alpha = 0.2)+
  theme_classic()

短径の追加

終わりに

文字や直線、短径を追加することで、グラフの解釈がしやすくなります。
今まで論文を書く際にはイラストレータなどの他のソフトウェアで注釈を追加していましたが、Rで一気にできるのであれば作業の自動化が進みそうです。
グラフのデータ点に注釈をつける際にはまた別の方法がありますので、そちらはそちらで紹介します。


  1. Sepalはがく、Petalは花びら (正確な言葉遣いではないです) をさしています。

多重比較をRで実施する方法

はじめに

Rは細々続けています。

2群以上実験や複数水準のデータを取得したとき、多重比較をしたくなることがあると思います。

総当たりの多重比較は実験デザインとしてあまりよくないとはわかっていつつも、そういう場面があることもしばしばです。ということで、今回は多重比較をRでやる方法についてです。 参考は以下のページです。(ANOVAの部分をすっ飛ばしていますが)

ANOVA in R: The Ultimate Guide - Datanovia

最初に使用するパッケージを読み込んでおきます。

library(tidyverse)
library(ggpubr)
library(rstatix)
library(gt)
library(datarium)
library(emmeans)

前提

以下の解析はtidyverseのパッケージを用いています。 実際の解析で自分のデータを使う際にはtidy dataに変換してから実行が必要です。 tidy dataについてはたくさん解説があるので、そちらを参照 (例えばこちら) 私たちのR - 17  整然データ構造

一元配置デザインの場合

要するに3群以上の平均値について比較したい場合です。教科書的には (?) よく、ANOVAしてpost hoc みたいな言い方をしますが、別にANOVAは必須ではないと思われます。

使用するサンプルデータはplant growthで、3条件 (control, treatment1, treatment2) での植物の成長を見ています。 組み込まれているので、

data("PlantGrowth")

で読み込めます。

まずはグラフをチェックしてみます。(ggpubrだとこういうとき楽)

ggboxplot(PlantGrowth, x = "group", y = "weight", color = "group")

箱ヒゲ図

なんとなく群間に差がありそうな気がします。

Tukey-HSD

よく使うTukey-HSDはいかで実行可能です。

PlantGrowth %>% tukey_hsd(weight ~ group)

ちなみにgtパッケージを使うと以下のような整った表形式で表示可能です。

PlantGrowth %>% tukey_hsd(weight ~ group) %>% gt()

tukey

Tukey-HSDは論文でもよくみますが、独立、正規性、等分散性の3条件が必要なので、適応可能な範囲は比較的狭いと思われます。

Bonferroniとその変法たち

Bonferroniもよく見かける手法です。これも解説はいくらでもネットにあるので飛ばしますが、検定の回数で有意水準を割るというのが簡単な解説です。(逆にp-valueに検定の回数をかけて表示することも) そこから派生したHolmやShafferの方法は割る数の方を少しずつ減らしたりしている方法になります。 計算方法からわかるように、Bonferroniやその変法はかなり適応範囲が広いです。

pairwise_t_test(data=PlantGrowth,
                weight~group, pool.sd = TRUE,
                p.adjust.method = "holm") %>% gt()

ここではHolmをやってみて、gt()で表形式で出してみました。

holm

pairwise_t_testの引数は最初はデータ、2つ目はfomula形式になっています。 fomula形式は 独立変数 ~ 従属変数 という形式での記載になっております。この記載が変数間の関係性が明示的で私は好みです。 3つ目はpool.sdで、(間違っているかもしれませんが) 検定統計量を計算する際の分散としてプールされた分散を使用しているという意味だと思っております。 最後のp.adjust.methodで調整方法を選べて、ここではholmにしていますが、他にも"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"が選択可能です。(あまり使わないので、bonferroniとholm以外わかりませんが)

ちなみに違いを見るためにBonferroniも出してみます。

bonf

最もp-valueが小さい比較のp-valueは同じですが、他で違いがあります。holmはp-valueを小さい順に並べて、1つずつp-valueにかける値を減らします。そのため最初に評価される、最も小さいp-valueは同一になります、

二元配置デザインの場合

二元配置の場合、多重比較をやりたくなる場面としては交互作用があって、全群比較をしたくなった場合だと思います。 ここでは上記と異なるデータで、jobsatisfactionというデータを用います。

データの見た目
データはid列とgender列、education列とscore列からなります。 まずは簡単に図示してみます。

ggboxplot(
  jobsatisfaction, x = "gender", y = "score",
  color = "education_level"
)

jobsatisfaction

今回はANOVAの結果も載せておきます。

jobsatisfaction %>% anova_test(score ~ gender * education_level) %>% gt()

2 way ANOVA

interactionが有意になっています。 この場合、全群での比較をしたいわけですが、1元配置のときのpairwise_t_testはそのままだと使えません。 一度データをグループに分けてから使用する必要があります。

jobsatisfaction %>% 
  tukey_hsd(score~education_level + gender,
                  pool.sd = TRUE,
                  p.adjust.method = "holm")

multiple comp

三元以上の場合

おそらく同様にできると思います。 ただ三元配置以上は交互作用の解釈がかなり難しくなるので、可能であれば実験デザインの段階で避けるべきでしょうね(少なくとも私がやっている生物系の実験の場合には)

感想

Rでの多重比較について書いてみました。 以前はSPSSを使っていたのですが、SPSSではGUIだったのでRとはだいぶん感覚が違います。Rの方が手法を勉強する気にはなるので良いと思っております。

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) も出してみています。

感想

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