理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

独立している2群を比較するグラフ ggplot2

データの準備

dat <- c(c(rnorm(15)),c(rnorm(15,3,1.2)),c(rnorm(15,10,2)),c(rnorm(15,17,1.5)))
pre_post <- c(rep("前", 30),rep("後", 30))
treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15))
dataf <- data.frame(dat, pre_post, treat)
head(dataf)
          dat pre_post treat
1 -0.27386197       前     A
2 -0.03189246       前     A
3 -0.33039790       前     A
4 -0.82833272       前     A
5  0.33606814       前     A
6  0.51831329       前     A

パッケージのインストールを忘れずに

install.packages("ggplot2")
library(ggplot2)

stripchartで簡単に作図してみます

stripchart(dat~pre_post, data = dataf, 
           vertical=TRUE,
           pch=16,
           xlim=c(0.5,2.5))

f:id:yoshida931:20190815173750p:plain:w400

ggplot2で2群を比較するグラフの編集に挑戦してみます

(g1 <- ggplot(dataf, aes(x = pre_post, y = dat)) + geom_point())

f:id:yoshida931:20190815174001p:plain:w400

x軸の「前」「後」の順を変えてみます
前後を数字(0,1)にした列pre_post2を追加します(前=0、後=1で設定します)

# dplyrパッケージのmutateを使用します
library(dplyr)
dataf2 <- mutate(dataf, pre_post2 = if_else(pre_post == "前" , true = 0, false = 1))
head(dataf2, n=3)

          dat pre_post treat pre_post2
1 -0.27386197       前     A         0
2 -0.03189246       前     A         0
3 -0.33039790       前     A         0

#pre_post2という新しい列が追加されました
#データセット名もdataf2に変更しました

reorder で前後の順を変更します

(g2 <- ggplot(dataf2, aes(x = reorder(pre_post, pre_post2), y = dat)) + geom_point())

f:id:yoshida931:20190815174633p:plain:w400

このグラフに回帰直線を描き加えてみます

fit <- lm(dat~pre_post2 , data = dataf2)
(g2 + geom_abline(slope=fit$coef[2], intercept=fit$coef[1]))   

切片がずれる!!!
f:id:yoshida931:20190815174905p:plain:w400

仕方がないので、曝露を数字「0、1」に戻してみます

(g3 <- ggplot(dataf2, aes(x = pre_post2, y = dat)) +
    geom_point() +
    geom_abline(slope=fit$coef[2], intercept=fit$coef[1])
  )    

f:id:yoshida931:20190815175702p:plain:w400

バックの色や軸を編集

(g3 +
    xlim(-0.5,1.5) +
    theme(panel.background = element_blank(),
          axis.line = element_line(colour="black"))
)

f:id:yoshida931:20190815180055p:plain:w400

軸を文字にしたい場合にはパワーポイントで編集するしかないか・・・?

最後にt検定で確認してみます

t.test(dat ~ pre_post2, data = dataf2)

    Welch Two Sample t-test

data:  dat by pre_post2
t = -13.911, df = 45.549, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.751126  -9.526776
sample estimates:
mean in group 0 mean in group 1 
       1.899348       13.038299 

fit$coef[1]  #非曝露群の平均
(Intercept) 
   1.899348 
fit$coef[2]+fit$coef[1]  #曝露群の平均
pre_post2 
  13.0383