myscratchpad

# homalsも勉強してみる
# http://www.jstatsoft.org/v31/i04/paper
# http://www.amstat.org/sections/srms/Proceedings/y2008/Files/301376.pdf
# http://www.jstatsoft.org/v32/i09/paper

a1  <- factor(c(1,1,2,3,1,2,2,1,3,1,3,2,2,1,3,2,1,2,3,2))
a2  <- factor(c(1,3,1,4,4,3,1,3,2,1,1,3,1,3,2,4,2,1,1,2))
a3  <- factor(c(1,2,1,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1))
a4  <- factor(c(2,1,1,2,1,2,1,2,2,1,1,2,2,2,1,1,1,1,2,1))
a5  <- factor(c(1,2,1,2,1,2,1,2,3,2,3,1,2,1,3,3,1,3,1,2))
dat <- data.frame(a1,a2,a3,a4,a5)

library(homals)
hres <- homals(dat)
hres # 固有値を10倍するとspssの結果にほぼ一致する
summary(hres)

# カテゴリースコア
hres$catscores

# オブジェクトスコア
hres$objscores

# プロット
 # オブジェクトプロット
 plot(hres, plot.type="objplot")
 # バイプロット
 plot(hres, plot.type="jointplot") # デフォルト plot(hres) でもいい

# カテゴリープロットは自分でやる
x <- hres$catscore # リストオブジェクト
cl <- rep(1:length(x), time=unlist(lapply(x, nrow)))
d1 <- unlist(lapply(x, function(a) a[,1]))
d2 <- unlist(lapply(x, function(a) a[,2]))
plot(d1, d2, col=cl, pch=cl)
text(d1, d2+0.005, names(d1))