1 fca <- function(xtab) {
3 # Correspondence analysis of principal table. List returned with
4 # projections, correlations, and contributions of rows (observations),
5 # and columns (attributes). Eigenvalues are output to display device.
10 fI <- apply(fIJ, 1, sum)
11 fJ <- apply(fIJ, 2, sum)
12 fJsupI <- sweep(fIJ, 1, fI, FUN="/")
13 #fIsupJ <- sweep(fIJ, 2, fJ, FUN="/")
14 s <- as.matrix(t(fJsupI)) %*% as.matrix(fIJ)
15 s1 <- sweep(s, 1, sqrt(fJ), FUN="/")
16 s2 <- sweep(s1, 2, sqrt(fJ), FUN="/")
17 # In following s2 is symmetric. However due to precision S-Plus didn't
18 # find it to be symmetric. And function eigen in S-Plus uses a different
19 # normalization for the non-symmetric case (in the case of some data)!
20 sres <- eigen(s2,symmetric=T)
21 sres$values[sres$values < 1.0e-8] <- 0.0
22 factexpl<-sres$values[-1]
23 #tot <- sum(sres$values[-1])
24 #factexplp<-100*sres$values[-1]/tot
25 # Eigenvectors divided rowwise by sqrt(fJ):
26 evectors <- sweep(sres$vectors, 1, sqrt(fJ), FUN="/")
28 # PROJECTIONS ON FACTORS OF ROWS AND COLUMNS
29 rproj <- as.matrix(fJsupI) %*% evectors
30 #temp <- as.matrix(s2) %*% sres$vectors
31 # Following divides rowwise by sqrt(fJ) and columnwise by sqrt(eigenvalues):
32 # Note: first column of cproj is trivially 1-valued.
33 # NOTE: VBESxFACTORS. READ PROJS WITH FACTORS 1,2,... FROM COLS 2,3,...
34 #cproj <- sweep(sweep(temp,1,sqrt(fJ),FUN="/"),2,sqrt(sres$values),FUN="/")
36 # CONTRIBUTIONS TO FACTORS BY ROWS AND COLUMNS
37 # Contributions: mass times projection distance squared.
38 #temp <- sweep( rproj^2, 1, fI, FUN="*")
39 # Normalize such that sum of contributions for a factor equals 1.
40 #sumCtrF <- apply(temp, 2, sum)
41 # Note: Obs. x factors. Read cntrs. with factors 1,2,... from cols. 2,3,...
42 #rcntr <- sweep(temp, 2, sumCtrF, FUN="/")
43 #temp <- sweep( cproj^2, 1, fJ, FUN="*")
44 #sumCtrF <- apply(temp, 2, sum)
45 # Note: Vbs. x factors. Read cntrs. with factors 1,2,... from cols. 2,3,...
46 #ccntr <- sweep(temp, 2, sumCtrF, FUN="/")
48 # CORRELATIONS WITH FACTORS BY ROWS AND COLUMNS
49 # dstsq(i) = sum_j 1/fj (fj^i - fj)^2
50 #temp <- sweep(fJsupI, 2, fJ, "-")
51 #dstsq <- apply( sweep( temp^2, 2, fJ, "/"), 1, sum)
52 # NOTE: Obs. x factors. Read corrs. with factors 1,2,... from cols. 2,3,...
53 #rcorr <- sweep(rproj^2, 1, dstsq, FUN="/")
54 #temp <- sweep(fIsupJ, 1, fI, "-")
55 #dstsq <- apply( sweep( temp^2, 1, fI, "/"), 2, sum)
56 # NOTE: Vbs. x factors. Read corrs. with factors 1,2,... from cols. 2,3,...
57 #ccorr <- sweep(cproj^2, 1, dstsq, "/")
59 # Value of this function on return: list containing projections, correlations,
60 # and contributions for rows (observations), and for columns (variables).
61 # In all cases, allow for first trivial first eigenvector.
62 list(rproj=rproj[,-1],factexpl=factexpl)#,rcorr=rcorr[,-1], rcntr=rcntr[,-1],factexpl=factexpl,factexplp=factexplp,
63 #cproj=cproj[,-1], ccorr=ccorr[,-1], ccntr=ccntr[,-1])