#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 1200
######################################################################
## Copyright 2024 Carl F. Falk
##
## This program is free software: you can redistribute it and/or
## modify it under the terms of the GNU General Public License as
## published by the Free Software Foundation, either version 3 of
## the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
## <http://www.gnu.org/licenses/>
library(shiny)
library(xtable)
library(gtools)
library(dplyr)
library(knitr)
library(kableExtra)
library(effectsize)
#library(report) # auto-generates APA-style report...
#library(shinyjs) # not used
library(ggplot2) # not used
# Define server logic
server <- function(input, output, session) {
# create dataset and descriptives
dat <- eventReactive({input$newdat},{
# get options
J <- input$ngroups
opts <- input$genoptions
#"Normality"
#"Equal variances"
#"Ensure equal sample sizes"
#"Null hypothesis true"
# Sample size
# N size hard coded between 2 and 6
if("Ensure equal sample sizes" %in% opts){
Ng <- rep(input$nsize, J)
} else {
alph <- rep(1, length(2:6))
alph[input$nsize-1] <- 5
probs <- as.vector(rdirichlet(1, alph))
Ng <- sample(2:6, J, replace=TRUE, prob = probs)
}
# generate population parameters
if ("Normality" %in% opts){
if("Null hypothesis true" %in% opts){
mu <- rep(rnorm(1, 5, 2),J)
} else {
mu <- rnorm(J, 5, 2)
}
if("Equal variances" %in% opts){
sd <- rep(abs(rnorm(1, 2, 1.5)), J)
} else {
sd <- abs(rnorm(J, 2, 1.5))
}
} else {
# flip a coin to do: log-normal or contaminated normal
mix <- .5
whichdist <- rbinom(1, 1, prob=c(.5,.5))
if("Null hypothesis true" %in% opts){
if(whichdist == 0 ){
# log-normal
mu <- rep(rnorm(1, 0, 1),J)
} else {
# contaminated normal
mu <- rep(rnorm(1, 5, 2),J)
}
} else {
if(whichdist == 0){
# log-normal
mu <- rnorm(J, 0, 1)
} else {
# contaminated normal
mu <- rnorm(J, 5, 2)
}
}
if("Equal variances" %in% opts){
if(whichdist == 0){
# log-normal
sd <- rep(abs(rnorm(1, 1, 1.5)), J)
} else {
# contaminated normal
sd <- rep(abs(rnorm(1, 1, 1.5)), J)
}
} else {
if(whichdist == 0){
# log-normal
sd <- abs(rnorm(J, 1, 1.5))
} else {
# contaminated normal
sd <- abs(rnorm(J, 1, 1.5))
}
}
}
# generate data
Y<-vector("numeric")
X<-vector("character")
for(j in 1:J){
X<-c(X, rep(paste("Group",j), Ng[j]))
if ("Normality" %in% opts){
# normal
Y<-c(Y, round(rnorm(Ng[j], mu[j], sd[j]),2))
} else {
if(whichdist == 0){
# log-normal
Y<-c(Y, round(rlnorm(Ng[j], mu[j], sd[j]),2))
} else {
# contaminated normal? did not check carefully yet...
# mixing proportion hard-coded for now
contam <- rnorm(Ng[j], mu[j], 2*sd[j])
Y<-c(Y, round(.5*rnorm(Ng[j], mu[j], sd[j]/2)+.5*contam,2))
}
}
}
data <- data.frame(X=X,Y=Y)
out <- list(data = data
)
out
})
# display of raw data
output$rawData <- renderDataTable({
data <- dat()
data$data}
)
# display of descriptive statistics
output$desc <- function(){
d <- dat()
d$data %>% group_by(X) %>%
summarise(N = n(),
mean = mean(Y),
sd = sd(Y),
var = var(Y)) %>%
knitr::kable(format="html", digits=3) %>%
kableExtra::kable_styling("striped", full_width=T)
}
# plots
output$plots <- renderPlot({
data <- dat()
ggplot(data$data, aes(x=X, y=Y)) + geom_boxplot()
})
# F test
output$results <- renderPrint({
data <- dat()
data <- data$data
alpha <- as.numeric(input$alpha)
N <- nrow(data)
k <- input$ngroups
# df
dfM <- k-1
dfR <- N-k
dfT <- N-1
# SS
Xbar <- round(mean(data$Y), 3)
SST <- sum(round((data$Y-Xbar)^2,3)) %>% round(digits=3)
SSM <- data %>% group_by(X) %>%
summarise(SS = n()*(mean(Y)-Xbar)^2) %>%
summarise(sum(SS)) %>% unlist() %>% round(digits=3)
SSR <- data %>% group_by(X) %>%
summarise(SS = sum((Y-mean(Y))^2)) %>%
summarise(sum(SS)) %>% unlist() %>% round(digits=3)
# MS
MSM <- round(SSM/dfM, digits=3)
MSR <- round(SSR/dfR, digits=3)
# F
Fratio <- round(MSM/MSR, digits=3)
Fcrit <- round(qf(alpha, dfM, dfR, lower.tail=FALSE), digits=3)
reject <- ifelse(Fratio > Fcrit, TRUE, FALSE)
# create output
out <- ""
out <- paste0(out, "alpha: ", input$alpha, "\n")
out <- paste0(out, "Grand mean: ", Xbar, "\n")
out <- paste0(out, "SSM: ", SSM, "\n")
out <- paste0(out, "SSR: ", SSR, "\n")
out <- paste0(out, "SST: ", SST, "\n")
out <- paste0(out, "dfM: ", dfM, "\n")
out <- paste0(out, "dfR: ", dfR, "\n")
out <- paste0(out, "dfT: ", dfT, "\n")
out <- paste0(out, "MSM: ", MSM, "\n")
out <- paste0(out, "MSR: ", MSR, "\n")
out <- paste0(out, "F: ", Fratio, "\n")
out <- paste0(out, "\n")
out <- paste0(out, "Critical F: ", Fcrit, "\n")
out <- paste0(out, "\n")
out <- paste0(out, "reject H0? ", ifelse(reject, "Yes", "No"), "\n")
cat(out)
})
}
ui <- fluidPage(
# Application title
#titlePanel(""),
# Layout of UI
sidebarLayout(
sidebarPanel(strong("Data Generation Options"),
numericInput("ngroups",
"Number of groups:",
value = 3,
step = 1,
min = 2,
max = 5),
numericInput("nsize",
"Typical sample size per group",
value = 3,
step = 1,
min = 2,
max = 12),
checkboxGroupInput("genoptions",
NULL,
choices = c("Normality",
"Equal variances",
"Ensure equal sample sizes",
"Null hypothesis true"),
selected= c("Normality",
"Equal variances",
"Ensure equal sample sizes",
"Null hypothesis true")),
actionButton("newdat",
"Generate New Dataset"),
br(),
br(),
strong("Analysis options"),
selectInput("alpha",
label = "alpha:",
choices = c(.05,.005,.01,.02, .1, .2, .3, .4, .5)
)
),
# Main Panel
mainPanel(
tabsetPanel(
tabPanel("Raw Data",
dataTableOutput("rawData")),
tabPanel("Descriptive Stats",
tableOutput("desc")),
tabPanel("Results",
verbatimTextOutput("results")
),
tabPanel("Visualization",
plotOutput("plots")
)
)
)
) # end sidebarLayout
)
shinyApp(ui = ui, server = server)