library(shiny)

par2eta <- function(beta, means, S) {
  Smarg <- cov(means) + S
  betaz <- colMeans(
    (t(beta) - c(colMeans(means), rep(0, 9))) / sqrt(diag(Smarg))
  )
  s2margz <- sum(cov2cor(Smarg))/9
  s2meanz <- 1/4 * 2*(betaz["mu"]^2 + (betaz["mu"] + betaz["b2"])^2)
  sigmaz <- sqrt(s2margz - s2meanz)  # s2marg = s2means + s2
  (betaz["ab22"]/4)^2 / ((betaz["ab22"]/4)^2 + sigmaz^2)
}

ui <- fluidPage(
  withMathJax(),
  titlePanel("From model parameters to partial eta squared"),
  p("Challenge: Find a set of plausible (!) parameter values so that partial
    eta squared (\\(\\eta_p^2\\)) becomes 0.02."),
  fluidRow(
    column(3, ""),
    column(1, "\\(\\mu\\)"),
    column(1, "\\(\\alpha_2\\)"),
    column(1, "\\(\\beta_2\\)"),
    column(1, "\\((\\alpha \\beta)_{22}\\)"),
    column(1, "\\(\\sigma\\)")
  ),
  fluidRow(
    column(3, HTML("<br>Redwood (m)")),
    column(1, numericInput("muR",    "", value = 0)),
    column(1, numericInput("a2R",    "", value = 0)),
    column(1, numericInput("b2R",    "", value = 0)),
    column(1, numericInput("ab22R",  "", value = 0)),
    column(1, numericInput("sigmaR", "", value = 0))
  ),
  fluidRow(
    column(3, HTML("<br>United nations (n)")),
    column(1, numericInput("muU",    "", value = 0)),
    column(1, numericInput("a2U",    "", value = 0)),
    column(1, numericInput("b2U",    "", value = 0)),
    column(1, numericInput("ab22U",  "", value = 0)),
    column(1, numericInput("sigmaU", "", value = 0))
  ),
  fluidRow(
    column(3, HTML("<br>Cat speed (km/h)")),
    column(1, numericInput("muC",    "", value = 0)),
    column(1, numericInput("a2C",    "", value = 0)),
    column(1, numericInput("b2C",    "", value = 0)),
    column(1, numericInput("ab22C",  "", value = 0)),
    column(1, numericInput("sigmaC", "", value = 0))
  ),
  fluidRow(
    column(3, HTML("<br>Within-group correlation")),
    column(1, numericInput("rho",    "", value = 0,
                           min = 0, max = 1, step = 0.1))
  ),
  fluidRow(
    column(3, "Partial eta squared (\\(\\eta_p^2\\))"),
    column(1, textOutput("eta")),
  ),
  br(),
  plotOutput("plot")
)

server <- function(input, output) {
  n <- 388
  dat <- data.frame(
    A = factor(rep(1:2, each = n/2), labels = c("low", "high")),
    B = factor(rep(rep(1:2, each = n/4), 2), labels = c("without", "with"))
  )
  X <- model.matrix(~ A*B, dat)

  r <- reactive({
    beta <- cbind(
      redwood = c(mu = input$muR, a2 = input$a2R, b2 = input$b2R, ab22 = input$ab22R),
      unation = c(mu = input$muU, a2 = input$a2U, b2 = input$b2U, ab22 = input$ab22U),
      catspee = c(mu = input$muC, a2 = input$a2C, b2 = input$b2C, ab22 = input$ab22C)
    )
    means <- X %*% beta
    s <- c(input$sigmaR, input$sigmaU, input$sigmaC); r <- input$rho
    S <- r * s %o% s; diag(S) <- s^2  # within-group covariance
    list(beta = beta, means = means, S = S)
  })

  output$eta <- renderText(
    par2eta(r()$beta, r()$means, r()$S)  # partial eta squared from parameters
  )

  output$plot <- renderPlot({
    y <- r()$means + MASS::mvrnorm(n, mu = c(0, 0, 0), Sigma = r()$S)
   
    lattice::xyplot(redwood + unation + catspee + z ~ B,
      cbind(dat, y, z = rowMeans(scale(y))), groups = A,
      type = c("g", "a", "p"), scales = list(y = list(relation = "free")),
      xlab = "Time pressure (B)", ylab = "Numerical estimate",
      auto.key = list(title = "Anchor (A)", cex.title = 1, lines.title = 1.2,
        background = "lightgray", x = -0.01, y = 1.13, columns = 2,
        lines = TRUE))
  })
}

shinyApp(ui = ui, server = server)