-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbeta_paired_plots.R
127 lines (88 loc) · 4.59 KB
/
beta_paired_plots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
####PLOT ALL BETA
pdf("beta.rep.diff.all.pdf", width=20)
par(mfrow=c(1,6))
for(i in 1:6) {
pre<-top6.beta[i,rep1]
pre.case<-top6.beta[i,rep1 & phe$grp2=="case"]
pre.ctrl<-top6.beta[i,rep1 & phe$grp2=="ctrl"]
post<-top6.beta[i,rep2]
post.case<-top6.beta[i,rep2 & phe$grp2=="case"]
post.ctrl<-top6.beta[i,rep2 & phe$grp2=="ctrl"]
s<-seq(length(pre))
par(bty="l")
boxplot(pre,post,main=paste(top6[i], "(",top6.names[i],")"),ylab="Beta",names=c("baseline","follow-up"))
stripchart(list(pre,post),vertical=T,pch=16,cex=0.5,add=T)
segments(rep(1,length(pre.case))[s],pre.case[s],rep(2,length(pre.case))[s],post.case[s],col="red",lwd=0.5)
segments(rep(1,length(pre.ctrl))[s],pre.ctrl[s],rep(2,length(pre.ctrl))[s],post.ctrl[s],col="blue",lwd=0.5)
}
dev.off()
####PLOT ALL BETA VERSION 2
pdf("beta.rep.diff.all.2.pdf", width=14)
par(mfrow=c(1,6))
for(i in 1:6) {
pre<-top6.beta[i,rep1]
pre.case<-top6.beta[i,rep1 & phe$grp2=="case"]
pre.ctrl<-top6.beta[i,rep1 & phe$grp2=="ctrl"]
post<-top6.beta[i,rep2]
post.case<-top6.beta[i,rep2 & phe$grp2=="case"]
post.ctrl<-top6.beta[i,rep2 & phe$grp2=="ctrl"]
res.case<-wilcox.test(post.case,pre.case,paired=T,conf.int=T)
res.ctrl<-wilcox.test(post.ctrl,pre.ctrl,paired=T,conf.int=T)
stripchart(list(post.case-pre.case,post.ctrl-pre.ctrl),vertical=T,pch=16,method="jitter",main=paste(top6[i], "(",top6.names[i],")"),ylab="Follow-up - Baseline",xlab="Median+/-95%CI", col=c("red","blue"))
points(1,res.case$estimate,col="red",pch=16,cex=2)
arrows(1,res.case$conf.int[1],1,res.case$conf.int[2],col="red",code=3,lwd=3,angle=90)
points(2,res.ctrl$estimate,col="blue",pch=16,cex=2)
arrows(2,res.ctrl$conf.int[1],2,res.ctrl$conf.int[2],col="blue",code=3,lwd=3,angle=90)
abline(h=0,lty=2)#Zero-effectline
}
dev.off()
####PLOT ALL BETA VERSION 3
pdf("beta.rep.diff.all.3.pdf", width=21)
par(mfrow=c(1,6))
for(i in 1:6) {
pre<-top6.beta[i,rep1]
pre.case<-top6.beta[i,rep1 & phe$grp2=="case"]
pre.ctrl<-top6.beta[i,rep1 & phe$grp2=="ctrl"]
post<-top6.beta[i,rep2]
post.case<-top6.beta[i,rep2 & phe$grp2=="case"]
post.ctrl<-top6.beta[i,rep2 & phe$grp2=="ctrl"]
s<-seq(length(pre))
par(bty="l")
boxplot(pre.case,pre.ctrl,post.case,post.ctrl,xlim = c(0, 8), at=c(1,2,7,8),main=paste(top6[i], "(",top6.names[i],")"),ylab="Beta",names=c("baseline","","follow-up", ""), col=c("palevioletred","lightblue", "palevioletred","lightblue"))
stripchart(list(pre,post),vertical=T,pch=16,cex=0.5,add=T)
segments(rep(1,length(pre.case))[s],pre.case[s],rep(7,length(pre.case))[s],post.case[s],col="red",lwd=0.5)
segments(rep(2,length(pre.ctrl))[s],pre.ctrl[s],rep(8,length(pre.ctrl))[s],post.ctrl[s],col="blue",lwd=0.5)
}
dev.off()
#### PLOT ALL + DIFFERENCES
pdf("beta.rep.diff.all.1.pdf", width=21)
#PLOT
pre<-top6.beta[i,rep1]
pre.case<-top6.beta[i,rep1 & phe$grp2=="case"]
pre.ctrl<-top6.beta[i,rep1 & phe$grp2=="ctrl"]
post<-top6.beta[i,rep2]
post.case<-top6.beta[i,rep2 & phe$grp2=="case"]
post.ctrl<-top6.beta[i,rep2 & phe$grp2=="ctrl"]
#Settinguptwoscreens
par(mfrow=c(1,2))
#FirstGraph
s<-seq(length(pre))
par(bty="l")
boxplot(pre,post,main=paste("Res at",top6[i], "(",top6.names[i],")"),xlab=paste("pre/post p =",signif(p.time[i,3],3),", ","case/control p =",signif(p.time[i,2],3),"\n","interaction p =",signif(p.time[i,1],3)),ylab="Res",names=c("pre","post"),col=c("lightblue","lightgreen"))
stripchart(list(pre,post),vertical=T,pch=16,method="jitter",cex=0.5,add=T)
segments(rep(0.95,length(pre.case))[s],pre.case[s],rep(2,length(pre.case))[s],post.case[s],col="red",lwd=0.5)
segments(rep(0.95,length(pre.ctrl))[s],pre.ctrl[s],rep(2,length(pre.ctrl))[s],post.ctrl[s],col="blue",lwd=0.5)
#Secondgraph
#Confidenceintervals eitherparametric (t.test) or non-parametric (wilcox.text)
#res<-t.test(post,prä,paired=T,conf.int=T)
res.case<-wilcox.test(post.case,pre.case,paired=T,conf.int=T)
res.ctrl<-wilcox.test(post.ctrl,pre.ctrl,paired=T,conf.int=T)
stripchart(post-pre,vertical=T,pch=16,method="jitter",main="Difference in Res",ylab="Difference:Post–Pre",xlab="Median+/-95%CI", , col="white")
stripchart(post.case-pre.case,vertical=T,pch=16,method="jitter",main="Difference in Beta",ylab="Difference:Post–Pre",xlab="Median+/-95%CI", , col="red",add=T)
stripchart(post.ctrl-pre.ctrl,vertical=T,pch=16,method="jitter",main="Difference in Beta",ylab="Difference:Post–Pre",xlab="Median+/-95%CI", add=T, col="blue")
points(1,res.case$estimate,col="red",pch=16,cex=2)
arrows(1,res.case$conf.int[1],1,res.case$conf.int[2],col="red",code=3,lwd=3,angle=90)
points(1,res.ctrl$estimate,col="blue",pch=16,cex=2)
arrows(1,res.ctrl$conf.int[1],1,res.ctrl$conf.int[2],col="blue",code=3,lwd=3,angle=90)
abline(h=0,lty=2)#Zero-effectline
dev.off()