library(tidyr)
d1<-read.csv("sbp.csv")
# group명 변경
library(car)
## Loading required package: carData
d1$group <-recode(d1$group, "1='P'")
d1$group <-recode(d1$group, "2='K1'")
d1$group <-recode(d1$group, "3='K2'")
# wide -> long
d1_long <- gather(d1, time_point, measurement, before_sbp:sbp_8min, factor_key=TRUE)
head(d1_long)
## record_id group time_point measurement
## 1 1 K2 before_sbp 146
## 2 2 K2 before_sbp 132
## 3 3 K1 before_sbp 128
## 4 4 P before_sbp 116
## 5 5 K1 before_sbp 119
## 6 6 P before_sbp 107
# stat으로 변형
library(plyr)
cdata <- ddply(d1_long, c("group", "time_point"), summarise,
N = length(measurement),
mean = mean(measurement,na.rm=T),
sd = sd(measurement,na.rm=T),
se = sd / sqrt(N) )
cdata
## group time_point N mean sd se
## 1 K1 before_sbp 40 126.2250 18.87542 2.984466
## 2 K1 sbp_30s 40 118.6667 19.02538 3.008176
## 3 K1 sbp_2min 40 120.6410 17.88926 2.828540
## 4 K1 sbp_4min 40 122.5128 17.85301 2.822809
## 5 K1 sbp_6min 40 119.5758 17.64942 2.790618
## 6 K1 sbp_8min 40 123.1250 19.01885 3.007144
## 7 K2 before_sbp 40 123.5000 17.05647 2.696864
## 8 K2 sbp_30s 40 114.5500 14.61813 2.311329
## 9 K2 sbp_2min 40 121.9500 15.86796 2.508946
## 10 K2 sbp_4min 40 128.2250 15.96228 2.523857
## 11 K2 sbp_6min 40 130.7222 17.99780 2.845701
## 12 K2 sbp_8min 40 134.0526 17.17063 2.714914
## 13 P before_sbp 40 125.6250 20.44215 3.232187
## 14 P sbp_30s 40 114.7000 18.09689 2.861370
## 15 P sbp_2min 40 104.5000 22.69757 3.588800
## 16 P sbp_4min 40 109.3333 15.83855 2.504295
## 17 P sbp_6min 40 111.5882 16.16790 2.556370
## 18 P sbp_8min 40 106.4583 14.74413 2.331252
# CI까지 가능한 함수
library(Rmisc)
## Loading required package: lattice
d2_long <- summarySE(d1_long, measurevar="measurement", groupvars=c("group", "time_point"), na.rm=TRUE)
d2_long
## group time_point N measurement sd se ci
## 1 K1 before_sbp 40 126.2250 18.87542 2.984466 6.036652
## 2 K1 sbp_30s 39 118.6667 19.02538 3.046498 6.167313
## 3 K1 sbp_2min 39 120.6410 17.88926 2.864574 5.799027
## 4 K1 sbp_4min 39 122.5128 17.85301 2.858770 5.787278
## 5 K1 sbp_6min 33 119.5758 17.64942 3.072369 6.258211
## 6 K1 sbp_8min 16 123.1250 19.01885 4.754713 10.134430
## 7 K2 before_sbp 40 123.5000 17.05647 2.696864 5.454923
## 8 K2 sbp_30s 40 114.5500 14.61813 2.311329 4.675105
## 9 K2 sbp_2min 40 121.9500 15.86796 2.508946 5.074821
## 10 K2 sbp_4min 40 128.2250 15.96228 2.523857 5.104983
## 11 K2 sbp_6min 36 130.7222 17.99780 2.999633 6.089578
## 12 K2 sbp_8min 19 134.0526 17.17063 3.939212 8.275977
## 13 P before_sbp 40 125.6250 20.44215 3.232187 6.537716
## 14 P sbp_30s 40 114.7000 18.09689 2.861370 5.787666
## 15 P sbp_2min 40 104.5000 22.69757 3.588800 7.259034
## 16 P sbp_4min 39 109.3333 15.83855 2.536198 5.134264
## 17 P sbp_6min 34 111.5882 16.16790 2.772772 5.641248
## 18 P sbp_8min 24 106.4583 14.74413 3.009634 6.225902
# plot
library(ggplot2)
pd <- position_dodge(0.2)
ggplot(d2_long, aes(x=time_point, y=measurement, colour=group, group=group)) +
geom_errorbar(aes(ymin=measurement-ci, ymax=measurement+ci), width=.1, position=pd) +
geom_line(position=pd) +
geom_point(aes(shape=group),size =3, position=pd) +
coord_cartesian(xlim = c(1,6), ylim = c(80,150))+
labs(y=expression(paste("Heart rate (beats/min)")),x =expression(paste("Time point"))) +
theme_classic() +
theme(legend.position = c(0.8,0.9)) +
theme(legend.title = element_blank())
