inst/apps/leistungstest_variable/design_aid/server_design_aid_plot.R

In stephanGit/leistungstests: Supplies Functions to calculate 'leistungstests'

# bootstrap data --------------
bootstrap_data <- reactive({

  if(ostit$sample_size > 1){
    messwerte <- ostit$valid_data$Messwert
  }else{
    return(NA)
  }

  boot_obj <-
    boot::boot(
      messwerte,
      statistic = function(data, indices){c(mean(data[indices]), sd(data[indices]))},
      R = 1000,
      sim = input$bootstrap_type,
      ran.gen = function(data, mle){rnorm(n = length(data), mean = mle[1], sd = mle[2])},
      mle = c(mean(messwerte), sd(messwerte))
    )

  boot_data <-
    data.frame(
      boot_mean = boot_obj$t[,1],
      boot_sd = boot_obj$t[,2]
    )

  #add test-result to bootstrap-statistics
  lsl <- ostit$spec_type == "lsl"
  spec <- ostit$spec
  kc <- ostit$kc

  boot_data <-
    boot_data %>%
    mutate(lsl = lsl,
           spec = spec,
           kc = kc) %>%
    mutate(test_result = ifelse(lsl,
                                boot_mean - kc * boot_sd > spec,
                                boot_mean + kc * boot_sd < spec)
    ) %>%
    mutate(test_result = factor(test_result,
                                levels = c(TRUE, FALSE), labels=c("bestanden","nicht bestanden"))
    ) %>%
    select(-c(lsl, spec, kc))

  return(boot_data)
}) %>% debounce(500)

# scales design-aid-plot-----------------
plot_scales_desgin_aid <- reactive({

  if (input$scaling_design_aid_plot == "manual"){

    mean_min <- input$mean_min_aid
    mean_max <- input$mean_max_aid
    sd_min <- input$sd_min_aid
    sd_max <- input$sd_max_aid

  }else # automatic
  {
    mean_min <- round(min(round(ostit$sample_mean * 0.5), ostit$spec, bootstrap_data()$boot_mean))
    mean_max <- round(max(round(ostit$sample_mean * 1.5), ostit$spec, bootstrap_data()$boot_mean))
    sd_min <- 0
    sd_max <- round(max(round(ostit$sample_sd * 1.5), bootstrap_data()$boot_sd))
    updateNumericInput(session, "mean_min_aid", value = mean_min)
    updateNumericInput(session, "mean_max_aid", value = mean_max)
    updateNumericInput(session, "sd_min_aid", value = sd_min)
    updateNumericInput(session, "sd_max_aid", value = sd_max)
  }

  return(list(
    mean = c(mean_min, mean_max),
    sd = c(sd_min, sd_max)
  ))

})

# design-aid-plot -----------------
output$design_aid_plot <- renderPlot({

  validate(
    need({ifelse(ostit$spec_type == "lsl",
                 ostit$sample_mean > ostit$spec,
                 ostit$sample_mean < ostit$spec)},
         "Stichprobenmittelwert außerhalb der Spezifikation"),
    need(!is.na(ostit$sample_mean), "Stichprobenmittelwert unbekannt"),
    need(!is.na(ostit$sample_sd),"Stichproben-Variabilität unbekannt. (Zur Bestimmung müssen mindestens zwei Messwerte vorliegen)")
  )

  # Create a Progress object
  {
    progress <- shiny::Progress$new()
    # Make sure it closes when we exit this reactive, even if there's an error
    on.exit(progress$close())
    progress$set(message = "Please wait - calculating data", value = 0)
  }

  design_aid_data <-
    expand.grid(
      process_mean = seq(plot_scales_desgin_aid()$mean[1], plot_scales_desgin_aid()$mean[2],
                         by = (plot_scales_desgin_aid()$mean[2] - plot_scales_desgin_aid()$mean[1])/100),
      process_sd = seq(plot_scales_desgin_aid()$sd[1], plot_scales_desgin_aid()$sd[2],
                       by = (plot_scales_desgin_aid()$sd[2] - plot_scales_desgin_aid()$sd[1])/100)
    ) %>%
    mutate(
      p_pass = (
        1 - p_reject_ostit(process_mean = process_mean,
                           process_sd = process_sd,
                           spec_limit = ostit$spec,
                           spec_type = ostit$spec_type,
                           p_min = ostit$p_min,
                           sample_size = input$sample_size_desired,
                           power = ostit$power)
      )*100
    ) %>%
    mutate(p_pass = if_else(p_pass == 100, NA_real_, p_pass))

  ggplot_object <-
    design_aid_data %>%
    # setup plot
    ggplot(aes(x = process_sd, y = process_mean)) +
    scale_x_continuous(limits = c(plot_scales_desgin_aid()$sd[1], plot_scales_desgin_aid()$sd[2])) +
    scale_y_continuous(limits = c(plot_scales_desgin_aid()$mean[1], plot_scales_desgin_aid()$mean[2])) +
    # add probability-lines
    geom_contour(aes(z = p_pass, colour = ..level..), binwidth = 1) +
    scale_color_gradient(paste("Wahrscheinlichkeit den Test mit Stichprobengröße n=",input$sample_size_desired,"\nbei dem geschätzten Prozess zu bestehen"),
                         high = "green", low = "red",
                         breaks = c(0, 25, 50, 75, 100),
                         labels = paste(c(0, 25, 50, 75, 100), "%", sep = ""),
                         limits = c(0,100))


  if (input$display_bootstrap_samples == "show"){
    ggplot_object <- ggplot_object +
      geom_point(data = bootstrap_data(), aes(x = boot_sd, y = boot_mean, shape = test_result),
                 position = ifelse(input$bootstrap_type == "ordinary", "jitter", "identity"),
                 color = "blue", alpha = 0.5, size = 4) +
      #geom_density2d(data = bootstrap_data(), aes(x = boot_sd, y = boot_mean, alpha = ..level..), color = "blue") +
      scale_alpha(guide = 'none') +
      guides(shape = guide_legend(title = paste("Ergebnis sim. Tests mit\nStichprobengröße n=", ostit$sample_size, sep=""),
                                  title.position = "top", label.position = "bottom", direction = "horizontal"))
  }

  ggplot_object +
    # add Prozess-parameter-schätzung
    geom_point(x = ostit$sample_sd, y = ostit$sample_mean, size = 10, color = "black", shape = 4) +
    # add line that separates parameter-space in i.o. and n.i.o processes
    switch(ostit$spec_type,
           lsl = geom_line(aes(x = process_sd, y = ostit$spec + qnorm(ostit$p_min) * process_sd)),
           usl = geom_line(aes(x = process_sd, y = ostit$spec - qnorm(ostit$p_min) * process_sd))
    ) +
    # add annotations to plot
    labs(
      x = "Geschätzte Produktionsprozess-Standardabweichung bzw. Variabiltät",
      y = "Geschätzter Produktions-Mittelwert",
      title = paste("Schwarzes Kreuz: Geschätzter Prozess-Mittelwert (",round(ostit$sample_mean, 2), ostit$variable_unit,
                    ") und Variabilität (",round(ostit$sample_sd, 2), ostit$variable_unit,") aus der Vorversuchsreihe"),
      subtitle = paste(
                       switch(input$display_bootstrap_samples,
                              show = paste("Blau: Mittelwert und Varianz von 1000 simulierten Leistungstests (mit n=", ostit$sample_size, ") an einem Prozess mit den geschätzten Parametern."),
                              dont_show = "Ausschluß einzelner einzelner Messwerte zeigt den Einfluß auf die geschätzten Prozessparameter."
                       ), sep = "")
    ) +
    # plot appearance
    guides(color = guide_legend(title.position = "left", label.position = "bottom", direction = "horizontal", order = 1)) +
    theme(legend.position = c(0,1), legend.justification = c(0,1))

})

# VERBATIM HINTS ------------------

# EMPFEHLUNG SAMPLE-SIZE
observe({

  sample_size_hint <-
    sample_size_hint_ostit(
      p_success = ostit$p_pass_desired,
      process_mean = ostit$sample_mean,
      process_sd = ostit$sample_sd,
      spec_limit = ostit$spec,
      spec_type = ostit$spec_type,
      p_min = ostit$p_min,
      power = ostit$power
    )

  output$sample_size_recommendation <- renderUI({tags$b(paste("Empfehlung(!) für n: ",sample_size_hint))})
})

# BERECHNETE ERFOLGSWAHRSCHEINLICHKEIT
observe({

  p_reject <-
    p_reject_ostit(
      process_mean = ostit$sample_mean,
      process_sd = ostit$sample_sd,
      sample_size = input$sample_size_desired,
      spec_limit = ostit$spec,
      spec_type = ostit$spec_type,
      p_min = ostit$p_min,
      power = ostit$power
    )

  output$p_pass_calculated <- renderUI({tags$b(paste("Geschätzte Erfolgswahrscheinlichkeit:",round((1 - p_reject)*100), "%"))})

})