rm(list = ls())
workdir = "G:/Dropbox/idx/topic/20200930T145416"
setwd(workdir)
df <- read.csv(file = 'data.csv',
header = TRUE, sep = ",", dec = ".")
head(df)
## ID X F1 M Y F2 F3
## 1 1 0 F 124.92947 -1.7235400 up H
## 2 2 0 F 136.35355 -2.1410510 up H
## 3 3 0 F 118.97132 -1.8739778 up H
## 4 4 0 F 118.38304 -0.9486887 up H
## 5 5 0 F 85.96932 -1.7548541 up H
## 6 6 0 F 133.42065 -2.3422194 up H
df[which(df$X==0),"F3"] <- "H"
df[which(df$X==60 | df$X==-60),"F3"] <- "M"
df <- na.omit(df)
library(corrplot)
library(ggpubr)
library(psych)
Boxplots with signif test
lim.y <- c(-4,2)
BoxPlotByTests <- function(df,x,y,f1,f2,m,testmethod){
ggpubr::ggboxplot(df, x = x, y = y,
#palette = "npg",
add = "jitter",
color = m,
fill = x,
outlier.shape=NA,
notch=TRUE) +
facet_grid(df[,f1]~df[,f2])+
scale_colour_gradient(high = "black",low = "grey80")+
geom_text(aes(label=paste("n=",..count..,sep="")),
y=-3.8,
stat='count', size=4) +
ggpubr::stat_compare_means(comparisons = list(c("H", "M")),
aes(label = paste(..method.., ..p.format..)),
paired=FALSE,
# label="{p.signif}{method}",
method = testmethod,
# p.adjust.methods="bonferroni",
label.x = 1.4,
label.y = 1.2
)+
stat_summary(fun = median, geom = "label",
#color=color_small_dots,
size = 4,
hjust = -0.5,
# alpha=0.5,
aes(label=gsub("\\.",".",as.character(round(..y..,1)))))+
#theme(strip.background = element_rect(fill = "white"))+
theme_bw()+
xlab("")+
ggtitle(testmethod)+
ylab("")+
ylim(lim.y)
}
p.wilcox <- BoxPlotByTests(df,"F3","Y","F2","F1","M","wilcox.test")
p.t <- BoxPlotByTests(df,"F3","Y","F2","F1","M","t.test")
ggpubr::ggarrange(p.wilcox,p.t,common.legend = TRUE)
Boxplots with scatterplots
p.scatter<-ggplot(df,aes(M,Y, color = F3))+
geom_point()+
# scale_colour_gradient(high = "black",low = "grey80")+
facet_grid(F2~F1)+
ylab("")+
xlab("")+
stat_smooth(mapping=aes(M,Y,
colour=as.factor(F3)
),
method = "lm",
se = FALSE,
formula = y ~ x) +
ggpmisc::stat_poly_eq(formula = y~x,
mapping=aes(M,Y,
label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "~~~"),
colour=as.factor(F3)),
parse = TRUE) +
# theme(strip.background = element_rect(fill = "white"))+
theme_bw()+
# scale_color_grey(start=0.6, end=0)+
# ggtitle("Wilcoxon Rank Sum and Signed Rank Tests")+
# stat_ellipse()+
xlim(c(30,220))+
ylim(lim.y)
ggpubr::ggarrange(p.wilcox,p.scatter,align="hv",labels=c("A","B"),common.legend = TRUE)
Scatterplot with boxplots (ggscatterhist)
p.scatter1.FU <- ggpubr::ggscatterhist(subset(df,F1=="F" & F2=="up"),
x="M",
y="Y",
color = "F3",
margin.params = list(fill = "F3",
color = "black",
size = 0.2),
margin.plot = "boxplot")
p.scatter1.FU$sp <- p.scatter1.FU$sp +
geom_smooth(mapping=aes(M,Y,colour=as.factor(F3)),
method='lm', se=FALSE, formula = y ~ x)+
ggpmisc::stat_poly_eq(formula = y~x,
mapping=aes(M,Y,
label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "~~~"),
colour=F3),
parse = TRUE)+
stat_ellipse(mapping=aes(M,Y,colour=as.factor(F3)))
p.scatter1.FU
Scatterplot with boxplots (ggMarginal)
MyScatterMarginal <- function(df,t){
p <- ggplot(df,
aes(M,
Y,
color = F3))+
geom_point()+
stat_smooth(formula = y ~ x,
method = "lm",
se = T) +
ggpmisc::stat_poly_eq(formula = y~x,
mapping=aes(M,Y,
label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "~~~"),
colour=F3),
parse = TRUE)+
xlim(c(30,220))+
ylim(c(-4.1,1.1))+
ggtitle(t)+
theme_classic()+
theme(legend.position = c(0.8, 0.8))
p <- ggExtra::ggMarginal(p,
type = "boxplot",
margins = "both",
xparams = list(width=2),
yparams = list(width=2),
groupColour = TRUE,
groupFill = TRUE)
}
p.FU <- MyScatterMarginal(subset(df,F1=="F" & F2=="up"),"FU")
p.CU <- MyScatterMarginal(subset(df,F1=="C" & F2=="up"),"CU")
p.FD <- MyScatterMarginal(subset(df,F1=="F" & F2=="down"),"FD")
p.CD <- MyScatterMarginal(subset(df,F1=="C" & F2=="down"),"CD")
ggpubr::ggarrange(p.FU,p.CU,p.FD,p.CD,ncol=2,nrow=2,
labels=c("A","B","C","D"))
Scatterplot with boxplots and signif test
x.lim <- c(min(df$M)-0.001*min(df$M),max(df$M)+0.001*max(df$M))
y.lim <- c(min(df$Y)-0.001*min(df$Y),max(df$Y)+0.2*max(df$Y))
MyScatterBoxplot <- function(df,t){
pt <- ggplot(df,
aes(x=F3,
y=M,
fill=F3)) +
# geom_violin(alpha=0.2)+
geom_boxplot(notch = TRUE)+
stat_summary(fun = median, geom = "label",
#color=color_small_dots,
size = 4,
hjust = -0.1,
# alpha=0.5,
aes(label=gsub("\\.",".",as.character(round(..y..,1)))))+
# stat_compare_means(paired=FALSE,
# comparisons =list(c("H","M")),
# label="p.signif",
# method = "t.test",
# vjust=10)+
# xlab("")+
ylab("")+
ggtitle(t)+
ylim(x.lim)+
theme(axis.line.x =element_blank(),
axis.line.y.right = element_blank(),
axis.line.y.left = element_line(size = rel(1)),
axis.text.x = element_blank(),
panel.background =element_blank())+
# theme_classic()+
coord_flip()
pr <- ggplot(df,
aes(x=F3,
y=Y,
fill=F3),show.legend = FALSE) +
# geom_violin(alpha=0.2)+
# geom_jitter()+
#geom_histogram()+
geom_hline(yintercept=0, linetype="dashed",
size=0.5)+
geom_boxplot(notch = TRUE)+
# geom_text(aes(label=paste("n=",..count..,sep="")),
# y=-3.8, stat='count', size=4) +
stat_compare_means(paired=FALSE,
label.y = y.lim[2],
comparisons =list(c("H","M")),
label="p.signif",
method = "wilcox.test"
# label.x = 1.25
)+
geom_text(aes(label=paste("n=",..count..,sep="")),
y=-3.7,
stat='count', size=4) +
stat_summary(fun = median, geom = "label",
#color=color_small_dots,
size = 4,
hjust = 0.3,
# alpha=0.5,
aes(label=gsub("\\.",".",as.character(round(..y..,1)))))+
ylab("")+
ylim(y.lim)+
theme(axis.line.y =element_blank(),
axis.line.x.top = element_blank(),
axis.line.x.bottom = element_line(size = rel(1)),
axis.text.y = element_blank(),
panel.background =element_blank())
ps <- ggplot(df,
aes(x=M,
y=Y,
color=F3)) +
geom_point(size = rel(1.5)) +
stat_smooth(formula = y ~ x,
method = "lm",
se = T) +
# ggpmisc::stat_poly_eq(formula = y~x,
# mapping=aes(M,Y,
# label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "~~~"),
# colour=F3),
# parse = TRUE)+
stat_cor(p.digits = 3, p.accuracy = 0.001)+
ylim(y.lim)+
geom_hline(yintercept=0, linetype="dashed",
size=0.5)+
xlim(x.lim)+
ylab(t)+
theme_pubr()
# theme(axis.line=element_line(),
# axis.line.y.right=element_line())
blank <- ggplot() +
geom_point(aes(1,1),
colour="white") +
theme(axis.ticks=element_blank(),
panel.background=element_blank(),
panel.grid=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
ggarrange(pt, blank, ps, pr, ncol=2, nrow=2,
common.legend = TRUE,
legend="none",
align = "hv",
widths=c(4, 1),
heights=c(1, 4))
# gridExtra::grid.arrange(pt, blank, ps, pr, ncol=2, nrow=2,
# widths=c(4, 1), heights=c(1, 4))
}
p.FU<-MyScatterBoxplot(subset(df,F1=="F" & F2=="up"),"FU")
p.FU
p.CU<-MyScatterBoxplot(subset(df,F1=="C" & F2=="up"),"CU")
p.FD<-MyScatterBoxplot(subset(df,F1=="F" & F2=="down"),"FD")
p.CD<-MyScatterBoxplot(subset(df,F1=="C" & F2=="down"),"CD")
windows(width=15,height=15)
# ggarrange(p.FU,p.CU,p.FD,p.CD,common.legend = TRUE,align="hv"
# labels=c("A","B","C","D")
# )
Correlation (lower)/Partial correlation (upper) Matrix
CorPcorPlot <- function(xtitle,x){
# graphics.off()
# windows()
x.cor <- corr.test(x)
x.pcor <- corr.p(partial.r(x),n=x.cor$n)
x.lowup <- lowerUpper(lower=x.cor$r,upper=x.pcor$r,diff=FALSE)
x.lowup_pvalue <- lowerUpper(lower=x.cor$p,upper=x.pcor$p,diff=FALSE)
corrplot::corrplot(x.lowup,
method="ellipse",
# add = FALSE,
diag=FALSE,
# title=xtitle,
# label=TRUE,
p.mat = x.lowup_pvalue,
# sig.level = c(0.001,0.01,0.05),
# sig.level = 0.01,
# insig="label_sig",
insig="blank",
# tl.cex = 1,
addCoef.col = "black",
# addCoefasPercent = TRUE,
outline=FALSE,
pch.cex = 1.5,
#cl.length = 21,
# pch= "*"
)
# savePlot(filename = "1.png",type="png")
# windows()
}
p.Corr.FU <- CorPcorPlot("Cor/PartialCor-Plot",subset(df,F1=="F" & F2=="up",select=c(2,4,5)))
ggpairs matrix
dff <- df
dff$X <- as.factor(dff$X)
GGally::ggpairs(subset(dff,F1 == "F" & F2=="up"),columns=c("Y","M","X"),aes(color=F3))
Mediation Analysis
myMediation <- function(var_names, var_values){
data <- c(0, var_values[2], 0,
0, 0, 0,
var_values[3], var_values[1], 0)
M <- matrix (data=data,
nrow=3, ncol=3, byrow = TRUE)
p <- diagram::plotmat(M,
pos=c(1,2),
name= var_names,
box.type = "rect",
box.size = 0.12,
arr.width = 0.1,
arr.type = "triangle",
arr.length = 0.2,
box.prop=0.5,
curve=0)
}
# label signif stars for mediation diagram
p2star <- function(x){
if(x > 0.05){
print(" ")
} else if (x > 0.01) {
print("*")
} else if (x > 0.001) {
print("**")
} else if (x <= 0.001) {
print("***")
}
}
myMedCalc <- function(x,y,m,xymLabels){
df <- c()
df$x <- x
df$y <- y
df$m <- m
df <- as.data.frame(df)
model.c <- lm(y ~ x, data = df)
model.a <- lm(m ~ x, data = df)
model.b <- lm(y ~ x + m, data = df)
# standardized
model.cs <- lm(scale(y) ~ scale(x), data = df)
model.as <- lm(scale(m) ~ scale(x), data = df)
model.bs <- lm(scale(y) ~ scale(x) + scale(m), data = df)
# number of digits shown on mediation diagram
rd <- 2
# labels for mediation diagram
b <- c()
b[1]<- paste("'",as.character(round(model.cs$coefficients[2],digits=rd)),
p2star(summary(model.cs)$coefficients[,4][2]),
" (",as.character(round(model.bs$coefficients[2],digits=rd)),
p2star(summary(model.bs)$coefficients[,4][2]),
")'",sep="")
b[2]<- paste("'",as.character(round(model.as$coefficients[2],digits=rd)),
p2star(summary(model.as)$coefficients[,4][2]),"'",
sep="")
b[3]<- paste("'",as.character(round(model.bs$coefficients[3],digits=rd)),
p2star(summary(model.bs)$coefficients[,4][3]),"'",
sep="")
# mediation diagram with signif stars
myMediation(c(xymLabels[3],xymLabels[1],xymLabels[2]),b)
# compare my myMediation with psych::mediate
psych::mediate(y ~ x + (m),
data=df,
n.iter=500,
main = "std=F",
std=FALSE,
zero=FALSE)
psych::mediate(y ~ x + (m),
data=df,
main = "std=T",
n.iter=500,
std=TRUE,
zero=FALSE)
# bootstrapping
rb <- mediation::mediate(model.a,
model.b,
treat="x",
mediator="m",
boot=TRUE,
sims=500)
# summary(rb)
plot(rb)
# print(m1,short=FALSE)
}
myMedCalc(df$X,df$Y,df$M,
c("X","Y","Mediator"))
## [1] "*"
## [1] "***"
## [1] "*"
## [1] "***"
