compatibility R 4.0
[iramuteq] / Rscripts / afc.R
1 fca <- function(xtab) {
2         
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.
6 # FM, 2003/12.
7         
8         tot <- sum(xtab)
9         fIJ <- xtab/tot
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="/")
27         
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="/")
35         
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="/")
47         
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, "/")
58         
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])
64         
65 }