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

統計学備忘録 since2016

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

エラーバー付きのグラフ ggplot2

使用するパッケージ

install.packages("ggplot2")
library(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))
da01 <- data.frame(dat, pre_post, treat);head(dataf)

            dat pre_post treat
1  0.0004020798       前     A
2  0.3236960969       前     A
3  1.0026311892       前     A
4  0.3526106439       前     A
5 -0.1307876208       前     A
6 -1.5476749466       前     A

列pre_postの配列を確認します

levels(da01$pre_post)
[1] "後" "前"

このままではx軸が後→前の配列になるので入れ替えます

da01 <- transform(dataf, pre_post= factor(pre_post, levels = c("前", "後")))
levels(da01$pre_post)
[1] "前" "後"

二群比較

ggplot(da01, aes(x = pre_post, y = dat)) + geom_point()+theme_bw()

f:id:yoshida931:20190825213846p:plain:w400

標準誤差をエラーバーとしたグラフを描きます

#position_jitterdodgeで図の重なり具合を調整しておきます
pos_1 <- position_jitterdodge(
  jitter.width  = 0,#横にずらす
  jitter.height = 0,#縦にずらす
  dodge.width   = 0.5#カテゴリで分ける
)

(g <- ggplot(data = da02,
       aes(
         x      = pre_post,
         y      = dat,
         colour = treat,
       )) +
  geom_jitter(alpha = 0.4, position = pos_1) +  #alpha重なりの濃さ
  stat_summary(aes(x = as.numeric(pre_post) + 0.07), fun.y = "mean", geom = "point", size = 3, position = pos_1)+#平均値のプロット
  stat_summary(aes(x = as.numeric(pre_post) + 0.07),fun.data = "mean_se",#mean_seで標準誤差、mean_cl_normalで95%信頼区間
                  # 標準偏差   fun.ymin = function(x) mean(x) - sd(x),fun.ymax = function(x) mean(x) + sd(x),
               geom = "errorbar",
               width = 0.1,#ひげの横幅
               lwd = 1,#ひげの棒の幅
               #fun.args = list(conf.int = 0.95),#信頼区間のときに使用する
               position = pos_1
  ) +
  theme_bw()) #背景の設定

f:id:yoshida931:20190825195621p:plain:w400

単回帰分析

fitA <- lm(dat ~ pre_post , data =da01[da01$treat=="A",])
summary(fitA)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.2681     0.5873  -0.456    0.652    
pre_post後    4.0898     0.8306   4.924 3.42e-05 ***

fitB <- lm(dat ~ pre_post , data =da01[da01$treat=="B",])
summary(fitB) 

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.0459     0.6852   4.445 0.000126 ***
pre_post後    1.0633     0.9690   1.097 0.281873  

単回帰の結果を書き込んでみます

(g2 <- g + 
geom_signif(stat = "identity",
            data = data.frame(x = c(0.87, 1.12),
                              xend = c(1.87, 2.12),
                              y = c(9, 12),
                              annotation = c("**", "N.S.")),
            aes(x = x, xend = xend, y = y, yend = y, annotation = annotation),
            col = c("red","#20B2AA")) ) +
  ylim (-5,13)

f:id:yoshida931:20190825223325p:plain:w400

重回帰分析

fit <- lm(dat ~ pre_post + treat + pre_post * treat, data =da01)
summary(fit)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.2681     0.6381  -0.420 0.676005    
pre_post後          4.0898     0.9025   4.532 3.12e-05 ***
treatB              3.3140     0.9025   3.672 0.000539 ***
pre_post後:treatB  -3.0265     1.2763  -2.371 0.021189 *  

#分散分析表
anova(fit)

Response: dat
               Df Sum Sq Mean Sq F value    Pr(>F)    
pre_post        1  99.58  99.579 16.3016 0.0001658 ***
treat           1  48.64  48.640  7.9627 0.0065965 ** 
pre_post:treat  1  34.35  34.349  5.6231 0.0211892 *  
Residuals      56 342.08   6.109  

この記事で使用したデータフレーム
コピーしてRにペースト
da01 <- read.table("clipboard", header=T)
これで再現できます

> da01
    X         dat pre_post treat
1   1 -0.29733004       前     A
2   2 -2.63812280       前     A
3   3 -0.90097072       前     A
4   4  1.06843016       前     A
5   5 -3.03846846       前     A
6   6  2.14097694       前     A
7   7  2.47494865       前     A
8   8 -0.02154341       前     A
9   9  0.56411223       前     A
10 10  2.81826063       前     A
11 11 -1.22557353       前     A
12 12 -1.96770187       前     A
13 13 -2.41952982       前     A
14 14 -2.04788520       前     A
15 15 -1.49396856       前     A
16 16  4.42978119       前     B
17 17  1.14121396       前     B
18 18  2.81408173       前     B
19 19  3.88617124       前     B
20 20  0.27836700       前     B
21 21  7.04775580       前     B
22 22  6.76956940       前     B
23 23 -1.01243943       前     B
24 24  2.28982618       前     B
25 25  2.56697055       前     B
26 26  3.92271476       前     B
27 27  0.56754099       前     B
28 28  2.55950816       前     B
29 29  5.45781910       前     B
30 30  5.97792826       前     B
31 31  3.57321119       後     A
32 32  3.84212737       後     A
33 33  1.38192000       後     A
34 34  1.89198329       後     A
35 35  5.90437940       後     A
36 36  3.16262753       後     A
37 37  1.29383105       後     A
38 38  7.08585415       後     A
39 39  1.53787689       後     A
40 40  5.82484000       後     A
41 41  4.45300261       後     A
42 42  3.94823267       後     A
43 43  2.22314781       後     A
44 44  5.21597500       後     A
45 45 -1.10466187       後     A
46 46  1.03430878       後     B
47 47  3.26798971       後     B
48 48  4.46817078       後     B
49 49  4.30615587       後     B
50 50  3.13762549       後     B
51 51  3.12041534       後     B
52 52  5.40843117       後     B
53 53  2.01223446       後     B
54 54 -2.57111800       後     B
55 55  4.01156402       後     B
56 56  6.58882849       後     B
57 57  6.94067008       後     B
58 58  3.12482569       後     B
59 59  3.81695957       後     B
60 60  8.53560520       後     B