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

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

580 15 Correlation

recsep logwl recsep logwl recsep logwl

2.9370 0.0694 2.9284 0.0433 2.9212 0.0331

2.9196 0.0288 2.9149 0.0221 2.9106 0.0121

2.9047 0.0080 2.9047 0.0069 2.9031 0.0049

2.9320 0.0714 2.9154 0.0427 2.9165 0.0336

2.9085 0.0292 2.9090 0.0240 2.9047 0.0197

2.9058 0.0052 2.9025 0.0104 2.9025 0.0087

2.8976 0.0078 2.8971 0.0065 2.8976 0.0062

(a) Test the hypothesis that the population coefﬁcient of correlation for

the transformed measures is

ρ = 0.96 against the alternative ρ < 0.96. Use

α =0.05.

(b) Find a 95% conﬁdence interval for

ρ.

From the MATLAB code below, we see that r

=0.9246 and the p-value for

the test is 0.0832. Thus, at a 5% signiﬁcance level H

: ρ =0.96 is not rejected.

The 95% conﬁdence interval for

ρ in this case is [0.8203,0.9694].

%nanoprism.m

load ’nanoprism.dat’

recsep = nanoprism(:,1);

logwl = nanoprism(:,2);

r = corr(recsep, logwl) %0.9246

n = length(logwl);

fisherz = @(x) atanh(x); invfisherz = @(x) tanh(x);

%fisher z and inverse transformations as pure functions

%Forming z for testing H0: rho = rho0:

z = (fisherz(r) - fisherz(0.96))/(1/sqrt(n-3)) %-1.3840

pval = normcdf(z) %0.0832

%95% confidence interval

invfisherz([fisherz(r) - norminv(0.975)/sqrt(n-3) ,...

fisherz(r) + norminv(0.975)/sqrt(n-3)])

%0.8203 0.9694

15.2.1.3 Test for the Equality of Two Correlation Coefﬁcients

In some testing scenarios we might be interested to know if the correlation

coefﬁcients from two bivariate populations,

and ρ

, are equal. From the

ﬁrst population the pairs (X

(1)

), i = 1,.. . , n

are observed. Analogously,

from the second population the pairs (X

(2)

), i =1, ..., n

are observed and

sample correlations r

and r

are calculated. The populations are assumed

normal and independent, but components X and Y within each population

might be correlated.

The test statistic for testing H

: ρ

= ρ

versus H

: ρ

>,6=,< ρ

is ex-

pressed in terms of Fisher’s z-transformations of the sample correlations r

and r

15.2 The Pearson Coefﬁcient of Correlation 581

z =

−w

−3

where w

log

1+r

1−r

, i = 1,2 and n

and n

are the number of pairs in the

ﬁrst and second sample, respectively.

Alternative α-level rejection region p-value

: ρ

>ρ

1−α

,∞) 1-normcdf(z)

: ρ

6=ρ

(−∞, z

α/2

] ∪[z

1−α/2

,∞) 2

normcdf(-abs(z))

: ρ

<ρ

(−∞, z

] normcdf(z)

If the H

is not rejected, then one may be interested in pooling the two

sample estimators r

and r

. This is done in the domain of z-transformed

values as

−3)w

+(n

−3)w

−6

and inverting w

to r

via

1 −exp{2w

}

1 +exp{2w

}

Example 15.5. Swallowtail Butterﬂies. The following data were extracted

from a larger study by Brower (1959) on a speciation in a group of swallowtail

butterﬂies (Fig. 15.3). Morphological measurements are in millimeters coded

× 8 (X – length of eighth tergile, Y – length of superuncus).

(a) (b)

Fig. 15.3 (a) Papilio multicaudatus. (b) Papilio rutulus.

582 15 Correlation

Species X Y X Y X Y X Y

Papilio 24 14 21 15 20 17.5 21.5 16.5

multicaudatus 21.5 16 25.5 16 25.5 17.5 28.5 16.5

23.5 15 22 15.5 22.5 17.5 20.5 19

21 13.5 19.5 19 26 18 23 17

21 18 21 17 20.5 16 22.5 15.5

Papilio 20 11.5 21.5 11 18.5 10 20 11

rutulus 19 11 20.5 11 19.5 11 19 10.5

21.5 11 20 11.5 21.5 10 20.5 12

20 10.5 21.5 12.5 17.5 12 21 12.5

21 11.5 21 12 19 10.5 19 11

18 11.5 21.5 10.5 23 11 22.5 11.5

19 13 22.5 14 21 12.5 19.5 12.5

The observed correlation coefﬁcients are r

= −0.1120 (for P. multicauda-

tus) and r

= 0.1757 (for P. rutulus). We are interested if the corresponding

population correlation coefﬁcients

and ρ

are signiﬁcantly different.

The Fisher z-transformations of r

and r

are w

= −0.1125 and w

0.1776. The test statistic is z =

−0.1125−0.1776

1/17+1/25

= −0.9228. For this value of z

the p-value against the two-sided alternative is 0.3561, and the null hypothe-

sis of the equality of population correlations is not rejected. Here is a MATLAB

session for the above exercise.

PapilioM=[24, 14; 21, 15; 20, 17.5; 21.5, 16.5; ...

21.5, 16; 25.5, 16; 25.5, 17.5; 28.5, 16.5; ...

23.5, 15; 22, 15.5; 22.5, 17.5; 20.5, 19; ...

21, 13.5; 19.5, 19; 26, 18; 23, 17; ...

21, 18; 21, 17; 20.5, 16; 22.5, 15.5];

PapilioR=[20, 11.5; 21.5, 11; 18.5, 10; 20, 11; ...

19, 11; 20.5, 11; 19.5, 11; 19, 10.5; ...

21.5, 11; 20, 11.5; 21.5, 10; 20.5, 12; ...

20, 10.5; 21.5, 12.5; 17.5, 12; 21, 12.5; ...

21, 11.5; 21, 12; 19, 10.5; 19, 11; ...

18, 11.5; 21.5, 10.5; 23, 11; 22.5, 11.5; ...

19, 13; 22.5, 14; 21, 12.5; 19.5, 12.5];

PapilioMX=PapilioM(:,1); % X

PapilioMY=PapilioM(:,2); % Y

PapilioRX=PapilioR(:,1); % X

PapilioRY=PapilioR(:,2); % Y

n1=length(PapilioMX);

n2=length(PapilioRX);

r1=corr(PapilioMX, PapilioMY); % -0.1120

r2=corr(PapilioRX, PapilioRY); % 0.1757

%test for rho1 = 0

pval1 = 2

tcdf(-abs(r1

sqrt(n1-3)/sqrt(1-r1^2)), n1-3);

15.2 The Pearson Coefﬁcient of Correlation 583

% 0.6480

%test for rho2 = 0

pval2 = 2

tcdf(-abs(r2

sqrt(n2-3)/sqrt(1-r2^2)), n2-3);

%0.3806

fisherz = @(x) 1/2

log( (1+x)/(1-x) );

%Fisher z transformation as pure function

w1 = fisherz(r1); %-0.1125

w2 = fisherz(r2); %0.1776

%test for rho1 = rho2 vs rho1 ~= rho2

z = (w1 - w2)/sqrt(1/(n1-3) + 1/(n2-3)) %-0.9228

pval = 2

normcdf(-abs(z)) %0.3561

15.2.1.4 Testing the Equality of Several Correlation Coefﬁcients

We are interested in testing H

: ρ

=ρ

=···=ρ

Let r

be the correlation coefﬁcient based on n

pairs, i = 1,..., k, and let

log

1+r

1−r

be its Fisher transformation. Deﬁne

= (n

−3) +(n

−3) +···+(n

−3) =

i=1

−3k, and

w =

−3)w

+(n

−3)w

+···+(n

−3)w

Then the statistic

=(n

−3)w

+(n

−3)w

+···+(n

−3)w

−N(w)

has a χ

-distribution with k −1 degrees of freedom.

Example 15.6. In testing whether the correlations between sepal and petal

lengths differ for the species: setosa, versicolor and virginica, the p-value

found was close to 0.

load fisheriris

%Correlations between sepal and petal lengths

%for species setosa, versicolor and virginica

n1=50; n2=50; n3=50; k=3; N=n1+n2+n3-3

r =[ corr(meas(1:50, 1), meas(1:50, 3)), ...

corr(meas(51:100, 1), meas(51:100, 3)), ...

corr(meas(101:150, 1), meas(101:150, 3)) ]

%0.2672 0.7540 0.8642

fisherz = @(r) 1/2

log( (1+r)./(1-r) );

w=fisherz(r)

%0.2738 0.9823 1.3098

wbar = [n1-3 n2-3 n3-3]/N

w’ %0.8553

584 15 Correlation

chi2 = [n1-3, n2-3, n3-3]

(w.^2)’- N

wbar^2 %26.3581

pval = 1-chi2cdf(chi2, k-1) %1.8897e-006

15.2.1.5 Power and Sample Size in Inference About Correlations

Power and sample size computations for testing correlations use the fact that

Fisher’s

-transformed sample correlation coefﬁcients have an approximately

normal distribution, as we have seen before.

The statistic t

= r

n−2

1−r

, which has a t-distribution with d f = n−2 degrees

of freedom, determines the critical points of the rejection region as

α,n−2

−2,1−α

+(n −2)

where

α is replaced by α/2 for the two-sided alternative, and where the sign of

the square root is negative if H

: ρ <0.

Let z-transformations of r and r

α,n−2

be w and w

, respectively. Then the

power, as a function of w, is

−β =Φ

(w −w

)

n −3

The sample size needed to achieve a power of 1

−β, in a test of level α

against the one-sided alternative, is

1−α

1−β

+3.

Here w

is z-transformation of ρ

from the speciﬁc alternative, H

: ρ = ρ

and z

1−α

and z

1−β

are quantiles of a standard normal distribution. When the

alternative is two-sided, z

1−α

is replaced by z

1−α/2

Example 15.7. One wishes to determine the size of a sample sufﬁcient to reject

: ρ = 0 with a power of 1 −β = 0.90 in a test of level α = 0.05 whenever

ρ ≥ 0.4. Here we take H

: ρ = 0.4 and ﬁnd that a sample of size 51 will be

necessary:

fisherz = @(x) atanh(x);

w1 = fisherz(0.4) % 0.4236

n =((norminv(0.95)+norminv(0.90))/w1)^2 + 3 %50.7152

If H

: ρ =0 is to be rejected whenever |ρ|≥0.4, then a sample of size 62 is

needed:

n =((norminv(0.975)+norminv(0.90))/w1)^2 + 3 %61.5442

15.2 The Pearson Coefﬁcient of Correlation 585

15.2.1.6 Multiple Correlation Coefﬁcient

Consider three variables X

, X

, and Y . The correlation between Y and the

pair X

, X

is measured by the multiple correlation coefﬁcient

y.x

−2r

·r

1 −r

This correlation is signiﬁcant if

y.x

1 −R

y.x

/(n −3)

is large. Here n is the number of triplets (X

, X

,Y ), and statistic F has an

F-distribution with 2, n

−3 degrees of freedom.

The general case R

= R

y.x

...x

will be covered in the multiple regression

section, p. 605; in fact, R

is called the coefﬁcient of determination and testing

its signiﬁcance is equivalent to testing the signiﬁcance of the multiple regres-

sion of Y on X

,... , X

15.2.2 Bayesian Inference for Correlation Coefﬁcients

To conduct a Bayesian inference on a correlation coefﬁcient, a bivariate nor-

mal distribution for the data is assumed and Wishart’s prior is placed on

the inverse of the covariance matrix. Recall that the Wishart distribution is

a multivariate counterpart of a gamma distribution (more precisely, of a

distribution) and the model is in fact a multivariate analog of a gamma prior

on the normal precision parameter.

To illustrate this Bayesian model, a bivariate normal sample of size n

=46

is generated. For this sample Pearson’s coefﬁcient of correlation was found to

be r

=0.9908, with sample variances of s

=8.1088 and s

=35.0266.

WinBUGS code

corr.odc, as given bellow, is run and Bayes estimators

for population correlation

ρ and component variances σ

and σ

are obtained

as 0.9901, 8.114, and 34.99. These values are close to the classical estimators

since the priors are noninformative. The hyperparameters of Wishart’s prior

are matrix

W and degrees of freedom df, and low degrees of freedom, df=3, make

this prior “vague.”

As an exercise, compute a classical 95% conﬁdence interval for

ρ (p. 577)

and compare it with the 95% credible set [0.9823,0.9952]. Are the intervals

similar?

586 15 Correlation

model{

for( i in 1:nn){

y[i,1:2] ~ dmnorm( mu[i,], Tau[,] )

mu[i,1] ~ dnorm(mu.x, tau1)

mu[i,2] ~ dnorm(mu.y, tau2)

}

mu.x ~ dnorm(0, 0.0001)

mu.y ~ dnorm(0, 0.0001)

tau1 ~ dgamma(0.001, 0.001)

tau2 ~ dgamma(0.001, 0.001)

Tau[1:2,1:2] ~ dwish( W[,], df )

df <- 3

Sigma[1:2, 1:2] <- inverse(Tau[,])

rho <- Sigma[1,2]/sqrt(Sigma[1,1]

Sigma[2,2])

}

DATA

list( nn=46, y = structure(.Data = c( 0.5674, -1.6458,

-0.4656, -4.8531,

1.5253 , -2.1200,

...

8.3435, 14.7693,

11.0151, 18.8988,

8.9949, 15.3797,

10.5287, 18.1455), .Dim=c(46,2)) ,

W = structure(.Data = c(1,0,0,1),.Dim=c(2,2) ) )

INIT

list(Tau = structure(.Data = c(1,1,1,1), .Dim = c(2,2)),

mu.x = 1, mu.y = 1, tau1=1, tau2=1)

mean sd MC error val2.5pc median val97.5pc start sample

Sigma[1,1] 8.114 1.757 0.01462 5.361 7.881 12.21 1001 100000

Sigma[1,2] 16.68 3.625 0.03035 11.0 16.2 25.13 1001 100000

Sigma[2,1] 16.68 3.625 0.03035 11.0 16.2 25.13 1001 100000

Sigma[2,2] 34.99 7.570 0.06331 23.15 33.98 52.63 1001 100000

rho 0.9901 0.0033 3.287E-5 0.9823 0.9906 0.9952 1001 100000

15.3 Spearman’s Coefﬁcient of Correlation

Charles Edward Spearman (Fig. 15.4) was a late bloomer, academically speak-

ing. He received his Ph.D. at the age of 48, after having served as an ofﬁcer

in the British army for 15 years. He is most famous in the ﬁeld of psychology,

where he theorized that “general intelligence” was a function of a compre-

hensive mental competence rather than a collection of multifaceted mental

abilities. His theories eventually led to the development of factor analysis.

15.3 Spearman’s Coefﬁcient of Correlation 587

Fig. 15.4 Charles Edward Spearman (1863–1945).

Spearman (1904) proposed the rank correlation coefﬁcient long before

statistics became a scientiﬁc discipline. For bivariate data, an observation

has two coupled components (X,Y ) that may or may not be related to each

other. Let

ρ=Corr(X ,Y ) represent the unknown correlation between two com-

ponents. In a sample of n, let R

,... , R

denote the ranks for the ﬁrst compo-

nent X and S

,... , S

denote the ranks for Y . For example, if x

= x

(n)

is the

largest value from x

,... , x

and y

= y

(1)

is the smallest value from y

,... , y

then (R

, S

) = (n, 1). Corresponding to Pearson’s (parametric) coefﬁcient of

correlation, the Spearman coefﬁcient of correlation is deﬁned as

ρ =

−R)(S

−S)

−R)

−S)

. (15.1)

This expression can be simpliﬁed. From (15.1),

R = S = (n +1)/2 and

−

−S)

= nV ar(R

) = n(n

−1)/12. Deﬁne D as the difference between

ranks, i.e., D

−S

. With R =S, we can see that

=(R

−R) −(S

−S)

and

i=1

−R)

i=1

−S)

−2

i=1

−R)(S

−S),

i.e.,

i=1

−R)(S

−S) =

n(n

−1)

−

i=1

588 15 Correlation

By dividing both sides of the equation by

−R)

−S)

−R)

= n(n

−1)/12, we obtain

ρ =1 −

n(n

−1)

. (15.2)

Consistent with Pearson’s coefﬁcient of correlation (the standard paramet-

ric measure of covariance), Spearman’s coefﬁcient of correlation ranges be-

tween

−1 and 1. If there is perfect agreement, i.e., all the differences are 0,

then

ρ =1. The scenario that maximizes

occurs when ranks are perfectly

opposite: R

= n −S

+1.

If the sample is large enough, then Spearman’s statistic can be approxi-

mated using the normal distribution. It was shown that if n

>10, then

ρ −ρ)

n −1 ∼N (0,1).

Example 15.8. Stichler et al. (1953) list tread wear for tires, each tire mea-

sured by two methods based on (a) weight loss and (b) groove wear.

Weight Groove

45.9 35.7 41.9 39.2

37.5 31.1

33.4 28.1

31.0 24.0

30.5 28.7

30.9 25.9

31.9 23.3

30.4 23.1

27.3 23.7

20.4 20.9

24.5 16.1

20.9 19.9

18.9 15.2

13.7 11.5

11.4 11.2

For this data,

ρ = 0.9265. Note that if we opt for the parametric measure

of correlation, the Pearson coefﬁcient is 0.948.

Ties in the Data: The statistics in (15.1) and (15.2) are not designed for

paired data that include tied measurements. If ties exist in the data, a simple

adjustment should be made. Deﬁne u

u(u

−1)/12 and v

v(v

−1)/12

where the us and vs are the ranks for X and Y adjusted (e.g., averaged) for

ties. Then

n(n

−1) −6

−6(u

)

{[n(n

−1) −12u

][n(n

−1) −12v

]}

1/2

and it holds that, for large n,

−ρ

)

n −1 ∼N (0,1).

15.4 Kendall’s Tau 589

The MATLAB function corr(x,y,’type’,’Spearman’) computes the Spear-

man correlation coefﬁcient for column vectors x and y.

15.4 Kendall’s Tau

M. G. Kendall (Fig. 15.5) formalized an alternative measure of dependence

(originally proposed and used in the nineteenth century) by ﬁnding out how

many pairs in a bivariate sample are “concordant,” which means that the signs

between X and Y agree in the pairs. Pairs for which one sign is plus and the

other is minus are “discordant.” From (X

), i = 1,.. . , n one can choose

different pairs. The pair (X

),(X

) is concordant if either X

≤ X

and

≤Y

or X

≥ X

and Y

≥Y

. The pair is called discordant if either X

≤ X

and Y

≥Y

or X

≥ X

and Y

≤Y

. For example, the pairs (2,4) and (1,−1) are

concordant, while the pairs (

−2,4) and (1,−1) are discordant.

Fig. 15.5 Sir Maurice George Kendall (1907–1983).

Kendall’s

τ-statistic (Kendall, 1938) is deﬁned as

τ =

n(n −1)

, S

i=1

j=i+1

sign{r

−r

where r

s are deﬁned via ranks of the second sample corresponding to the

ordered ranks of the ﬁrst sample,

{1,2,.. . , n}, i.e.,

1 2 . . . n

... r

In this notation

from Spearman’s coefﬁcient of correlation becomes

−i)

. In terms of the number of concordant (n

) and discordant (n

n −n

) pairs,

τ =

2(n

−n

)

n(n −1)