Note: You may work with a partner but you must writeup and submit your own assignment. Submit your Rmarkdown (.Rmd) and the compiled pdf on Gradescope in a zip file. Include any addition files (e.g. scanned handwritten solutions) in zip file with the pdf.

Problem 1. Cancer Research in Laboratory Mice

As a reminder from homework 2, a laboratory is estimating the rate of tumorigenesis (the formation of tumors) in two strains of mice, A and B. They have tumor count data for 10 mice in strain A and 13 mice in strain B. Type A mice have been well studied, and information from other laboratories suggests that type A mice have tumor counts that are approximately Poisson-distributed. Tumor count rates for type B mice are unknown, but type B mice are related to type A mice. Assuming a Poisson sampling distribution for each group with rates \(\theta_A\) and \(\theta_B\). We assume \(\theta_A \sim \text{gamma}(120, 10)\) and \(\theta_B\sim\text{gamma}(12, 1)\). We observe \(y_A = (12,9,12,14,13,13,15,8,15,6)\) and \(y_B = (11,11,10,9,9,8,7,10,6,8,8,9,7)\). Now we will actually investigate evidence that Type A mice are have higher rates of tumor formation than Type B mice.

a.(20 points) For \(n_0 \in \{1, 2, ... , 50 \}\), obtain \(Pr(\theta_B < \theta_A \mid y_A, y_B)\) via Monte Carlo sampling for \(\theta_A \sim \text{gamma}(120, 10)\) and \(\theta_B \sim \text{gamma}(12 \times n_0, n_0)\). Make a line plot of \(Pr(\theta_B < \theta_A \mid y_A, y_B\)) vs \(n_0\). Describe how sensitive the conclusions about the event \(\{\theta_B < \theta_A\}\) are to the prior distribution on \(\theta_B\).

yA <- c(12, 9, 12, 14, 13, 13, 15, 8, 15, 6)
yB <- c(11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)

n_sim <- 10000

n0_range <- 1:50

set.seed(123)
# calculate the probability Pr(θB < θA | yA, yB) for each n0
calc_probability <- function(n0) {
  thetaA_samples <- rgamma(n_sim, shape=120, rate=10)
  thetaB_samples <- rgamma(n_sim, shape=12*n0, rate=n0)
  
  # probability that θB < θA
  mean(thetaB_samples < thetaA_samples)
}

# probabilities for each n0
probabilities <- sapply(n0_range, calc_probability)

# data frame for plotting
plot_data <- data.frame(n0 = n0_range, Probability = probabilities)

# plotting
ggplot(plot_data, aes(x = n0, y = Probability)) +
  geom_line() +
  labs(title = "Probability that θB < θA given yA, yB", x = "n0", y = "Probability")

The graph fluctuates as values of n0 vary, which means the priors strength is very sensitive. It is confined to a relatively narrow range and does not show a consistent trend as n0 increases.

b.(20 points) Repeat the previous part replacing the event \(\{\theta_B < \theta_A \}\) with the event \(\{\tilde Y_B < \tilde Y_A\}\), where \(\tilde Y_A\) and \(\tilde Y_B\) are samples from the posterior predictive distribution.

y_A = c(12, 9, 12, 14, 13, 13, 15, 8, 15, 6)
y_B = c(11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)

n_sim <- 10000

n0_range <- 1:50

set.seed(113)
# calculate the probability that the average predicted ŶB < ŶA
calc_predictive_probability <- function(n0) {
  # smpling from the posterior distributions
  postA_samples <- rgamma(n_sim, shape = 120 + sum(yA), rate = 10 + length(yA))
  postB_samples <- rgamma(n_sim, shape = 12 * n0 + sum(yB), rate = n0 + length(yB))
  
  # poisson distributed tumor counts based on the posterior theta samples
  yA_sim <- rpois(n_sim, lambda = postA_samples)
  yB_sim <- rpois(n_sim, lambda = postB_samples)
  
  # probability that average ŶB < average ŶA
  mean(yB_sim < yA_sim)
}

# probabilities for each n0
predictive_probabilities <- sapply(n0_range, calc_predictive_probability)

# data frame for plotting
plot_data <- data.frame(n0 = n0_range, Probability = predictive_probabilities)

# plotting
ggplot(plot_data, aes(x = n0, y = Probability)) +
  geom_line(color = "red") +
  labs(title = "Probability that average ŶB < ŶA given yA, yB", x = "n0", y = "Probability")

This graph shows a downward trend as n0 increases. It is confined to a much larger range and is much less sensitive than the previous graph. It can be concluded that the prior strength has a significant factor on the posterior predictive outcomes in this question.

  1. In the context of this problem, describe the meaning of the events \(\{\theta_B < \theta_A\}\) and \(\{\tilde Y_B < \tilde Y_A\}\). How are they different?

    \(\{\theta_B < \theta_A\}\) describe the rate parameters for tumor formation which are simulated using the prior gamma distributions given. In this context, we are looking for the probability that the rate parameter B is less than the rate parameter A. \(\{\tilde Y_B < \tilde Y_A\}\) describes the predicted counts of tumors observed by using the posterior gamma distribution. We first gain access to the updated rate parameters by using the posterior gamma distribution, then use these updated rate parameters to simulate the counts of tumors using the Poisson distribution. In this context, we are trying to find the probability that the counts of the sample B, provided by the poisson(posterior_B), are less than the counts of A, provided by poisson(posterior_A).

2. Posterior Predictive Model Checking

Model checking and refinement is an essential part of Bayesian data analysis. Let’s investigate the adequacy of the Poisson model for the tumor count data. Consider strain A mice only for now, and generate posterior predictive datasets \(y_A^{(1)}, ..., y_A^{(1000)}\). Each \(y_A^{(s)}\) is a sample of size \(n_A = 10\) from the Poisson distribution with parameter \(\theta_A^{(s)}\), \(\theta_A^{(s)}\) is itself a sample from the posterior distribution \(p(\theta_A \mid y_A)\) and \(y_A\) is the observed data. For each \(s\), let \(t^{(s)}\) be the sample average divided by the sample variance of \(y_A^{(s)}\).

a.(10 points) If the Poisson model was a reasonable one, what would a “typical” value \(t^{(s)}\) be? Why?

set.seed(212)

y_A = c(12, 9, 12, 14, 13, 13, 15, 8, 15, 6)

num_samples <- 1000

generate_and_compute_ts <- function() {
  
  # sample theta_A from posterior gamma
  theta_A <- rgamma(1, shape = 237, rate = 20)
  
  # generate poisson
  yA_s <- rpois(10, lambda = theta_A)
  
  # calculate mean and variance
  mean_yA_s <- mean(yA_s)
  var_yA_s <- var(yA_s)
  
  # calculate t(s)
  if (var_yA_s > 0) {
    return(mean_yA_s / var_yA_s)
  } else {
    return(NA) # avoids division by zero
  }
}

# generate dataset
ts_values <- replicate(num_samples, generate_and_compute_ts())

# remove any NA
ts_values <- na.omit(ts_values)

head(ts_values)
## [1] 1.2756827 1.1034871 0.7110061 0.6712527 1.4238026 2.0374332
mean(ts_values)
## [1] 1.286646

If the Poisson model was a reasonable one, a typical value of \(t^{s}\) would be close to one. This is because the mean and variance are equal in a reasonable Poisson model. The mean value of \(t^{s}\) for our model was around 1.29, which means this model is somewhat reasonable.

b.(10 points) In any given experiment, the realized value of \(t^{s}\) will not be exactly the “typical value” due to sampling variability. Make a histogram of \(t^{(s)}\) and compare to the observed value of this statistic, \(\frac{\text{mean}(y_A)}{\text{var}(y_A)}\). Based on this statistic, make a comment on if the Poisson model seems reasonable for these data (at least by this one metric).

yA <- c(12, 9, 12, 14, 13, 13, 15, 8, 15, 6)
mean_yA <- mean(yA)
var_yA <- var(yA)
t_observed <- mean_yA / var_yA

data <- data.frame(ts_values)


ggplot(data, aes(x = ts_values)) +
  geom_histogram(bins = 30, fill = "lightblue", color = "black", alpha = 0.7) +
  geom_vline(aes(xintercept = t_observed), color = "red", linetype = "dashed", size = 1) +
  labs(title = "Histogram of t^{(s)} and Observed t Value",
       x = "t^{(s)} Value",
       y = "Frequency") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

t_observed
## [1] 1.252081

Although the observed t value did not line up perfectly with our most frequently observed model t values, it was still very close to the mean and was not a complete outlier. Therefore, this Poisson model seems reasonable for the data that we were given.

c.(10 points) Repeat the part b) above for strain B mice, using \(Y_B\) and \(n_B = 13\) to generate the samples. Assume the prior distribution \(p(\theta_B) \sim \text{Gamma}(12, 1)\). Again make a comment on the Poisson model fit.

set.seed(210)

yB = c(11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)

num_samples <- 1000

generate_and_compute_ts <- function() {
  
  # sample theta_A from posterior gamma
  theta_B <- rgamma(1, shape = 125, rate = 14)
  
  # generate poisson
  yB_s <- rpois(10, lambda = theta_B)
  
  # calculate mean and variance
  mean_yB_s <- mean(yB_s)
  var_yB_s <- var(yB_s)
  
  # calculate t(s)
  if (var_yB_s > 0) {
    return(mean_yB_s / var_yB_s)
  } else {
    return(NA) # avoids division by zero
  }
}

# generate dataset
ts_values2 <- replicate(num_samples, generate_and_compute_ts())

# remove any NA
ts_values2 <- na.omit(ts_values2)

head(ts_values2)
## [1] 1.2558140 0.9949749 0.8861004 0.5153268 1.1177977 0.8575393
mean(ts_values2)
## [1] 1.294823
mean_yB <- mean(yB)
var_yB <- var(yB)
t_observed2 <- mean_yB / var_yB

data2 <- data.frame(ts_values2)


ggplot(data2, aes(x = ts_values2)) +
  geom_histogram(bins = 30, fill = "lightblue", color = "black", alpha = 0.7) +
  geom_vline(aes(xintercept = t_observed2), color = "red", linetype = "dashed", size = 1) +
  labs(title = "Histogram of t^{(s)} and Observed t Value",
       x = "t^{(s)} Value",
       y = "Frequency") +
  theme_minimal()

t_observed2
## [1] 3.625668

The observed t value for this data is nowhere close to the most common t value for the generated data. Given this information, we can assume that this model does not properly fit our data.

3. Interval estimation with rejection sampling.

a.(15 points) Use rejection sampling to sample from the following density:

\[p(x)= \frac{1}{4}|\text{sin}(x)|\times I\{x \in [0, 2\pi]\}\]

Use a proposal density which is uniform from 0 to \(2\pi\) and generate at least 1000 true samples from p(x). Compute and report the Monte Carlo estimate of the upper and lower bound for the 50% quantile interval using the quantile function on your samples. Compare this to the 50% HPD region calculated on the samples. What are the bounds on the HPD region? Report the length of the quantile interval and the total length of the HPD region. What explains the difference? Hint: to compute the HPD use the hdi function from the HDInterval package. As the first argument pass in density(samples), where samples is the name of your vector of true samples from the density. Set the allowSplit argument to true and use the credMass argument to set the total probability mass in the HPD region to 50%.

library(HDInterval)

# target density function
p <- function(x) (1/4) * abs(sin(x))

# uniform from 0 to 2*pi
proposal <- function(n) runif(n, min = 0, max = 2*pi)

# rejection sampling function
rejection_sampling <- function(target, proposal, n) {
  samples <- numeric(n) 
  count <- 1
  while(count <= n) {
    x_proposed <- proposal(1)
    accept_prob <- target(x_proposed) / (1/(2*pi))
    if(runif(1) < accept_prob) {
      samples[count] <- x_proposed
      count <- count + 1
    }
  }
  samples
}

# generate samples
set.seed(123)
samples <- rejection_sampling(p, proposal, 1000)

# compute 50% quantile interval
quantile_interval <- quantile(samples, probs = c(0.25, 0.75))

# compute 50% HPD interval
hpd_interval <- hdi(samples, credMass = 0.5, allowSplit = TRUE)

# Output the results
list(
  Quantile_Interval = quantile_interval,
  HPD_Interval = hpd_interval
)
## $Quantile_Interval
##      25%      75% 
## 1.657015 4.683583 
## 
## $HPD_Interval
##    lower    upper 
## 3.305142 6.017216 
## attr(,"credMass")
## [1] 0.5

b.(20 points) Plot \(p(x)\) using the curve function (base plotting) or stat_function (ggplot). Add lines corresponding to the intervals / probability regions computed in the previous part to your plot using them segments function. To ensure that the lines don’t overlap visually, for the HPD region set y0 and y1 to 0 and for the quantile interval set y0 and y1 to 0.01. Make the segments for HPD region and the segment for quantile interval different colors. Report the length of the quantile interval and the total length of the HPD region, verifying that indeed the HPD region is smaller.

set.seed(123)
  # Rejection sampling and interval construction
rejection_sampling <- function(target, proposal, n) {
  samples <- numeric(n) 
  count <- 1
  while(count <= n) {
    x_proposed <- proposal(1)
    accept_prob <- target(x_proposed) / (1/(2*pi))
    if(runif(1) < accept_prob) {
      samples[count] <- x_proposed
      count <- count + 1
    }
  }
  samples
}

  # hd_region is the result of calling hdi function
  hd_region <- hdi(samples, credMass = 0.5, allowSplit = TRUE)
  print(hd_region)
##    lower    upper 
## 3.305142 6.017216 
## attr(,"credMass")
## [1] 0.5
  print(sprintf("Total region length: %.02f", sum(hd_region[2] - hd_region[1])))
## [1] "Total region length: 2.71"
  quantile_interval <- quantile(samples, probs = c(0.25, 0.75)) # YOUR CODE HERE
  print(quantile_interval)
##      25%      75% 
## 1.657015 4.683583
  print(sprintf("Total region length: %.02f", quantile_interval[2] - quantile_interval[1]))
## [1] "Total region length: 3.03"
ggplot(data = data.frame(x = c(0, 2*pi)), aes(x = x)) +
  stat_function(fun = p, geom = "line", color = "blue") +
  geom_segment(aes(x = hd_region[1], y = 0, xend = hd_region[2], yend = 0),
               color = "red", size = 1.5) +
  geom_segment(aes(x = quantile_interval[1], y = 0.01, xend = quantile_interval[2], yend = 0.01),
               color = "green", size = 1.5) +
  labs(title = "Density of p(x) with Quantile and HPD Intervals",
       x = "x", y = "Density") +
  theme_minimal()
## Warning in geom_segment(aes(x = hd_region[1], y = 0, xend = hd_region[2], : All aesthetics have length 1, but the data has 2 rows.
## ℹ Did you mean to use `annotate()`?
## Warning in geom_segment(aes(x = quantile_interval[1], y = 0.01, xend = quantile_interval[2], : All aesthetics have length 1, but the data has 2 rows.
## ℹ Did you mean to use `annotate()`?