library(car)
library(BayesFactor)
main_effect_finder <- function(df,VAL = NULL,F1 = NULL, F2 = NULL, p_val, BF){
if(is.null(VAL) == TRUE | (is.null(F1) == TRUE & is.null(F2) == TRUE)){
return(NULL)
}
#Computing the results
if(is.null(F1) == FALSE & is.null(F2) == FALSE){
working.data <- df[c(VAL,F1,F2)]
working.data <- as.data.frame(working.data)
names(working.data) <- c("VAL","F1","F2")
ANOVA_P <- Anova(lm(as.formula("VAL ~ F1*F2"), data=working.data, contrasts=list(F1=contr.poly, F2=contr.poly)), type=3)
working.data[,2] <- as.factor(working.data[[2]])
working.data[,3] <- as.factor(working.data[[3]])
ANOVA_bf = anovaBF(as.formula("VAL ~ F1*F2"), data=working.data)
#print(ANOVA_bf)
return.data <- as.data.frame(matrix(rep(NA, times = (length(c(F1,F2))+1)*5),nrow = (length(c(F1,F2))+1),ncol = 5))
}
if(is.null(F1) == FALSE & is.null(F2) == TRUE){
working.data <- df[c(VAL,F1)]
names(working.data) <- c("VAL","F1")
working.data <- as.data.frame(working.data)
ANOVA_P <- Anova(lm(as.formula("VAL ~ F1"), data=working.data, contrasts=list(F1=contr.poly)), type=3)
working.data[,2] <- as.factor(working.data[,2])
ANOVA_bf = anovaBF(as.formula("VAL ~ F1"), data=working.data)
#print(ANOVA_bf)
return.data <- as.data.frame(matrix(rep(NA, times = (length(c(F1,F2)))*5),nrow = (length(c(F1,F2))),ncol = 5))
}
if(is.null(F1) == TRUE & is.null(F2) == FALSE){
working.data <- df[c(VAL,F2)]
names(working.data) <- c("VAL","F2")
working.data <- as.data.frame(working.data)
ANOVA_P <- Anova(lm(as.formula("VAL ~ F2"), data=working.data, contrasts=list(F2=contr.poly)), type=3)
working.data[,2] <- as.factor(working.data[,2])
ANOVA_bf = anovaBF(as.formula("VAL ~ F2"), data=working.data)
#print(ANOVA_bf)
return.data <- as.data.frame(matrix(rep(NA, times = (length(c(F1,F2)))*5),nrow = (length(c(F1,F2))),ncol = 5))
}
#Creating return DF
names(return.data) <- c("Factor","P-value","ANOVA - parametric vs crit", "BF","Bayes factor vs crit")
if(is.null(F1) == FALSE & is.null(F2) == FALSE){
return.data[,1] <- c(F1,F2,paste(F1,":",F2))
}
if(is.null(F1) == FALSE & is.null(F2) == TRUE){
return.data[,1] <- c(F1)
}
if(is.null(F1) == TRUE & is.null(F2) == FALSE){
return.data[,1] <- c(F2)
}
#FILLING IN parametric analysis results
for(i in c(1:dim(return.data)[1])){
return.data[i,2] <- ANOVA_P$"Pr(>F)"[i+1]
if(i <= 2){
return.data[i,4] <- exp(1)**ANOVA_bf@bayesFactor[,"bf"][i]
}
if(i > 2){
return.data[i,4] <- exp(1)**ANOVA_bf@bayesFactor[,"bf"][1+i]
}
if(ANOVA_P$"Pr(>F)"[i+1] < p_val){
return.data[i,3] <- "significant"
}
if(ANOVA_P$"Pr(>F)"[i+1] >= p_val){
return.data[i,3] <- "N.S"
}
if(i <= 2){
if(exp(1)**ANOVA_bf@bayesFactor[,"bf"][i] >= BF){
return.data[i,5] <- "significant"
}
if(exp(1)**ANOVA_bf@bayesFactor[,"bf"][i] < BF){
return.data[i,5] <- "N.S"
}
}
if(i > 2){
if(exp(1)**ANOVA_bf@bayesFactor[,"bf"][i+1] >= BF){
return.data[i,5] <- "significant"
}
if(exp(1)**ANOVA_bf@bayesFactor[,"bf"][i+1] < BF){
return.data[i,5] <- "N.S"
}
}
}
return(return.data)
}
library(shiny)
library(shinythemes)
library(dplyr)
library(markdown)
library(ggplot2)
library(ggthemes)
library(extrafont)
library(scales)
library(reshape2)
library(stringr)
library(stargazer)
library(Cairo)
library(datasets)
library(DT)
library(car)
library(agricolae)
library(BayesFactor)
library(cowplot)
library(BEST)
source("worker_functionsAppHDI.R")
source("main_effects_computer.R")
#HELPERS FOR PREDICTOR AND OUTCOM SLECTION
Predictor_col_finder <- function(df,down.crit = 2, up.crit = 8){
col_names <- names(df)
retur_results <- c()
for(i in col_names){
if(class(df[,i]) == "character" & length(unique(df[,i])) >= down.crit & length(unique(df[,i])) <= up.crit){
retur_results[length(retur_results)+1] <- i
}
}
return(retur_results)
}
Outcom_col_finder <- function(df){
col_names <- names(df)
retur_results <- c()
for(i in col_names){
if((class(df[,i]) == "numeric" | class(df[,i]) == "integer" )& length(df[,i]) > sum(is.na(df[,i]))){
retur_results[length(retur_results)+1] <- i
}
}
return(retur_results)
}
#Pre-provided data factors
provided_data_factors <- list()
provided_data_outcomes <- list()
files <- dir(paste(getwd(),"/Data",sep = ""))
if(length(files) > 0){
for(i in files){
temp_data <- read.csv(paste(getwd(),"/Data/",i,sep = ""),header = TRUE,stringsAsFactors = FALSE)
provided_data_factors[[i]] <- Predictor_col_finder(temp_data)
provided_data_outcomes[[i]] <- Outcom_col_finder(temp_data)
}
}
#Shiny server function
shinyServer(function(input,output,session){
#DEFAULT UI
output$data_selector <- renderUI({selectizeInput("data_selector_in","Select dataset to use:",
choices = dir(paste(getwd(),"/Data",sep = "")),
multiple = FALSE,
options = list(maxItems = 1),
selected = dir(paste(getwd(),"/Data",sep = ""))[1])})
output$post_hoc_selector <- renderUI({selectizeInput("post_hoc_selector_in","Select parametric post hoc to use:",
choices = c("LSD.test", "HSD.test", "duncan.test"),
multiple = FALSE,
options = list(maxItems = 1),
selected = "LSD")})
output$FAC_selector <- renderUI({selectizeInput("FAC_selector_in","Select independent variables:",
choices = provided_data_factors[as.character(unlist(input$data_selector_in))],
multiple = TRUE,
options = list(maxItems = 2),
selected = NULL)})
output$VAL_selector <- renderUI({selectizeInput("VAL_selector_in","Select dependent variable:",
choices = provided_data_outcomes[as.character(unlist(input$data_selector_in))],
multiple = FALSE,
options = list(maxItems = 1),
selected = NULL)})
#DATA UPLOAD
input_file <- reactive({
input_file <- input$OwnData
if (is.null(input_file)) {
return(NULL)
}
read.csv(input_file$datapath,header = TRUE,stringsAsFactors = FALSE)
})
user.data <- observeEvent(input$upload,{
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait! Adding data.", value = 0)
if(is.null(input$OwnData) == TRUE){
return(NULL)
}
if(is.null(input$OwnData) == FALSE){
file_name <- paste(gsub(x=input$user.dataset.name,pattern = ".",replacement = "",fixed = TRUE),".csv",sep = "")
path <- paste(getwd(),"/Data/",file_name,sep = "")
write.csv(input_file(),path,row.names = FALSE)
}
})
user.data.factors <- eventReactive(input$upload,{})
user.data.factors <- eventReactive(input$upload,{})
#Updating UI elements
observeEvent(input$upload,{
updateSelectizeInput(session, "data_selector_in", label = "Select dataset to use:", choices = dir(paste(getwd(),"/Data",sep = "")),
selected = dir(paste(getwd(),"/Data",sep = ""))[1], options = list(maxItems = 1))
})
observeEvent(input$upload2,{
provided_data_factors <- list()
files <- dir(paste(getwd(),"/Data",sep = ""))
if(length(files) > 0){
for(i in files){
temp_data <- read.csv(paste(getwd(),"/Data/",i,sep = ""),header = TRUE,stringsAsFactors = FALSE)
provided_data_factors[[i]] <- Predictor_col_finder(temp_data)
}
}
factors <- provided_data_factors[[as.character(unlist(input$data_selector_in))]]
provided_data_outcomes <- list()
files <- dir(paste(getwd(),"/Data",sep = ""))
if(length(files) > 0){
for(i in files){
temp_data <- read.csv(paste(getwd(),"/Data/",i,sep = ""),header = TRUE,stringsAsFactors = FALSE)
provided_data_outcomes[[i]] <- Outcom_col_finder(temp_data)
}
}
vals <- provided_data_outcomes[[as.character(unlist(input$data_selector_in))]]
updateSelectizeInput(session, "VAL_selector_in", label = "Select dependent variable:", choices = vals,
selected = "ext", options = list(maxItems = 1))
updateSelectizeInput(session, "FAC_selector_in", label = "Select independent variables:", choices = factors,
selected = "Study", options = list(maxItems = 2))
})
dfc <- eventReactive(input$go,{
file_name <- input$data_selector_in
data <- read.csv(paste(getwd(),"/Data/",file_name,sep = ""),header = TRUE,stringsAsFactors = FALSE)
return(data)
})
f1c <- eventReactive(input$go,{
if(length(input$FAC_selector_in)==1){
return(as.character(input$FAC_selector_in[1]))
}else if(length(input$FAC_selector_in)==2){
return(as.character(input$FAC_selector_in[1]))
}else{
return(NULL)
}
})
f2c <- eventReactive(input$go,{
if(length(input$FAC_selector_in)==2){
return(as.character(input$FAC_selector_in[2]))
}
else if(length(input$FAC_selector_in) == 1){
return(NULL)
}
})
varc <- eventReactive(input$go,{return(as.character(input$VAL_selector_in))})
output$tbl = DT::renderDataTable(dfc(),
server = TRUE,
options = list(dom = 'C<"clear">lfrtip'))
output$debug1 <- renderText({dim(dfc())})
output$debug2 <- renderText({(f1c())})
output$debug3 <- renderText({(is.null(f2c()))})
output$debug4 <- renderText({(varc())})
aaa <- eventReactive(input$go,{
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
ac <- interaction.maker(df = dfc() ,factor1 = f1c(), factor2 = f2c(), value = varc())
aac <- post_hoc(ac, post.hoc.type = input$post_hoc_selector_in)
aaa <- post_hoc_plot(post.hoc.object = aac,interaction.object = ac, p.val.criteria = input$pvalue,BF.criteria = input$BF,y.axis.length.factor = input$height)
return(aaa)
})
meffect <- eventReactive(input$go,{
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
if(is.null(f1c()) == FALSE & is.null(varc()) == FALSE){
meffect <- main_effect_finder(df = dfc(),VAL = varc(),F1 = f1c(), F2 = f2c(), p_val = input$pvalue, BF = input$BF)
return(meffect)
}
else(return(NULL))
})
output$MEout <- renderUI({
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
if(is.null(meffect()) == FALSE){(HTML(stargazer(meffect(), type = "html", summary = FALSE,
title = "The main effects ",
notes="The singificance and non-significance (N.S) is based on the specified p-value and Bayes Factor.",
nobs = FALSE, font.size = "Huge",column.sep.width = "10pt",min.max= FALSE)))}
else(NULL)
})
output$F1resultsout <- renderUI({
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
if(is.null(aaa()) == FALSE){
(HTML(stargazer(aaa()@plot.comp.f1[,c(1,2,5)],title = "Comparisons between the levels of factor 1 ",
notes="Comparisons based on specified p-value and Bayes Factor", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
}
else(NULL)
})
output$F2resultsout <- renderUI({
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
if(is.null(aaa()) == FALSE){
if(length(aaa()@plot.comp.f2) != 0){
(HTML(stargazer(aaa()@plot.comp.f2[,c(1,2,5)],title = "Comparisons between the levels of factor 2 ",
notes="Comparisons based on specified p-value and Bayes Factor", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
}else(NULL)
}
})
output$INTresults <- renderUI({
progress <- shiny::Progress$new()
on.exit(progress$close())
progress$set(message = "Please wait, computing the results.", value = 0)
if(is.null(aaa()) == FALSE){
if(length(aaa()@plot.comp.inter)!=0){
(HTML(stargazer(aaa()@plot.comp.inter[,c(1,2,5)],title = "Comparisons between the interactions of factor 1 and factor 2 levels ",
notes="Comparisons based on specified p-value and Bayes Factor", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
}else(NULL)
}
})
output$F1plotout <- renderPlot(
if(is.null(aaa()) == FALSE){
{aaa()@plot.f1}
}else(NULL)
)
output$F2plotout <- renderPlot(
if(is.null(aaa()) == FALSE){
{aaa()@plot.f2}
}else(NULL)
)
output$Interplotout <- renderPlot(
if(is.null(aaa()) == FALSE){
{aaa()@plot.inter}
}else(NULL)
)
})
library(shiny)
library(shinythemes)
shinyUI(fluidPage(theme = shinytheme("flatly"),
titlePanel("Comparison of Frequentist and Bayesian analysis in ANOVA post-hoc settings"),
includeScript("../../../Matomo-tquant.js"),
sidebarLayout(
sidebarPanel("Please specify analysis settings",width = 3,
uiOutput("data_selector"),
p(""),
uiOutput("post_hoc_selector"),
p(""),
sliderInput("pvalue", "P-value:",
min = 0.01, max = 1, value = 0.05
),
p(""),
sliderInput("BF", "BF in absolute terms:",
min = 0, max = 100, value = 10
),
p(""),
actionButton("upload2","After data upload - update variable selection"),
p(""),
uiOutput("VAL_selector"),
p(""),
uiOutput("FAC_selector"),
p(""),
textInput("user.dataset.name", "Give name for dataset",value = "UserDataset"),
p(""),
fileInput("OwnData", "Upload your own data...",
accept=c("text/csv", "text/comma-separated-values,text/plain")),
p(""),
actionButton("upload","Add uploaded dataset to selection"),
p(""),
actionButton("go","Compute the results"),
sliderInput("height", "Adjust plot height:",
min = 1.1, max = 2, value = 1.45
),
p("")
),
mainPanel(
tags$style(type="text/css",
".shiny-output-error { visibility: hidden; }",
".shiny-output-error:before { visibility: hidden; }"),
tabsetPanel(
tabPanel(title = "Overview",textOutput("test"),
fluidRow(includeMarkdown("documentation.md"))),
tabPanel("Explore the data",
DT::dataTableOutput('tbl')),
tabPanel("Main effects",
uiOutput("MEout")),
tabPanel("Post hocs in table",
uiOutput("F1resultsout"),
uiOutput("F2resultsout"),
uiOutput("INTresults")),
tabPanel("Post hocs plotted on a graph",
plotOutput("F1plotout"),
plotOutput("F2plotout"),
plotOutput("Interplotout")
),
# tabPanel("DEBUG PRINT PANEL",
# textOutput("debug1"),
# textOutput("debug2"),
# textOutput("debug3"),
# textOutput("debug4"),
# textOutput("debug5"),
# ("debug6")
# ),
tabPanel(title = "About",textOutput("test2"),
fluidRow(includeMarkdown("About.md")))
)
))))
#workerk_function without HDI implementation
library(ggplot2)
library(car)
library(agricolae)
library(BayesFactor)
library(cowplot)
library(ggthemes)
library(stringr)
#Helper functions
s.err <- function(x) sd(x,na.rm = TRUE)/sqrt(sum(!is.na(x)))
RobustMax <- function(x) {if (sum(!is.na(x))>0) max(x,na.rm = TRUE) else NA}
RobustMin <- function(x) {if (sum(!is.na(x))>0) min(x,na.rm = TRUE) else NA}
setClass("interaction",
representation(df = "data.frame",
f1.label = "character",
f2.label = "character",
value.label = "character",
f.levels1 = "character",
f.levels2 = "character",
f.levels.int = "character",
f.levels1.null = "logical",
f.levels2.null = "logical",
f1.averages.se = "data.frame",
f2.averages.se = "data.frame",
inter.averages.se = "data.frame"
))
interaction.maker <- function(df,factor1=NULL,factor2=NULL,value){
return.results <- new("interaction")
df[,"interaction"] <- rep(NA, times = dim(df)[1])
if(is.null(factor1) == TRUE & is.null(factor2) == TRUE){
return.results@f.levels1.null <- TRUE
return.results@f.levels2.null <- TRUE
return(return.results)
}
if(is.null(factor1) == FALSE & is.null(factor2) == TRUE){
return.results@f1.label <- factor1
return.results@value.label <- value
df[,factor1] <- as.character(df[,factor1])
return.results@df <- df
df[,factor1] <- as.factor(df[,factor1])
return.results@f.levels1 <- levels(df[,factor1])
return.results@f.levels1.null <- FALSE
return.results@f.levels2.null <- TRUE
#Averages and SE computations
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels1)*3),
nrow = length(return.results@f.levels1),ncol = 3))
results[,1] <- return.results@f.levels1
names(results) <- c("Factor", "Mean", "SE")
for (i in results[,"Factor"]){
subset <- df[df[,factor1] == i,]
results[results[,"Factor"] == i,c("Mean", "SE")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]))
}
return.results@f1.averages.se <- results
#Returning results
return(return.results)
}
if(is.null(factor1) == TRUE & is.null(factor2) == FALSE){
return.results@f2.label <- factor2
return.results@value.label <- value
df[,factor2] <- as.character(df[,factor2])
return.results@df <- df
df[,factor2] <- as.factor(df[,factor2])
return.results@f.levels2 <- levels(df[,factor2])
return.results@f.levels1.null <- TRUE
return.results@f.levels2.null <- FALSE
#Averages and SE computations
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels2)*3),
nrow = length(return.results@f.levels2),ncol = 3))
results[,1] <- return.results@f.levels2
names(results) <- c("Factor", "Mean", "SE")
for (i in results[,"Factor"]){
subset <- df[df[,factor2] == i,]
results[results[,"Factor"] == i,c("Mean", "SE")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]))
}
return.results@f2.averages.se <- results
return(return.results)
}
if(is.null(factor1) == FALSE & is.null(factor2) == FALSE){
return.results@f.levels1.null <- FALSE
return.results@f.levels2.null <- FALSE
return.results@f1.label <- factor1
return.results@f2.label <- factor2
return.results@value.label <- value
#Making interaction term
factor1_vec <- unique(df[,factor1])
factor2_vec <- unique(df[,factor2])
combined_vec <- c()
n = 0
for (i in c(1:length(factor1_vec))){
for (m in c(1:length(factor2_vec))){
combined_vec[n+1] <- paste(factor1_vec[i],"###",factor2_vec[m], sep = "")
n <- n+1
}
}
#Adding interaction term to df
for(i in c(1:dim(df)[1])){
df[i,"interaction"] <- paste(df[i,factor1],"###",df[i,factor2], sep = "")
}
df[,"interaction"] <- factor(df[,"interaction"], levels = combined_vec)
df[,"interaction"] <- as.character(df[,"interaction"])
return.results@df <- df
df[,factor1] <- as.factor(df[,factor1])
df[,factor2] <- as.factor(df[,factor2])
return.results@f.levels1 <- levels(df[,factor1])
return.results@f.levels2 <- levels(df[,factor2])
return.results@f.levels.int <- combined_vec
#Computing the averages
#First factor
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels1)*3),
nrow = length(return.results@f.levels1),ncol = 3))
results[,1] <- return.results@f.levels1
names(results) <- c("Factor", "Mean", "SE")
for (i in results[,"Factor"]){
subset <- df[df[,factor1] == i,]
results[results[,"Factor"] == i,c("Mean", "SE")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]))
}
return.results@f1.averages.se <- results
#Seccond factor
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels2)*3),
nrow = length(return.results@f.levels2),ncol = 3))
results[,1] <- return.results@f.levels2
names(results) <- c("Factor", "Mean", "SE")
for (i in results[,"Factor"]){
subset <- df[df[,factor2] == i,]
results[results[,"Factor"] == i,c("Mean", "SE")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]))
}
return.results@f2.averages.se <- results
#Interaction
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels.int)*3),
nrow = length(return.results@f.levels.int),ncol = 3))
results[,1] <- return.results@f.levels.int
names(results) <- c("Factor", "Mean", "SE")
for (i in results[,"Factor"]){
subset <- df[df[,"interaction"] == i,]
results[results[,"Factor"] == i,c("Mean", "SE")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]))
}
return.results@inter.averages.se <- results
return(return.results)
}
}
setClass("post.hoc",
representation(df = "data.frame",
f1.comps = "data.frame",
f2.comps = "data.frame",
interaction.comps = "data.frame",
f.levels1.null = "logical",
f.levels2.null = "logical",
post.hoc.type = "character"
))
post_hoc <- function(interaction.object,
post.hoc.type = "LSD.test"
){
return.results <- new("post.hoc")
return.results@df <- interaction.object@df
return.results@post.hoc.type <- post.hoc.type
return.results@f.levels1.null <- interaction.object@f.levels1.null
return.results@f.levels2.null <- interaction.object@f.levels2.null
if(interaction.object@f.levels1.null == FALSE & interaction.object@f.levels2.null == FALSE){
factor1 <- interaction.object@df[,interaction.object@f1.label]
factor2 <- interaction.object@df[,interaction.object@f2.label]
inter <- interaction.object@df[,"interaction"]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f1 <- LSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- LSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- LSD.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f1 <- HSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- HSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- HSD.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f1 <- duncan.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- duncan.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- duncan.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
comp.f1[,2] <- comp.f1[,1]
comp.f1[,1] <- row.names(comp.f1)
names(comp.f1) <- c("comp", "pvalue")
comp.f1$group1 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",1))
comp.f1$group2 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",2))
rownames(comp.f1) <- c(1:dim(comp.f1)[1])
comp.f2[,2] <- comp.f2[,1]
comp.f2[,1] <- row.names(comp.f2)
names(comp.f2) <- c("comp", "pvalue")
comp.f2$group1 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",1))
comp.f2$group2 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",2))
rownames(comp.f2) <- c(1:dim(comp.f2)[1])
comp.inter[,2] <- comp.inter[,1]
comp.inter[,1] <- row.names(comp.inter)
names(comp.inter) <- c("comp", "pvalue")
comp.inter$group1 <- unlist(lapply(strsplit(comp.inter[,1],split = " - "),"[",1))
comp.inter$group2 <- unlist(lapply(strsplit(comp.inter[,1],split = " - "),"[",2))
rownames(comp.inter) <- c(1:dim(comp.inter)[1])
#BF computations F1
for(i in c(1:dim(comp.f1)[1])){
crit <- interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group1"] | interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f1.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f1[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
#BF computations F2
for(i in c(1:dim(comp.f2)[1])){
crit <- interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group1"] | interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f2.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f2[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
#BF computations inter
for(i in c(1:dim(comp.inter)[1])){
crit <- interaction.object@df[,"interaction"] == comp.inter[i,"group1"] | interaction.object@df[,"interaction"] == comp.inter[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", "interaction",sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.inter[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f1.comps <- comp.f1
return.results@f2.comps <- comp.f2
return.results@interaction.comps <- comp.inter
return(return.results)
}
if(interaction.object@f.levels1.null == FALSE & interaction.object@f.levels2.null == TRUE){
factor1 <- interaction.object@df[,interaction.object@f1.label]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f1 <- LSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f1 <- HSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f1 <- duncan.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
comp.f1[,2] <- comp.f1[,1]
comp.f1[,1] <- row.names(comp.f1)
names(comp.f1) <- c("comp", "pvalue")
comp.f1$group1 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",1))
comp.f1$group2 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",2))
rownames(comp.f1) <- c(1:dim(comp.f1)[1])
#F1 BF computations
for(i in c(1:dim(comp.f1)[1])){
crit <- interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group1"] | interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f1.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f1[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f1.comps <- comp.f1
return(return.results)
}
if(interaction.object@f.levels1.null == TRUE & interaction.object@f.levels2.null == FALSE){
factor2 <- interaction.object@df[,interaction.object@f2.label]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f2 <- LSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f2 <- HSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f2 <- duncan.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
comp.f2[,2] <- comp.f2[,1]
comp.f2[,1] <- row.names(comp.f2)
names(comp.f2) <- c("comp", "pvalue")
comp.f2$group1 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",1))
comp.f2$group2 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",2))
rownames(comp.f2) <- c(1:dim(comp.f2)[1])
#F2 BF computations
for(i in c(1:dim(comp.f2)[1])){
crit <- interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group1"] | interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f2.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f2[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f2.comps <- comp.f2
return(return.results)
}
}
setOldClass(c("gg"))
setClass("post.hoc.plot",
representation(plot.data.f1 = "data.frame",
plot.data.f2 = "data.frame",
plot.data.inter = "data.frame",
plot.comp.f1 = "data.frame",
plot.comp.f2 = "data.frame",
plot.comp.inter = "data.frame",
plot.f1 = "gg",
plot.f2 = "gg",
plot.inter = "gg"
))
post_hoc_plot <- function(post.hoc.object,
interaction.object,
plot.title.f1.par = "Parametric Comparisons for F1",
plot.title.f2.par ="Parametric Comparisons for F2",
plot.title.inter.par = "Parametric Comparisons for Interactions",
plot.title.f1.bayes ="Bayes Factor for F1",
plot.title.f2.bayes ="Bayes Factor for F2",
plot.title.inter.bayes = "Bayes Factor for Interactions",
x.axis.labs = c("Factor no 1", "Factor no 2", "Interactions"),
y.axis.lab = "Mean +/- 1 SE",
BF.criteria = 5,
p.val.criteria = 0.1,
y.axis.length.factor = 1.5,
str.width = 15
){
return_results <- new("post.hoc.plot")
#Factor I
if(post.hoc.object@f.levels1.null == FALSE){
w.data.f1.comps <- post.hoc.object@f1.comps
w.data.f1.vals <- interaction.object@f1.averages.se
factor.levels <- interaction.object@f.levels1
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.f1.vals[,"Factor"] <- str_wrap(w.data.f1.vals[,"Factor"], width = str.width)
w.data.f1.comps[,"group1"] <- str_wrap(w.data.f1.comps[,"group1"], width = str.width)
w.data.f1.comps[,"group2"] <- str_wrap(w.data.f1.comps[,"group2"], width = str.width)
w.data.f1.vals[,"Factor"] <- factor(w.data.f1.vals[,"Factor"],levels = factor.levels)
w.data.f1.comps[,"group1"] <- factor(w.data.f1.comps[,"group1"], levels = factor.levels)
w.data.f1.comps[,"group2"] <- factor(w.data.f1.comps[,"group2"], levels = factor.levels)
return_results@plot.data.f1 <- w.data.f1.vals
return_results@plot.comp.f1 <- w.data.f1.comps
names(return_results@plot.comp.f1) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
y_axis_max <- (RobustMax(w.data.f1.vals[,"Mean"]) + 3*RobustMax(w.data.f1.vals[,"SE"]))
y_axis_min <- (RobustMin(w.data.f1.vals[,"Mean"]) - RobustMin(w.data.f1.vals[,"SE"]))
y_axis_max_p <- (RobustMax(w.data.f1.vals[,"Mean"]) + 3*RobustMax(w.data.f1.vals[,"SE"]) + 0.025*dim(w.data.f1.comps)[1])
f1.plot <- ggplot(data = w.data.f1.vals, aes(x= Factor, y= Mean))
f1.plot <- f1.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor))
f1.plot <- f1.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
f1.plot <- f1.plot + theme_fivethirtyeight() + scale_colour_economist()
f1.plot <- f1.plot + geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
#Annotations based on p-vals
base.length <- 1.05
f1.plot.par <- f1.plot
for (i in c(1:dim(w.data.f1.comps)[1])){
if(w.data.f1.comps[i,"pvalue"] < p.val.criteria){
f1.plot.par <- f1.plot.par + annotate("segment", x = w.data.f1.comps[i,"group1"],
xend = w.data.f1.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
f1.plot.BF <- f1.plot
for (i in c(1:dim(w.data.f1.comps)[1])){
if(w.data.f1.comps[i,"BF"] > BF.criteria){
f1.plot.BF <- f1.plot.BF + annotate("segment", x = w.data.f1.comps[i,"group1"],
xend = w.data.f1.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
f1.plot.BF <- f1.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
f1.plot.par <- f1.plot.par + ggtitle(plot.title.f1.par) + theme(plot.title = element_text(size = 10, face = "bold"))
f1.plot.BF <- f1.plot.BF + ggtitle(plot.title.f1.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))
f1.plot <- ggdraw() +
draw_plot(f1.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(f1.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.f1 <- f1.plot
}
#Factor II
if(post.hoc.object@f.levels2.null == FALSE){
w.data.f2.comps <- post.hoc.object@f2.comps
w.data.f2.vals <- interaction.object@f2.averages.se
factor.levels <- interaction.object@f.levels2
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.f2.vals[,"Factor"] <- str_wrap(w.data.f2.vals[,"Factor"], width = str.width)
w.data.f2.comps[,"group1"] <- str_wrap(w.data.f2.comps[,"group1"], width = str.width)
w.data.f2.comps[,"group2"] <- str_wrap(w.data.f2.comps[,"group2"], width = str.width)
w.data.f2.vals[,"Factor"] <- factor(w.data.f2.vals[,"Factor"],levels = factor.levels)
w.data.f2.comps[,"group1"] <- factor(w.data.f2.comps[,"group1"], levels = factor.levels)
w.data.f2.comps[,"group2"] <- factor(w.data.f2.comps[,"group2"], levels = factor.levels)
return_results@plot.data.f2 <- w.data.f2.vals
return_results@plot.comp.f2 <- w.data.f2.comps
names(return_results@plot.comp.f2) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
y_axis_max <- (RobustMax(w.data.f2.vals[,"Mean"]) + 3*RobustMax(w.data.f2.vals[,"SE"]))
y_axis_min <- (RobustMin(w.data.f2.vals[,"Mean"]) - RobustMin(w.data.f2.vals[,"SE"]))
y_axis_max_p <- (RobustMax(w.data.f2.vals[,"Mean"]) + 3*RobustMax(w.data.f2.vals[,"SE"]) + 0.025*dim(w.data.f2.comps)[1])
f2.plot <- ggplot(data = w.data.f2.vals, aes(x= Factor, y= Mean))
f2.plot <- f2.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor))
f2.plot <- f2.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
f2.plot <- f2.plot + theme_fivethirtyeight() + scale_colour_economist()
f2.plot <- f2.plot + geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
#Annotations based on p-vals
base.length <- 1.05
f2.plot.par <- f2.plot
for (i in c(1:dim(w.data.f2.comps)[1])){
if(w.data.f2.comps[i,"pvalue"] < p.val.criteria){
f2.plot.par <- f2.plot.par + annotate("segment", x = w.data.f2.comps[i,"group1"],
xend = w.data.f2.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
f2.plot.BF <- f2.plot
for (i in c(1:dim(w.data.f2.comps)[1])){
if(w.data.f2.comps[i,"BF"] > BF.criteria){
f2.plot.BF <- f2.plot.BF + annotate("segment", x = w.data.f2.comps[i,"group1"],
xend = w.data.f2.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
f2.plot.BF <- f2.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
f2.plot.par <- f2.plot.par + ggtitle(plot.title.f2.par) + theme(plot.title = element_text(size = 10, face = "bold"))
f2.plot.BF <- f2.plot.BF + ggtitle(plot.title.f2.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))
f2.plot <- ggdraw() +
draw_plot(f2.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(f2.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.f2 <- f2.plot
}
#Interaction plot
if(post.hoc.object@f.levels1.null == FALSE & post.hoc.object@f.levels2.null == FALSE){
w.data.inter.comps <- post.hoc.object@interaction.comps
w.data.inter.vals <- interaction.object@inter.averages.se
factor.levels <- interaction.object@f.levels.int
w.data.inter.vals[,"Factor"] <- factor(w.data.inter.vals[,"Factor"],levels = factor.levels)
w.data.inter.comps[,"group1"] <- factor(w.data.inter.comps[,"group1"], levels = factor.levels)
w.data.inter.comps[,"group2"] <- factor(w.data.inter.comps[,"group2"], levels = factor.levels)
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.inter.vals[,"Factor"] <- str_wrap(w.data.inter.vals[,"Factor"], width = str.width)
w.data.inter.comps[,"group1"] <- str_wrap(w.data.inter.comps[,"group1"], width = str.width)
w.data.inter.comps[,"group2"] <- str_wrap(w.data.inter.comps[,"group2"], width = str.width)
return_results@plot.data.inter <- w.data.inter.vals
return_results@plot.comp.inter <- w.data.inter.comps
names(return_results@plot.comp.inter) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
return_results@plot.comp.inter[,1] <- gsub("###", "***",(return_results@plot.comp.inter[,1]))
return_results@plot.comp.inter[,1] <- gsub("-", "vs",(return_results@plot.comp.inter[,1]))
y_axis_max <- (RobustMax(w.data.inter.vals[,"Mean"]) + 3*RobustMax(w.data.inter.vals[,"SE"]))
y_axis_min <- (RobustMin(w.data.inter.vals[,"Mean"]) - RobustMin(w.data.inter.vals[,"SE"]))
y_axis_max_p <- (RobustMax(w.data.inter.vals[,"Mean"]) + 3*RobustMax(w.data.inter.vals[,"SE"]) + 0.025*dim(w.data.inter.comps)[1])
inter.plot <- ggplot(data = w.data.inter.vals, aes(x= Factor, y= Mean))
inter.plot <- inter.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor*1.5))
inter.plot <- inter.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
inter.plot <- inter.plot + theme_fivethirtyeight() + scale_colour_economist()
inter.plot <- inter.plot + geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
inter.plot <- inter.plot + theme(axis.text.x = element_text(angle = 90, hjust = 1))
#Annotations based on p-vals
base.length <- 1.05
inter.plot.par <- inter.plot
for (i in c(1:dim(w.data.inter.comps)[1])){
if(w.data.inter.comps[i,"pvalue"] < p.val.criteria){
inter.plot.par <- inter.plot.par + annotate("segment", x = w.data.inter.comps[i,"group1"],
xend = w.data.inter.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
inter.plot.BF <- inter.plot
for (i in c(1:dim(w.data.inter.comps)[1])){
if(w.data.inter.comps[i,"BF"] > BF.criteria){
inter.plot.BF <- inter.plot.BF + annotate("segment", x = w.data.inter.comps[i,"group1"],
xend = w.data.inter.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
inter.plot.BF <- inter.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
inter.plot.par <- inter.plot.par + ggtitle(plot.title.inter.par) + theme(plot.title = element_text(size = 10, face = "bold"))
inter.plot.BF <- inter.plot.BF + ggtitle(plot.title.inter.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))
inter.plot <- ggdraw() +
draw_plot(inter.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(inter.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.inter <- inter.plot
}
return(return_results)
}
library(ggplot2)
library(car)
library(agricolae)
library(BayesFactor)
library(cowplot)
library(ggthemes)
library(stringr)
library(BEST)
#Helper functions
s.err <- function(x) sd(x,na.rm = TRUE)/sqrt(sum(!is.na(x)))
RobustMax <- function(x) {if (sum(!is.na(x))>0) max(x,na.rm = TRUE) else NA}
RobustMin <- function(x) {if (sum(!is.na(x))>0) min(x,na.rm = TRUE) else NA}
setClass("interaction",
representation(df = "data.frame",
f1.label = "character",
f2.label = "character",
value.label = "character",
f.levels1 = "character",
f.levels2 = "character",
f.levels.int = "character",
f.levels1.null = "logical",
f.levels2.null = "logical",
f1.averages.se = "data.frame",
f2.averages.se = "data.frame",
inter.averages.se = "data.frame"
))
interaction.maker <- function(df,factor1=NULL,factor2=NULL,value, credvalue = 0.68){
return.results <- new("interaction")
df[,"interaction"] <- rep(NA, times = dim(df)[1])
if(is.null(factor1) == TRUE & is.null(factor2) == TRUE){
return.results@f.levels1.null <- TRUE
return.results@f.levels2.null <- TRUE
return(return.results)
}
if(is.null(factor1) == FALSE & is.null(factor2) == TRUE){
return.results@f1.label <- factor1
return.results@value.label <- value
df[,factor1] <- as.character(df[,factor1])
return.results@df <- df
df[,factor1] <- as.factor(df[,factor1])
return.results@f.levels1 <- levels(df[,factor1])
return.results@f.levels1.null <- FALSE
return.results@f.levels2.null <- TRUE
#Averages and SE computations
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels1)*3),
nrow = length(return.results@f.levels1),ncol = 3))
results[,1] <- return.results@f.levels1
results$HDIlow <- NA
results$HDIhigh <- NA
names(results) <- c("Factor", "Mean", "SE", "HDIlow", "HDIhigh")
for (i in results[,"Factor"]){
subset <- df[df[,factor1] == i,]
results[results[,"Factor"] == i,c("Mean", "SE", "HDIlow", "HDIhigh")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]),
hdi(subset[,value], credMass = credvalue)[1], hdi(subset[,value], credMass = credvalue)[2])
}
return.results@f1.averages.se <- results
#Returning results
return(return.results)
}
if(is.null(factor1) == TRUE & is.null(factor2) == FALSE){
return.results@f2.label <- factor2
return.results@value.label <- value
df[,factor2] <- as.character(df[,factor2])
return.results@df <- df
df[,factor2] <- as.factor(df[,factor2])
return.results@f.levels2 <- levels(df[,factor2])
return.results@f.levels1.null <- TRUE
return.results@f.levels2.null <- FALSE
#Averages and SE computations
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels2)*3),
nrow = length(return.results@f.levels2),ncol = 3))
results[,1] <- return.results@f.levels2
results$HDIlow <- NA
results$HDIhigh <- NA
names(results) <- c("Factor", "Mean", "SE", "HDIlow", "HDIhigh")
for (i in results[,"Factor"]){
subset <- df[df[,factor2] == i,]
results[results[,"Factor"] == i,c("Mean", "SE", "HDIlow", "HDIhigh")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]),
hdi(subset[,value], credMass = credvalue)[1], hdi(subset[,value], credMass = credvalue)[2])
}
return.results@f2.averages.se <- results
return(return.results)
}
if(is.null(factor1) == FALSE & is.null(factor2) == FALSE){
return.results@f.levels1.null <- FALSE
return.results@f.levels2.null <- FALSE
return.results@f1.label <- factor1
return.results@f2.label <- factor2
return.results@value.label <- value
#Making interaction term
factor1_vec <- unique(df[,factor1])
factor2_vec <- unique(df[,factor2])
combined_vec <- c()
n = 0
for (i in c(1:length(factor1_vec))){
for (m in c(1:length(factor2_vec))){
combined_vec[n+1] <- paste(factor1_vec[i],"###",factor2_vec[m], sep = "")
n <- n+1
}
}
#Adding interaction term to df
for(i in c(1:dim(df)[1])){
df[i,"interaction"] <- paste(df[i,factor1],"###",df[i,factor2], sep = "")
}
df[,"interaction"] <- factor(df[,"interaction"], levels = combined_vec)
df[,"interaction"] <- as.character(df[,"interaction"])
return.results@df <- df
df[,factor1] <- as.factor(df[,factor1])
df[,factor2] <- as.factor(df[,factor2])
return.results@f.levels1 <- levels(df[,factor1])
return.results@f.levels2 <- levels(df[,factor2])
return.results@f.levels.int <- combined_vec
#Computing the averages
#First factor
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels1)*3),
nrow = length(return.results@f.levels1),ncol = 3))
results[,1] <- return.results@f.levels1
results$HDIlow <- NA
results$HDIhigh <- NA
names(results) <- c("Factor", "Mean", "SE", "HDIlow", "HDIhigh")
for (i in results[,"Factor"]){
subset <- df[df[,factor1] == i,]
results[results[,"Factor"] == i,c("Mean", "SE", "HDIlow", "HDIhigh")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]),
hdi(subset[,value], credMass = credvalue)[1], hdi(subset[,value], credMass = credvalue)[2])
}
return.results@f1.averages.se <- results
#Seccond factor
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels2)*3),
nrow = length(return.results@f.levels2),ncol = 3))
results[,1] <- return.results@f.levels2
results$HDIlow <- NA
results$HDIhigh <- NA
names(results) <- c("Factor", "Mean", "SE", "HDIlow", "HDIhigh")
for (i in results[,"Factor"]){
subset <- df[df[,factor2] == i,]
results[results[,"Factor"] == i,c("Mean", "SE", "HDIlow", "HDIhigh")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]),
hdi(subset[,value], credMass = credvalue)[1], hdi(subset[,value], credMass = credvalue)[2])
}
return.results@f2.averages.se <- results
#Interaction
results <- as.data.frame(matrix(rep(NA, times = length(return.results@f.levels.int)*3),
nrow = length(return.results@f.levels.int),ncol = 3))
results[,1] <- return.results@f.levels.int
results$HDIlow <- NA
results$HDIhigh <- NA
names(results) <- c("Factor", "Mean", "SE", "HDIlow", "HDIhigh")
for (i in results[,"Factor"]){
subset <- df[df[,"interaction"] == i,]
results[results[,"Factor"] == i,c("Mean", "SE", "HDIlow", "HDIhigh")] <- c(mean(subset[,value],na.rm = TRUE), s.err(subset[,value]),
hdi(subset[,value], credMass = credvalue)[1], hdi(subset[,value], credMass = credvalue)[2])
}
return.results@inter.averages.se <- results
return(return.results)
}
}
setClass("post.hoc",
representation(df = "data.frame",
f1.comps = "data.frame",
f2.comps = "data.frame",
interaction.comps = "data.frame",
f.levels1.null = "logical",
f.levels2.null = "logical",
post.hoc.type = "character"
))
post_hoc <- function(interaction.object,
post.hoc.type = "LSD.test"
){
return.results <- new("post.hoc")
return.results@df <- interaction.object@df
return.results@post.hoc.type <- post.hoc.type
return.results@f.levels1.null <- interaction.object@f.levels1.null
return.results@f.levels2.null <- interaction.object@f.levels2.null
if(interaction.object@f.levels1.null == FALSE & interaction.object@f.levels2.null == FALSE){
factor1 <- interaction.object@df[,interaction.object@f1.label]
factor2 <- interaction.object@df[,interaction.object@f2.label]
inter <- interaction.object@df[,"interaction"]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f1 <- LSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- LSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- LSD.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f1 <- HSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- HSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- HSD.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f1 <- duncan.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
comp.f2 <- duncan.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
comp.inter <- duncan.test((aov(lm(value ~ inter, data=interaction.object@df))),
"inter",console=FALSE,group = F)$comparison[2]
}
comp.f1[,2] <- comp.f1[,1]
comp.f1[,1] <- row.names(comp.f1)
names(comp.f1) <- c("comp", "pvalue")
comp.f1$group1 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",1))
comp.f1$group2 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",2))
rownames(comp.f1) <- c(1:dim(comp.f1)[1])
comp.f2[,2] <- comp.f2[,1]
comp.f2[,1] <- row.names(comp.f2)
names(comp.f2) <- c("comp", "pvalue")
comp.f2$group1 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",1))
comp.f2$group2 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",2))
rownames(comp.f2) <- c(1:dim(comp.f2)[1])
comp.inter[,2] <- comp.inter[,1]
comp.inter[,1] <- row.names(comp.inter)
names(comp.inter) <- c("comp", "pvalue")
comp.inter$group1 <- unlist(lapply(strsplit(comp.inter[,1],split = " - "),"[",1))
comp.inter$group2 <- unlist(lapply(strsplit(comp.inter[,1],split = " - "),"[",2))
rownames(comp.inter) <- c(1:dim(comp.inter)[1])
#BF computations F1
for(i in c(1:dim(comp.f1)[1])){
crit <- interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group1"] | interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f1.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f1[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
#BF computations F2
for(i in c(1:dim(comp.f2)[1])){
crit <- interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group1"] | interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f2.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f2[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
#BF computations inter
for(i in c(1:dim(comp.inter)[1])){
crit <- interaction.object@df[,"interaction"] == comp.inter[i,"group1"] | interaction.object@df[,"interaction"] == comp.inter[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", "interaction",sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.inter[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f1.comps <- comp.f1
return.results@f2.comps <- comp.f2
return.results@interaction.comps <- comp.inter
return(return.results)
}
if(interaction.object@f.levels1.null == FALSE & interaction.object@f.levels2.null == TRUE){
factor1 <- interaction.object@df[,interaction.object@f1.label]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f1 <- LSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f1 <- HSD.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f1 <- duncan.test((aov(lm(value ~ factor1, data=interaction.object@df))),
"factor1",console=FALSE,group = F)$comparison[2]
}
comp.f1[,2] <- comp.f1[,1]
comp.f1[,1] <- row.names(comp.f1)
names(comp.f1) <- c("comp", "pvalue")
comp.f1$group1 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",1))
comp.f1$group2 <- unlist(lapply(strsplit(comp.f1[,1],split = " - "),"[",2))
rownames(comp.f1) <- c(1:dim(comp.f1)[1])
#F1 BF computations
for(i in c(1:dim(comp.f1)[1])){
crit <- interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group1"] | interaction.object@df[,interaction.object@f1.label] == comp.f1[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f1.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f1[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f1.comps <- comp.f1
return(return.results)
}
if(interaction.object@f.levels1.null == TRUE & interaction.object@f.levels2.null == FALSE){
factor2 <- interaction.object@df[,interaction.object@f2.label]
value <- interaction.object@df[,interaction.object@value.label]
if(post.hoc.type == "LSD.test"){
comp.f2 <- LSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "HSD.test"){
comp.f2 <- HSD.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
if(post.hoc.type == "duncan.test"){
comp.f2 <- duncan.test((aov(lm(value ~ factor2, data=interaction.object@df))),
"factor2",console=FALSE,group = F)$comparison[2]
}
comp.f2[,2] <- comp.f2[,1]
comp.f2[,1] <- row.names(comp.f2)
names(comp.f2) <- c("comp", "pvalue")
comp.f2$group1 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",1))
comp.f2$group2 <- unlist(lapply(strsplit(comp.f2[,1],split = " - "),"[",2))
rownames(comp.f2) <- c(1:dim(comp.f2)[1])
#F2 BF computations
for(i in c(1:dim(comp.f2)[1])){
crit <- interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group1"] | interaction.object@df[,interaction.object@f2.label] == comp.f2[i,"group2"]
subset <- interaction.object@df[crit,]
comp.formula <- paste(interaction.object@value.label," ~ ", interaction.object@f2.label,sep = "")
BF.res = ttestBF(formula = as.formula(comp.formula), data = subset)
comp.f2[i,"BF"] <- abs(BF.res@bayesFactor[,"bf"])
}
return.results@f2.comps <- comp.f2
return(return.results)
}
}
setOldClass(c("gg"))
setClass("post.hoc.plot",
representation(plot.data.f1 = "data.frame",
plot.data.f2 = "data.frame",
plot.data.inter = "data.frame",
plot.comp.f1 = "data.frame",
plot.comp.f2 = "data.frame",
plot.comp.inter = "data.frame",
plot.f1 = "gg",
plot.f2 = "gg",
plot.inter = "gg"
))
post_hoc_plot <- function(post.hoc.object,
interaction.object,
plot.title.f1.par = "Parametric Comparisons for F1",
plot.title.f2.par ="Parametric Comparisons for F2",
plot.title.inter.par = "Parametric Comparisons for Interactions",
plot.title.f1.bayes ="Bayes Factor for F1",
plot.title.f2.bayes ="Bayes Factor for F2",
plot.title.inter.bayes = "Bayes Factor for Interactions",
x.axis.labs = c("Factor no 1", "Factor no 2", "Interactions"),
y.axis.lab = "Mean +/- 1 SE",
BF.criteria = 5,
p.val.criteria = 0.1,
y.axis.length.factor = 1.5,
str.width = 15
){
return_results <- new("post.hoc.plot")
if(interaction.object@f.levels1.null == TRUE & interaction.object@f.levels2.null == TRUE){
return(NULL)
}
#Factor I
if(post.hoc.object@f.levels1.null == FALSE){
w.data.f1.comps <- post.hoc.object@f1.comps
w.data.f1.vals <- interaction.object@f1.averages.se
factor.levels <- interaction.object@f.levels1
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.f1.vals[,"Factor"] <- str_wrap(w.data.f1.vals[,"Factor"], width = str.width)
w.data.f1.comps[,"group1"] <- str_wrap(w.data.f1.comps[,"group1"], width = str.width)
w.data.f1.comps[,"group2"] <- str_wrap(w.data.f1.comps[,"group2"], width = str.width)
w.data.f1.vals[,"Factor"] <- factor(w.data.f1.vals[,"Factor"],levels = factor.levels)
w.data.f1.comps[,"group1"] <- factor(w.data.f1.comps[,"group1"], levels = factor.levels)
w.data.f1.comps[,"group2"] <- factor(w.data.f1.comps[,"group2"], levels = factor.levels)
return_results@plot.data.f1 <- w.data.f1.vals
return_results@plot.comp.f1 <- w.data.f1.comps
names(return_results@plot.comp.f1) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
y_axis_max <- max((RobustMax(w.data.f1.vals[,"Mean"])), RobustMax(w.data.f1.vals[,"HDIhigh"]), RobustMax(w.data.f1.vals[,"SE"]))
y_axis_min <- min(RobustMin(w.data.f1.vals[,"Mean"]), RobustMin(w.data.f1.vals[,"HDIlow"]))
y_axis_max_p <- (RobustMax(w.data.f1.vals[,"Mean"]) + 0.5*RobustMax(w.data.f1.vals[,"HDIhigh"]) + 0.025*dim(w.data.f1.comps)[1])
f1.plot <- ggplot(data = w.data.f1.vals, aes(x= Factor, y= Mean))
f1.plot <- f1.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor))
f1.plot <- f1.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
f1.plot <- f1.plot + theme_fivethirtyeight() + scale_colour_economist()
#Annotations based on p-vals
base.length <- 1.05
f1.plot.par <- f1.plot
for (i in c(1:dim(w.data.f1.comps)[1])){
if(w.data.f1.comps[i,"pvalue"] < p.val.criteria){
f1.plot.par <- f1.plot.par + annotate("segment", x = w.data.f1.comps[i,"group1"],
xend = w.data.f1.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
f1.plot.BF <- f1.plot
for (i in c(1:dim(w.data.f1.comps)[1])){
if(w.data.f1.comps[i,"BF"] > BF.criteria){
f1.plot.BF <- f1.plot.BF + annotate("segment", x = w.data.f1.comps[i,"group1"],
xend = w.data.f1.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
f1.plot.BF <- f1.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
f1.plot.par <- f1.plot.par + ggtitle(plot.title.f1.par) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
f1.plot.BF <- f1.plot.BF + ggtitle(plot.title.f1.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=HDIlow, ymax=HDIhigh),
width=.2,
position=position_dodge(.9), size = 1)
f1.plot <- ggdraw() +
draw_plot(f1.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(f1.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.f1 <- f1.plot
}
#Factor II
if(post.hoc.object@f.levels2.null == FALSE){
w.data.f2.comps <- post.hoc.object@f2.comps
w.data.f2.vals <- interaction.object@f2.averages.se
factor.levels <- interaction.object@f.levels2
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.f2.vals[,"Factor"] <- str_wrap(w.data.f2.vals[,"Factor"], width = str.width)
w.data.f2.comps[,"group1"] <- str_wrap(w.data.f2.comps[,"group1"], width = str.width)
w.data.f2.comps[,"group2"] <- str_wrap(w.data.f2.comps[,"group2"], width = str.width)
w.data.f2.vals[,"Factor"] <- factor(w.data.f2.vals[,"Factor"],levels = factor.levels)
w.data.f2.comps[,"group1"] <- factor(w.data.f2.comps[,"group1"], levels = factor.levels)
w.data.f2.comps[,"group2"] <- factor(w.data.f2.comps[,"group2"], levels = factor.levels)
return_results@plot.data.f2 <- w.data.f2.vals
return_results@plot.comp.f2 <- w.data.f2.comps
names(return_results@plot.comp.f2) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
y_axis_max <- max((RobustMax(w.data.f2.vals[,"Mean"])), RobustMax(w.data.f2.vals[,"HDIhigh"]),RobustMax(w.data.f2.vals[,"SE"]))
y_axis_min <- min(RobustMin(w.data.f2.vals[,"Mean"]), RobustMin(w.data.f2.vals[,"HDIlow"]))
y_axis_max_p <- (RobustMax(w.data.f2.vals[,"Mean"]) + 0.5*RobustMax(w.data.f2.vals[,"HDIhigh"]) + 0.025*dim(w.data.f2.comps)[1])
f2.plot <- ggplot(data = w.data.f2.vals, aes(x= Factor, y= Mean))
f2.plot <- f2.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor))
f2.plot <- f2.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
f2.plot <- f2.plot + theme_fivethirtyeight() + scale_colour_economist()
#Annotations based on p-vals
base.length <- 1.05
f2.plot.par <- f2.plot
for (i in c(1:dim(w.data.f2.comps)[1])){
if(w.data.f2.comps[i,"pvalue"] < p.val.criteria){
f2.plot.par <- f2.plot.par + annotate("segment", x = w.data.f2.comps[i,"group1"],
xend = w.data.f2.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
f2.plot.BF <- f2.plot
for (i in c(1:dim(w.data.f2.comps)[1])){
if(w.data.f2.comps[i,"BF"] > BF.criteria){
f2.plot.BF <- f2.plot.BF + annotate("segment", x = w.data.f2.comps[i,"group1"],
xend = w.data.f2.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
f2.plot.BF <- f2.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
f2.plot.par <- f2.plot.par + ggtitle(plot.title.f2.par) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
f2.plot.BF <- f2.plot.BF + ggtitle(plot.title.f2.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=HDIlow, ymax=HDIhigh),
width=.2,
position=position_dodge(.9), size = 1)
f2.plot <- ggdraw() +
draw_plot(f2.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(f2.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.f2 <- f2.plot
}
#Interaction plot
if(post.hoc.object@f.levels1.null == FALSE & post.hoc.object@f.levels2.null == FALSE){
w.data.inter.comps <- post.hoc.object@interaction.comps
w.data.inter.vals <- interaction.object@inter.averages.se
factor.levels <- interaction.object@f.levels.int
w.data.inter.vals[,"Factor"] <- factor(w.data.inter.vals[,"Factor"],levels = factor.levels)
w.data.inter.comps[,"group1"] <- factor(w.data.inter.comps[,"group1"], levels = factor.levels)
w.data.inter.comps[,"group2"] <- factor(w.data.inter.comps[,"group2"], levels = factor.levels)
factor.levels <- str_wrap(factor.levels, width = str.width)
w.data.inter.vals[,"Factor"] <- str_wrap(w.data.inter.vals[,"Factor"], width = str.width)
w.data.inter.comps[,"group1"] <- str_wrap(w.data.inter.comps[,"group1"], width = str.width)
w.data.inter.comps[,"group2"] <- str_wrap(w.data.inter.comps[,"group2"], width = str.width)
return_results@plot.data.inter <- w.data.inter.vals
return_results@plot.comp.inter <- w.data.inter.comps
names(return_results@plot.comp.inter) <- c("Comparisons","P-value", "Group1", "Group2", "Bayes Factor")
return_results@plot.comp.inter[,1] <- gsub("###", "***",(return_results@plot.comp.inter[,1]))
return_results@plot.comp.inter[,1] <- gsub("-", "vs",(return_results@plot.comp.inter[,1]))
y_axis_max <- max(RobustMax(w.data.inter.vals[,"Mean"]), RobustMax(w.data.inter.vals[,"HDIhigh"]),RobustMax(w.data.inter.vals[,"SE"]))
y_axis_min <- min(RobustMin(w.data.inter.vals[,"Mean"]), RobustMin(w.data.inter.vals[,"HDIlow"]))
y_axis_max_p <- (RobustMax(w.data.inter.vals[,"Mean"]) + 0.5*RobustMax(w.data.inter.vals[,"HDIhigh"]) + 0.025*dim(w.data.inter.comps)[1])
inter.plot <- ggplot(data = w.data.inter.vals, aes(x= Factor, y= Mean))
inter.plot <- inter.plot + coord_cartesian(ylim = c(y_axis_min/y.axis.length.factor,y_axis_max_p*y.axis.length.factor*1.5))
inter.plot <- inter.plot + geom_bar(position=position_dodge(), stat="identity", fill = "#56B4E9")
inter.plot <- inter.plot + theme_fivethirtyeight() + scale_colour_economist()
inter.plot <- inter.plot + theme(axis.text.x = element_text(angle = 90, hjust = 1))
#Annotations based on p-vals
base.length <- 1.05
inter.plot.par <- inter.plot
for (i in c(1:dim(w.data.inter.comps)[1])){
if(w.data.inter.comps[i,"pvalue"] < p.val.criteria){
inter.plot.par <- inter.plot.par + annotate("segment", x = w.data.inter.comps[i,"group1"],
xend = w.data.inter.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
#Annotations based on BF-s
base.length <- 1.05
inter.plot.BF <- inter.plot
for (i in c(1:dim(w.data.inter.comps)[1])){
if(w.data.inter.comps[i,"BF"] > BF.criteria){
inter.plot.BF <- inter.plot.BF + annotate("segment", x = w.data.inter.comps[i,"group1"],
xend = w.data.inter.comps[i,"group2"] ,
y = y_axis_max*base.length,
yend = y_axis_max*base.length, colour = "red", size = 1)
}
base.length <- base.length + 0.025
}
inter.plot.BF <- inter.plot.BF + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#F1 return plot
inter.plot.par <- inter.plot.par + ggtitle(plot.title.inter.par) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=Mean-1*SE, ymax=Mean+1*SE),
width=.2,
position=position_dodge(.9), size = 1)
inter.plot.BF <- inter.plot.BF + ggtitle(plot.title.inter.bayes) + theme(plot.title = element_text(size = 10, face = "bold"))+
geom_errorbar(aes(ymin=HDIlow, ymax=HDIhigh),
width=.2,
position=position_dodge(.9), size = 1)
inter.plot <- ggdraw() +
draw_plot(inter.plot.par, x = 0,y = 0,height = 1,width = 0.55) +
draw_plot(inter.plot.BF, x = 0.55,y = 0,height = 1,width = 0.45)
return_results@plot.inter <- inter.plot
}
return(return_results)
}