# Load the shiny and longpower packages
library(shiny)
library(shinythemes)
library(longpower)
 
# UI ----------------------------------------------------------------------
ui <- fluidPage(
  # Custom theme
  theme = shinytheme("flatly"),
  # Add title
  titlePanel("Linear Mixed Model Sample Size Calculations"),
   
  # Add layout
  sidebarLayout(
    # Input
    sidebarPanel(
      selectInput(
        "typeInput",
        "Use effect size (Cohen's d) or mean difference:",
        choices = c("Effect size", "Mean difference")
      ),
      conditionalPanel(
        condition = "input.typeInput == 'Effect size'",
        numericInput(
          "numDelta",
          "Cohen's d:",
          min = 0,
          step = 0.1,
          value = 0.5
        )
      ),
      conditionalPanel(
        condition = "input.typeInput == 'Mean difference'",
        numericInput(
          "numDelta",
          "Mean difference:",
          min = 0,
          step = 0.1,
          value = 0.5
        ),
        numericInput("sdA",
                     "SD for Group A:",
                     min = 0,
                     value = 1
        ),
        numericInput("sdB",
                     "SD for Group B:",
                     min = 0,
                     value = 1
        )
      ),
      numericInput(
        "numLambda",
        "Group allocation ratio (N1:N2):",
        value = 1,
        min = 0.01,
        step = 0.1
      ),
      numericInput(
        "numTime",
        "Number of time-points:",
        min = 2,
        max = 5,
        value = 2
      ),
      numericInput(
        "corTime",
        "Correlation between time-points:",
        min = 0.1,
        max = 0.9,
        step = 0.1,
        value = 0.5
      ),
      numericInput(
        "numAlpha",
        "Alpha:",
        min = 0,
        max = 1,
        step = 0.01,
        value = 0.05
      ),
      selectInput(
        "sideInput",
        "One- or two-sided test:",
        choices = c("Two-sided", "One-sided")
      ),
      numericInput(
        "numPower",
        "Power:",
        min = 0.8,
        max = 1,
        step = 0.1,
        value = 0.8
      ),
      # Attrition
      strong("Attrition"),
      p("Specify attrition rate between the time-points as a proportion."),
      conditionalPanel(
        condition = "input.numTime >= 2",
        numericInput(
          "numAttr1",
          "Time 1 to 2:",
          min = 0,
          max = 0.99,
          step = 0.1,
          value = 0
        )
      ),
      conditionalPanel(
        condition = "input.numTime >= 3",
        numericInput(
          "numAttr2",
          "Time 1 to 3:",
          min = 0,
          max = 0.99,
          step = 0.1,
          value = 0
        )
      ),
      conditionalPanel(
        condition = "input.numTime >= 4",
        numericInput(
          "numAttr3",
          "Time 1 to 4:",
          min = 0,
          max = 0.99,
          step = 0.1,
          value = 0
        )
      ),
      conditionalPanel(
        condition = "input.numTime == 5",
        numericInput(
          "numAttr4",
          "Time 1 to 5:",
          min = 0,
          max = 1,
          step = 0.1,
          value = 0
        )
      )
    ),
     
    # Output
    mainPanel(
      h2("Overview"),
      tags$p(
        "This function calulates the sample size for a mixed model of repeated measures with a compound symmetry correlation structure."
      ),
      tags$p(
        "It uses the",
        tags$code("power.mmrm"),
        "function from the",
        tags$a(tags$code("longpower"), href = "https://cran.r-project.org/web/packages/longpower/index.html"),
        "package. See Lu, Luo, and Chen (2008) for more information."),
      tags$p("The function has been modified to round the sample size estimate up to the nearest integer."),
      h3("Important Information"),
      tags$p("The calculations assume that the correlation structure and attrition rate are the same for both groups."),
      h2("Results"),
      verbatimTextOutput("console"),
      h2("Sample Size Paragraph for Studies or Grants"),
      textOutput("samp_paragraph"),
      h3("Reference"),
      tags$p(
        "Lu, K., Luo, X., Chen, P.-Y. (2008). Sample size estimation for repeated measures analysis in randomized clinical trials with missing data.",
        tags$i("International Journal of Biostatistics"), ", ", tags$i("4"), "(1)"
      ),
      h2("Disclaimer"),
      tags$p(tags$b("Notice to any user of this website")),
      p("The author of this website makes no warranty, expressed or implied, as to the results obtained from the use of the information on the website."),
      p("The author of this website shall have no liability for the accuracy of the information and cannot be held liable for any third-party claims or losses of any damages."),
      p("The author of this website may, at any time, revise the information on this website without notice and makes no commitment to update this information.")
    ),
  )
)
 
# Server ------------------------------------------------------------------
server <- function(input, output, session) {
  # Updates the attrition so it's equal to the previous timepoint
  # Done before user specifies the value
   
  observe({
    numAttr1 <- as.numeric(input$numAttr1)
    numAttr2 <- as.numeric(input$numAttr2)
    numAttr3 <- as.numeric(input$numAttr3)
    numAttr4 <- as.numeric(input$numAttr4)
     
    if (numAttr2 < numAttr1) {
      updateNumericInput(session, "numAttr2", value = numAttr1)
    }
     
    if (numAttr3 < numAttr2) {
      updateNumericInput(session, "numAttr3", value = numAttr2)
    }
     
    if (numAttr4 < numAttr3) {
      updateNumericInput(session, "numAttr4", value = numAttr3)
    }
  })
   
  powerout <- reactive({
    rho <- as.numeric(input$corTime)
    k <- as.numeric(input$numTime)
    lambda <- as.numeric(input$numLambda)
     
    # Correlation matrix
    Ra <- matrix(rho,
                 nrow = k,
                 ncol = k
    )
     
    diag(Ra) <- 1
     
    # Delta
    delta <- as.numeric(input$numDelta)
     
    # Standard deviations
    if (input$typeInput == "Effect size") {
      sigmaa <- 1
      sigmab <- 1
    } else {
      sigmaa <- as.numeric(input$sdA)
      sigmab <- as.numeric(input$sdB)
    }
     
    power <- as.numeric(input$numPower)
    sig.level <- as.numeric(input$numAlpha)
    alternative <- ifelse(input$sideInput == "Two-sided",
                          "two.sided",
                          "one.sided"
    )
     
    # Attrition vector
    ra <- vector("numeric", length = k)
     
    for (i in 1:k) {
      if (i == 1) {
        ra[i] <- 1
      } else {
        ra[i] <- 1 - input[[paste0("numAttr", i - 1)]]
      }
    }
     
    samp <- power.mmrm(
      N = NULL,
      Ra = Ra,
      Rb = Ra,
      ra = ra,
      rb = ra,
      delta = delta,
      sigmaa = sigmaa,
      sigmab = sigmab,
      power = power,
      sig.level = sig.level,
      lambda = lambda,
      alternative = alternative
    )
     
    # Round up the sample sizes to nearest integer
    samp$n1 <- ceiling(samp$n1)
    samp$n2 <- ceiling(samp$n2)
     
    return(
      list(
        samp = samp, 
        n1 = samp$n1, 
        n2 = samp$n2, 
        attr = round(100 * (1 - ra), 2)
      )
    )
     
  })
   
  # Display output from powr.mrmm
  output$console <- renderPrint({
    powerout()$samp
  })
   
  output$samp_paragraph <- renderPrint({
     
    # Helper function for collapsing character vectors
    fPaste <- function(vec) sub(",\\s+([^,]+)$", " and \\1", toString(vec))
     
    tplParagraph <- "A sample size of %s will be required to test the primary outcome of interest. The calculation was based on a mixed model of repeated measures with a general correlation structure (Lu et al., 2008). The calculations assumed a group allocation of %s, alpha of %s, %s%% power, %s assessment points, a between-groups %s at post-intervention, a compound symmetry correlation matrix (rho = %s), and %s."
     
    # Number of time-points in written form
    num_string <- c("one", "two", "three", "four", "five")
     
    # Group allocation
    if (input$numLambda != 1) {
       
      n_text <- paste0(powerout()$n1, " and ", powerout()$n2, " in each respective group (", sum(powerout()$n1, powerout()$n2), " in total)")
       
      n_alloc <- paste0(input$numLambda, ":1")
       
    } else {
       
      n_text <- paste0(powerout()$n1, " per group (", sum(powerout()$n1, powerout()$n2), " in total)")
       
      n_alloc <- "1:1"
    }
     
    # Mean difference or effect size
    if (input$typeInput == "Mean difference") {
       
      es_text <- paste0("mean difference of ", input$numDelta, " (SD Group A = ", input$sdA, " and SD Group B = ", input$sdB, ")")
       
    } else {
       
      es_text <- paste0("effect size of Cohen's d = ", input$numDelta)
    }
     
    # Attrition Text
    if (all(powerout()$attr == 0)) {
       
      att_text <- "no attrition between the assessment points"
       
    } else {
       
      att_text <- paste0("attrition rates of ", fPaste(paste0(powerout()$attr[-1], "%")), " between assessment points")
       
    }
     
    cat(
      "\n",
      sprintf(
        tplParagraph,
        n_text,
        n_alloc,
        input$numAlpha,
        round(100 * input$numPower, 1),
        num_string[input$numTime],
        es_text,
       input$corTime,
       att_text
      ),
      "\n",
      sep = ""
    )
     
  })
}
 
shinyApp(ui, server)