# Keith Lohse
# 20022-07-20
library("tidyverse"); library("patchwork")


# Figure 1: Example NHST/MET -----
labels <- factor(c("A", "B", "C", "D", "E"))
means <- c(0.01, 1.5, 2.2, 1, 2.5)
SEs <- c(0.25, 1.3, 1, 0.25, 0.5)


# NHST
2*pt(abs(means/SEs), 1000, lower.tail = FALSE) #abs() to make all ts positive
NHST <- 2*pt(abs(means/SEs), 1000, lower.tail = FALSE) # 2* for two-tailed p

# MET
pt((means-1)/SEs, 1000, lower.tail = FALSE) # -1 for test value of positive 1
MET <- pt((means-1)/SEs, 1000, lower.tail = FALSE)
example_data <- data.frame(labels, means, SEs, NHST, MET)


ggplot(data=example_data, aes(x=labels, y=means))+
  geom_hline(yintercept=c(0,1), lty=2, col=c("black", "grey30"))+
  geom_errorbar(aes(group=labels, ymin=means-1.96*SEs, ymax=means+1.96*SEs),
                col="black", width=0.2)+
  geom_point(aes(fill=factor(NHST<0.05)), 
             col="black", stroke=1.5, size=3, shape=21)+
  scale_x_discrete(name = NULL, limits=rev(levels(labels))) +
  scale_y_continuous(name = "Mean Difference", breaks=c(-1:6)) +
  scale_fill_manual(values = c("black", "white"))+
  #scale_shape_manual(values =c(21,24))+
  labs(fill = "NHST p<0.05")+ 
  theme_classic()+
  coord_flip()+
  theme(axis.text=element_text(size=12, color="black", face="bold"), 
        legend.text=element_text(size=12, color="black"),
        legend.title=element_text(size=12, face="bold"),
        axis.title=element_text(size=12, face="bold"),
        plot.title=element_text(size=12, face="bold", hjust=0.5),
        panel.grid.minor = element_blank(),
        strip.text = element_text(size=12, face="bold"),
        legend.position = "none")

print(NHST)
print(MET)

# Figure 2: Sampling Distributions ----
# Simulation for Independent Samples T-Test

ES<-0.5
s <- 1
n_group <- 64
NoOfTrials <- 5000


# # Example before looping through script ----
# set.seed(17)
# SampleA<-rnorm(n=n_group, mean=0, sd=s)
# SampleB<-rnorm(n=n_group, mean=ES, sd=s) 
# t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE)$estimate[1]
# t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE)$estimate[2]
# mean(SampleA)
# mean(SampleB)
# var(SampleA)
# var(SampleB)
# as.numeric(t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE)$p.value)
# 
# str(t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE))


# Simulating a Data Set ----

index<-c(seq(from =1, to=NoOfTrials, by=1))
sim_data<-data.frame(index, ES, n_group)

head(sim_data)

set.seed(17)

for (i in 1:length(index)){
    SampleA<-rnorm(n=n_group, mean=0, sd=s)
    SampleB<-rnorm(n=n_group, mean=ES, sd=s) 
    
    sim_data$groupA_mean[i] <- mean(SampleA)
    sim_data$groupB_mean[i] <- mean(SampleB) 
    sim_data$groupA_var[i] <- var(SampleA) 
    sim_data$groupB_var[i] <- var(SampleB)
    sim_data$mean_diff[i] <- mean(SampleB) - mean(SampleA)
    sim_data$SE_diff[i] <- sqrt((var(SampleA)/n_group)+(var(SampleB)/n_group))
    sim_data$t_value[i] <- as.numeric(t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE)$statistic)
    sim_data$p_value[i] <- as.numeric(t.test(SampleA, SampleB, var.equal=TRUE, paired=FALSE)$p.value)
    
}

head(sim_data)
summary(sim_data$p_value)
sum(sim_data$p_value<0.05)/5000
hist(sim_data$p_value, breaks = c(seq(from=0, to=1, by=0.05)))

sim_data$p_group <- factor(cut(sim_data$p_value, breaks=c(seq(from=0, to=1, by=0.1))))
sim_data$cohen_d <- with(sim_data, mean_diff/sqrt(((groupA_var+groupB_var)/2))                                        )
sim_data$d_group <- factor(cut(sim_data$cohen_d, breaks=c(seq(from=-2, to=2, by=0.2))))

# Effect Size Plot
z_plot <- ggplot(data=sim_data, aes(x=cohen_d))+
  geom_histogram(aes(fill=p_group), col="black", 
                 breaks=c(seq(from=-2, to=2, by=0.1)))+
  scale_x_continuous(name = "Effect Size (d)") +
  scale_y_continuous(name = "Count") +
  #labs(fill = "Group", col="Group", shape="Estimation", lty="Estimation")+ 
  theme_bw()+
  theme(axis.text=element_text(size=10, color="black"), 
        legend.text=element_text(size=12, color="black"),
        legend.title=element_text(size=12, face="bold"),
        axis.title=element_text(size=12, face="bold"),
        plot.title=element_text(size=12, face="bold", hjust=0.5),
        panel.grid.minor = element_blank(),
        strip.text = element_text(size=12, face="bold"),
        legend.position = "none")

# P-Values Plot
p_plot <- ggplot(data=sim_data, aes(x=p_value))+
  geom_histogram(aes(fill=p_group), col="black", 
                 breaks=c(seq(from=0, to=1, by=0.05)))+
  scale_x_continuous(name = "P-Value") +
  scale_y_continuous(name = "Count") +
  #labs(fill = "Group", col="Group", shape="Estimation", lty="Estimation")+ 
  theme_bw()+
  theme(axis.text=element_text(size=10, color="black"), 
        legend.text=element_text(size=12, color="black"),
        legend.title=element_text(size=12, face="bold"),
        axis.title=element_text(size=12, face="bold"),
        plot.title=element_text(size=12, face="bold", hjust=0.5),
        panel.grid.minor = element_blank(),
        strip.text = element_text(size=12, face="bold"),
        legend.position = "none")


z_plot/p_plot