rm(list = ls()) # clear the memory ########################### ## Example 10.1 in text # Enter or compute covariances: # Example 10.1: S11 = cbind(c(1,.4), c(.4,1)) S12 = cbind(c(.5, .3), c(.6, .4)) S21 = t(S12) S22 = cbind(c(1,.2), c(.2, 1)) p = nrow(S11) q = nrow(S22) # Get symmetric square roots: s1 = eigen(S11) rootS11 = s1\$vectors%*%diag(sqrt(s1\$values))%*%t(s1\$vectors) s2 = eigen(S22) rootS22 = s2\$vectors%*%diag(sqrt(s2\$values))%*%t(s2\$vectors) S12tilde = solve(rootS11, S12)%*%solve(rootS22) s = svd(S12tilde, nu = p, nv = q) G = s\$u H = s\$v R = cbind(diag(s\$d),matrix(nrow=p, data = rep(0, p*(q-p)))) #check: round(S12tilde - G%*%R%*%t(H),5) # Compute matrices Atilde and Btilde of (unnormalized) canonical coefficients Atilde = t(solve(rootS11,G)) Btilde = t(solve(rootS22,H)) Btilde = Btilde[1:p,1:q] # Normalize these so their rows have unit norm: Anorms = sqrt(diag(Atilde%*%t(Atilde))) A = solve(diag(Anorms),Atilde) Bnorms = sqrt(diag(Btilde%*%t(Btilde))) B = solve(diag(Bnorms),Btilde) round(cbind(A, B),2) # coefficients of the canonical pairs; unit norms round(1/diag(Atilde%*%t(Atilde)),2) round(1/diag(Btilde%*%t(Btilde)),2) # Canonical correlations: round(diag(R),2) ################### job satisfaction ############ # Triangular data in E10-5.DAT L = matrix(0, 12, 12) L[outer(1:12, 1:12, '<=')] <- scan("http://www.stat.ualberta.ca/~wiens/stat575/datasets/E10-5.DAT") S = L + t(L) - diag(diag(L)) rownames(S) = NULL colnames(S) = NULL p = 5 q = 7 S11 = S[1:p, 1:p] S12 = S[1:p, (p+1):(p+q)] S21 = t(S12) S22 = S[(p+1):(p+q),(p+1):(p+q)] # Get symmetric square roots: s1 = eigen(S11) rootS11 = s1\$vectors%*%diag(sqrt(s1\$values))%*%t(s1\$vectors) s2 = eigen(S22) rootS22 = s2\$vectors%*%diag(sqrt(s2\$values))%*%t(s2\$vectors) S12tilde = solve(rootS11, S12)%*%solve(rootS22) s = svd(S12tilde, nu = p, nv = q) G = -s\$u H = -s\$v R = cbind(diag(s\$d),matrix(nrow=p, data = rep(0, p*(q-p)))) # Compute matrices Atilde and Btilde of (unnormalized) canonical coefficients Atilde = t(solve(rootS11,G)) Btilde = t(solve(rootS22,H)) Btilde = Btilde[1:p,1:q] # Normalize: Anorms = sqrt(diag(Atilde%*%t(Atilde))) A = solve(diag(Anorms),Atilde) Bnorms = sqrt(diag(Btilde%*%t(Btilde))) B = solve(diag(Bnorms),Btilde) round(cbind(A, B),2) # coefficients; unit norms # Canonical correlations: round(diag(R),2) corUX = Atilde%*%S11/sqrt(diag(S11)) corUY = Atilde%*%S12/sqrt(diag(S22)) corVX = Btilde%*%S21/sqrt(diag(S11)) corVY = Btilde%*%S22/sqrt(diag(S22)) corU1X = corUX[1,] corV1X = corVX[1,] corU1Y = corUY[1,] corV1Y = corVY[1,] round(corU1X,2) round(corV1X,2) round(corU1Y,2) round(corV1Y,2) # LR tests pval = vector(length = p) chisq = vector(length = p) df = vector(length = p) rho = diag(R) for (k in 0:(p-1)) { n = 784 deg = (p-k)*(q-k) chisquare = -(n-1-(p+q+1)/2)*log(prod((1-rho^2)[(k+1):p])) pvalue = 1 - pchisq(chisquare, df = deg) pval[k+1] = pvalue df[k+1] = deg chisq[k+1] = chisquare } num.sig = 0:(p-1) round(cbind(num.sig,chisq, df, pval),2)