Project 3

Deadline: April 27, 23:59.

Oral presentations: Thursday, April 24 25, 16:15-18:00 in R60

Project report upload link

Oral presentation slides upload link

Plan for oral presentations

Zoom link to oral presentations

Suppose that \(X_1,X_2,\dots,X_n\) is an iid sample from an exponential distribution with scale parameter \(\beta\).

a) Show that the pivotal quantity \(2\sum_{i=1}^n X_i/\beta\) is chi-square with \(2n\) degrees of freedom. Use this to derive an exact \((1-\alpha)\) confidence interval for \(\beta\).

b) Suppose that we instead were to use parameteric bootstrapping and constructed a bootstrap confidence interval for \(\beta\) using the percentile method (Givens & Hotings, section 9.3.1), that is, using the empirical \(\alpha/2\) and \(1-\alpha/2\) quantiles of the distribution of bootstrap replicates \(\hat\beta^*\) of \(\hat\beta\) as confidence limits, where \(\hat\beta\) is the MLE of \(\beta\). When \(\hat\beta^*\) are based on bootstrap samples from \(F(x;\hat\beta)\), what is the exact distribution of \(\hat\beta^*\)? Find analytic formulas for the resulting confidence limits as functions of \(\hat\beta\) for a given sample size \(n\) assuming that the Monte-Carlo error is negligible (that is, in the limit \(B\rightarrow\infty\)).

c) Find an expression for the exact coverage of the parametric bootstrap percentile interval in point b) in terms of the cdf and quantile function of the chi-square distribution. Compute the exact coverage for \(n=5, 10, 20, 50, 100\) for \(\alpha=0.05\).

The following vector contains a iid sample of \(n=50\) observations from a Cauchy distribution with probability density function \[ f(x) = \frac1{s\pi(1+(\frac{x-\mu}s)^2)} \] where \(\mu\) and \(s\) are unknown location and scale parameters.

 x <- c(12.205, 14.708, 1.415, 9.407, 11.469, 9.339, 9.649, 6.384, 
 5.343, 10.393, 11.511, 11.239, 6.997, 15.248, 8.235, 286.695, 
 7.545, 9.949, 15.053, 8.318, 9.584, 11.573, 6.125, 10.833, 12.229, 
 15.349, 10.084, 15.164, 9.132, 13.647, 45.488, 3.816, 108.554, 
 11.325, 8.795, 6.581, 8.491, 10.705, 7.639, 16.988, 8.739, 5.984, 
 8.376, -1.892, -11.359, 8.442, 10.147, 37.911, 7.764, 7.111)
 

The objective of this exercise is to construct a 90% confidence interval for \(\mu\).

a) Find functions \(\hat\mu\) and \(\hat s\) of the data that can be used as estimators of \(\mu\) and \(s\). The functions don't need to be available in closed form. Implement a function in R that computes \(\hat\mu\) and \(\hat s\) from a vector x containing the data. Note that the likelihood may have multiple optima so maximising the likelihood numerically (e.g. using optim) may require carefully chosen initial parameter values. Alternative estimators can perhaps be obtained by equating empirical and theoretical quantiles.

b) Consider \[ R = \frac{\hat\mu - \mu}{\hat s}. \] For your choice of \(\hat\mu\) and \(\hat s\), prove or disprove that \(R\) is exactly pivotal. Hint: A Cauchy random variable \(X\) with location \(\mu\) and scale \(s\) can be represented as \(X=\mu + s V\) where \(V\) is standard Cauchy. Using parametric bootstrapping from the estimated model, generate \(B=99\) bootstrap replicates \[ R_b^*=\frac{\hat\mu^{b*}-\hat\mu}{\hat s^{b*}}, \] \(b=1,2,\dots,B\) of \(R\) and estimate the upper and lower quantiles of \(R^*\) by the empirical 0.05 and 0.95-quantiles of the bootstrap replicates using the R-function quantile with the type=6 argument. Relying on the assumption that \(R\) is pivotal, derive a 90% confidence interval for \(\mu\).

c) Given that the quantiles are only estimated, is the confidence interval computed using the method in point b) exact or only approximate, that is, does it have the exact nominal coverage of 90% (if we were to repeat the whole procedure)? Hint: When the quantiles are estimated by the type=6 estimator, the above empirical quantiles are equal to the 5th and 95th order statistics of \(R_1^*, R_2^*, ...,R_{99}^*\), that is, \(R_{(5)}^*\) and \(R_{(95)}^*\) (see ?quantile and Hyndman and Fan (1996) for details). Provided that \(R\) is exactly pivotal, what can you say about the probability of the event \(R_{(5)}^* < R < R_{(95)}^*\)?

d) Using Monte-Carlo integration, compute estimates of the expected length of the confidence interval for \(B=19, 39, 99, 199, 399, 999\) assuming that the original data gives an estimate of \(s\) equal to \(\hat s=1\) and an original sample size of \(n=50\) as above. Briefly comment on what is gained by using a higher number of bootstrap samples \(B\). Update April 8: Is the expected length always finite?

Let \(x_1,x_2,\dots,x_n \overset{\text{iid}}{\sim} N(\mu_0, \sigma^2)\) and independent of \(y_1,y_2,\dots,y_n \overset{\text{iid}}{\sim} N(\mu_1, \sigma^2)\). Suppose we don't observe each \(x_i\) and \(y_i\) directly but only \[ z_i =\max(x_i,y_i) \] and \[ u_i = I(x_i<y_i) \] for \(i=1,2,\dots,n\). The aim of this task is to develop an EM-algorithm for computing the MLEs of \(\mu_0,\mu_1,\sigma\) based on the observed data.

a) Write an R function that simulates a realization of \(\mathbf{u},\mathbf{z}\) given \(n\) and the model parameters.

b) Write down the "complete data" log likelihood, \(\ln f(\mathbf{x},\mathbf{y}|\mu_0,\mu_1,\sigma^2)\).

c) Explain why, conditional on \(u_i=0\) and \(z_i\), \(x_i=z_i\) and \(y_i\) has a truncated normal distribution with the reverse applying when \(u_i=1\). Let \(m_i^{(t)}\) and \(v_i^{(t)}\) denote means and variances of these truncated normals (either \(x_i\) or \(y_i\)) and develop formulas for these as functions of the current parameter values \(\mu_0^{(t)},\mu_1^{(t)},\sigma^{2(t)}\) and the observed data \(\mathbf{z},\mathbf{u}\). You may build on the formulas given in wikipedia.

d) Show that \begin{align} Q(\mu_0,\mu_1,\sigma^2|\mu_0^{(t)},\mu_0^{(t)},\sigma^{2(t)}) &=E(\ln f(\mathbf{x},\mathbf{y}|\mu_0,\mu_1,\sigma^2)|\mathbf{u},\mathbf{z},\mu_0^{(t)},\mu_0^{(t)},\sigma^{2(t)}) \\&=-n\ln(2\pi)-n\ln\sigma^2 \\&-\frac1{2\sigma^2}\sum_{i=1}^n \left\{ v_i^{(t)} + (1-u_i)\Big[(z_i-\mu_0)^2 + (m_i^{(t)}-\mu_1)^2\Big] + u_i\Big[(z_i-\mu_1)^2 + (m_i^{(t)}-\mu_0)^2\Big]\right\} \end{align}

e) Find formulas for the values of \(\mu_0,\mu_1,\sigma^2\) that maximises \(Q\).

f) Implement an EM-algorithm based on the above as a function in R and test if the parameter estimates are close to the true values for a simulated data set of size \(n=1000\).

Consider a diploid organism that has two genetic variants (alleles) \(A\) and \(a\) at a given position (locus) of a the genome and two other alleles \(B\) and \(b\) at another nearby linked locus. Each diploid individual of a particular generation is then formed by the union of haploid gametes of four different types \(AB\), \(Ab\), \(aB\) and \(aa\). Since we don't care about which gamete have been inherited maternally and paternally (e.g., the multilocus genotype \(Ab/AB\) is considered the same as \(AB/Ab\)), this results in \(10\) different multilocus genotypes given by \begin{align} AB/AB,AB/Ab,AB/aB,AB/ab \\ Ab/Ab,Ab/aB,Ab,ab \\ aB/aB,aB/ab \\ ab/ab. \tag{*} \end{align}

We will assume that individuals have been formed by random mating among parents and that the probabilities of haploid gamete types produced by the parents are given \(p_{AB},p_{Ab},p_{aB},p_{ab}\). Our aim is to estimate the so-called the amount of linkage disequilibrium between the two loci defined as \(D_{AB}=p_{AB}-p_Ap_B\) where \(p_A\) and \(p_B\) are the probabilities that a gamete carries alleles \(A\) and \(B\) respectively. Linkage disequilibrium \(D_{AB}\) quantifies the degree of non-independence which alleles are present at different loci in a given individual sampled from the population. Such non-independence is of interest in how it is affected by linkage, recombination and various other evolutionary forces such as natural selection and migration.

a) Given the above assumptions, what are the probabilities of each of the 10 multilocus genotypes?

b) Assuming that corresponding counts of these multilocus genotypes in a sample of \(n\) individuals are available and follow a multinomial distribution, show that the "complete data" MLE of \(p_{AB}\) is \[ \hat p_{AB}=\frac1{2n}(2n_{AB/AB} + n_{AB/Ab} + n_{AB/aB} + n_{AB/ab}), \] that is, the MLE of \(p_{AB}\) is equal to the proportion of that gamete type in the sample.

Derive similar formulas for the MLEs of \(p_{Ab}\) and \(p_{aB}\).

How does \(p_A\) and \(p_B\) relate to the gamete type probabilities \(p_{AB},p_{Ab},p_{aB},p_{ab}\)? Use this and functional invariance of MLEs to derive the MLE also of \(D_{AB}\).

c) In practice, with commonly used genetic methods, it is not possible to distinguish between the two types of double heterozygotes \(Ab/aB\) and \(AB/ab\), that is, double heterozygotes formed by the union of gametes of type \(Ab\) and \(aB\) and by union of \(AB\) and \(ab\), respectively. Thus, only \(n_{AaBb}=n_{AB/ab}+n_{Ab/aB}\) can be observed.

Derive an Expectation-Maximisation algorithm for finding the MLE of the model parameters (the gamete probabilities) as well as \(D_{AB}\). In particular, show that the conditional expectation of the complete data log likelihood has the same form as the ordinary multinomial likelihood considered in point b) but with some of the missing counts replaced by their conditional expected values. What is the conditional distribution of \(n_{AB/ab}\) conditional on \(n_{AaBb}\) and given \(\theta^{(t)}\). Use this to find formulas for \(n_{AB/ab}^*=E(n_{AB/ab}|n_{AaBb},\theta^{(t)})\) and \(n_{Ab/aB}^*=E(n_{Ab/aB}|n_{AaBb},\theta^{(t)})\).

Implement the algorithm in R and compute the MLEs of the gamete type probabilities and \(D_{AB}\) when the data (in the same order as the genotypes listed in (*)) for a sample of \(n=500\) individuals is given by

 n <- c(17,22,86,NA,5,NA,38,58,124,51)

and \(n_{AaBb}=99\).