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

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

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

更新日2020.3.13

使用するパッケージ

install.packages("ggplot2")
install.packages("ggsignif")

データセットを作成します

dat <- c(-0.29733004,-2.63812280,-0.90097072,1.06843016,-3.03846846,2.14097694,2.47494865,-0.02154341,0.56411223,2.81826063,-1.22557353,-1.96770187,-2.41952982,-2.04788520,
-1.49396856,4.42978119,1.14121396,2.81408173,3.88617124,0.27836700,7.04775580,6.76956940,-1.01243943,2.28982618,2.56697055,3.92271476,0.56754099,2.55950816,5.45781910,
5.97792826,3.57321119,3.84212737,1.38192000,1.89198329,5.90437940,3.16262753,1.29383105,7.08585415,1.53787689,5.82484000,4.45300261,3.94823267,2.22314781,5.21597500,
-1.10466187,1.03430878,3.26798971,4.46817078,4.30615587,3.13762549,3.12041534,5.40843117,2.01223446,-2.57111800,4.01156402,6.58882849,6.94067008,3.12482569,3.81695957,
8.53560520)
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.2973300       前     A
2 -2.6381228       前     A
3 -0.9009707       前     A
4  1.0684302       前     A
5 -3.0384685       前     A
6  2.1409769       前     A

列pre_postの配列を確認します

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

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

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

二群比較
パッケージggplot2を使用

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

f:id:yoshida931:20200313164433p:plain:w500

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

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

(g <- ggplot(data = da01,
       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 = "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:20200313164648p:plain:w500

単回帰分析

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

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.4656     0.5305  -0.878    0.388    
pre_post後    3.8146     0.7503   5.084  2.2e-05 ***

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

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.2465     0.6520   4.979 2.94e-05 ***
pre_post後    0.5671     0.9221   0.615    0.544    

単回帰の結果を書き込んでみます
パッケージ ggsignifを使用

library("ggsignif")
(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:20200313165330p:plain:w500

重回帰分析

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

                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.4656     0.5944  -0.783  0.43671    
pre_post後          3.8146     0.8406   4.538 3.05e-05 ***
treatB              3.7121     0.8406   4.416 4.64e-05 ***
pre_post後:treatB  -3.2475     1.1888  -2.732  0.00841 ** 

#分散分析表
anova(fit)

               Df  Sum Sq Mean Sq F value    Pr(>F)    
pre_post        1  71.995  71.995 13.5855 0.0005159 ***
treat           1  65.416  65.416 12.3440 0.0008831 ***
pre_post:treat  1  39.549  39.549  7.4629 0.0084092 ** 
Residuals      56 296.767   5.299