Homework 03

Instructions

For this homework, you may choose to work with another classmate enrolled in the course. Choose only one member of the pair to submit at SPOC. Otherwise, submit your solutions in the usual way.
Submit your solutions in PDF format to our SPOC website. The file name should be of the form IDNumber1_IDNumber2_HW02.pdf, where IDNumber1 and IDNumber2 are the corresponding university ID numbers of the pair who worked on this homework. If you worked alone, then use the file name IDNumber_HW03.pdf. The file size should be less than 20 MB.
The deadline is on Apr 9, 2023 at 1700 Beijing time. This deadline is also indicated in the submission form at SPOC. Late submissions will not be accepted, regardless of the validity of the reason. You may still choose to have your solutions checked, but it gets zero credit if submitted late.
Make sure to click on the Submit button after you have uploaded your assignment, so that it does not stay in Draft mode.
You are allowed to discuss the exercises with other classmates, but you have to write up your own solutions. This means you are not allowed to directly copy or even look at solutions from any other source.
Some of the solutions may be chosen for presentation in class. If you or your partner are unable to explain your answers in class, then both of you will get zero credit for this homework.
If you do not provide the table in Exercise A, you automatically will get zero credit for this homework. If you do provide the table in Exercise A, at least you get a 1 point bonus right away.

Exercise A: ID numbers and contributions

Create a table which lists every classmate involved in the development of your solutions. Only put down ID numbers. You will also indicate the contributions made by each member of the pair. You will also indicate what you have discussed with other classmates, along with their ID numbers. If you work alone, then there is no need to be explicit about the contributions. But if you discussed with another classmates, make sure to indicate what you have discussed with them.

The table has the following form:

ID number	Write YES if ID number belongs to submitter.	Contribution

Here is an example, which is not very complete, but conveys the general idea.

ID number	Write YES if ID number belongs to submitter.	Contribution
ID number 1	YES	Worked on R code.
ID number 2	YES	Wrote up the solutions.
ID number 3		Discussed Exercise B Item 3.

If you are working alone, your table might look like:

ID number	Write YES if ID number belongs to submitter.	Contribution
ID number 1	YES
ID number 2		Discussed Exercise B Item 3.
ID number 3		Discussed Exercise A Item 4.

Exercise B: Chips Ahoy, Part 3

You have personally observed that the number of chocolate chips on a Chips Ahoy cookie is not the same across all cookies. In other words, there is variability in the number of chocolate chips even though these cookies are mass-produced and have quality control standards. Chips Ahoy, Part 1 gave you a chance to work through what could be a suitable model for the number of chocolate chips on a cookie. Chips Ahoy, Part 2 gave you a chance to start from the model developed in Part 1 and to produce a finite-sample confidence interval based on Chebyshev’s inequality that is quite conservative given the Monte Carlo simulation.

In this exercise, you will be constructing a confidence interval based on likelihood theory. Afterwards, you will be making comparisons of the coverage properties of the finite-sample confidence interval based on Chebyshev’s inequality and a likelihood-based confidence interval.

Recall once again that the assumed parametric model is that \(Y_1,\ldots,Y_J \sim\ \mathsf{IID}\ \mathsf{Poi}\left(\lambda\right)\). Under “regularity” conditions which are satisfied for this model, the MLE \(\widehat{\lambda}_{\mathsf{MLE}}\) has the following asymptotic distribution as \(n\to\infty\) (or, if you wish, the MLE obeys the following central limit theorem): \[\frac{\widehat{\lambda}_{\mathsf{MLE}}-\lambda}{\sqrt{1/I\left(\widehat{\lambda}_{\mathsf{MLE}}\right)}} \overset{d}{\to} N\left(0,1\right)\]

(2 points) Derive the log-likelihood function for \(\lambda\), treating the \(Y_i\)’s as random. Show the details of your work.
(2 points) Derive the score function for \(\lambda\). Show directly without using the proof presented in the notes, that the expected value is equal to zero.
(2 points) Derive the Fisher information \(I\left(\lambda\right)\). Show the details of your work. Do you notice something?
(2 points) Use the asymptotic distribution of \(\widehat{\lambda}_{\mathsf{MLE}}\) to construct a 95% confidence interval for \(\lambda\). Make sure to simplify \(L\) and \(U\).
(3 points) Adapt the code found in the computational examples I shared with you, and use nlm() to report an estimate of \(\lambda\), a standard error for \(\lambda\), and a 95% confidence interval for \(\lambda\) using the cookie dataset we have collected together in class. Compare with the finite-sample 95% confidence interval you already obtained using Chebyshev’s inequality. Comment.

# Load cookie dataset
cookie <- 
# Assign to y the relevant column of the cookie dataset
y <- 
# Setup the negative log-likelihood
mlnl <- 
# Use nlm() to find optima, you can choose to experiment with different starting points
output <- nlm(, hessian = TRUE, print.level = 2)
# Report estimate of lambda, its standard errors, and a 95% confidence interval for lambda

(4 points) Conduct a Monte Carlo simulation where you are to compare the coverage properties and the average lengths of the two procedures for constructing a confidence interval: (a) finite-sample confidence interval based on Chebyshev’s inequality (b) large-sample confidence interval based on likelihood theory. Coverage is how often you can “capture” the true value and length is \(U-L\). The idea is to extend the R code used for Homework 02. You have to fill up the missing parts accordingly. I have fixed the value of \(\lambda\) and \(J\) for this item. You should explore the properties of a 75% and a 95% interval. Discuss your findings.

# Set a value for lambda to generate artificial datasets following a Poisson distribution
lambda <- 10
# Set c for the desired 1-1/c^2  
c <- 
# Set up negative log-likelihood
# Make sure to point to the data 
mlnl <- function(par, y)
{
  
}
# Create a function to record confidence intervals
# Argument is the sample size J
cons.ci <- function(J)
{
  # Random draws from Poisson
  y <- rpois(J, lambda)
  # Calculate the sample mean (which is really a method of moments estimator of lambda)
  ybar <- mean(y)
  # Calculate the MLE, include hessian = TRUE, and point to the data
  output <- 
  # Finite-sample confidence interval
  ci.finite <- c(2*ybar + c^2/J - sqrt(4*ybar*c^2/J + c^4/J^2), 2*ybar + c^2/J + sqrt(4*ybar*c^2/J + c^4/J^2))/2
  # Likelihood-based confidence interval
  # Make sure to use results from nlm() output
  ci.large <- c( , )
  # Collect four numbers lower and upper limits of the first interval
  # and the lower and upper limits of the second interval
  return(c(ci.finite, ci.large))
}
# Set the sample size for yourself 
J <- 20
# Construct interval 10000 times
results <- replicate(10^4, cons.ci(J))
# Calculate rate of capture

# Calculate the average length of the intervals

Exercise C: Inference for the binomial parameter

This exercise should force you to start working on Sections 5.3 and 6.3 of the book, especially the IID Bernoulli case. The context of this exercise was based on Diez, Barr, and Cetinkaya-Rundel (2014).

Assisted Reproductive Technology (ART) is a collection of techniques that help facilitate pregnancy (e.g. in vitro fertilization). A 2008 report by the Centers for Disease Control and Prevention estimated that ART has been successful in leading to a live birth in 31% of cases. A new fertility clinic claims that their success rate is higher than average. A random sample of 30 of their patients yielded a success rate of 40%. You are to evaluate the claim that this new fertility clinic’s success rate is no different from what was reported in 2008.

We can assume that the we have IID random variables \(Y_1,\ldots, Y_{30}\) are IID \(\mathsf{Ber}\left(p\right)\), where \(p\) is the unknown estimand. The context described above is very similar to Examples 6.2.1 and 6.2.2. The difference from a statistical view is the distribution of the random variables. Your task is to calculate a \(p\)-value in different ways.¹ Let \(\widehat{p}=\overline{Y}\). which is really just a sample proportion. Notice that the book uses \(\widehat{p}\) in two senses: as an estimator and an estimate.

In the current context, the required \(p\)-value is given by \[\mathbb{P}\left(\widehat{p} \geq 0.4 | p = 0.31\right).\]

(2 points) Why do you think this is an appropriate \(p\)-value to calculate in the given context?
(2 points) What is the exact distribution of \(\displaystyle\sum_{i=1}^n Y_i\)? Show that you can simply rewrite the required \(p\)-value as \[\mathbb{P}\left(\widehat{p} \geq 0.4 | p = 0.31\right)=\mathbb{P}\left(\sum_{i=1}^{30} Y_i \geq 12 \bigg| p =0.31\right).\] Find this value using the appropriate commands in R. Be careful here and read the documentation to be sure.
(3 points) Consider R code based on LM Example 6.2.2 with modifications presented below. Fill in the missing parts of the code and use it to calculate a simulated \(p\)-value. Report the code you used and your result.

nsim <- 10^4
# Create artificial data 
art.mat <- replicate(nsim, )
# Calculate the simulated success rate
sim.success <- colMeans(art.mat)
# Obtain a simulated p-value
mean()

(4 points) Yet another way to calculate the required \(p\)-value is to either apply LM Theorem 4.3.1 or to also include a continuity correction (see LM page 240). Find these two values by hand and verify using the appropriate commands in R.
(1 point) Given what you have seen so far, do you think the continuity correction is needed? Why?
(1 point) Based on the \(p\)-values calculated in different ways, do you think there is evidence to support the claim that this new fertility clinic’s success rate is no different from what was reported in 2008? What is your reasoning?

Consider a different approach to evaluating the claim. Someone makes the following suggestion: Since we really do not know \(p\), we need a way to control how often, in repeated samples, we make the mistake of rejecting \(p=0.31\) even if \(p=0.31\) is true. Let us control this rate at 5%. Based on this target rate, we should reject the claim whenever \(\displaystyle\sum_{i=1}^{30} Y_i \geq 14\) (or if you wish \(\widehat{p} \geq 14/30 \approx 0.47\)).

If we adopt this suggestion, we will have to reject the claim because we actually observe a sum equal to 12 (or if you wish we observe a proportion equal to 0.4). You are now to going to explore the properties of the suggestion in repeated samples. In what sense is the suggestion useful?

(2 points) Given that you know the exact distribution of \(\displaystyle\sum_{i=1}^{30} Y_i\), write down (do not calculate!!) an expression for \[\mathbb{P}\left(\sum_{i=1}^{30} Y_i \geq 14\bigg| p=0.31\right).\] Use R to calculate the exact value of this probability.
(2 points) Because the previous probability is the same as \[\mathbb{P}\left(\widehat{p} \geq \frac{14}{30} \bigg| p=0.31\right),\] you could use a Monte Carlo simulation to understand what this probability actually means. Below you will see pretty much the same code as in Item 3. Fill up the missing parts of the R code.

nsim <- 10^4
# Create artificial data again
art.mat <- replicate(nsim, )
# Calculate the simulated total successes
sim.success <- colMeans(art.mat)
# Estimate the desired probability using simulation

(3 points) Let us see if you understood what you found in Item 8. What distribution was used to generate the 10000 artificial datasets? In how many out of the 10000 artificial datasets did you reject the claim? What do you think was achieved by applying the suggestion?

Footnotes

Don’t confuse the \(p\) in \(p\)-value with the estimand \(p\).↩︎