Coursehero >> Utah >> BYU >> LAW 666

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

Multivariate III. Statistical Inference Why use a multivariate approach when conducting tests on p variables? 1. Type I error protection EX p = 10 univariate tests at = .05 If variables are independent, Pr{at least one rejection} = 1 Pr{all 10 tests accept } = 1 (.95)10 .40 = In practice, the overall ( experimentwise ) Type I error rate will fall in what range?? 2. Power Multivariate test is more powerful in my cases. EX All p univariate tests fail to reject, but multivariate test is signi cant due to combination of small e ects on some variables. 3. Understanding variables acting in combination. 1 A.i. Hotelling s T 2 First consider univariate test of H0 : = 0 vs. H1 : = 0 when is known. (Consider only two-sided tests, since one-sided don t readily generalize for p > 1) Test statistic using r.s. (x1 , . . . , xn ): z= or x 0 N (0, 1) under H0 / n z2 = square of standardized distance x 0 2 =n 2 under H0 1 2 Multivariate generalization ( is known): 01 02 H0 : = 0 = . vs. H1 : = 0 . . At least one i 0p is not equal to 0i Test statistic using r.s. ( x1 , . . . , xn ): p 1 z 2 = n( 0 ) 1 ( 0 ) 2 x x p under H0 3 More frequently in practice, is unknown. Univariate test statistic using r.s. (x1 , . . . , xn ): t =n = 2 n( 0 )(s2 ) 1 ( 0 ) n x x 2 x 0 s 2 t2 n 1 under H0 = N1 0, scaled 2 df 1 N1 0, 2 4 Multivariate generalization (Hotelling s T 2 ): 2 T 2 = n( 0 ) S 1 ( 0 ) Tp,n 1 x x unbiased estimate of 1 under H0 = ( 0 ) x x and S are indep. since they are based on a r.s. from MVN S n ( 0 ) x characteristic form inverse sample cov. matrix for x n i=1 (xi = n( 0 ) x Np (0, ) random vector x)(xi x) n 1 1 n( 0 ) x Wp (n 1, ) random matrix divided by d.f. Np (0, ) random vector 5 6 Important properties of T 2 2 1. Sometimes we refer to the subscripts for Tp, distribution as 2 dimension and df (e.g., Tdim,df ) 2. Must have n > p Otherwise S is singular and T 2 cannot be computed. 3. Degrees of freedom for T 2 is same as for analogous univariate t-test: = n 1 for one-sample test = n1 + n2 2 for two-sample test 4. Alternative hypothesis is 2-sided (no such thing as H1 : > 0 ) Critical region is one-tailed (reject for large values) since test statistic is squared distance 7 5. p+1 2 p Tp,v = Fp, p+1 So, p-value for T 2 test is p-value = Pr Fp, p+1 > Critical value for T 2 test is 2 T ,p, = p+1 2 T p (n 1)p p F ,p, p+1 or F ,p,n p when = n 1 p+1 n p p 1 p p p 1 6. T 2 invariant under transformations of the form x = C where C is nonsingular x + d, 8 7. T 2 is the likelihood ratio test (LRT) of H0 : = 0 Under H0 the likelihood is 1 (2 )np/2 | |n/2 1 (2 )np/2 | |n/2 1 exp 2 n L( 0 , ) = = (xi 0 ) 1 (xi 0 ) i=1 n 1 exp tr 1 ( (xi 0 )(xi 0 ) ) 2 i=1 Using Result 4.10 (again), we obtain max L( 0 , ) = 1 (2 ) np 2 | 0 | n 2 exp np 2 where 0 = 1 n n i=1 (xi 0 )(xi 0 ) 9 Recall max , L( , ) = 1 (2 ) np 2 | | n 2 exp np 2 where = 1 n n i=1 (xi x)(xi x) Likelihood Ratio: max L( 0 , ) max ( , ) L( , ) 2 = = | | | 0 | n 2 < c = n = | | 1 = 1 1 + n 1 T 2 | 0 | Wilks Lambda is rejected for small or large T 2 * 2 ln 2 0 where = # of unrestricted parameters and 0 = # of parameters under H0 ex Turnips * T 2 = (n 1) || 0 | (n 1) | 10 A.ii. Con dence Regions Con dence region R(X): Set of possible values of in based on X R(X) is 100(1 )% C.R. if, before the sample is selected Pr{R(X) will cover the true } = 1 . C.R. for [100(1 )%] {all or {all 2 n( ) S 1 ( ) T ,p, } x x squared mult. distance from x x n( ) S 1 ( ) x p p+1 F ,p, p+1 } 11 Axes of the ellipsoid (based on eigenvalues 1 , . . . , p and eigenvectors e1 , . . . , ep of S): i n 2 T ,p, along ei 1 2 Elongation of ellipsoid: 12 Interest in C.I. s for individual components of x or linear combination a x. De ne z = a x 2 z N1 (a , a a) = N1 ( z , z ) Sample statistics: z=ax s2 = a Sa z Note: a1 = [0, 1, 0, . . . , 0] will yield a1 x = x2 and a2 = [1, 1, 0, . . . , 0] implies that a2 x = x1 x2 , etc. 100(1 )% C.I. for z is a x t ,n 1 2 a Sa n t-interval 13 Experimentwise Type I error rate (EER) Pr {at least one C.I. wrong } = 1 Pr{no C.I. s are wrong} = 1 (1 )p ex assuming independence of C.I. s = .05: EER for p = 10 is 1 (.95)10 .40 = Rewrite t-interval as all a x (a ( ))2 t2 n n 1 a Sa Is there a bound c2 which can replace t2 n 1 and de nes a C.R. that simultaneously contains a for all a?? 14 Preliminary result (2 50, JW) For B p.d. and x = 0 p p max x=0 (x d)2 = d B 1 d x Bx with maximum attained when x = cB 1 d, c = 0 So, max a=0 x (a ( )) a Sa = n ( ) S 1 ( ) = T 2 x x with maximum at a = c S 1 ( ) , x = c=0 discriminant function 15 = Simultaneously for all a, the interval ax or ax p a Sa F ,p, p+1 p+1 n 2 T ,p, a Sa n T 2 interval or when = n 1 ax (n 1)p F ,p,n p a Sa n(n p) will contain a with probability 1 . More conservative (wider) than t-interval Preserve EER Allows data-snooping 16 If we re willing to specify a few linear combinations a1 , . . . , ak before collecting the data, we might consider using intervals based on the Bonferroni inequality which are narrower than T 2 intervals but still protect EER for a nite set of l.c. s. Given C.I. s for k l.c. s a1 , . . . , ak , Ei : event that ith interval contains ai c P {Ei } = i Pr{all Ei } = = c 1 Pr{at least one Ei } c c c 1 Pr{E1 E2 Ek } k c Pr{Ei } i=1 1 = 1 i 17 Usually, specify i = So, k a x t 2k ,n 1 a Sa n Bonferroni Interval 18 Critical values for 95% C.I. s for 1 , . . . , p 19 Notes: Often useful to examine the discriminant function a = S 1 ( ) x in n(a ( ))2 x max = T2 a=0 a Sa a indicates the relative contribution of the x s to the separation of x and 0 Comparisons of a1 , . . . , ap only informative when x s are commensurate (i.e., measured on the same scale with comparable variances) 20 If x s are not commensurate, consider coe cients a , . . . , a that p 1 are applicable to standard variables. Discriminant function in terms of standardized variables z = a 1 instead of z = a1 x1 + + ap xp OR a = D a 1 2 x1 x1 xp x1 + + a p s1 sp s 11 where D = 0 0 .. . spp ex Turnips 21 III.B. Comparison of Several Mean Vectors III.B.i. Paired Observations Let x1i and x2i be 2 p-variate responses for observation i (i = 1, . . . , n) ex LaVerl s SAT pre-class test grades and post-class grades Pre-class grades: x1i = (Quant = 640, Analyt = 610, Verbal = 490) Post-class grades: x2i = (Quant = 680, Analyt = 620, Verbal = 560) 22 1. Calculate di = x1i x2i 2. Calculate 1 d= n and 1 Sd = n 1 q 2 3. T 2 = nd S 1 d Tp,n 1 d n di i=1 n (di d)(di d) i=1 (n 1)p n p Fp, n p p+1 Note: is shorthand for the equivalence of the quantiles of two dist ns q Same follow-up analyses as in one-sample T 2 test/intervals apply here Con dence regions/intervals Discriminant functions 23 Alternatively, think of each observation x1i xi = 2p 1 x2i x1 x = 2p 1 x2 S11 S = 2p 2p S21 pre-tests post-tests S12 S22 0 1 .. . 1 0 .. . 1 Interest is in Cxi , where 1 1 C = p 2p 0 0 1 24 Note di = Cxi d = C x Sd = CSC 2 and T 2 = n C (CSC ) 1 C Tp,n 1 x x (n 1)p Fp,n 1 p+1 (n 1 p + 1) q (n 1)p Fp,n p n p q 25 An extension to a comparison of p treatments given to each subject over time x i1 evaluation after day 1 dosage xi2 evaluation after day 2 dosage xi = . . . xip evaluation after day p dosage Interest may lie in comparisons of treatment means 1 0 C = . . (p 1) p . 0 1 0 0 0 . . . 1 1 0 1 2 0 2 3 2 = . . . . . . . . . 1 p p p 1 i = 1, . . . , n 1 1 . . . . . . 0 0 26 2 T 2 = n(C ) (CSC ) 1 C Tp 1,n 1 x x (n 1)(p 1) F(p 1),n 1 (p 1)+1 (n 1 (p 1) + 1) (n 1)(p 1) Fp 1,n p+1 n p+1 e.g., if comparing 3 days, we might use 1 C = 2 3 1 0 2 1 linear increase/decrease in response 1 quadratic e ect on response e.g., if comparing 4 days, we might use 3 C = 1 3 4 1 1 linear 1 1 1 quadradic 3 3 1 cubic 1 3 27 B.ii. Two-Sample Comparisons Interest in 1 2 (di erence in two population means). Assumptions: x11 , x12 , . . . , x1n1 is a r.s. from Np ( 1 , ) x21 , x22 , . . . , x2n2 is a r.s. from Np ( 2 , ) Note that 1 = 2 = The two samples are independent In practice, we can relax these assumptions somewhat for large n. Let 1 xi = ni ni xij , j=1 ni i = 1, 2 1 Si = ni 1 (xij xi )(xij xi ) , j=1 i = 1, 2 28 Since and (n1 1)S1 Wp (n1 1, ) (n2 1)S2 Wp (n2 1, ) (n1 1)S1 + (n2 1)S2 Wp (n1 + n2 2, ) =(n1 +n2 2)Sp = E{Sp } = Since the two samples are independent ( 1 x2 ) Np 1 2 , x and x T = [ 1 x2 ( 1 2 )] 2 2 2 Tp, = Tp,n1 +n2 2 1 1 + n1 n2 1 1 1 + n1 n2 Sp [ 1 x2 ( 1 2 )] x q (n1 + n2 2)p Fp,(n1 +n2 2) p+1 (n1 + n2 2) p + 1 29 100(1 )% C.R. for 1 2 = : all 2 T 2 T ,p, = n1 + n2 2 where T 2 is the squared mult. distance between x1 and x2 or all T2 (n1 + n2 2)p F ,p,(n1 +n2 2) p+1 (n1 + n2 2) p + 1 Follow-up analyses t-interval : a x1 a x2 t ,n1 +n2 2 2 1 1 + n1 n2 a Sp a 30 Bonferroni interval : a x1 a x2 t 2k ,n1 +n2 2 1 1 + n1 n2 a Sp a k is # of contrasts of interest ex want intervals for each of p variables Then, [a1 , . . . , ap ] = Ip and k = p T 2 -interval a x1 a x2 2 where T ,p, 2 T ,p, 1 1 + n1 n2 a Sp a (n1 +n2 2)p (n1 +n2 2) p+1 F ,p,(n1 +n2 2) p+1 Examine discriminant function a = S 1 ( 1 x2 ) x p for indication of contribution of the variables to separation of the groups 31 If x s are not commensurate consider standardized coe cients 2 a = Dp a 1 where Dp s 11,p = diag{Sp } = 0 0 .. . spp,p ex Duchenne muscular dystrophy Individual tests using t /2 , t /2p , 2 T ,p, as critical values Test H0 : 1 = 2 using x3 , x4 , x5 , &x6 Examine discriminant function coe . Standardized coe cients 32 Testing 1 = 2 when 1 = 2 Univariate case ( Behrens-Fisher Problem ): t = where = s2 1 n1 s2 1 n1 2 x1 x2 n1 s2 1 approx + n2 s2 2 t + + s2 2 n2 s2 2 n2 2 (Welch, 1937, 1947) n1 +1 n2 +1 Hsu (1938) and Sche e (1959) argue that signi cance level for usual t-test is preserved when n1 = n2 33 Multivariate case: 2 T = ( 1 x2 ) x S1 S2 + n1 n2 1 ( 1 x2 ) 2 x p as (n1 p) , (n2 p) Signi cance level preserved for usual T 2 test when n1 = n2 and n1 and n2 are very large (Ito and Schull, 1964) If sample sizes are equal the signi cance level [of usual T 2 test] is not a ected (Carter, Khatri, and Srivastava, 1979) ? But do these properties hold with small to moderate sample sizes ? 34 Simulation Study in Christensen & Rencher (1997) For matrices of form 2 = k 1 , equality of sample sizes (n1 = n2 ) is less able to protect Type I error rate as p increases (Study considered small to moderate n1 , n2 (2p, 10p)) 35 For multivariate Behrens-Fisher problem, consider T = ( 1 x2 ) S 1 ( 1 x2 ) x x e 2 as a statistic, where Se = S1 S2 + n1 n2 36 There are several tests for 1 = 2 when 1 = 2 , and many of these use T For example: Yao (1965) test uses 1 1 = 2 )2 (T 2 2 approx 2 Tp, i=1 1 Si ( 1 x2 ) S 1 S 1 ( 1 x2 ) x x e ni 1 ni e 2 Note: this is a multivariate extension of Welch s approach to univariate problem Nel and Van der Merwe (1986) test uses = tr S2 + (tr{Se }) e 2 2 i=1 1 tr ni 1 Si ni 2 + tr Si ni 2 37 Simulation study: Nel and Van der Merwe (1986) and Kim (1992) have highest power among tests with unin ated Type I error rate ex Muscular Dystrophy 38 Tests for additional information Let x1i = y1i p 1 and x2i , i = 1, . . . , n1 be a r.s. from Np+q ( 1 , ) z1i q 1 y2i p 1 = , i = 1, . . . , n2 be a r.s. from Np+q ( 2 , ) z2i q 1 Start with y measurements Will the q 1 subvector z measured in addition to y signi cantly increase the separation of the two samples (or is z redundant in presence of y?) y1 y2 Sample means: x1 = and x2 = z1 z2 Syy Syz Common sample covariance matrix: Sp = Szy Szz 39 If y and z are independent: 2 2 2 Tp+q = Tp + Tq 2 2 If not independent: Compare Tp+q with Tp 2 Tp+q = n1 n2 ( 1 x2 ) S 1 ( 1 x2 ) x x p n1 + n2 2 Tp = n1 n2 ( y2 ) S 1 ( 1 y2 ) y yy y n1 + n2 1 Then, we can show that 2 Tadd 2 2 Tp+q Tp 2 = ( p) Tq, p 2 + Tp or ( p)q Fq, p q+1 p q+1 40 Fadd = p q+1 q 2 2 Tp+q Tp Fq, p q+1 2 + Tp where = n1 + n2 2 If just checking the addition of one variable: 2 Tadd F1, p ex Duchenne muscular dystrophy x3 and x4 are relatively inexpensive to measure compared to x5 and x6 . Are x5 and x6 important above and beyond x3 and x4 x3 , x4 , x5 , x6 may depend on age and season. Are x1 = age and x2 = season important? 41 B.iii. MANOVA (one-way) Comparing means from g groups Sample from population 1: x11 , x12 , . . . , x1n1 Sample from population 2: x21 , x22 , . . . , x2n2 . . . Sample from population g: xg1 , xg2 , . . . , xgng x j independent random samples N ( , ) is the common covariance matrix 42 Instead of testing H0 : 1 = 2 = = g vs. H1 : at least two s are unequal we usually reparameterize = + Thus x j treatment e ect N ( + , ) and H0 : 1 = 2 = = g Our model: x j = + + e j, = 1, . . . , g, j = 1, . . . , n For uniqueness (identi ability), we impose the constraint g n = 0 =1 43 Decomposition of sample: observed x j = overall sample mean x + ( x) + (x j x ) x estimated treatment e ect residual ej 44 Multivariate analog of total (corrected) sum of squares is g n g g n (x j x)( ) = =1 j=1 total corrected sum of squares and cross products matrix =1 n ( x)( ) + x =1 j=1 =H Between groups matrix (x j x )( ) =E Within groups matrix g = =1 (n 1)S Notes: Assuming no linear dependencies, rank{H} = min(p, H ) S is the covariance matrix for the E ( g =1 th sample. So, = 1 E n ) g where rank{E} = min(p, E ) 45 MANOVA TABLE (one-way) Source Treatment Error Total (corrected) Wilks The likelihood ratio test of H0 : 1 = 2 = = g rejects H0 when = |E| ,p, H , E |E + H| SS Matrix H E H+E ( d.f. H = g 1 E = ( g =1 n ) g g =1 n ) 1 Note: Reject for small values of . As in univariate anova F -test, we accept when total SS (E + H) is dominated by error (E). Note: We sometimes refer to the subscripts of the p, H , E distribution as dimension, numerator df, and denominator df (e.g., dim,dfnum ,dfden ) 46 Properties of Wilk s : 1. For statistic to be obtained, we need E p. 2. Degrees of freedom H and E are the same as in analogous univariate case; e.g., one-way model: H = g 1 and g E = =1 n g 3. Let 1 , . . . , s be the s non-zero eigenvalues of E 1 H, where s 1 s = min(p, H ). Then = i=1 1+ i . 4. Critical value ,p, H , E decreases as p increases. Thus, adding variables decreases power unless variables contribute to separation. 47 5. When H = 1 or H = 2 or p = 1 or p = 2, can be transformed to follow an F distribution. If H = 1 E p + 1 1 Fp, E p+1 p If H = 2 E p + 1 1 F2p,2( E p+1) p If p = 1 E H If p = 2 ( E 1) 1 F2 H ,2( E 1) H 1 F H , E 48 6. Approximate tests For p > 2 or H > 2 and n large 2 = E 1 approx (p H + 1) ln 2 H p 2 1 3 2 Approximately valid when p2 + H E 1 2 (p H + 1) More correct approximate distribution for (exact when H or p is 1 or 2): F= 1 1/t df2 df1 1/t approx Fdf1 ,df2 df1 = p H 1 df2 = wt 2 (p H 2) w = E + H 1 (p + H + 1) 2 2 4 p2 H 2 for p2 + H 5 > 0 2 p2 + H 5 t= 2 1 for p2 + H 5 0 (or p + H > 3) (or p + H 3) 49 Other MANOVA Tests Let ( 1 , . . . , s ) be the ordered eigenvalues of E 1 H, where s = min(p, H ) = rank of H Roy s Largest Root: = 1 Note: SAS and most authors denote Roy s Largest Root as 1 (the largest root of E 1 H). ACR de nes Roy s Largest Root as 1 1 = 1+ 1 , which is the largest root of (E + H) 1 H. 50 Approximate F -statistic (used by SAS): ( E d + H ) 1 d is an upper bound for true F which is distributed F = Fd , E d + H where (d = max (p, H )) Thus, F -test is anti-conservative (yields lower bound on p-value) The eigenvector a1 corresponding to 1 comprises the discriminant function coe cients. For programs unable to obtain eigenvalues of nonsymmetric matrices, we can use the fact that 1 is a solution to both (E 1 H I)a = 0 and (E 2 HE 2 I) 1 1 E2 a e vector of E 2 HE 2 1 1 1 =0 symmetric 51 Pillai s Trace: s V= i=1 i 1 + i 1 s = tr (E + H) H= i=1 i 1 where 1 , . . . , s are the s ordered e vals of (E + H) Note 1: E 1 SS H is analagous to betweenSS within + H) 1 H is analagous to between SS total SS H (E Large Ratio = Reject H0 Note 2: i = i i and i = 1 + i 1 i 52 Approximate F -statistic (used in SAS): FV = where s = min( H , p) m= N= 1 2 (| H p| 1) 1 2 ( E p 1) (2N + s + 1) (2m + s + 1) V s V Fs(2m+s+1),s(2N +s+1) 53 Lawley-Hotelling Trace s U= i=1 i = tr{E 1 H} Approximate F -statistic (used in SAS): Fu = 2(sN + 1) U Fs(2m+s+1),2(sN +1) s2 (2m + s + 1) Also known as Hotelling s generalized T 2 54 Why four test statistics? All 4 are exact tests (i.e., have size ), but when H0 not true they have di erent power For p = 1, 1 , . . . , k can be ordered along 1 dimension (line) and F -test is U.M.P. For p > 1, 1 , . . . , k are points in s = min(p, H ) dimensions. But means may in fact occupy only a subspace of the s dimensions; e.g., they may lie close to a line (1-D) or a plane (2-D). 55 56 ex Visual memory task x1 = % correct on positive stimulus questions x2 = % correct on negative stimulus questions g = 3 (One healthy group and two impaired groups) p=2 s = 2 out of p = 2 eigenvalues should be > 0 57 ex Egyptian skulls x1 = maximum breadth of skull (mm) x2 = basibregmatic height of skull (mm) x3 = basialveolar length of skull (mm) x4 = nasal height of skull (mm) g = 3 (4000 B.C., 3300 B.C., 1850 B.C.) p=4 s = 2 out of p = 4 e vals > 0 ex Rootstock x1 = trunk girth at 4 years (mm 100) x2 = extension growth at 4 years (m) x3 = trunk girth at 15 years (mm 100) x4 = weight of tree above ground at 15 years (lb 1000) g=6 p=4 s = 4 out of p = 4 eigenvalues > 0 58 Follow-up analyses Although only multivariate tests could detect group di erences above, we still are often interested in follow-up analyses after conducting a multivariate analysis. Univariate hypothesis (F ) tests Multivariate contrasts Con dence intervals/tests for ij kj (treatment di erences for j th variable) Analysis of discriminant function 59 Univariate F -tests Often interested in univariate ANOVA for testing H0,i : 1i mean of ith var. for 1st group = 2i = = gi , i = 1, . . . , p Some advocate a protected univariate test approach: 1. Conduct overall size test of H0 : 1 = = g using multivariate test (e.g. ) 2. Test each of H0,i (i = 1, ..., p) at level only if multivariate test in step 1 rejects. [That is, when H0 is accepted this approach automatically accepts H0,1 , H0,2 , . . . , H0,p .] 60 De ning our experiment by the p tests in step 2, the overall EER is (for independent variables when H0 is true): Pr{at least one H0,i rejects} = ( ) (1 (1 ) ) p What about properties of individual tests when H0 is false?? Suppose: i1 12 0 i2 1 = . = + . and i = . = , i = 2, . . . , g . p 1 . . . . . ip 0 1p 11 Let be some value such that our test of H0 : 1 = = p using has power = .50. 61 Consider the partial experiment de ned by the p 1 tests of H0,2 , H0,3 , . . . , H0,p . The partial EER for this scenario (assuming independence) is Pr {at least one rejection among H0,2 , . . . , H0,p } A = = Pr { A | rejects} Pr { rejects} 1 (1 ) p 1 (.50) Thus, the partial EER can be dramatically larger than ex p = 10, = .05 partial EER .20 = Conclusion: Protected F test approach protects overall EER, but may have poor properties for other inferences Consider tests at p level 62 Contrasts (multivariate) Already considered contrasts of the type C for testing q pp 1 H0 : C = 0, where each row of C sums to 0 ex Linear trend among 4 observations? day1 day2 H0 : 3 1 1 3 day3 day4 Here consider contrasts of the type = c1 1 + c2 2 + + cg g = Mc where M = 1 2 g p g 63 = c1 x1 + c2 x2 + + cg xg ci = ni i=1 2 g g g var{ } = i=1 ci 2 ni var{ } = where Sp = 1 E E and E = i=1 g i=1 (ni ci 2 ni 1). Sp So, our test is based on T 2 = var or de ne H1 = ex Rootstock 1 g i=1 ni c2 i 1 2 Tp, E and = |E| p,1, E |E + H1 | 64 Con dence intervals for treatment di erences ( ij kj ) Interest in components of di erence between groups i k = i k Speci cally, interested in the jth component of this di erence vector ij kj = ij kj which is estimated by xij xkj Because we often want to obtain con dence intervals for all g(g 1)/2 pairwise comparisons for each of p variables simultaneously, we use a Bonferroni adjustment to protect overall EER. 65 100(1 )% (Simultaneous) Con dence Interval for ij kj is: 1 1 + ni nk ( ij xkj ) t x [ / (pg(g 1))] ,[ divided by # of comparisons = pg(g 1)/2 2 g i=1 (ni 1)] sp ,jj where spl, is the j th diagonal element of Sp = E/ ( g i=1 (ni 1)) 66 Warning for SAS implementation SAS uses upper distribution g(g 1) quantile instead of upper pg(g 1) quantile of t (Bonferroni intervals are part of univariate output.) Adjust by specifying ALPHA in MEANS statement ex p = 3, g = 5, desired overall EER = .05 proc glm; class group; model y1 y2 y3 = group; means group/bon alpha = .016667 run; ex Rootstock .05 p; 67 Analysis of discriminant function (More detail to come in Section V of the course) g = 2 case: Choose a to maximize (for a = 0): x [a ( 1 x2 )] x x a ( 1 x2 ) ( 1 x2 ) a = a Sp a a Sp a = a = S 1 ( 1 x2 ) x p g > 2 case: Choose a to maximize (for a = 0): a Ha a Ea = 1 = largest e value of E 1 H and a1 is corresp. e vec 1 = Relative importance of 1st disc fcn = 1 s i=1 i 2 68 ex Rootstock Girth4 Growth 2.91 .2627 = ai Girth15 11.97 .6532 so .6532 4.2908 .0738 1.712] Weight 12.16 .0738] Univariate F 1.93 a = [.4703 Recall that a = ai sp i a = = ,ii 1 E E, 1 [.4703 .3200 42 [.0411 .1413 .2627 12.1428 .2088 .0149] 69 Tests for additional information y j p 1 = , j = 1, . . . , n z j q 1 th Let x j be (p + q)-variate observations from the group Wish to determine if z makes a signi cant contribution beyond y in detecting separation of groups Calculate: Eyy = E (p+q) (p+q) Ezy Eyz Hyy and = H (p+q) (p+q) Ezz Hzy yz = |E| |E + H| Hyz Hzz y = |Eyy | |Eyy + Hyy | 70 Test of additional info: yz q , H , E p y # of vars in z z|y = partial statistic # of vars in y 71 Two-Way MANOVA ( xed-e ects) Model: xijk = + i + j + ij + eijk i = 1, . . . , a j = 1, . . . , b k = 1, . . . , n a i=1 (for simplicity, assume nij = n i, j) b j=1 i = i = a i=1 ij = b j=1 ij = 0 Assume eijk Np (0, ) 72 xi. = average over ith level of factor A x.j = average over jth level of factor B xij = average over ith level of A and jth level of B 73 Test for A, B, and AB (interaction): |E| p,a 1,ab(n 1) |E + HA | |E| B = p,b 1,ab(n 1) |E + HB | |E| AB = p,(a 1)(b 1),ab(n 1) |E + HAB | A = 74 Follow-up analyses Individual F -tests (univariate anova s) As You Might Expect (AYME) Contrasts AYME C.I. s for treatment e ects AYME Analysis of discriminant function AYME For analyzing contribution of p variables to separation of levels of A use rst dicrim. function (e vector) of E 1 HA Analyzing levels of B use E 1 HB Analyzing levels of AB use E 1 HAB 75 An interlude about interactions... Suppose we have two levels of factor A and two levels of factor B: yijk = + i + j + ij + ijk , i = 1, 2, j = 1, 2 76 Is Main E ect for A interpretable in Scenarios II, III, and IV? Yes, if signi cance simply refers to size of e ect 1 2 (i.e., e ect of A averaged over levels of B ). Signi cant doesn t mean one level is best Signi cance of Main E ect for A is a ected by number of levels of B and sample size for each level 77 Mixed Model MANOVA (Split-plot) ex 4 3 3 temperatures days metal alloys (Ti , i = 1, . . . , 4) (Dj , j = 1, . . . , 3) t=4 d=3 (Mk , k = 1, . . . , 3) m = 3 xijk is a p-variate response of metal strength 78 79 III.B.iv. Pro le Analysis p-variate response consists of tests, questions, etc. measured on members of g groups. ex Guinea pigs on three diets Weights measured at ends of week 1, 3, 4, 5, 6, & 7 Break hypothesis: H0 : 1 = 2 = = g into three more speci c hypotheses: H01 : The g pro les are parallel ex H0,1 true might yield a pro le plot like: 80 H02 : The g pro les are at same level ex H0,1 and H0,2 true might yield a pro le plot like: H03 : The g pro les are at ex H0,1 , H0,2 , and H0,3 yields the pro le plot: 81 Formalizing the null hypotheses Parallelism : Di erence in responses between any time points is the same for all groups. H01 : 1j 1(j 1) = 2j 2(j 1) = = gj g(j 1) for j = 2, . . . , p OR C 1 = C 2 = = C g 1 1 0 0 1 1 where C = . . . . (p 1) p . . . . . 0 0 0 0 0 . . . 1 0 0 . . . 1 or C can be any other full row rank (p 1) p matrix such that C1 = 0 82 Same level : Total (or average) response (over time) is the same for all groups. H02 : 1 1 = 1 2 = = 1 g Note: If H01 holds, we can also refer to H02 as the hypothesis of coincident pro les and H02 can be written: H02 : 1j = 2j = = gj for j = 1, . . . , p 83 Flatness : No change in response (over time) for the pro les responses at each time (averaged across groups) are the same. H03 : OR C 11 + 21 + = g1 1p + 2p + = gp = = g g 1 + + g g = O or C = O average pro le Note: If H01 and H02 hold, H03 can also be written H03 : 11 = 12 = = 1p = 21 = = gp that is, all pg response means are equal. 84 Tests Test for H01 : = |CEC | p 1, H , E |C(E + H)C | Test for H02 : = 1 E1 1, H , E 1 E1 + 1 H1 1 E F H , E H Test for H03 : g T= =1 2 n (C ) x 1 CEC e 1 2 C Tp 1, E x E (p 1) + 1 2 T Fp 1, E (p 1)+1 E (p 1) 85 ex Guinea Pigs H01 : parallel? H02 : same level? H03 : at? 86 III.B.v. Repeated Measures Similarities to pro le analysis Each subject measured under several treatments or time points Comparing means of treatments applied to each subject: within-subjects tests Comparing levels of factors assigned to groups of subjects: between-subjects tests 87 Structure of g-groups R.M. experiment 88 Univariate model: (split-plot): xijr = + Bi + S(i)j + Ar + BAir + ijr between within interaction ANOVA Table 89 Standard univariate assumption: var{xij } = = 2 Ip i, j Univariate F -tests still valid as long as (p 1) p orthonormal contrast matrix sphericity C C = 2 I This condition is often called sphericity (but we ll say generalized sphericity for clarity) ex For p = 4, we could use 3/ 12 1/ 12 1/ 12 1/ 12 C= 0 1/ 6 1/ 6 2/ 6 1/ 2 0 0 1/ 2 90 Special case of C C = 2 I: 1 p = 2 . . . p p 1 .. .. . p p 2 . = [(1 p)I + p11 ] . . 1 . p This case is called compound symmetry 91 Univariate strategies 1. Assume generalized sphericity Fehlberg (1980): Use H0 : C C = 2 I preliminary test using = .40. [This test to be discussed later in the course] If hypothesis is accepted, use standard F -tests . . . M SA . . . for A: F = Fp 1,g(n 1)(p 1) M SE M SAB . . . for AB: F = F(g 1)(p 1),g(n 1)(p 1) M BUT, SE even mild departures from C C = 2 I can seriously in ate Type I error (Boik, Psychometrika, 1981). 2. Conservative test: M SA F1,g(n 1) M SE M SAB . . . for AB: F = F(g 1),g(n 1) M SE Too conservative . . . for A: F = 92 3. Adjusted F -tests A compromise between approaches 1 and 2 when sphericity violated. Greenhouse and Geisser (1959) recommend approximate F -tests involving within-subjects factor which reduce numerator and denominator d.f. by a factor of 1 tr p 11 2 = (p 1)tr 1 p 11 2 SAS: G G E To estimate , use = E F -tests . . . M SA . . . for A: F = M SE F (p 1), g(n 1)(p 1) M SAB . . . for AB: F = M SE F (g 1)(p 1), g(n 1)(p 1) 93 1 and ( p 1 , 1) sphericity general holds (non-spherical) Approach is generally too conservative, especially for small n Huynh and Feldt (1976) give another expression SAS: H F for Less conservative H F can exceed 1 set equal to 1 94 Multivariate Model: xij = + i + ij Notes about i i is a p-vector of main e ects for group i Tests on factor A (within subjects) and AB interaction constructed with contrasts of i (as in pro le analysis) Standard multivariate assumption: var{xij } = i, j Note: is completely unrestricted (no sphericity requirement, etc.) Several similarities with g groups pro le analysis 95 Contrast Matrices in SAS Proc GLM ( repeated statement) Assume p=5 times/variables contrast or contrast(5): 1 0 0 0 1 0 1 0 0 1 0 0 1 0 1 0 0 0 1 1 1 0 0 0 contrast(2): 1 0 0 0 1 1 0 0 1 0 1 0 1 0 0 1 2 2 1 1 polynomial 1 2 4 0 1 2 1 2 0 6 1 2 2 1 4 1 repeated time 5 (1 2 5 10 20) polynomial : .43 .36 .17 .15 .80 .43 .21 .33 .71 .39 .43 .14 .73 .51 .08 .49 .78 .37 .09 .01 96 helmert: 1 .25 .25 .25 .25 0 1 .33 .33 .33 0 0 1 .5 .5 0 0 0 1 1 mean or mean(5) The helmert contrast matrix identi es the time at which the treatments cease to change or plateau 1 .25 .25 .25 .25 .25 1 .25 .25 .25 .25 .25 1 .25 .25 .25 .25 .25 1 .25 1 1 0 1 0 0 00 profile: 0 1 1 0 0 0 0 0 1 0 1 1 97 Test for A (within subjects) Analogous to atness test in pro le analysis Want to compare means for x1 , . . . , xp averaged across levels of B g Let = i=1 i /g = ( 1 , . . . , p ) H0 : 1 = . . . = p or C = O 1 0 C= . . . 0 1 0 0 0 . . . 0 0 or some similar . contrast matrix . . ex 1 1 . . . . . . 0 0 1 1 98 Test statistic for A: 2 T 2 = N (C x ) (C Sp C ) 1 ( C x) Tp 1, E g i=1ni grand mean E E (p 1) p E (p 1) + 1 2 T Fp 1, E (p 1)+1 E (p 1) 99 Test for B (between subjects) Analogous to same level test in pro le analysis Want to compare group means (averaging over p levels of A) 1 1 H0 : 1 1 = = 1 g or 1 1 = = 1 g p p SAS That is, we can just conduct one-way ANOVA on zij = 1 xij , so the test statistic for B is = 1 E1 1, H , E 1 E1 + 1 H1 Multivariate between-subject test results generally agree with univariate test results 100 Test for AB interaction Analogous to parallelism in pro le analysis H0 : C 1 = = C g Test statistic for AB: = |CEC | p 1, H , E |C(E + H)C | 101 ex Wear of fabrics Measured in 3 periods (within-subjects factor) 1st 1000 revolutions 2nd 1000 revolutions 3rd 1000 revolutions 2 abrasive surfaces (between subjects factor #1) 2 llers (between subjects factor #2) 3 levels of proportion of ller (between subjects factor #3) 25% ller 50% ller 75% ller ? Linear or Quadratic trend in proportion of ller? ? Linear or Quadratic trend in periods ? How do univariate and multivariate tests compare? 102 Repeated Measures with 2 Within-Subjects Factors 103 Model: xij = + i + ij i is a p-vector of main e ects for group i e ects for A, B, AB, AC, BC, ABC assessed with contrasts Denote a = number of levels for factor A Denote b = number of levels for factor B To test factors A, B, and AB, specify contrast matrices with (a 1), (b 1), and (a 1)(b 1) linearly independent rows, respectively. 104 ex Blood data Compare 4 di erent reagents used in blood testing. (Reagent 1 is standard and reagents 2, 3, 4 are inexpensive alternatives.) Measuring 3 blood counts (white blood, red blood, hemoglobin) 2 groups of 10 subjects with potentially di erent blood properties each subject s sample has 12 measures 105 a1 1 12 R1 vs. R2 a2 R1 vs. R3 a3 R1 vs. R4 b1 1 12 white vs. hemoglobin (or linear in bc s) (or quadratic in bc s) b2 red vs. white+hemo 2 a1 b 1 a b 1 2 G = . where is an element-wise product . 6 12 . a3 b 2 (So, rst row of G is [1 0 -1 -1 0 1 0 0 0 0 0 0]) 106 Test for A (Reagents): 2 T 2 = N (A x ) (A Sp A ) 1 A Ta 1, E x grand mean 1 E E 1) when only one between subjects factor is used g i=1 (ni ex OR = 2 2 TReagent T3,18 where H = N xx is from the partitioning g n1 |AEA | |A(E+H )A | a 1,1, E xij xij = E + H + N xx i=1 j=1 107 Test for B (Blood counts): 2 T 2 = N (B ) (BSp B ) 1 B Tb 1, E x x OR = |BEB | |B(E+H )B | b 1,1, E Test of AB interaction: 2 T 2 = N (G ) (GSp G ) 1 G T(a 1)(b 1), E x x OR = |GEG | |G(E+H )G | (a 1)(b 1),1, E 108 Test for C (groups): Conduct ANOVA test (F -test) using zij = 1 xij , i = 1, . . . , g, j = 1, . . . , ni OR = 1 E1 1 E1+1 H1 1, H , E Tests for AC, BC, ABC interactions: = = = |AEA | |A(E+H)A | |BEB | |B(E+H)B | |GEG | |G(E+H)G | a 1, H , E b 1, H , E (a 1)(b 1), H , E Note: Between subjects e ects (e.g., C) and associated interactions (e.g., AC, BC, ABC) use H (not H ) ex Blood data in SAS 109 III.B.vii. Tests on Covariance Matrices (Reference: Rencher, 1995, Ch. 7) H0 : = 0 vs. H1 : = 0 (assuming MVN) u = ln | 0 | ln |S| + tr{S 1 } p 0 is a modi cation of the likelihood ratio with = degrees of freedom for S. large: u 2 p(p+1) 1 2 small to moderate: 1 1 6 1 2p + 1 2 p+1 u 2 p(p+1) 1 2 110 H0 : = 2 I ( sphericity . . . assuming MVN) Likelihood ratio test: |S| = (tr{S}/p) 2 ln = n ln u where u = 2 n n 2 pp pp |S| = = p (tr{S}) ( p i=1 p i=1 i p i ) and 1 , . . . , p are the e vals of S large: n ln u 2 p(p+1) 1 1 2 small to moderate: 2p2 +p+2 6p ln u 2 p(p+1) 1 1 2 111 Note: Testing C C = 2 I use CSC in place of S in the test, i.e., n ln where ex p=4 3/ 12 C= 0 0 1/ 12 1/ 12 1/ 6 2/ 6 0 1/ 2 1/ 12 1/ 6 1/ 2 (p 1) p (p 1)p 1 |CSC | (tr{CSC }) p 1 2 (p 1)(p) 1 1 2 C has orthonormal contrasts as its rows Often called Mauchly s test Calculated by SAS with PRINTE option of REPEATED statement in PROC GLM. Fehlberg (1980) recommends a preliminary test of = 2 I at = .40 before using standard univariate F -tests in r.m. analysis. 112 H0 : 1 = 2 = = g (assuming MVN for all groups) Box s M = |S1 | 1 2 |S2 | 2 2 |Sp | i i2 |Sg | g 2 where i = ni 1, i = 1, . . . , g and Sp = M near 0 reject H0 M near 1 accept H0 g g i=1 i Si g i=1 i Note: M = i=1 |Si | |Sp | i 2 . . . is maximized at 1 when S1 = = Sg . . . approaches 0 when one or more |Si | is very small (with other |Si | large) u = 2(1 c1 ) ln M 21 (g 1)p(p+1) ] [2 where c1 = g 1 i=1 i g i=1 1 i 2p2 +3p 1 6(p+1)(g 1) 113 Note: M -test not recommended pre-test before T 2 or MANOVA tests Sensitive to nonnormality (often of little concern) and innocuous forms of heterogeneity (e.g., varying amounts of kurtosis) Note: A better approximation is u Fa1 ,a2 . See ACR for details. 114 Multivariate Multiple Regression A short review of vec and Kronecker notation a 1 an m-vector a 1 a2 Let A = . = m n . . am a 1 a 2 vec A = . . . a n a 2 a n = (aij ) Splus: c(A) gives vec A 115 Let B = (bij ) p q a11 B a12 B a21 B a22 B A B = . . . m n p q . . . am1 B am2 B .. . a2n B . . . amn B a1n B mp nq Splus: kronecker(A, B) gives A B 116 Some properties (without proof) Assuming that all dimensions are appropriate for matrix multiplication. . . (a) (A B)(C D) = (AC) (BD) (b) vec (ABC) = (C A)vec B (c) tr{AB} = (vec A ) vec B = (vec A) vec B (d) tr{ABCD} = (vec A ) (D B) vec C = (vec A) (B D ) vec C (e) (A B) = A B (f) (A B) 1 = A 1 B 1 117 Univariate Multiple Regression: n 1 y=X +e n r r 1 n 1 Assume E{e} = 0 and var{e} = 2 In . Then O.L.S. estimator = (X X) 1 X y is B.L.U.E. for . Note: we ll use q to denote the # of xs and r = q + 1 to denote the # of columns in the X matrix when using an intercept 118 Multivariate Multiple Regression: Y = XB + where y p y 1 . Y = . = y 1 . yn B = 1 2 p e 1 . = . = e 1 . en e p 119 Note that y j n 1 = X j + e j n r r 1 p p n 1 Assume E{ } = 0, var{ei } = , and cov{ei , ek } = 0 for all p p i=k Question: Is B = (X X) 1 X Y a B.L.U.E.? r p 120 Rewrite model: vec Y = vec (XB) + vec ( ) = (Ip X) rank = pr when rank(X)=r vec B + vec ( ) pr 1 e np 1 Note: E{e} = 0 I e 1 11 n 21 In . . = . var{e} = var . . . e p p1 In and np 1 12 In 22 In . . . p2 In .. . 2p In . = p p In . . pp In 1p In Since var{ e } does not take the form 2 Inp , the B.L.U.E. for will be the G.L.S. estimator for (which depends on the unknown ) BUT. . . 121 = Ip X n p (Ip X ) 1 ( In ) 1 Ip X n p (Ip X) ( In ) 1 vec Y ( 1 In) 1 = 1 X (Ip X) 1 (Ip X ) 1 In vec Y [by prop s (a),(e),(f)] = 1 (X X) = (X X) 1 1 X vec Y vec Y [by prop (a)] [by prop (f)] [by prop (a)] 1 X vec Y = Ip (X X) 1 X B = (X X) 1 X Y [by prop (b)] O.L.S. = G.L.S. is BLUE! (Even when is unknown) 122 Despite the fact that the p variables yi1,...,yip are correlated, all the info needed to estimate i is found in y i only. That is, r 1 multivariate regression coe cient matrix B can be formed by r p pasting together the p columns from p separate univariate regressions (as long as each regression uses the same predictors X ) n r But all ij in B are intercorrelated . . . must take multivariate approach to inference 123 Assumptions for Multivariate Multiple Regression: Model is: Y = X B + or vec Y = (I X) vec B +vec n p n r r p n p = Assumptions: 1. E{Y} = XB or E{ } = 0 2. var{vec Y} = var{vec } = In (That is, var{yi } = for all i = 1, . . . , n and cov{yi , yj } = 0 for all i = j) p p 124 Some properties of B = (X X) 1 X Y 1. B is called the least squares estimator because it minimizes E = = (Y XB) (Y XB) (where E is an error matrix p p analogous to the E matrix in MANOVA). Matrix is minimized in several senses: (a) Let B be some other estimate of B. Then, (Y XB) (Y XB) = (Y XB) (Y XB) + A where A is a positive de nite matrix (b) B = B minimizes tr{(Y XB) (Y XB) (c) B = B minimizes |(Y XB) (Y XB)| 125 2. Let Y = XB = X(X X) 1 X Y be predicted values and = Y Y = (I X(X X) 1 X )Y be residuals Then (a) Residuals are perpendicular to the columns of X X = X (I X(X X) 1 X )Y = 0 (b) Residuals are perpendicular to the columns of Y Y = B X (I X(X X) 1 X )Y = 0 p p r p (c) Total sum of squares and cross products ( Total SS and CP ) can be partitioned as: Y Y = (Y + ) (Y + ) YY = YY + total SS&CP matrix predicted SS&CP matrix error SS&CP matrix 126 3. B is B.L.U.E. for B Minimum variance estimator among all unbiased estimators If columns of are normal, B is B.U.E. 4. Elements of B are intercorrelated 01 11 B = . . r p . q1 02 12 . . . q2 .. . 0p 1p . . . qp s in each row are correlated due to correlation in y s in each column are correlated due to correlation in x 127 5. Unbiased estimate of var(yi ) = var(ei ) = . E (Y XB) (Y XB) = = n q 1 n q 1 n q 1 1 = (Y Y B X Y) n q 1 S= Proof: E{ } = E{Y XB} n p = E{(In X(X X) 1 X )Y} = (In X(X X) 1 X ) E{XB + } = (In X(X X) 1 X ) E{ } =0 n p 128 E{ i e j } = E e =E (In X(X X) 1 X )y i (In X(X X) 1 X )e i (In X(X X) 1 X )y j (In X(X X) 1 X )e j In X(X X) 1 X e j } = E{e i In X(X X) 1 X = E tr In X(X X) 1 X In X(X X) 1 X e j e i ij In = tr{ In X(X X) 1 X E {e j e i }} = ij tr In X(X X) 1 X = ij tr {In } tr X X(X X) 1 = ij (n (q + 1)) E 1 n q 1 = = ( ij ) 129 Note: If X is not full rank, we can obtain similar results based on B = (X X) X Y . . . we ll leave that discussion for linear models ! Another note: var and E 1 n q 1 E np 1 e = 1 n q 1 E In I n = In 6. Variance of (i.e., var{vec B}) var = var (I X) (I X) 1 1 1 1 rp 1 (I X) vec Y 1 = Ip (X X) = Ip (X X) = Ip (X X) = (X X) = (X X) X var{vec Y} Ip (X X) X var{vec } Ip (X X) X ( In ) Ip X (X X) X In X (X X) 1 X X 1 1 1 1 130 Notes: S (a) var{ } = 1 1 E (X X) n q 1 (b) cov{ i , j } = ij (X X) 1 (c) cov{ i , e j } = r 1 n 1 cov (X X) 1 X y i , In X (X X) 1 1 X 1 y j X e j = cov (X X) = (X X) 1 X e i , In X (X X) 1 X ij In In X (X X) 1 X 1 = ij (X X) =0 r n X (X X) 1 X X (X X) X 131 (d) Estimating mean of x0 B 1 r r p x0 B is an unbiased estimator of x0 B var{x0 B} = (x0 (X X) 1 x0 ) scalar (e) Estimating a new observation y0 using x0 y0 = x0 B + e0 x0 B is an unbiased estimator of y0 forecast error variance var{y0 x0 B} Note that cov{y0 , x0 B} = cov e0 , x0 (X X) 1 X (XB + ) = cov e0 , x0 (X X) 1 X e1 . since e0 is indep. of = . . en =0 p p 132 So var{y0 x0 B} = var{y0 } + var{x0 B} 2 cov{y0 , x0 B} = + x0 (X X) 1 x0 + 0 = 1 + x0 (X X) 1 x0 133 7. MLE s of B and r p Thus far, we have assumed E{e } = 0 and var{e} = In vec If we assume: np 1 e Nnp (0, In ) then MLE s of B and are B = (X X) 1 X Y and 1 1 = = E p p n n E Wp (n q 1, ) Proof: omitted. where 134 8. Model Corrected for Means Rewrite Y = X B+ n p as Yc n p and = Xc Bc + where q = # of predictors = r 1 n q q p y11 y 1 y12 y 2 y1p y p . . . . Yc = . . yn1 y 1 yn2 y 2 ynp y p x x 1 x1q x q 11 . . . . Xc = . . xn1 x 1 xnq x q 135 Then Bc = S 1 Sxy xx Syy where S = Sxy Y = y 1 1n Syx Sxx y p 1n + Xc Bc is the sample covariance matrix of the p + q variables (y1 , . . . , yp , x1 , . . . , xq ) n q q p 136 Hypothesis Tests (assuming e Nnp {0, Ir }) H0 : B1 = 0 q p (Test of overall regression) 0 1 p vector of intercepts where B = r p B1 q p Partition the total SS and CP matrix: pnnp Y Y = Y XB Y XB +B X Y =Y Y B X Y = E To avoid inclusion of 0 = 0 as part of the null hypothesis, we subtract n y : y Y Y n y y corrected total SS & CP = Y Y B X Y + B X Y n y y =E pp =H p p 137 = |E| |Y Y B X Y| = p, |E + H| |Y Y n y | y q ,n q 1 r 1 n r H is large when ij s are large The 4 MANOVA statistics can be calculated as functions of the eigenvalues of E 1 H, ( 1 , . . . , s ): Wilks : = Roy s: = 1 Pillai s Trace: V = s i i=1 1+ i s i=1 s 1 i=1 1+ i Lawley-Hotelling Trace: U = i Critical values (and p-values) based on approximate F -distributions given on the MANOVA pages on these notes . . . use: s = min(p, q) m = 1 (|q p| 1) 2 N = 1 (n q p 2) 2 138 Essential dimensionality of E 1 H is the essential dimensionality of p p B1 . q p For example, a single non-zero eigenvalue (i.e., rank of B1 is 1) could be due to several causes: 1. B1 has only one nonzero row only one of the x s predicts the y s 2. B1 has only one nonzero column only one of the y s is predicted by the x s 3. All of the rows of B1 are linear combinations of each other x s act alike in predicting y s [or, in other words] All of the columns of B1 are linear combinations of each other only one dimension in the y s as they relate to x s 139 Essential dimensionality of E 1 H is number of substantially non-zero eigenvalues and takes value less than or equal to s = min(p, q) can also be calculated from the partitioned sample covariance matrix of (y1 , . . . , yp , x1 , . . . , xq ) S = p p (p+q) (p+q) Sxy q p Syy Syx q q p q Sxx using |S| |Sxx ||Syy | which is essentially a test of independence between y and x since Independence of y and x |S| = |Syy ||Sxx | 140 H0 : Ba = 0 Tests on a subset of the x s Hypothesis states that the y s do not depend on the last h of the x s. That is, H0 : Badd = 0 where Bred (r h) p B = r p Badd h p Compare SS and CP matrix for full and reduced models: Hdi = B X Y Br Xr Y p p di erence in regression SS and CP and Efull = Y Y B X Y p p E matrix based on full model 141 Then xq h+1 , ,xq |x1 , ,xq h = = |Efull | |Efull | = |Efull + Hdi | |Ered | |Y Y B X Y| |Y Y Bred Xred Y| p,h,n q 1 # of xs Note: xq h+1 , ,xq |x1 , ,xq h = = |Y Y B X Y| |Y Y n y | y |Y Y Bred Xred Y| |Y Y n y | y full red makes full vs. reduced testing simple to carry out 142 Note: , V , and U can be calculated from eigenvalues of E 1 Hdi full with s = min(p, h) m= N= 1 (|h p| 1) 2 1 (n h p 2) 2 143 Subset Selection Finding a subset of the x s to include in a model Forward Selection (Step 1) Start with 01 Bi = i1 02 i2 0p ip and calculate xi = |Y Y Bi X Y| p,1,n 2 |Y Y n y | y for i = 1, . . . , q. Add the xi that minimizes xi (as long as xi < ,p,1,n 2 stop otherwise) 144 (Step j + 1, j = 1, 2, . . .) Let x1 , . . . , xj be the variables added in previous steps. Calculate xi |x1 ,...,xj p,1,n j 1 for all xi among the q j remaining candidate variables. For the xi that minimizes xi |x1 ,...,xj : add xi if xi |x1 ,...,xj < ,p,1,n j 1 stop the procedure if xi |x1 ,...,xj > ,p,1,n j 1 Backward Elimination Start with all x s and delete one at a time until the least valuable remaining x is signi cant. For the m remaining x s after a given step, nd the xi maximizing xi |x1 ,...,xi 1 ,xi+1 ,...,xm p,1,n m 1 drop xi if xi |x1 ,...,xi 1 ,xi+1 ,...,xm > ,p,1,n m 1 stop the procedure if xi if xi |x1 ,...,xi 1 ,xi+1 ,...,xm < ,p,1,n m 1 145 Stepwise Add most signi cant candidate xi if partial is less than critical value Then, remove least signi cant selected xi if partial is greater than critical value Best Subsets Choose best subset of size , for = 1, . . . , q, with respect to some criterion (e.g., a multivariate extension of Mallow s Cp , or tr{S}, etc.) After selecting a subset of the x s, subset of y s may be selected using stepwise discriminant approach . . . to be discussed later. 146 ex chemical reaction data How are the responses (y1 , y2 , and y3 ) a ected by the inputs (x1 , x2 , and x3 )? y1 = % of unchanged starting material y2 = % converted to the desired product y3 = % of unwanted by-product x1 = temperature x2 = concentration x3 = time 0 Regress y on x to obtain B and test B1 = 0, where B = . B1 Determine what the eigenvalues of E 1 H reveal about the essential rank of B1 and the power of the 4 ( MANOVA ) statistics. NOTE: MTEST/PRINT DETAILS gives eigenvalues of i i (E + H) 1 H ( 1 , . . . , s ) and i = 1+ i , and i = 1 i 147 Check the signi cance of x1 x2 , x1 x3 , x2 x3 , x2 , x2 , and x2 1 2 3 adjusted for x1 , x2 , and x3 148 III.C.ii Seemingly Unrelated Regressions (SUR) Standard multivariate regression: Y = XB + or X y 1 . . = . 0 y p vec Y .. . e 1 1 . . . + . . . X p e p 0 vec B vec p p I X n p Note: each y 1 uses the same regressors X vec B = (I (X X) 1 X )vec Y n r or B = (X X) 1 X Y is B.L.U.E. even though var{vec Y} = I What if each y j uses di erent regressors Xj ? n rj 149 SUR Model y j n 1 = X j j + ej n rj rj 1 n 1 j = 1, . . . , p ex y j is the jth economic outcome for n regions and Xj is the matrix of economic indicators (unemployment, housing starts, etc.) and the indicators used in the model are potentially di erent for each outcome that is, Xj = Xj 150 Model assumptions: X1 n r1 y 1 y 2 = . . . y p 0 np 1 0 X2 n r2 .. . Xp n rp np ( p j=1 ( 1 2 . . . p p j=1 e 1 e2 + . . . ep np 1 rj ) 1 rj ) OR y = X + e cov{ei , ej } = ij In p p independent observations np np var{e } = In = 151 OLS estimator OLS (X X ) X1 y 1 11 1,OLS 1 . . . = X X X y = = . . . 1 X y p X Xp 1 p p p,OLS is not BLUE in general 152 Zyskind condition states that the OLS estimator is BLUE only if there exists a matrix Q such that X = X Q In standard multivariate regression, = In and X = Ip X and n r p p X = p p In Ip X n r = X = Ip X X n r p p In Q Therefore, OLS is BLUE under standard multivariate regression assumptions In SUR case, there is no simple way of writing X , and in general, there exists no Q satisfying X = X Q Use SUR = X 1 X 1 X 1 y 153 III.C.iii Canonical Correlation Analysis Objective: summarize the linear relationship between two groups of variables y = (y1 , . . . , yp ) and x = (x1 , . . . , xq ) Neither x nor y considered dependent Multivariate extension of the squared multiple correlation coe cient (used to relate a single response y with x). Consider a single random variable y and a random vector x Recall y Syy var = S = x Sxy and y Ryy corr = R = x Rxy Syx p=1 if = 1 1 Sxx sxy syy syx Sxx q q 1 p=1 if = Rxx rxy q q Ryx ryx Rxx 154 Squared multiple correlation between y and (x1 , . . . , xq ) is 2 Ry|x = syx S 1 sxy xx = ryx R 1 rxy xx syy 2 where Ry|x is themaximum correlation between y and a linear combination of the x s Extending to the case with y = (y1 , . . . , yp ) and x = (x1 , . . . , xq ), a measure of association between y and x is s 2 RM = |S 1 Syx S 1 Sxy | yy xx = i=1 2 ri (s = min(p, q)) 2 2 where r1 , . . . , rs are the eigenvalues of S 1 Syx S 1 Sxy . yy xx 2 RM too small and too heavily dependent on smallest eigenvalues. 2 2 Instead, work directly with r1 , . . . , rs called the squared canonical correlations 155 2 Largest squared canonical correlation r1 is the maximum squared correlation between a linear combination of x and a linear combination of y. 2 r1 = corr{a1 y , b1 x } = max corr{a y, b x} a,b u1 v1 u1 = a1 y and v1 = b1 x are the rst canonical variates First s eigenvalues of S 1 Sxy S 1 Syx are same as rst s eigenvalues xx yy q q of S 1 Syx S 1 Sxy , yy xx p p but eigenvectors are di erent. (S 1 Syx S 1 Sxy r2 Ip )a = 0 yy xx If q < p, only q of the eigenvectors a are meaningful (S 1 Sxy S 1 Syx r2 Iq )b = 0 xx yy If p < q, only p of the eigenvectors b are meaningful 156 The canonical correlations r1 , . . . , rs respond to the s pairs of canonical variates: the s nonredundant u1 = a1 y and v1 = b1 x dimensions of the relationship u2 = a2 y and v2 = b2 x (s = min(p, q)) . . . . . . . . . us = as y and vs = bs x ui s are uncorrelated (so are vi s) ui uncorrelated with vj for i = j If software requires a symmetric matrix to obtain eigenvalues and eigenvectors, use 1/2 S 1/2 Sxy S 1 Syx Sxx xx yy 2 2 which has eigenvalues r1 , . . . , rs and eigenvectors Sxx bi and 1/2 S 1/2 Syx S 1 Sxy Syy yy xx 2 2 which has eigenvalues r1 , . . . , rs and eigenvectors Syy ai 1/2 1/2 157 Importance of the relationship between ui and vi (that is, 2 2 importance of ri ) can be judged by the relative size of ri : s 2 j=1 rj 2 ri ai = p 1 1 1 ri Syy Syx bi p p p q q 1 1 1 ri Sxx Sxy ai q q q p bi q 1 = 158 When interpreting the canonical variates, we prefer to use standardized coe cient vectors 0 0 s y1 0 0 sy2 ci = . . . ai , . .. . . . . . . 0 0 syp and s x1 0 di = . . . 0 where syi = 0 sx2 . . . 0 .. . 0 0 . . . sxq bi var{yi } and sxi = var{xi }. More simply, conduct analysis using R 1 Ryx R 1 Rxy and yy xx 1 1 Rxx Rxy Ryy Ryx which have eigenvectors ci and di , respectively, 2 2 and have eigenvalues r1 , . . . , rs . 159 Properties of canonical correlations: 2 ri invariant to change of scale on y s or x s r1 exceeds the absolute value of the correlation between any y and any or all of the x s. 2 2 2 2 ri = Rui |x = Rvi |y (where Rui |x is the squared multiple correlation between ui and (x1 , . . . , xq )) 160 Statistical Inference: H0 : no linear relationship between y s and x s or H0 : B1 = 0 q p or H0 : independence of y and x Test statistic: |S| |R| = p,q,n 1 q |Syy ||Sxx | |Ryy ||Rxx | q,p,n 1 p since p,q,n 1 q = q,p,n 1 p [An exact F test exists for s = min(p, q) 2, see -to-F conversions in the MANOVA section] q q 161 2 2 2 = i=1 (1 ri ) is a function of r1 , . . . , rs , as are Pillai s (V ), Lawley-Hotelling (U ), and Roy s Largest root ( ) 2 As strength of relationship between x and y increases, ri s increase and decreases s Testing if the s canonical correlations (combined) are signi cant 2 Note: If p = 1, = 1 Ry|x 2 If test rejects, next consider how many ri s are signi cant H0 : The canonical correlations rm , . . . , rs are non-signi cant s m = 2 (1 ri ) p m+1,q m+1,n m q i=m or q m+1,p m+1,n m p [since p, H , E = H ,p, H + E p ] An approach: Check 2 , . . . , s to determine number of signi cant 2 ri values. q 162 Interpretation of canonical variates (e.g., u1 = a1 y and v1 = b1 x) Wish to assess the contribution of each variable to the canonical 2 correlation ri . standardized coe cients correlations between y1 and uj = aj y Standardized coe cients Use ci and di to account for di erences in scaling among the variables Absolute values of coe cients ci show contribution of each yi in the presence of the other yi s. Add or remove yi s ci changes We want this property in multivariate analysis!! 163 Correlations between yi and uj = aj y (and between xi and vj = bj x) More frequently used and widely claimed to yield more valid interpretation of canonical variates (a.k.a. structure coe cients ) corr{yi , uj } is stable (not dramatically di erent) if we add or remove yi s . . . sounds nice, but it s not! In fact, these correlations provide no information about the multivariate contribution of the variable yi to the correlation 2 structure. (Analogous to Tp -test vs. p univariate t-tests.) Rencher (1988, 1992) showed that s 2 2 (corr{yi , uj })2 rj = Ryi |x j=1 2 where Ryi |x is the multiple correlation between yi and (x1 , . . . , xq ) 164 Although corr{yi , uj } might seem to quantify the importance of yi in a multivariate relationship with x in the presence of the other y variables, it summarizes only a univariate relationship. ex chemical reaction data 2 ? How many ri s are necessary? ? Interpret u1 and v1 . Recall 2 Rxi |y = 3 22 j=1 (corr{xi , v1 }) rj . First canonical correlation is mostly due to relationship of: temperature and concentration (suppressed by x1 x2 , and to a lesser degree x1 x3 and x2 ) with: % changed. 1 165 INPUTS x1 = temperature x2 = concentration x3 = time x1 x2 x1 x3 x2 x3 x2 1 x2 2 x2 3 YIELDS y1 = % unchanged y2 = % converted y3 = % by-product ** * ** * ** ** d1i 5.01 5.86 1.65 -3.92 -2.30 0.53 -2.67 -1.23 0.57 corr{xi , v1 } .69 .23 .45 .41 .54 .45 .69 .23 .42 corr{yi , u1 } -.996 .64 .85 166 2 Rxi |y .69 .24 .51 .43 .58 .48 .69 .23 .47 2 Ryi |x c1i -1.54 -0.21 -0.47 .987 .92 .91

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: