サイトの紹介です．
日本社会心理学会広報委員会では，このたびの災害に対して被災者や非被災者，
あるいはさまざまな組織がいかに対応するべきかについて，社会心理学を研究する
個人あるいは関連団体が自身を含むこれまでの研究成果にもとづいて提言を
おこなっているサイトの情報を集めたリンク集を作成・公開することになりました．
東北地方太平洋沖地震関連ページ（日本社会心理学会）
このページは被災者やその支援者向けというだけでなく，地震の直接的影響のない一般市民が地震後の生活の中で起こる不安や疑問にも答えるものとなっています．
サイトは整備中で，まだまだ内容は充実していないかもしれませんが，ぜひご参考にしていただければと思います
Rでコレスポンデンス分析を行う話の続き．カテゴリー・ウェイト行列では直感的な理解とならないので，やはりこういった分析の時にはプロットをしましょう，というのが今回のお題．
> postscript(file="plot_mca.eps",family="Japan1GothicBBB")# 日本語を表示
> plot(out.mca$cs[,1], out.mca$cs[,2], pch=".",
+ xlab="成分1", ylab="成分2",
+ main="タイトル")
> text(out.mca$cs[,1], out.mca$cs[,2], rownames(out.mca$cs)) # アイテム-カテゴリー名を表示
> dev.off()
当方はlinux上で分析しているので，図の出力はEPSにした．他のOSの方は適当な出力デバイスを開いてください．family="Japan1GothicBBB"という引数を追加しているのは，自分の環境では日本語フォントが化けないようにするため．
デフォルトで2成分しか値を求めていないなら関数biplot()が簡単でいいです．今回は3成分なので少し変える．
結果のオブジェクトはデータフレームになっていて，列データを表すcsもデータフレームである．なので，変数名や列番号を駆使して，プロットに必要なデータを指示してあげればよい．成分1と成分3のプロットなども，これを応用するだけである．
ここでは，関数plot()は，xlab, ylab, mainなどの引数を入力して，それぞれ横軸ラベル，縦軸ラベル，図のタイトルを指定した(上の関数plot()は3行で1つの関数となっている．カッコの対応関係に注意．関数は数行にわたって書いてもいいので，長ったらしい場合は複数行に分けて書くことを積極的に薦める．見やすさが全然違うから，後で再利用する際にも便利)．残った引数pchでは，プロットに用いる記号を指定している．今回のように pch="." としておくと，プロット記号が目立たなくなるので，変数ラベルを表示させたい時には便利だ．覚えておくとよい引数の一つである．
変数ラベルを表示するために，関数text()を追加した．3つの引数は，順に，ラベル表示位置のx座標，ラベル表示位置のy座標，表示するテキスト，である．
関数plot()の引数fileで，出力はファイルに保存するようにしたが，この引数がなければ別ウィンドウが開いて図が表示されるはずである．最後はdev.off()とし，出力デバイスを忘れずに閉じておいたほうがいいようだ．
この前，コレスポンデンス分析のやり方について尋ねられた．以前に書いた準備編の続きを放っておいたのをはっと思い出した．いろいろやらなければいけないことをおいとおいて，続きをまとめておこうと思う(実際のところ，Rによる多変量解析の解説本とたいして違う話が書けるわけではないのだが)．
今回は関数mca()を使うことする．分析のための下準備は，以前に書いた準備編を参考のこと．関数corresp()もデータの形式がちがうだけで同じようにできるでしょう．
関数mca()を使う前に，MASSパッケージを呼び出しておく．関数mca()はデフォルトで2成分の結果を返す．それでもいいが，今回は必要な成分数をこちらで決めるため，最大の成分数で解いてみる．最大の成分数はアイテム-カテゴリー数-1となる．これは自分で数えて引数nfに入れる．
> out.mca <- mca(dat, nf=67) # デフォルトは2成分抽出．3成分以上を抽出するときには
> # nf=で設定．まずはアイテム数の成分を求め寄与率から
> # 必要な成分数を探る
> summary(out.mca)
Length Class Mode
rs 1716 -none- numeric
cs 204 -none- numeric
fs 1716 -none- numeric
d 3 -none- numeric
p 1 -none- numeric
call 3 -none- call
> out.mca$d # 正準相関係数
[1] 4.921493e-01 2.769077e-01 2.518887e-01 2.189185e-01 2.047861e-01
(以下，略)
> mca.eig <- out.mca$d^2 ; print(mca.eig)# 固有値を求める(正準相関係数の2乗)
[1] 2.422109e-01 7.667789e-02 6.344793e-02 4.792532e-02 4.193733e-02
(以下，略)
> kiyoritsu <- round(100 * mca.eig / sum(mca.eig) , 2); print(kiyoritsu)
[1] 20.26 6.41 5.31 4.01 3.51 3.34 3.27 3.14 2.96 2.78 2.68 2.62
(以下，略)
Rでは，分析結果は一旦適当な名前のオブジェクトに代入しておくと良い．どのような結果が返ってくるかは，関数summary()で知ることができる．rsは行スコア(各回答者別のスコア)，csは列スコア(各アイテム-カテゴリーごとのスコア)，dは正準相関係数である．正準相関係数の2乗が固有値となるので，これをもとに固有値と寄与率を求めた．
今回は，3成分で解くことにする．
> # 3成分抽出とする
> out.mca <- mca(dat04, nf=3)
> out.mca$cs # 各アイテム-カテゴリーごとのウェイト
1 2 3
Q1.男性 2.319587e-04 -2.347945e-03 7.938858e-04
Q1.女性 -6.569684e-05 8.633984e-04 -3.224514e-04
(以下，略)
回転とかはないので，固有値などの値は前のものと変わらない．
あと一息だけど，やはり長くなったので，続く．
> dat[[1]] <- factor(dat[[1]], levels=1:2)データフレーム内の変数は変数名を使って，たとえばdat$sexなどと表記してもよい．関数factor()がファクタへの変換を宣言し，関数levels()で変数の水準に対応したラベルをふっていく．
> levels(dat[[1]]) <- c("男性", "女性")
> # リッカート尺度は,forループを用いて,簡潔に書くあらかじめ，繰り返す列番号を確認しておき，その数字を利用する．上の場合だと12番目から45番目までの変数に対して同じ変換を繰り返すことになる．
> for (i in 12:45) {
> dat[[i]] <- factor(dat[[i]], levels=1:4)
> levels(dat[[i]]) <- c("全くない", "それほどない", "そう思う", "強く思う")
> }
> # リッカート尺度は,forループを用いて,簡潔に書く．レベルが減っていくことに注意
> for (i in 6:32) {
> levels(dat[[i]])[1] <- "思わない"
> levels(dat[[i]])[2] <- "思わない"
> levels(dat[[i]])[3] <- "そう思う" # レベル1と2が統合されることによって最大値が3になっている
> }
上の例では，4件法を「そう思う」「思わない」の2水準にまとめている．
注意したいのは，水準1と2をまとめてしまうと，この時点で水準は3になっている点である．最後は levels(dat[[i]])[4] <- "そう思う" と書きたくなってしまうが，それではミスとなる．
これだけ準備してようやくコレスポンデンス分析にとりかかることができる．が，ちょっと長い記事になったので，今日はここまで．
でも，実際の分析については，いろいろ書籍やウェブみれば分かるので，実はここまでの作業のほうが貴重な情報なのかもしれない．
> dat <- read.csv("xxx.csv, header=T)これは標準化残差のようだ．
> table(dat$var02, dat$var02)
A1 A2 A3 A4
B1 7 11 10 16
B2 52 55 41 48
B3 43 43 70 67
B4 16 18 33 42
> ct <- chisq.test(dat$var02, dat$var01)
ピアソンのカイ二乗検定（連続性補正なし）
データ: dat$var02 と dat$var01
カイ二乗値 = 21.6707, 自由度 = 9, P値 = 0.009983
> summary(ct)
Length Class Mode
statistic 1 -none- numeric
parameter 1 -none- numeric
p.value 1 -none- numeric
method 1 -none- character
data.name 1 -none- character
observed 16 table numeric
expected 16 -none- numeric
residuals 16 table numeric
> ct$residuals
dat$var01
dat$var02 A1 A2 A3 A4
B1 -0.68936794 0.39377346 -0.53638819 0.73802881
B2 1.81898217 1.74062507 -1.62015922 -1.46502575
B3 -0.44282457 -0.92549477 1.28561709 -0.05428336
B4 -1.36779699 -1.26051514 0.67448868 1.57327321
Vij = { 1 - ( n.j / N ) } { 1 - ( ni. / N ) }n.j, ni. は，各行・列の周辺度数，Nは全度数である．
> c.phfreq <- margin.table(ct$observed, 1) # 観測度数表から周辺度数を計算求めたc.phfreq, r.phfreqは，ベクトルになっているから，これを行列とみなして，[4行1列]の行列と[1行4列]の行列積を計算すればいいんじゃないか．
> r.phfreq <- margin.table(ct$observed, 2)
> N.ct <- sum(ct$observed)
> V.ct <- (1-(c.phfreq/N.ct)) %*% (1-(t(r.phfreq/N.ct))) # 各セルの分散を求める行列の転置を行うt()関数を使っているのがミソ．行列積は %*% で表す．
> V.ct
dat$var01
dat$var02 A1 A2 A3 A4
B1 0.7326520 0.7181280 0.6745562 0.6438946
B2 0.5217370 0.5113942 0.4803658 0.4585310
B3 0.4842718 0.4746717 0.4458714 0.4256046
B4 0.6424581 0.6297221 0.5915143 0.5646272
> Adj.ct <- ct$residuals / sqrt(V.ct) # 調整化残差を求める
> Adj.ct
dat$var01
dat$var02 A1 A2 A3 A4
B1 -0.8053826 0.4646709 -0.6530856 0.9197418
B2 2.5182719 2.4340379 -2.3376079 -2.1635200
B3 -0.6363369 -1.3433132 1.9253372 -0.0832077
B4 -1.7064723 -1.5884502 0.8769848 2.0937423
