pacman::p_load(tidyverse, knitr,DiagrammeR, kableExtra)

Function to draw a probability Tree

bayes_probability_tree <- function(prior, true_positive, true_negative) {
  
  if (!all(c(prior, true_positive, true_negative) > 0) && !all(c(prior, true_positive, true_negative) < 1)) {
    stop("probabilities must be greater than 0 and less than 1.",
         call. = FALSE)
  }
  c_prior <- 1 - prior
  c_tp <- 1 - true_positive
  c_tn <- 1 - true_negative
  
  round4 <- purrr::partial(round, digits = 4)
  
  b1 <- round4(prior * true_positive)
  b2 <- round4(prior * c_tp)
  b3 <- round4(c_prior * c_tn)
  b4 <- round4(c_prior * true_negative)
  
  bp <-  round4(b1/(b1 + b3))
  
  labs <- c("Start", prior, c_prior, true_positive, c_tp, true_negative, c_tn, b1, b2, b4, b3)
  
  tree <-
    create_graph() %>%
    add_n_nodes(
      n = 11,
      type = "path",
      label = labs,
      node_aes = node_aes(
        shape = "circle",
        height = 1,
        width = 1,
        fillcolor = "steelblue",
        fontsize = 12,
        x = c(0, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8),
        y = c(0, 2, -2, 3, 1, -3, -1, 3, 1, -3, -1))) %>% 
    add_edge(
      from = 1,
      to = 2,
      edge_aes = edge_aes(
        label = "Has disease",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 1, 
      to = 3,
      edge_aes = edge_aes(
        label = "Does Not",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 2,
      to = 4,
      edge_aes = edge_aes(
        label = "Sensitivity (True +)",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 2,
      to = 5,
      edge_aes = edge_aes(
        label = "False Negative",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 3,
      to = 7,
      edge_aes = edge_aes(
        label = "False Positive",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 3,
      to = 6,
      edge_aes = edge_aes(
        label = "Specificity (True -)",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 4,
      to = 8,
      edge_aes = edge_aes(
        label = "=",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 5,
      to = 9,
      edge_aes = edge_aes(
        label = "=",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 7,
      to = 11,
      edge_aes = edge_aes(
        label = "=",
        penwidth = 2, fontsize = 12
      )
    ) %>% 
    add_edge(
      from = 6,
      to = 10,
      edge_aes = edge_aes(
        label = "=",
        penwidth = 2, fontsize = 12
      )
    ) 
  # message(glue::glue("The probability of having {prior} after testing positive is {bp}"))
  print(render_graph(tree))
  #invisible(tree)
}

1.

Suppose a medical test has a 96% chance of detecting a disease if the person has it and a 92% chance of correctly indicating that the disease is absent if the person really does not have the disease. Assume 12% of a particular population has the disease.

a. (4 pts) Create a probability tree diagram to show probabilities of outcomes. Indicate on your diagram with labels the following: sensitivity, specificity, false positive, false negative.

tree <- bayes_probability_tree(.12, .96, .92)
tree

b. What is the probability that a randomly chosen person will test positive? (Show calculations)

Pr{test +} = Pr{+ | has disease} + Pr{+ | does not}

p1 <- .12 * .96   
paste("Pr{+ | has disease} = ", p1)
## [1] "Pr{+ | has disease} =  0.1152"
p2 <- .88 * .08  
paste("Pr{+ | does not} = ", p2)
## [1] "Pr{+ | does not} =  0.0704"
p3 <- p1 + p2
paste("Pr{test +} =  ", p1, " + ", p2, " = ", p3)
## [1] "Pr{test +} =   0.1152  +  0.0704  =  0.1856"
  1. Suppose that a randomly chosen person does test positive. What is the probability that this person really has the disease?

Pr{has disease | test +} = Pr{+ | has disease} / Pr{test +} =

(p1 / p3) %>% round(4)
## [1] 0.6207

2.

The prevalence of mild myopia in adults over age 40 is 23% in the U.S. Let random variable Y denote the number with myopia out of a random sample of six.

a. Complete the table to the right (3 pts)

y = c(0, 1, 2, 3, 4, 5, 6)
notY = c(6, 5, 4, 3, 2, 1, 0)
p = 0.23
pc = 1 - p
bc <- c(1, 6, 15, 20, 15, 6, 1)
bx <- paste0(bc," * (",p, ")^", y,"^ * (",pc,")^",notY,"^")
pnj <- (bc * (p^y * pc^notY)) %>% 
       round(3)

tibble(y, notY, bx, pnj) %>% kable(caption = "Binomial table for p = .23, n = 6",
               col.names = c("# myopic", "# not myopic", "Binom Expansion", "Probability"),
               align = "c") %>% 
              column_spec (1:4,border_left = T, border_right = T) %>% 
              row_spec(0, bold = T, color = "steelblue") %>% 
              kable_styling(bootstrap_options = "striped", full_width = F, position = "center",
                            font_size = 16)
Binomial table for p = .23, n = 6
# myopic # not myopic Binom Expansion Probability
0 6 1 * (0.23)0 * (0.77)6 0.208
1 5 6 * (0.23)1 * (0.77)5 0.374
2 4 15 * (0.23)2 * (0.77)4 0.279
3 3 20 * (0.23)3 * (0.77)3 0.111
4 2 15 * (0.23)4 * (0.77)2 0.025
5 1 6 * (0.23)5 * (0.77)1 0.003
6 0 1 * (0.23)6 * (0.77)0 0.000 
paste("Total = ", sum(pnj))
## [1] "Total =  1"

b. Find Pr{Y >= 3}

paste("Pr{Y = 3} + Pr{Y = 4} + Pr{Y = 5} + Pr{Y = 6} = ", sum(pnj[4:7]))
## [1] "Pr{Y = 3} + Pr{Y = 4} + Pr{Y = 5} + Pr{Y = 6} =  0.139"

c. Find Pr{Y <= 2}

paste("1 - Pr{Y >= 3} = ", 1 - sum(pnj[4:7]))
## [1] "1 - Pr{Y >= 3} =  0.861"