Vidakovic B. Statistics for Bioengineering Sciences: With Matlab and WinBugs Support

Подождите немного. Документ загружается.

12.9 Exercises 499

are being recorded on videotape. Subsequently three clinical experts, un-

involved in the experimental treatment and not knowing which subject re-

ceived which treatment, independently view the videotapes and rate each

subject according to the degree of claustrophobic tendency shown in the test

situations. Each subject is rated by the experts on a scale of 0 (no claustro-

phobia) to 10 (an extreme claustrophobia). The following table shows the

average ratings for each subject in each of the two groups.

4.6 4.7 4.9 5.1 7.0 4.9

5.1 5.2 5.5 4.8 5.7 5.0

5.8 6.1 6.5 7.0 6.4 5.2

4.6 4.7 4.9 6.4 5.9 4.7

5.8 5.2 5.4 6.1 7.7 6.2

5.8 5.1 6.5 2.2 6.9 5.0

6.5 7.2 8.2 6.7

5.2 5.3 5.4 7.7 8.1 4.9

5.6 6.2 6.3 7.0 7.0 7.8

6.8 7.7 8.0 6.6 5.5 8.2

8.1 5.0

The investigators expected treatment A to prove more effective, and sure

enough it is group A that appears to show the lower mean level of claustro-

phobic tendency in the test situations.

Using Wilcoxon’s sum rank test, test the hypothesis H

that treatments A

and B are of the same effectiveness, versus the alternative that treatment A

is more effective.

12.10. Nonparametric Stats with Raynaud’s Phenomenon. (a) Refer to Ray-

naud’s phenomenon data from Exercise 11.21. Using Wilcoxon’s signed-

rank test, compare the responses (number of attacks) for drug/placebo ef-

fect. Ignore the

Period variable. Compare this result with the result from a

paired t-test.

(b) Compare the drug/placebo effect using Bayesian inference. Write a Win-

BUGS program that will read the drug/placebo info, take the difference

between the measurements,

d, and model it as a normal, with noninforma-

tive priors on the mean

mu and precision prec. Check if the credible set for

mu contains 0 and draw the appropriate conclusion.

12.11. Cotinine and Nicotine Oxide. Garrod et al. (1974) measured the nico-

tine metabolites, cotinine, and nicotine-1’-N-oxide in 24-hour urine collec-

tions from normal healthy male smokers and smokers affected with cancer

of the urinary bladder. The data, also discussed by Daniel (1978), show the

ratio of cotinine to nicotine-1’-N-oxide in the two groups of subjects.

Ratio

Cancer 5.0 8.3 6.7 3.0 2.5 12.5 2.4 5.5 5.2 21.3 5.1 1.6

patients

2.1 4.6 3.2 2.2 7.0 3.3 6.7 11.1 3.4 5.9 27.4

Control 2.3 1.9 3.6 2.5 0.75 2.5 2.1 1.1 2.3 2.2 3.5 1.8

subjects

2.3 1.4 2.1 2.0 2.3 2.4 3.6 2.6 1.5

Test the difference between the two populations using WSuRT, at α =0.05.

500 12 Distribution-Free Tests

12.12. Coagulation Times. In Example 11.1 a standard one-way ANOVA gave p-

value of 4.6585e-05. Repeat the test of equality of means using the Kruskal–

Wallis procedure. If the test turns out signiﬁcant, compare the means using

the Kruskal–Wallis pairwise comparisons.

12.13. Blocking by Rats. In Example 11.5 factor

Procedure was found signiﬁ-

cant with p-value of 8.7127e-04. Is this factor signiﬁcant according to Fried-

man’s test?

MATLAB FILES USED IN THIS CHAPTER

http://springer.bme.gatech.edu/Ch12.NP/

claustrophobia.m, friedmanGT.m, friedmanpairwise.m, grubbs.m, kw.m,

miammt.m, NPdemo.m, npexamples.m, ranks.m, signtst.m, walshnp.m,

wsirt.m, wsurt.m

CHAPTER REFERENCES

Bradley, J. V. (1968). Distribution Free Statistical Tests. Prentice Hall, Englewood Cliffs.

Conover, W. J. (1999). Practical Nonparametric Statistics, 3rd edn. Wiley, New York.

Daniel, W. W. (1978). Applied Nonparametric Statistics, Houghton Mifﬂin, Boston.

Friedman, M. (1937). The use of ranks to avoid the assumption of normality implicit in the

analysis of variance. J. Am. Stat. Assoc., 32, 675–701.

Garrod, J. W., Jenner, P., Keysell, G. R., Mikhael, B. R. (1974). Oxidative metabolism of

nicotine by cigarette smokers with cancer of urinary bladder. J. Nat. cancer Inst., 52,

1421–1924.

Iman, R. L. and Davenport, J. M. (1980). Approximations of the critical region of the Fried-

man statistic. Commun. Stat. A., 9, 571–595.

Kruskal, W. H. and Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. J.

Am. Stat. Assoc., 47, 583–621.

Kvam, P. and Vidakovic, B. (2007). Nonparametric Statistics with Applications to Science

and Engineering. Wiley, Hoboken.

Madansky, A. (1988). Prescriptions for Working Statisticians. Springer, Berlin Heidelberg

New York.

Scanlon, T. J., Luben, R. N., Scanlon, F. L., and Singleton, N. (1993). Is Friday the 13th bad

for your health? Br. Med. J., 307, 1584–1586.

Vagenakis, A. G., Rapoport, B., Azizi, F., Portnay, G. I., Braverman, L. E., and Ingbar, S.

H. (1974). Hyperresponse to thyrothropin-releasing hormone acompanying small de-

creases in serum thyroid concentrations. J. Clin. Invest., 54, 913–918.

Chapter References 501

Walsh, J. E. (1962). Handbook of Nonparametric Statistics I and II. Van Nostrand, Prince-

ton.

Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics, 1, 80–83.

Chapter 13

Goodness-of-Fit Tests

When schemes are laid in advance, it is surprising how often the circumstances ﬁt in

with them.

– William Osler, Canadian Physician (1849–1919)

WHAT IS COVERED IN THIS CHAPTER

• Quantile–Quantile Plots

• Pearson’s χ

-Test

• Kolmogorov–Smirnov Goodness-of-Fit Test

• Smirnov’s Two-Sample Test

• Moran’s Test

• Testing Departures from Normality

13.1 Introduction

Goodness-of-ﬁt tests are batteries of tests that test that the distribution of a

sample is equal to some ﬁxed-in-advance distribution. We already saw Q–Q

plots in Chap. 5 where the samples were compared to some theoretical distri-

butions but in a descriptive fashion, without formal inference. In this chapter

we discuss the celebrated Pearson’s

-test and the Kolmogorov–Smirnov (KS)

test.

B. Vidakovic, Statistics for Bioengineering Sciences: With MATLAB and WinBUGS Support, 503

Springer Texts in Statistics, DOI 10.1007/978-1-4614-0394-4_13,

504 13 Goodness-of-Fit Tests

13.2 Quantile–Quantile Plots

Quantile–quantile, or Q–Q, plots are a popular and informal diagnostic tool

for assessing the distribution of data and comparing two distributions. There

are two kinds of Q–Q plot, one that compares an empirical distribution with a

theoretical one and another that compares two empirical distributions.

We will explain the plots with an example. Suppose that we generated

some random sample X from a uniform

U (−10, 10) distribution of size n =200.

Suppose for a moment that we do not know the distribution of this data and

want to check for normality by a Q–Q plot. We generate quantiles of a normal

distribution for n equally spaced probabilities between 0 and 1, starting at

0.5/n and ending at 1

−0.5/n. The empirical quantiles for X are simply the ele-

ments of a sorted sample. If the empirical distribution matches the theoretical,

the Q–Q plot is close to a straight line.

n=200;

X = unifrnd(-10, 10, [1 n]);

q=((0.5/n):1/n:(1-0.5/n));

qY = norminv(q);

qZ = unifinv(q);

qX = sort(X);

figure(1); plot(qX, qY, ’

’)

figure(2); plot(qX, qZ, ’

’)

The Q–Q plot is given in Fig.13.1a. Note that the ends of the Q–Q plot

curve up and down from the straight line, indicating that the sample is not

normal.

−10 −5 0 5 10

−3

−2

−1

Empirical Quantiles

Standard Normal Quantiles

−10 −5 0 5 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Empirical Quantiles

Uniform (0,1) Quantiles

(a) (b)

Fig. 13.1 (a) Q–Q plot of a uniform sample against normal quantiles. (b) Q–Q plot of a

uniform sample against uniform quantiles.

13.2 Quantile–Quantile Plots 505

then the Q–Q plot indicates that the match is good and that the sample is uni-

form (Fig. 13.1b). Note that while the moments (mean and variance) of the

sample and theoretical distributions may differ, matching the family of distri-

butions is important. If the mean and variance differ, the straight line in the

plot is shifted and has a different slope. Sometimes the straight line passing

through the ﬁrst and third quartiles is added to the Q–Q plot for reference.

MATLAB has built-in functions

qqplot, normplot, probplot, and cdfplot for

a variety of visualization and diagnostic tasks.



In our example we used p

-quantiles for p

=(i −0.5)/n, i =1,..., n. In gen-

eral, one can use p

= (i −c)/(n −2c +1) for any c ∈[0, 1], and popular choices

are c

= 1/3, c = 1/2 (our choice), and c = 1. The choice c = 0.3175 is recom-

mended if the distributions are long-tailed since the extreme points on a Q–Q

plot are more stable. This choice of c corresponds to medians of order statistics

from a uniform distribution on (0,1) and has invariance property, F

−1

) are

medians of order statistics from F.

Multivariate extensions are possible. A simple procedure that can be used

to assess multivariate normality is described next.



Given multivariate data y

=(y

, y

,... , y

i p

), i =1, . .., n one ﬁnds

=(y

− y

) S

−1

− y

)

, (13.1)

where

is the sample mean and S =

n−1

− y

)

− y

)

is the sample covariance matrix. If the y

s are distributed as multivariate

normal, then the d

s have an approximately χ

-distribution, and plotting the

empirical quantiles of d

against the corresponding quantiles of χ

assesses

the goodness of a multivariate normal ﬁt.

Example 13.1. The data set

pearson.dat compiled by K. Pearson contains

the hereditary stature of 1078 fathers and their sons. This data set will be

revisited in Chap. 16 in the context of regression. We will visualize the agree-

ment of pairs of heights to a bivariate normal distribution by inspecting a Q–Q

plot. We will calculate and sort d

s in (13.1) and plot them against quantiles

of a

-distribution (Fig. 13.2).

%qqbivariatenormal.m

load ’pearson.dat’

S = cov(pearson)

[m n] = size(pearson) %[1078 2]

mp = repmat(mean(pearson), m, 1);

d2 = diag((pearson - mp)

inv(S)

(pearson - mp)’);

% chi-square (2 df) quantiles

p = 0.5/m:1/m:1;

quanch2 = chi2inv(p,2);

If, on the other hand, the statement qY=norminv(q) is replaced by qY=unifinv(q),

506 13 Goodness-of-Fit Tests

figure(1)

loglog(quanch2, sort(d2),’o’,’Markersize’,msize,...

’MarkerEdgeColor’,’k’, ’MarkerFaceColor’,’g’)

hold on

loglog([0.0002 50],[0.0002 50],’r’,’linewidth’,lw)

axis([10^-4, 100, 10^-4, 100])

xlabel(’Quantiles of chi^2

2 ’)

ylabel(’Sorted d^2’)

Fig. 13.2 Sorted values of d

plotted against quantiles of a χ

-distribution. The linearity of

the plot indicates a good ﬁt of Pearson’s data to a bivariate normal distribution.

Example 13.2. Poissonness Plots. Poissonness plots were developed by Hoaglin

(1980) and consist of assessing the goodness of ﬁt of observed frequencies to

frequencies corresponding to Poisson-distributed data. Suppose that frequen-

cies n

, i =1,..., k corresponding to realizations 0,1,..., k −1 are observed and

that for i

≥ k all n

=0. Let N =

Then, the frequency n

and the theoretical counterpart e

= N ×

exp{−λ}

should be close, if the data are consistent with a Poisson distribution,

∼ N ×

exp

{−λ}.

By taking the logarithms of both sides one gets

log n

+log(i!) −log N ∼logλ ·i −λ.

Thus, the plot of log n

+log(i!) −log N against integer i should be linear with

slope log

λ if the frequencies are consistent with a Poisson P oi(λ) distribution.

Such a plot is called a poissonness plot.

13.2 Quantile–Quantile Plots 507

The script poissonness.m illustrates this. N = 50,000 Poisson random

variables with

λ =4 are generated, observed frequencies are calculated, and a

Poisonness plot is made (Fig.13.3).

rand(’state’,2)

N = 50000; xx = poissrnd(4, 1, N); %simulate N Poisson(4) rv’s.

% observed frequencies

ni = [sum(xx==0) sum(xx==1) sum(xx==2) sum(xx==3) sum(xx==4) ...

sum(xx==5) sum(xx==6) sum(xx==7) sum(xx==8) sum(xx==9) ...

sum(xx==10) sum(xx==11) sum(xx==12) sum(xx==13) sum(xx==14)]’ ;

i=(0:14)’

poissonness = log(ni) + log(factorial(i)) - log(N);

plot(i, poissonness,’o’)

xlabel(’i’)

ylabel(’log(n

i) + log(i!) - log(N)’)

hold on

plot(i, log(4)

i - 4,’r:’) %theoretical line with lambda=4

[beta] = regress(poissonness, [ones(size(i)) i])

% beta =

% -4.0591

% 1.4006

% Slope in the linear fit is 1.4006

% lambda can be estimated as exp(1.4006) = 4.0576,

% also, as negative intercept, 4.0509

0 2 4 6 8 10 12 14

−4

−2

log(n

) + log(i!) − log(N)

Fig. 13.3 Poissonness plot for simulated example involving N =50,000 Poisson P oi(4) ran-

dom variates. The slope of this plot estimates log

λ while the intercept estimates −λ. The

theoretical line, log(4)i

−4 (red), ﬁts the points (i,log(n

) +log(i!)−log(N)) well.

508 13 Goodness-of-Fit Tests

13.3 Pearson’s Chi-Square Test

Pearson’s χ

-test (Pearson, 1900) is the ﬁrst formally formulated testing proce-

dure in the history of statistical inference. In his 1984 Science article entitled

“Trial by number,” Ian Hacking says that the goodness-of-ﬁt chi-square test

introduced by Karl Pearson (Fig. 13.4) ushered in a new kind of decision mak-

ing and gives it a place among the top 20 discoveries since 1900 considering

all branches of science and technology.

Fig. 13.4 Karl Pearson (1857–1936).

Suppose that X

, X

,... , X

is a sample from an unknown distribution. We

are interested in testing that this distribution is equal to some speciﬁc distri-

bution F

, i.e., in testing the goodness-of-ﬁt hypothesis H

: F

(x) = F

(x).

Suppose that the domain (a, b) of the distribution F

is split into r nonover-

lapping intervals, I

=(a, x

], I

=(x

, x

] ... I

=(x

r−1

, b). Such intervals have

probabilities p

)−F

(a), p

)−F

), ..., p

(b)−F

r−1

under the assumption that H

is true. Of course it is understood that the obser-

vations belong to the domain of F

; if this is not the case, the null hypothesis

is automatically rejected.

Let n

, n

,... , n

be the observed frequencies of observations falling in the

intervals I

, I

,... , I

. In this notation, n

is the count of observations from the

sample X

,... , X

that fall in the interval I

. Of course, n

+···+n

= n since

the intervals partition the domain of the sample.

The discrepancy between observed frequencies n

and frequencies under

, np

is the rationale for forming the statistic

i=1

−n p

)

, (13.2)

which has a

-distribution with r −1 degrees of freedom. Alternative repre-

sentations include

13.3 Pearson’s Chi-Square Test 509

i=1

−n and χ

= n

i=1

−1

where

= n

/n.

In some cases, the distribution under H

is not fully speciﬁed; for example,

one might conjecture that the data is exponentially distributed without know-

ing the exact value of

λ. In this case, the unknown parameter can be estimated

from the sample.

If k parameters needed to fully specify F

are estimated from the sample,

the

-statistic in (13.2) has a χ

-distribution with r−k−1 degrees of freedom.

A degree of freedom in the test statistic is lost for each estimated parameter.

Firm recommendations on how to select the intervals or even the number of

intervals for a

-test do not exist. For a continuous distribution function one

may make the intervals approximately equal in probability. Practitioners may

want to arrange interval selection so that all np

>1 and that at least 80% of

the np

s exceed 5. The rule of thumb is n ≥10, r ≥3, n

/r ≥10, np

≥0.25.

Some statisticians refer to the

-test as Pearson’s “poorness-of-ﬁt” test.

This attribute is appropriate since the test can only conclude that the model

does not ﬁt the data well. If the model ﬁts the data well, then it may not be

unique and one could possibly ﬁnd other models that would pass the test as

well.

When some of the expected frequencies are less than 5, or when the number

of classes is small, then the

-test with Yates (Yates, 1934) corrections is

recommended,

i=1

(|n

−n p

|−0.5)

which, under H

, has a χ

-distribution with r −1 degrees of freedom.

Example 13.3. Weldon’s 26,306 Rolls of 12 Dice. Pearson (1900) discusses

Weldon’s data (table below) as the main illustration for his test. Raphael Wel-

don, an evolutionary biologist and founder of biometry, rolled 12 dice simulta-

neously and recorded the number of times a 5 or a 6 was rolled. In his letter to

Galton, dated 2 February 1894, Weldon provided the data and asked for Gal-

ton’s opinion about their validity. Three contemporary British statisticians,

Pearson, Edgeworth, and Yule, have also considered Weldon’s data. An inter-

esting historic account can be found in Kemp and Kemp (1991).

The results from 26,306 rolls of 12 dice are summarized in the following

table: