iramuteq.org Git - iramuteq/blob - 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 }