Lecture 27 Coverage and Level by Simulation

We first use the $z$-interval from example 26.3, with the numerical values $n=16$ and $\sigma=2$ (known) from that example, so that the $95\%$ interval is $\bar{X}\pm 1.96\cdot 0.5$. By definition 26.1 of the level $0.95$, the random interval contains the fixed value $\theta$ in $95\%$ of all samples. To verify this, we consider $\theta=10$ and compute the interval for $100$ samples. Before you execute the code, answer the following questions: (1) How many of the $100$ intervals will not cover the value $10$? (2) Will all intervals have the same width? (3) Will the misses be clustered or evenly distributed among the 100 samples?

As in lecture 12, covers is a vector of TRUE and FALSE, one for each interval, and the command sum can be used to count the TRUE values. In our example, $95$ out of $100$ intervals cover the truth, and the misses are for samples $1$, $29$, $34$, $63$ and $83$. The exact number $95$ is a coincidence of the seed: the number of intervals which cover the truth is binomially distributed with standard deviation $2.2$, so $92$ or $98$ would have been just as unremarkable.

Figure 27.1 shows each interval as a horizontal line, with the intervals which cover the truth in grey and the misses in black. The script R/S27-coverage.R shows how the command segments can be used to create this plot.

The intervals in the figure all have the same width, since $\sigma$ is known: the intervals move, while the dashed line stays fixed. The misses are scattered throughout the figure, since the samples are independent. If we had only observed sample $29$, the interval $[10.43,12.39]$ would have looked like any other, and there would have been no indication in the data that this was a miss. For this reason, the statement “$\theta$ is in the interval $[10.43,12.39]$ with probability $0.95$” is wrong: the value $10$ is either in the interval or it is not, and in this instance it is not. The probability refers to the procedure, not to the computed interval. In the script R/S27-coverage.R the experiment is repeated with $10000$ intervals, giving the coverage $0.947$, which is within simulation error of $0.95$. The interval is exact, and the simulation confirms this.

We now turn from an exact statement to an approximate one. Assume that the variance of the normal sample is not known, so that the model has the two parameters $\theta$ and $\sigma^{2}$, and consider the likelihood-ratio test of $H_{0}\colon\theta=\theta_{0}$, where $H_{0}$ leaves $\sigma^{2}$ free. By example 5.4, the log-likelihood is maximised over all parameters at $\hat{\sigma}^{2}=\frac{1}{n}\sum_{i}(X_{i}-\bar{X})^{2}$, and the same computation with $\theta_{0}$ in place of $\bar{X}$ gives the maximum $\hat{\sigma}_{0}^{2}=\frac{1}{n}\sum_{i}(X_{i}-\theta_{0})^{2}$ under $H_{0}$. In both cases the maximised log-likelihood equals $-\frac{n}{2}\log(2\pi\hat{\sigma}^{2})-\frac{n}{2}$, with $\hat{\sigma}_{0}^{2}$ in place of $\hat{\sigma}^{2}$ under $H_{0}$, and thus by (25.1) we find

By Wilks’ theorem (theorem 25.3), with $r=1$ as in lecture 25, the distribution of $W$ under $H_{0}$ converges to $\chi^{2}_{1}$ as $n\to\infty$. In the case of known variance from lecture 25 the chi-squared distribution was exact, but here it is only an approximation, and we need to check how good it is for the $n$ we have. We simulate $W$ under $H_{0}$ for $n=10$ and $n=100$, using $\theta_{0}=10$ and the value $\sigma=2$ from the first section, and we reject if $W$ exceeds the critical value $3.841$ of $\chi^{2}_{1}$, as if the theorem were exact. Before you run the code, predict the two rejection rates. Will the rate for $n=10$ be above or below $0.05$?

The rejection rate for $n=100$ is $0.052$, which is within simulation error of the nominal level. The rejection rate for $n=10$ is $0.075$: the test with nominal level $0.05$ has type I errors half again as often as it should. Figure 27.2 shows why. In the quantile plots of the bottom row, made in the same way as the normal quantile plots of lecture 12, the points for $n=10$ rise above the line from about $2$ onwards, and thus the upper tail of $W$ is heavier than that of the limit.

The script R/S27-wilks.R also computes the $0.95$-quantile of the simulated values of $W$, which is the threshold that does give size $0.05$: for $n=10$ the threshold needs to be about $4.65$ instead of $3.84$, while for $n=100$ the simulated quantile $3.90$ is close to the chi-squared value. The error carries over to intervals: by lemma 26.6, the confidence interval obtained by inverting the test with critical value $3.84$ covers the truth in only about $93\%$ of all samples of size $10$. For this model the exact distribution is known: exercise 25.1 uses the $t$ distribution of appendix A to compare the two tests, and lecture 31 builds the exact test. In general, an asymptotic level is only guaranteed in the limit. An exact test should be used if one is available, and if not, a simulation under $H_{0}$ can be used to determine the actual size at a given $n$.

Repeat the picture in figure 27.1 for the $t$-interval from example 26.4, with $n=5$, sigma replaced by sd(x) and qnorm(0.975) replaced by qt(0.975, df = 4). Estimate the long-run coverage from $10000$ samples. Then repeat the experiment with qnorm(0.975) in place of the $t$ quantile, and comment on your result.

Before the next session, you should type the commands used in this session, and should check that you get the same results as the text.