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