myscratchpad

# 参考URL。感謝
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/07_folder/da07_02.html


cr <- c(1,1,1,0,1,0,0,1,0,0)
ex <- c(4,5,3,2,4,3,2,4,3,1)
ky <- c(3,4,4,5,4,2,3,5,2,2)
ke <- c(2,3,2,2,4,3,4,5,3,1)
dat <- data.frame(cr, ex, ky, ke)


青木先生の関数
source("http://aoki2.si.gunma-u.ac.jp/R/src/disc.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/geneig.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/candis.R", encoding="euc-jp")
resd <- disc(dat[2:4], dat[1]) # 線形判別の関数
resc <- candis(dat[2:4], dat[1]) # 正準判別の関数

# MASSパッケージのlda関数
library(MASS)
ldares <- lda(cr~., data=dat)
ldares # 基本的に教科書の結果とは正負が逆になる。以下も同様

# 切片は自分で計算
apply(ldares$means %*% ldares$scaling, 2, mean)
  # 分解するとこんな感じ
  gm <- ldares$mean
  cf <- ldares$scaling
  x <- gm %*% cf
  apply(x, 2, mean)

# Wilks' lambda
 summary(manova(cbind(ex, ky, ke)~cr, data=dat), test="Wilks")
 # 複数の連続変数→カテゴリカル(2値) の一般線形モデルによる分析が判別分析で、カテゴリカル→複数の連続変数　の一般線形モデルの分析がmanovaにあたる

# 判別得点など
prres <- predict(ldares)
 # 判別結果
 prres$class
 cr # 元の群
 (x <- data.frame(prres$class, cr)) # 並べてみる
 t(table(x)) # 判別結果と所属群
 # 事後確率
  prres$posterior
 # 判別得点
 prres$x
 ## 表15.1 (B) p150
 data.frame("判別得点"=prres$x[,1], "判別結果"=prres$class)


# 交差妥当性
ldacv <- lda(cr~., data=dat, CV=T)
pcr <- predict(ldares)$class
cvpcr <- ldacv$class
table(pcr, cr)
 # 判別率・誤判別率
cprp(table(pcr, cr))$pt # prop.tableでもいい？
table(cvpcr, cr) # 交差妥当性の確認
 cprp(table(cvpcr, cr))$pt

# 太郎丸　博 (2005). 人文・社会科学のためのカテゴリカル・データ解析入門　ナカニシヤ出版

## p135
mat <- matrix(c(36,13,5,10,2,8), nr=2, dimnames=list(hus=c("jh",
"cu"), wif=c("wjh", "wc", "wu")))
xtb<- as.table(mat)
xtb
 library(MASS)
 lgres <- loglm(formula=~hus+wif, data=xtb)

 coef(lgres) # 係数
 fitted(lgres) # 期待値
 residuals(lgres) # 標準化残差
 (xtb-fitted(lgres))/ sqrt(fitted(lgres)) # ピアソンタイプの残差？
 xtb # 観測値

 lgres # 検定

# 組み込みのloglin関数
 linres <- loglin(xtb, margin=list(c(1), c(2)), param=T, fit=T)
 linres$para
 linres$fit
 linres$lrt # 尤度比
 linres$pearson # カイ二乗
 linres$df # 自由度
  # 検定
  1-pchisq(linres$lrt, linres$df)
  1-pchisq(linres$pearson, linres$df)

# 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))

library(FactoMineR)
data (poison)
dat <- poison[,-c(1,2)]

# 各変数の度数
sapply(dat, table)

MCres <- MCA(dat, graph=F)
MCres
MCres$eig # 固有値と説明率
MCres$var$coord[,1:2] # カテゴリースコア。spssだと重心座標。符号が入れ替わる
MCres$var$eta2 # 判別測定
 colSums(MCres$var$eta2) # 判別測定の合計
MCres$ind$coord[,1:2] # 個人スコア。spssとは結構違う。
 # カテゴリースコアで数量化された

 # 個人プロット
plot.MCA(MCres,invisible=c("var","quali.sup"))
 MCres$ind$coord[23,1:2]

 # カテゴリープロット
plot.MCA(MCres,invisible="ind")

 # バイプロット
plot.MCA(MCres)

# 個体スコアについて

x <- as.vector(as.matrix(dat[1,]))
cs <- MCres$var$coord[,1:2] # カテゴリースコア
dm <- colSums(MCres$var$eta2) # 判別測定の合計

# 個体1のカテゴリースコア
cs[x[1],]
cs[x[2],]
cs[x[3],]
 # 全部足す
(cs[x[1],]+cs[x[2],]+cs[x[3],])/dm[1:2] # spssではこれ
MCres$ind$coord[1,1:2]
 # 個体スコアの計算方法が異なることがわかる

# そろそろspssを気にしないことにしよう

参考。感謝
http://www.ec.kagawa-u.ac.jp/~hori/mva.html
http://www.ec.kagawa-u.ac.jp/~hori/chosadata/homalskom.htm

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(FactoMineR)
MCres <- MCA(dat, graph=F)
MCres
summary(MCres) # 何も出力されない
MCres$eig # 固有値と説明率
MCres$var$coord[,1:2] # カテゴリースコア
MCres$ind$coord[,1:2] # 個人スコア
plot(MCres)
 # 結果がみづらいな
 ## spssの結果とほぼ一致するのはこれ ##


library(homals)
hres <- homals(dat)
hres # 固有値を10倍するとspssの結果にほぼ一致する
summary(hres)
 # Category quantifications を10倍するとspssの結果にほぼ一致する
plot(hres)
 # プロットはみづらい


library(MASS)
Mres <- mca(dat)
Mres
Mres$cs
Mres$rs
 # いまいち見づらい、出力が少ない
plot(Mres) # みづらい。biplotに対応していない

library(ca)
mjres <- mjca(dat)
mjres
summary(mjres)
mjres$colcoord*mjres$sv # カテゴリースコア
mjres$rowcoord*mjres$sv # 個人スコア
par(mfrow=c(2,2))
plot(mjres)
plot(mjres,what=c("none","all")) # カテゴリーのプロット
plot(mjres,what=c("all","none")) # 個人のプロット
plot(mjres,what=c("none","all"),arrows=c(FALSE,TRUE)) # ベクトルのプロット
 # 結果オブジェクトの解釈にもっと勉強が必要

##############http://factominer.free.fr/index.html
# とりあえずFactoMineRを勉強してみよう
  # 本家HP: http://factominer.free.fr/index.html
   # 対応分析の例: http://factominer.free.fr/classical-methods/correspondence-analysis.html
   # 多重対応分析の例: http://factominer.free.fr/classical-methods/multiple-correspondence-analysis.html
  # Le, S., Josse, J., & Husson, F. (2008). FactoMineR: An R Package for Multivariate Analysis.  Journal of Statistical Software, 25, 1-18. http://factominer.free.fr/docs/article_FactoMineR.pdf

# 検索するとフランス語ばっかりでてくるな。どうしよう