dual <- function(tbl)
{
nr <- nrow(tbl)
nc <- ncol(tbl)
ft <- sum(tbl)
fr <- apply(tbl, 1, sum)
fc <- apply(tbl, 2, sum)
temp <- sqrt(outer(fc, fc))
w <- t(tbl)%*%(tbl/matrix(fr, nr, nc))/temp-temp/ft
res <- eigen(w)
ne <- length(res$values[res$values > 1e-5])
vec <- res$vectors[,1:ne]
val <- res$values[1:ne]
names(val) <- paste("Axis", 1:ne)
col.weight <- vec*matrix(sqrt(ft/fc),nc, ne)
row.weight <- tbl%*%col.weight/outer(fr, sqrt(val))
cont <- val/sum(diag(w))*100
chi.sq <- -(ft-1-(nr+nc-1)/2)*log(1-val)
df <- nr+nc-1-2*(1:ne)
result <- matrix(c(val, sqrt(val), cont, cumsum(cont), chi.sq, df, pchisq(chi.sq, df, lower=F)), byrow=TRUE, ncol=ne)
rownames(result) <- c("eta square", "correlation", "contribution", "cumulative contribution", "Chi square value", "degrees of freedom", "P value")
colnames(result) <- colnames(row.weight) <- colnames(col.weight) <- paste("Axis", 1:ne)
rownames(row.weight) <- paste("Row", 1:nr)
rownames(col.weight) <- paste("Col", 1:nc)
list(result=result, "row weight"=row.weight, "column weight"=col.weight, "row weight(weighted by eta square)"=row.weight*sqrt(matrix(val,byrow=TRUE,nr,ne)), "column weight(weighted by eta square)"=col.weight*sqrt(matrix(val,byrow=TRUE,nc,ne)))
}
# 使用例
tbl <- matrix(c(
2, 3, 5, 6,
5, 1, 7, 5,
5, 3, 4, 3
), byrow=TRUE, ncol=4)
dual(tbl)
$result
Axis 1 Axis 2
eta square 0.04905396 0.03765652
correlation 0.22148129 0.19405288
contribution 56.57212279 43.42787721
cumulative contribution 56.57212279 100.00000000
Chi square value 2.26340816 1.72727308
degrees of freedom 4.00000000 2.00000000
P value 0.68743886 0.42162603
$"row weight"
Axis 1 Axis 2
Row 1 1.4355847 -0.03995533
Row 2 -0.6650934 1.13131472
Row 3 -0.7331783 -1.31495864