忍者ブログ
9

×

[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。

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


PR

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

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

#dat <- Orange
for(i in 1:ncol(dat)) {
x <- paste(colnames(dat)[i], " <- c(", paste( dat[,i], collapse=","), ")", sep="")
print(x)
}

# あまり美しくない出力
xs <- vector()
#dat <- Orange
for(i in 1:ncol(dat)) {
x <- paste(" <- c(", paste( dat[,i], collapse=","), ")", sep="")
xs[i] <- x
}
names(xs) <- colnames(dat)
xd <- data.frame(xs)
print(xd, right=F)
Rで対応分析…の準備
CRAN Task Viewのまとめ http://cran.r-project.org/web/views/Psychometrics.html

# 対応パッケージと関数
 # MASSパッケージのcorresp, mca
  # 金先生の連載 http://mjin.doshisha.ac.jp/R/26/26.html
  # 多重対応分析mca関数の例 http://133.100.216.71/R_analysis/analysis_corresp02.html
  # 対応分析corresp関数の例 http://133.100.216.71/R_analysis/analysis_corresp01.html
  # Rで学ぶデータサイエンス1 カテゴリカルデータ解析 第9章 対応分析 p130もMASSパッケージの関数の解説


 # caパッケージのca, mjca
  # cran: http://cran.r-project.org/web/packages/ca/index.html マニュアルpdf http://cran.r-project.org/web/packages/ca/ca.pdf
  # Nenadic, O., & Greenacre, M. (2007). Correspondence Analysis in R, with Two- and Three-dimensional Graphics: The ca Package. Journal of Statistical Software, 20, 1-13.  http://www.jstatsoft.org/v20/a03/paper
  # Quick-Rの実行例: http://www.statmethods.net/advstats/ca.html
  # CARME-N.ORG: http://www.carme-n.org/?sec=ca
  # 中山 (2009). 対応分析によるデータ解析 http://www.kwansei.ac.jp/s_sociology/attached/6779_56340_ref.pdf


 # FactMineRパッケージのCA, MCA関数
  # cran: http://cran.r-project.org/web/packages/FactoMineR/index.html マニュアルpdf http://cran.r-project.org/web/packages/FactoMineR/FactoMineR.pdf
  # 本家HP: http://factominer.free.fr/book/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


 # homalsパッケージのhomals関数
  # cran: http://cran.r-project.org/web/packages/homals/ マニュアルpdf http://cran.r-project.org/web/packages/homals/homals.pdf
  # Leeuw, J., & Mair, P. (2009). Gifi methods for optimal scaling in R: The package homals. Journal of Statistical Software, 31, ... http://www.jstatsoft.org/v31/i04/paper
    # これも中身は同じ Leeuw, J., & Mair, P. (2008). Homogeneity analysis in R: The package homals. Journal of Statistical Software http://preprints.stat.ucla.edu/525/homalsR.pdf

# その他
 # 正準対応分析? http://preprints.stat.ucla.edu/531/anacor.pdf
 # 青木先生の数量化3類の関数 http://aoki2.si.gunma-u.ac.jp/R/qt3.html
# ダミー変数を生成するmake.dummy関数は大変便利。感謝。source("http://aoki2.si.gunma-u.ac.jp/R/src/make.dummy.R", encoding="euc-jp")

## 覚書
 # 足立浩平 (2006). 多変量データ解析法――心理・教育・社会系のための入門―― ナカニシヤ出版
 # 第13章 数量化分析 p125 より
  # 相対尺度法 (dual scaling) 、最適尺度法 (optimal scaling) 、等質性分析 (homogenuity analysis) 、対応分析 (corespondence analysis) 、多重対応分析・多重応答分析・同コレスポンデンス分析 (multiple correspondence analysis) 、数量化3類 (quantification method III, third method of quantification) 、主成分尺度分析 (principal components of scale analysis) 、質的データの要因分析 (factorial analysis of qualitative data) は本質的に同じ方法である。
  # 対応分析と多重対応分析は少し性質が異なる
  # 対応分析は分割表を直接解析対象にする。標準化比率にしたデータに主成分分析をすることに等しい
    # カテゴリカル主成分分析も同じようなもの (?) 。ダミー変数を使っているだけ (?)
    # カテゴリカル回帰は個々の因子変数を数量化した線形重回帰 (?)
  # 数量化III類 = カテゴリカル変数の主成分分析 = 対応分析
    # 数量化I類= 独立変数が連続変数で、従属変数がカテゴリカル変数。説明変数 (独立変数) がカテゴリカル変数のときの重回帰と同じ
    # 数量化II類 = 独立変数がカテゴリカル変数で、従属変数もカテゴリカル変数。説明変数 (独立変数) がカテゴリカル変数のときの判別分析と同じ
 # 多重対応分析と数量化3類は数値結果が異なる。たぶん等質性分析も異なる
  # 多重対応分析はYesを1、No を1のデータをつくる これをそのまま主成分分析しても同じ個体得点が得られる
  # 数量化はYesを1、Noを0のデータをつくる
 # 基本的な目的と解釈
  # 類似したパターンのサンプル (回答者) 、類似したパターンのカテゴリー (項目) を近くに移動する
  # 近くに位置づけられるカテゴリーどうし・個体どうしは互いに等質、似ている
    # 40代の人は健康に関心がある
  # 次元の解釈
    # 次元の両極に位置するカテゴリーを比較する

## Rでカテゴリカル分析のリンク
http://www.unt.edu/rss/class/Jon/R_SC/
  # かなりの充実度。海外はすごいな


プロフィール
HN:
tao
HP:
性別:
非公開
職業:
会社員
趣味:
アウトドア、自転車、ジョギング、英語学習
自己紹介:
・千葉在住のサラリーマンです。データ分析っぽいことが仕事。
・今年英検1級取得。今はTOEIC高得点を目指して勉強中。
・興味のあることは野球、アウトドア、英語学習、統計、プログラミング、PC関係などなど。
ブログ内検索
freead
順位表
プロ野球データFreak
セリーグ順位表
パリーグ順位表