Goodness-of-Fit for Composite Hypotheses

Section 11 Goodness-of-fit for composite hypotheses. Example. Let us consider a Matlab example. Let us generate 50 observations from N (1, 2): X=normrnd(1,2,50,1); Then, running a chi-squared goodness-of-fit test ’chi2gof’ [H,P,STATS]= chi2gof(X) outputs H = 0, P = 0.8793, STATS = chi2stat: 0.6742 df: 3 edges: [-3.7292 -0.9249 0.0099 0.9447 1.8795 2.8142 5.6186] O: [8 7 8 8 9 10] E: [8.7743 7.0639 8.7464 8.8284 7.2645 9.3226] The test accepts the hypothesis that the data is normal. Notice, however, that something is different. Matlab grouped the data into 6 intervals, so chi-squared test from previous lecture should have r − 1 = 6 − 1 = 5 degrees of freedom, but we have ’df: 3’ ! The difference is that now our hypothesis is not that the data comes from a particular given distribution but that the data comes from a family of distributions which is called a composite hypothesis. Running [H,P,STATS]= chi2gof(X,’cdf’,@(z)normcdf(z,mean(X),std(X,1))) would test a simple hypothesis that the data comes from a particular normal distribution N (µ̂, χ̂ 2 ) and the output H = 0, P = 0.9838 STATS = chi2stat: 0.6842 df: 5 edges: [-3.7292 -0.9249 0.0099 0.9447 1.8795 2.8142 5.6186] O: [8 7 8 8 9 10] E: [8.6525 7.0995 8.8282 8.9127 7.3053 9.2017] has ’df: 5.’ However, we can not use this test because we estimate the parameters µ̂ and χ̂ 2 of this distribution using the data so this is not a particular given distribution; in fact, this is the distribution that fits the data the best, so the T statistic in Pearson’s theorem will behave differently. Let us start with a discrete case when a random variable takes a finite number of values B1 , . . . , Br with probabilities p1 = P(X = B1 ), . . . , pr = P(X = Br ). We would like to test a hypothesis that this distribution comes from a family of distributions {P� : ν → �}. In other words, if we denote pj (ν) = P� (X = Bj ), we want to test H0 : pj = pj (ν) for all j � r for some ν → � H1 : otherwise. If we wanted to test H0 for one particular fixed ν we could use the statistic r � (�j − npj (ν))2 T = , npj (ν) j=1 and use a simple chi-squared goodness-of-fit test. The situation now is more complicated because we want to test if pj = pj (ν), j � r at least for some ν → � which means that we have many candidates for ν. One way to approach this problem is as follows. (Step 1) Assuming that hypothesis H0 holds, i.e. P = P� for some ν → �, we can find an estimate ν � of this unknown ν and then (Step 2) try to test if, indeed, the distribution P is equal to P�� by using the statistics r � (�j − npj (ν � ))2 T = npj (ν � ) j=1 in chi-squared goodness-of-fit test. This approach looks natural, the only question is what estimate ν � to use and how the fact that ν � also depends on the data will affect the convergence of T. It turns out that if we let ν � be the maximum likelihood estimate, i.e. ν that maximizes the likelihood function �(ν) = p1 (ν)�1 . . . pr (ν)�r then the statistic T = r � (�j − npj (ν � ))2 j=1 npj (ν � ) �d ϕ2r−s−1 (11.0.1) converges to ϕ2r−s−1 distribution with r − s − 1 degrees of freedom, where s is the dimension of the parameter set �. Of course, here we assume that s � r − 2 so that we have at least one degree of freedom. Very informally, by dimension we understand the number of free parameters that describe the set � � (p1 (ν), . . . , pr (ν)) : ν → � . Then the decision rule will be α= � H1 : T � c H2 : T > c where the threshold c is determined from the condition P(α = ≤ H0 |H0 ) = P(T > c|H0 ) � ϕ2r−s−1 (c, +≈) = � where � → [0, 1] is the level of sidnificance. Example 1. Suppose that a gene has two possible alleles A1 and A2 and the combina­ tions of these alleles define three genotypes A1 A1 , A1 A2 and A2 A2 . We want to test a theory that Probability to pass A1 to a child = ν Probability to pass A2 to a child = 1 − ν and that the probabilities of genotypes are given by p1 (ν) = P(A1 A1 ) = ν 2 p2 (ν) = P(A1 A2 ) = 2ν(1 − ν) p3 (ν) = P(A2 A2 ) = (1 − ν)2 . (11.0.2) Suppose that given a random sample X1 , . . . , Xn from the population the counts of each genotype are �1 , �2 and �3 . To test the theory we want to test the hypothesis H0 : p1 = p1 (ν), p2 = p2 (ν), p3 = p3 (ν) for some ν → [0, 1] H1 : otherwise. First of all, the dimension of the parameter set is s = 1 since the distributions are determined by one parameter ν. To find the MLE ν � we have to maximize the likelihood function p1 (ν)�1 p2 (ν)�2 p3 (ν)�3 or, equivalently, maximize the log-likelihood log p1 (ν)�1 p2 (ν)�2 p3 (ν)�3 = �1 log p1 (ν) + �2 log p2 (ν) + �3 log p3 (ν) = �1 log ν 2 + �2 log 2ν(1 − ν) + �3 log(1 − ν)2 . If we compute the critical point by setting the derivative equal to 0, we get ν� = 2�1 + �2 . 2n Therefore, under the null hypothesis H0 the statistic T = �d (�1 − np1 (ν � ))2 (�2 − np2 (ν � ))2 (�3 − np3 (ν � ))2 + + np1 (ν � ) np2 (ν � ) np3 (ν � ) ϕ2r−s−1 = ϕ23−1−1 = ϕ12 converges to ϕ21 -distribution with one degree of freedom. Therefore, in the decision rule � H1 : T � c α= H2 : T > c threshold c is determined by the condition ≤ H0 |H0 ) � ϕ21 (T > c) = �. P(α = For example, if � = 0.05 then c = 3.841. Example 2. A blood type O, A, B, AB is determined by a combination of two alleles out of A, B, O and allele O is dominated by A and B. Suppose that p, q and r = 1 − p − q are the population frequencies of alleles A, B and O correspondingly. If alleles are passed randomly from the parents then the probabilities of blood types will be Blood type Allele combinations Probabilities Counts O OO r2 �1 = 121 2 A AA, AO p + 2pr �2 = 120 �3 = 79 B BB, BO q 2 + 2pr AB AB 2pq �4 = 33 We would like to test this theory based on the counts of each blood type in a random sample of 353 people. We have four groups and two free parameters p and q, so the chi-squared statistics T under the null hypotheses will have ϕ24−2−1 = ϕ21 distribution with one degree of freedom. First, we have to find the MLE of parameters p and q. The log likelihood is �1 log r 2 + �2 log(p2 + 2pr) + �3 log(q 2 + 2qr) + �4 log(2pq) = 2�1 log(1 − p − q) + �2 log(2p − p2 − 2pq) + �3 log(2q − q 2 − 2pq) + �4 log(2pq). Unfortunately, if we set the derivatives with respect to p and q equal to zero, we get a system of two equations that is hard to solve explicitly. So instead we can minimize log likelihood numerically to get the MLE p̂ = 0.247 and q̂ = 0.173. Plugging these into formulas of blood type probabilities we get the estimated probabilities and estimated counts in each group p̂i np̂i O A B AB 0.3364 0.3475 0.2306 0.0855 118.7492 122.6777 81.4050 30.1681 We can now compute chi-squared statistic T � 0.44 and the p-value ϕ21 (T, ≈) = 0.5071. The data agrees very well with the above theory. We could also use a similar test when the distributions P� , ν → � are not necessarily supported by a finite number of points B1 , . . . , Br , for example, continuous distributions. In this case if we want to test the hypothesis H0 : P = P� for some ν → � we can group the data into r intervals I1 , . . . , Ir and test the hypothesis H0 : pj = pj (ν) = P� (X → Ij ) for all j � r for some ν. For example, if we discretize normal distribution by grouping the data into intervals I1 , . . . , Ir then the hypothesis will be H0� : pj = N(µ, χ 2 )(Ij ) for all j � r for some (�, χ 2 ). There are two free parameters µ and χ 2 that describe all these probabilities so in this case s = 2. Matlab function ’chi2gof’ tests for normality by grouping the data and computing statistic T in (11.0.1) - that is why it uses ϕ2r−s−1 distribution with r−s−1=r−2−1=r−3 degrees of freedom and, thus, ’df: 3’ in the example above. Example. Let us test if the data ’normtemp’ from normal body temperature dataset fits normal distribution. [H,P,STATS]= chi2gof(normtemp) gives H = 0, P = 0.0504 STATS = chi2stat: 9.4682 df: 4 edges: [1x8 double] O: [13 12 29 27 35 10 4] E: [9.9068 16.9874 27.6222 31.1769 24.4270 13.2839 6.5958] and we accept null hypothesis at the default level of significance � = 0.05 since p-value 0.0504 > � = 0.05. We have r = 7 groups and, therefore, r − s − 1 = 7 − 2 − 1 = 4 degrees of freedom. In the case when the distributions P� are continuous or, more generally, have infinite number of values that must be grouped in order to use chi-squared test (for example, normal or Poisson distribution), it can be a difficult numerical problem to maximize the “grouped“ likelihood function P� (I1 )�1 · . . . · P� (Ir )�r � max � ν � . � It is tempting to use a usual non-grouped MLE ν̂ of ν instead of the above ν � because it is often easier to compute, in fact, for many distributions we know explicit formulas for these MLEs. However, if we use ν̂ in the statistic T = r � (�j − npj (ν̂))2 j=1 npj (ν̂) (11.0.3) then it will no longer converge to ϕ2r−s−1 distribution. A famous result in [1] proves that typically this T will converge to a distribution ”in between” ϕ2r−s−1 and ϕ2r−1 . Intuitively this is easy to understand because ν � specifically fits the grouped data �1 , . . . , �r so the expected counts np1 (ν � ), . . . , npr (ν � ) should be a better fit compared to the expected counts np1 (ν̂), . . . , npr (ν̂). On the other hand, these last expected counts should be a better fit than simply using the true expected counts np1 (ν0 ), . . . , npr (ν0 ) since the MLE ν̂ fits the data better than the true distribution. So typically we would expect r r r � (�j − npj (ν � ))2 � (�j − npj (ν̂))2 � (�j − npj (ν0 ))2 � . � npj (ν � ) npj (ν0 ) npj (ν̂) j=1 j=1 j=1 But the left hand side converges to ϕ2r−s−1 and the right hand side converges to ϕ2r−1 . Thus, if the decision rule is based on the statistic (11.0.3): � H1 : T � c α= H2 : T > c then the threshold c can be determined conservatively from the tail of ϕ2r−1 distribution since P(α = ≤ H0 |H0 ) = P(T > c) � ϕ2r−1 (T > c) = �. References: [1] Chernoff, Herman; Lehmann, E. L. (1954) The use of maximum likelihood estimates in ϕ2 tests for goodness of fit. Ann. Math. Statistics 25, pp. 579-586.