[This article was first published on Analysis with Programming, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

More on Likelihood Ratio Test, the following problem is originally from Casella and Berger (2001), exercise 8.12.

Problem

For samples of size $n=1,4,16,64,100$ from a normal population with mean $mu$ and known variance $sigma^2$, plot the power function of the following LRTs (Likelihood Ratio Tests). Take $alpha = .05$.

1. $H_0:muleq 0$ versus $H_1:mu>0$
2. $H_0:mu=0$ versus $H_1:muneq 0$

Solution

1. The LRT statistic is given by $$ lambda(mathbf{x})=frac{displaystylesup_{muleq 0}mathcal{L}(mu|mathbf{x})}{displaystylesup_{-inftywhile the numerator is evaluated as follows: $$ begin{aligned} sup_{muleq 0}mathcal{L}(mu|mathbf{x})&=sup_{muleq 0}prod_{i=1}^{n}frac{1}{sqrt{2pi}sigma}expleft[-frac{(x_i-mu)^2}{2sigma^2}right]\ &=sup_{muleq 0}frac{1}{(2pisigma^2)^{1/n}}expleft[-displaystylesum_{i=1}^{n}frac{(x_i-mu)^2}{2sigma^2}right]. end{aligned} $$ Above expression will attain its maximum if the value inside the exponential function is small. And for negative values of $muin(-infty,0)$ the quantity $(x_i-mu)^2$ would be large, implies that the exponential term would become small. Therefore, the only value that will give us the supremum likelihood is $mu=mu_0=0$. Hence, $$ begin{aligned} sup_{muleq 0}mathcal{L}(mu|mathbf{x})&=frac{1}{(2pisigma^2)^{1/n}}expleft[-displaystylesum_{i=1}^{n}frac{(x_i-mu_0)^2}{2sigma^2}right]\ =frac{1}{(2pisigma^2)^{1/n}}&expleft[-displaystylesum_{i=1}^{n}frac{(x_i-bar{x}+bar{x}-mu_0)^2}{2sigma^2}right]\ =frac{1}{(2pisigma^2)^{1/n}}&expleft{-displaystylesum_{i=1}^{n}left[frac{(x_i-bar{x})^2+2(x_i-bar{x})(bar{x}-mu_0)+(bar{x}-mu_0)^2}{2sigma^2}right]right}\ =frac{1}{(2pisigma^2)^{1/n}}&expleft[-frac{(n-1)s^2+n(bar{x}-mu_0)^2}{2sigma^2}right], \ &text{since the middle term is 0.}\ =frac{1}{(2pisigma^2)^{1/n}}&expleft[-frac{(n-1)s^2+nbar{x}^2}{2sigma^2}right], text{since }mu_0=0.\ end{aligned} $$ So that $$ begin{equation} label{eq:lrtre} begin{aligned} lambda(mathbf{x})&=frac{frac{1}{(2pisigma^2)^{1/n}}expleft[-frac{(n-1)s^2+nbar{x}^2}{2sigma^2}right]}{frac{1}{(2pisigma^2)^{1/n}}expleft[-frac{(n-1)s^2}{2sigma^2}right]}\ &=expleft[-frac{nbar{x}^2}{2sigma^2}right].\ end{aligned} end{equation} $$ And we reject the null hypothesis if $lambda(mathbf{x})leq c$, that is $$ begin{aligned} expleft[-frac{nbar{x}^2}{2sigma^2}right]&leq c\ -frac{nbar{x}^2}{2sigma^2}&leq log c\ frac{lvertbar{x}rvert}{sigma/sqrt{n}}&geqsqrt{-2log c}=c’. end{aligned} $$

   Figure 1: Plot of Likelihood Ratio Test Statistic for $n = 4,sigma = 1$.

   
   Hence, rejecting the null hypothesis if $lambda(mathbf{x})leq c$, is equivalent to rejecting $H_0$ if $frac{bar{x}}{sigma/sqrt{n}}geq c’in[0,infty)$. Figure 1 depicts the plot of the LRT, the shaded region is on the positive side because that’s where the alternative region is, $H_1:mu>0$, in a sense that if the LRT is small enough to reject $H_0$, then it simply tells us that the plausibility of the parameter in the alternative in explaining the sample is higher compared to the null hypothesis. And if that’s the case, we expect the sample to come from the model proposed by $H_1$, so that the sample mean $bar{x}$, being an unbiased estimator of the population mean $mu$, a function of the LRT statistic, should fall on the side (shaded region) of the alternative.
   
   So that the power function, that is the probability of rejecting the null hypothesis given that it is true (the probability of Type I error) is, $$ begin{aligned} beta(mu)&=mathrm{P}left[frac{bar{x}-mu_0}{sigma/sqrt{n}}geq c’right],quadmu_0=0\ &=1-mathrm{P}left[frac{bar{x}+mu-mu-mu_0}{sigma/sqrt{n}}< c'right]\ &=1-mathrm{P}left[frac{bar{x}-mu}{sigma/sqrt{n}} + frac{mu-mu_0}{sigma/sqrt{n}}< c'right]\ &=1-mathrm{P}left[frac{bar{x}-mu}{sigma/sqrt{n}}< c'+ frac{mu_0-mu}{sigma/sqrt{n}}right]\ &=1-Phileft[c 1="frac{mu_0-mu}{sigma/sqrt{n}}right" language="'+"][/c]. end{aligned} $$ Values taken by $Phi$ are negative and so it decreases, but since we subtracted it to 1, then $beta(mu)$ is an increasing function. So that for $alpha=.05$, $$ begin{aligned} alpha&=sup_{muleq mu_0}beta(mu)\ .05&=beta(mu_0)Rightarrowbeta(mu_0)=1-Phi(c')\ .95&=Phi(c')Rightarrow c'=1.645. end{aligned} $$ Since, $$ begin{aligned} Phi(1.645)=int_{-infty}^{1.645}frac{1}{sqrt{2pi}}expleft[-frac{x^2}{2}right]operatorname{d}x=.9500151. end{aligned} $$ Therefore for $c'=1.645,mu_0=0,sigma=1$, the plot of the power function as a function of $mu$ for different sample size, $n$, is shown in Figure 2. For example, for $n=1$ we compute for the function begin{equation} label{eq:powcomp} begin{aligned} beta(mu)&=1-Phileft[c 1="frac{mu_0-mu}{sigma/sqrt{n}}right" language="'+"][/c]\ &=1-Phileft[1.645+ frac{0-mu}{1/sqrt{1}}right]\ &=1-int_{-infty}^{left(1.645+ frac{0-mu}{1/sqrt{1}}right)}frac{1}{sqrt{2pi}}expleft[-frac{x^2}{2}right]operatorname{d}x. end{aligned} end{equation} The obtained values would be the $y$. For $n = 64$ $$ begin{aligned} beta(mu)&=1-Phileft[c 1="frac{mu_0-mu}{sigma/sqrt{n}}right" language="'+"][/c]\ &=1-Phileft[1.645+ frac{0-mu}{1/sqrt{64}}right]\ &=1-int_{-infty}^{left(1.645+ frac{0-mu}{1/sqrt{64}}right)}frac{1}{sqrt{2pi}}expleft[-frac{x^2}{2}right]operatorname{d}x, end{aligned} $$ and so on.

   Figure 2: Power Function for Different Values of $n$.

   
   
2. The LRT statistic is given by $$ lambda(mathbf{x})=frac{displaystylesup_{mu= 0}mathcal{L}(mu|mathbf{x})}{displaystylesup_{-inftyFigure 3: Plot of Likelihood Ratio Test Statistic for $n = 4,sigma = 1$.
   So that the power function is, $$ begin{aligned} beta(mu)&=mathrm{P}left[frac{lvertbar{x}rvert}{sigma/sqrt{n}}geq c’right]\ &=1 – mathrm{P}left[frac{lvertbar{x}rvert}{sigma/sqrt{n}}< c'right]\ &=1 - mathrm{P}left[-c'Figure 4: Two-Sided Power Function for Different $n$.
   The points in the plot are computed by substituting values of $mu=0,sigma=1$ and $n$ to the power function just like we did in Equation (2).

---

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.