Rmanual 最近の更新履歴 緒賀研究室―個人別荘

心理学データ解析のための統計ソフト R のミニミニマニュアル

2007/1/11

目次

第１部：前書き... 2

R の初歩的操作...3

・データの入力...3

＜別のデータのやりとり方法＞... 4

person.data を利用して... 5

平均値の算出... 6

枝葉図とヒストグラムで分布を見よう... 7

要約統計量... 8

箱ひげ図を描いてみよう... 8

変数を選択してまとめよう...8

散布図を描こう... 9

相関係数...9

describe 関数を導入してみよう...9

by 関数を用いる... 10

subset 関数で分割したデータセットをつくる... 11

男女別データセットを作ろう... 12

群わけしたデータセットを作ろう・・変数の変換... 12

メモ：２つの変数から群わけするための変数をつくる方法... 12

t 検定... 13

ピアソンの相関係数... 13

相関係数の差の検定...14

散布図... 15

相関係数の検定... 15

クラスター分析...15

平均を示す図を描こう... 17

一要因の分散分析...17

多重比較！... 19

ホルムの多重比較...19

・回帰分析... 22

標準偏回帰係数について... 23

参考：...24

～～参考：統計処理の羅針盤～...27

第２部：データの入力から、尺度の構成までのミニミニヒント... 27

・質問紙を回収したら～...27

～～補足：エクセルのアドインを使っての因子分析... 28

尺度の構成へ... 29

せっかくですから R を使って因子分析... 29

尺度の信頼性の検討 ... 41

尺度得点の出し方...42

Rcommander を使ってみよう...45

第３部：クロス集計表の分析... 46

お勧め本！...46

比率の検定... 46

クロス集計の際のフィッシャーの直接確率検定...48

マクネマーの有意変化の検定... 50

心理学の実験や観察で、でてきる評点間一致について...50

＜おまけです＞... 51

分散分析について...51

リンク...51

一要因の分散分析... 52

～～ちょっと細かい話～～...55

２要因の分散分析...60

多重比較... 65

２要因の分散分析(混合計画）...72

～～～付記：重要な事柄～～...76

第４部：R でｔ検定...77

対応のある t 検定...81

参考リンク...84

第１部：前書き

おかしなところがあったらぜひ緒賀( s_oga@cc.gifu-u.ac.jp )までご連絡おねがいします。

もっとも書き散らかしていますので、おかしなところばかりですが・・・

フリーの統計ソフト R を使って心理学の統計分析をしていくやり方のヒント集です。

これは、

http://personality-project.org/r/

Using R for psychological research:A very simple guide to a very elegant package

にもっぱら基づいています。

==

統計の初歩は、以下のページを見て勉強してください。

アイスクリーム屋さんで学ぶ楽しい統計学──相関から因子分析まで── http://kogolab.jp/elearn/icecream/index.html

ハンバーガーショップで学ぶ楽しい統計学──平均から分散分析まで── http://kogolab.jp/elearn/hamburger/index.html

心理データ解析

http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/top_kaiseki.html

実際に論文に書いていくには、以下の本が参考になります。

小塩 真司 (著)　2005 研究事例で学ぶ SPSS と Amos による心理・調査データ解析　東京図書 田中 敏 (著) 　1996 実践心理データ解析―問題の発想・データ処理・論文の作成　新曜社 R については、以下の本と web をお勧めします。

データ解析環境「

_R 」―定番フリーソフトの基本操作からグラフィックス、統計解析まで I ・ O

BOOKS

 舟尾 暢男 (著), ^{高浪 洋平} (著)

http://cse.naro.affrc.go.jp/takezawa/r-tips/r2.html

↑　ほとんどこれを見れば分かります。

R の初歩的操作

・以下のページ（日本語です）を見て、R をインストールしよう。

http://www.okada.jp.org/RWiki/index.php?R%20%A4%CE%A5%A4%A5%F3%A5%B9%A5%C8%A1%BC%A5%EB

（　これは、

http://www.okada.jp.org/RWiki/index.php?RjpWiki の　＜R のインストール＞というと ころ）

・とりあえず立ち上げて、走らせてください。

> のプロンプト（入力を促す記号）がでてきます。このプロンプトの後に指示を入力していきます。

・終わるためには、quit の頭文字と()をあわせて、q()といれると終了します。

・困ったときには、help(関数名)を入れましょう。

・データの入力

　エクセルか、OpenOffice の Calc で入力します。最初の行には変数名、各列には被験者データを並べ て入れます。日本語も可能ですが、変数名は英語にしておいたほうがあとあと便利です。

　変数名も含めて、データ部分をコピーすると、clipboard にそのデータが入っています。 　そこで、

> x <- read.delim(“clipboard”)

　（左矢印にハイホン）とすると、x にデータが入ります。　ちなみに x はなんでもよく、mydata で も、好きな名前にしましょう。注意することは、R は、大文字と小文字は別物と扱います。

ここで、データの形式を確認します。

>class(x) とすると、 [1] "data.frame"

返事が返ってきます。data.frame 形式は、 数値や文字、因子などのデータをまとめたもの(object) です。

>names(x)

で、変数名が返ってきます。

>x

で、データをみることができます。

もし編集したければ、「編集」→「データエディタ」を選択しましょう。

・実際のデータを触りましょう。

＜別のデータのやりとり方法＞

http://androids.happy.nu/ja/rtips.html にある関数から

クリップボードから読み取ってデータフレームをつくる関数 read.cb <- function(header=TRUE, ...)

base::read.table(file="clipboard", header=header, ...) 変数名を含めてコピーして、

my.data <- read.cb()

で、my.data にデータが移る。

変数名がない場合は、read.cb(h=F)とする。

クリップボードにデータフレームを書き込む関数

write.cb <- function(data, sep="\t", header=TRUE, row.names=FALSE,

col.names=ifelse(header && row.names, NA, header), qmethod="double", ...) base::write.table(data, file="clipboard", sep=sep, row.names=row.names, col.names=col.names, qmethod=qmethod, ...)

これは、逆に R からクリップボードにデータフレームをコピーする関数です。 write.cb(my.data)で、クリップボードに入るので、貼り付ける。

person.data を利用して

以下を打ち込んでください。

datafilename <- "http://personality-project.org/r/datasets/maps.mixx.epi.bfi.data"

person.data <- read.table(datafilename,header=TRUE) #read the data file

#のマークの後は、メモです。R では読み飛ばされます。

datafilename という名前のオブジェクトに、web 上アドレスを読み込みました。それを使って、 person.data という名前のオブジェクトに、データを読み込んだわけです。

（このファイルの最後にデータを移しておきますので、それをコピーして、いったん表計算に貼り付 けて、clipboard にうつしてもかまいません。)

～～

もしコンピュータ上のファイルから読みたければ、以下の方法で可能です。 datafilename<-"Desktop/epi.big5.txt"

#now read the data file

person.data <-read.table(datafilename,header=TRUE) #read the data file

あるいは、CSV 形式のファイルの場合（エクセルから CSV 形式で保存が可能）は以下の方法でおこな う。

data <-read.table(datafilename,header=TRUE,sep=",") あるいは

data<-read.csv(datafilename)

～～～

（このマニュアルでは欠損値についてはとりあえず触れません。完全データとしてやっていきま す。）

～～～ 念のため

> class(person.data)

とすると、data.frame という答えが返ってきます。

>names(person.data) #list the names of the variables と問いかけると、次の変数名が返ってきます。

[1] "epiE" "epiS" "epiImp" "epilie" "epiNeur" "bfagree" [7] "bfcon" "bfext" "bfneur" "bfopen" "bdi" "traitanx" [13] "stateanx"

変数名から、内容は分かりますね？ epiE は、EPI 尺度の外向性得点。 epiS は、EPI 尺度の社交性得点。 epiImp は、EPI 尺度の衝動性得点。 epilie は、ライスケール得点。 epiNeur は、神経症傾向得点。

bfagree は、５因子性格の調和性得点。 bfcon は、５因子性格の生真面目さ得点。 bfext は、５因子性格の外向性得点。 bfneur は、５因子性格の神経症傾向得点。 bfopen は、５因子性格の開放性得点。

＜メモ＞

Extraversion (sometimes called Surgency). The broad dimension of Extraversion encompasses such more specific traits as talkative, energetic, and assertive.

Agreeableness. This dimension includes traits like sympathetic, kind, and affectionate. Conscientiousness. People high in Conscientiousness tend to be organized, thorough, and planful.

Neuroticism (sometimes reversed and called Emotional Stability). Neuroticism is characterized by traits like tense, moody, and anxious.

Openness to New Experiences (sometimes called Intellect or Culture). This dimension includes having wide interests, and being imaginative and insightful.

bdi は、ベックの抑うつ得点 traitanx は、特性不安得点。 stateanx は、状態不安得点。

平均値の算出

とりあえずは平均値を出してみましょう。平均は英語で mean ですね。

>mean(person.data$epiE) [1] 13.33333

とうつと、person.data の epiE得点の平均値が返ってきます。

通常、心理学では平均値は、小数点１桁ですので、round()という命令で丸めてみましょう。

>round(mean(person.data$epiE),1) [1] 13.3

小数点２桁にしたいならば、上を１を２にするだけです。round(変数名,桁数) さて、いちいち person.data$epiE と打つのはめんどくさいですね。

そういう時は、attach()を使います。

>attach(person.data)

これで、person.data の中の変数名を、ダイレクトに扱うことができます。

>mean(epiE)

で OK！　　（分析がすべて終わったら最後に、detach(person.data)としましょう。）

枝葉図とヒストグラムで分布を見よう

そうそう、平均値を見る前に、分布を確かめる必要があります。

>stem(epiE)

The decimal point is at the | 1 | 0

2 | 00 3 | 4 | 00 5 | 0000 6 | 00000000 7 | 00000 8 | 0000000 9 | 00000000

10 | 000000000000000 11 | 0000000000000000

12 | 000000000000000000000000 13 | 00000000000000000000000 14 | 00000000000000000000000000 15 | 000000000000000000000 16 | 0000000000000000 17 | 0000000000000000 18 | 0000000000 19 | 00000000000000 20 | 000000

21 | 0000 22 | 000

あるいは、ファンシーに figure を出すには、

>hist(epiE)

と打ってみましょう。

ところで、data.frame形式のオブジェクトは、そのまま使えることが多いです。

>mean(person.data)

とすると、すべての変数の平均値が出ますね。

>sd(person.data) では標準偏差です。

要約統計量

全部の変数の要約統計量をだすには、

>summary(person.data)

です。もちろん、( )内に変数を指定できます。 小数点２桁にしたいならば

>round(mean(person.data),2))

箱ひげ図を描いてみよう

箱ひげ図をえがきたいならば、

>boxplot(person.data)

・figure は、「メタファイル」形式でコピーをして、OpenOffice の draw にいったん貼り付けて、「切 り離し」をすると、自由に加工ができます！

・コマンドラインで、↑のキイを押すと、前に打ち込んだものが使えて重宝します。

ちょっと見にくい図ですね。

変数を選択してまとめよう

EPI 尺度だけをとりだしてみましょう。

>epi.df <- subset(person.data, select=epiE:epiNeur)

person.data の中から、epiE から epiNeur までの変数をとりだして、epi.df というオブジェクトにい れました。epi.df も data.frame形式です。.df というのは、その形式が分かるようにつけただけです ので、別になくてもかまいません。

>class(epi.df) で確かめましょう。

>boxplot(epi.df)

で EPI 尺度のみの箱ひげ図が完成。

同じように、ビッグファイブ尺度のみを取り出しましょう。

>bf.df <- subset(person.data,select=bfagree:bfopen) 箱ひげ図は

>boxplot(bf.df)

散布図を描こう

それぞれのテストの各下位尺度間についての散布図を見てみたいですね。

>pairs(epi.df) と

>pairs(bf.df)

で散布図がでました！

相関係数

それでは、ピアソンの相関係数は？

>cor(epi.df)

ちょっと見にくいですね。

>round(cor(epi.df),2) 見やすくなりました。 同じように、

>round(cor(bf.df),2)

~~

describe 関数を導入してみよう

さて、ここでネットにつなげられることができる人は、こうしましょう。 source("http://personality-project.org/r/useful.r")

を打ち込んでください。

そうすると、新たな関数が導入されます。

describe(person.data)

mean sd skew n median se epiE 13.33 4.14 -0.33 231 14 0.27 epiS 7.58 2.69 -0.57 231 8 0.18 epiImp 4.37 1.88 0.06 231 4 0.12 epilie 2.38 1.50 0.66 231 2 0.10 epiNeur 10.41 4.90 0.06 231 10 0.32 bfagree 125.00 18.14 -0.21 231 126 1.19 bfcon 113.25 21.88 -0.02 231 114 1.44 bfext 102.18 26.45 -0.41 231 104 1.74 bfneur 87.97 23.34 0.07 231 90 1.54 bfopen 123.43 20.51 -0.16 231 125 1.35 bdi 6.78 5.78 1.29 231 6 0.38 traitanx 39.01 9.52 0.67 231 38 0.63 stateanx 39.85 11.48 0.72 231 38 0.76

The describe function can be combined with the by function to provide evem more detailed tables. This example reports descriptive statistics for subjects with lie scores < 3 and those >= 3. The second element in the by command could be a categorical variable (e.g., sex).

by 関数を用いる

by 関数は、by(データ名,条件式,関数）　です。

以下は、虚偽尺度が３点より小さい人と、そうでない人を分けて、要約統計量を出す例です。 (条件式は、A>B, A>=B, A<=B, A==B, A!=B という形が使えます。最後のは not equal です。)

by(person.data,epilie<3,describe) epilie < 3: FALSE

mean sd skew n median se epiE 12.64 4.00 -0.68 90 13.0 0.42

epiS 7.61 2.81 -0.57 90 8.0 0.30 epiImp 3.97 1.67 0.15 90 4.0 0.18 epilie 3.89 1.08 1.02 90 4.0 0.11 epiNeur 9.33 5.20 0.03 90 9.0 0.55 bfagree 128.12 16.55 -0.08 90 129.0 1.74 bfcon 117.56 20.46 -0.01 90 118.0 2.16 bfext 100.88 25.24 -0.24 90 101.0 2.66 bfneur 82.22 22.80 0.27 90 81.5 2.40 bfopen 121.97 20.55 -0.14 90 121.0 2.17 bdi 5.77 4.71 1.27 90 5.0 0.50

traitanx 37.01 9.06 0.70 90 36.0 0.95 stateanx 38.41 11.36 0.59 90 36.5 1.20

--- ----

epilie < 3: TRUE

mean sd skew n median se epiE 13.77 4.17 -0.17 141 14 0.35

epiS 7.57 2.62 -0.56 141 8 0.22 epiImp 4.62 1.97 -0.08 141 5 0.17 epilie 1.41 0.73 -0.80 141 2 0.06 epiNeur 11.10 4.59 0.22 141 10 0.39 bfagree 123.00 18.88 -0.19 141 124 1.59 bfcon 110.50 22.38 0.03 141 111 1.88 bfext 103.01 27.25 -0.51 141 105 2.29 bfneur 91.64 23.01 -0.05 141 94 1.94 bfopen 124.36 20.50 -0.17 141 126 1.73 bdi 7.43 6.29 1.16 141 6 0.53 traitanx 40.28 9.62 0.65 141 39 0.81 stateanx 40.77 11.51 0.80 141 39 0.97 これでずいぶん便利になりましたね。

他にグラフィックで以下の関数がつかえるようになりました。下の２つを試してください。 multi.hist(person.data)

#draw the histograms for all variables

pairs.panels(person.data[1:5])

#draws scatterplots, histograms, and shows correlations

Tip! 図がでているところで、「記録」をチェックすると、戻って以前の図も見れるようになります。

～～

subset 関数で分割したデータセットをつくる

虚偽尺度３点以下の人を取り出して、新たなデータセットをつくってみましょう。subset 関数を使い ます。

>person.nolie.df<-subset(person.data,epilie<3) とりあえず作っただけです。あとで使いません。

>rm(person.nolie.df) で消去しましょう。

男女別データセットを作ろう

もし性別の変数があれば、同じようにして男だけのデータ、女だけのデータを取り出すことが可能で す。

たとえば、originaldata というなかの変数 sex にM と F がはいっていたとしたら次のようになりま す。

(data.frame では、文字列は自動的に因子形式として扱われます。）

>maledata.df<-subset(originaldata,sex==”M”) この場合のイコールは２つ続きです。

あるいは、女性だけを取り出すには、次のやりかたでも OK です。（注意！　,　を忘れないよう に。）

>femaledata.df<-originaldata[originaldata$sex==”F”,]

群わけしたデータセットを作ろう・・変数の変換

ところで、変数を他の変数に変えるには、次のようなやり方が可能です。 たとえば、３群にわけるとしましょう。

> x<-c(1:9) #x に１から９までの数字を入れたベクトルにする。combination の c です。

> x<-c(x,NA) #さいごに、欠損値をいれてみる。

> x

[1] 1 2 3 4 5 6 7 8 9 NA

> x.cat<-as.integer(x>3)+as.integer(x>6)+1

> x.cat

[1] 1 1 1 2 2 2 3 3 3 NA

( )内は条件式で、TRUE あるいは FALSE　を返してきます。

それを、as.integer を使うと、１　あるいは　０　という数字に変わります。 さいごに、１を足して、１、２、３　に変わりました。

epilie と epiNeur の関係は、

>cor(epilie,epiNeur) [1] -0.2508926

メモ：２つの変数から群わけするための変数をつくる方法

butterfly <- read.csv("butterfly.csv")

> attach(butterfly)

>

> group <- ifelse(day=="14" & genotype=="eng", 1, ifelse(day=="20" & 142

genotype=="eng", 2, ifelse(day=="14" & genotype=="swe", 3, ifelse (day=="20"

& genotype=="swe", 4, NA))))

>

> butterfly <- cbind(butterfly, group)

t 検定

ですが、群わけをして t 検定をして見ましょう。

>epi.cat <- as.integer(epilie>3)+1

>epi.cat <- factor(epi.cat) #因子形式のデータに変える。

>t.test(epiNeur~epi.cat, var.equal=T)

Two Sample t-test data: epiNeur by epi.cat

t = 2.5141, df = 229, p-value = 0.01262

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:

0.4339066 3.5790194 sample estimates:

mean in group 1 mean in group 2 10.810811 8.804348

このような形で、男女差をみることができますね。

ちなみに、var.equal=T は等分散を仮定しています。F とすると、Welch の検定になります。

F検定をするには、

>var.test(epiNeur~epi.cat) でみることができます。

ピアソンの相関係数

さて、EPI テストの下位尺度と、ビッグファイブの下位尺度の関係を知りたいですね。 ピアソンの相関係数を算出してみましょう。

>round(cor(epi.df,bf.df),2)

bfagree bfcon bfext bfneur bfopen epiE 0.18 -0.11 0.54 -0.09 0.14 epiS 0.20 0.05 0.58 -0.07 0.15 epiImp 0.08 -0.24 0.35 -0.09 0.07 epilie 0.17 0.23 -0.04 -0.22 -0.03 epiNeur -0.08 -0.13 -0.17 0.63 0.09

ふむふむ。こうですか！

やっぱり、EPI の神経症と、ビッグファイブの神経症は関連が強そうですね。

相関係数の差の検定

EPI の外向性(epiE)とビッグファイブの外向性(bfext)の相関係数は、0.54 で一番高いですが、 EPI の外向性とビッグファイブの調和性得点(bfagree)の相関係数の 0.18 より、統計的に差があるで しょうか？

その検定には、青木先生の HP から「同じサンプルからの相関係数の差」を使います。 http://aoki2.si.gunma-u.ac.jp/R/diff_r.html

使用法

diff.r(n, rxy, rvy, rxv)

引数

n サンプルサイズ

rxy, rvy 差を検定する二つの相関係数

rxv 先の二つの相関係数に関連するもう一つの相関係数

ソース

diff.r <- function(n, rxy, rvy, rxv) {

detR <- (1-rxy^2-rvy^2-rxv^2)+2*rxy*rxv*rvy

t0 <- (abs(rxy-rvy)*sqrt((n-1)*(1+rxv)))/sqrt(2*detR*(n-1)/(n- 3)+(rxy+rvy)^2*(1-rxv)^3/4)

df <- n-3

p <- pt(abs(t0), df, lower=FALSE)*2 c("t value"=t0, "df"=df, "P value"=p) }

ソースをコピーして、R に貼り付けると、diff.r という関数ができあがります。 人数である、データ数は　nrow 関数を使って確認します。

> nrow(person.data) [1] 231

次に、epiE と bfext の相関係数を確認

> cor(epiE,bfext) [1] 0.5434974

そして、epiE と bfagree の相関係数を確認

> cor(epiE,bfagree) [1] 0.1759432

それから、相関係数の差の検定をするために必要な、bfext と bfagree の相関を算出。

> cor(bfext,bfagree) [1] 0.4784819

>

さっきの関数、diff.r の関数に数値を入れます。 有意に差がある結果が確認できました。

> diff.r(291, 0.54, 0.18, 0.48)

t value df P value 7.099699e+00 2.880000e+02 9.787238e-12

散布図

散布図を見てみましょうか。

>plot(epiNeur,bfneur)

相関係数の検定

念のために、検定をして見ましょう。

>cor.test(epiNeur,bfneur)

しっかり有意な結果でますね。無相関ではなさそうですね。（注意：p が低くても、相関が強いという わけではありません。単に無相関でないということしか言えません！)

相関の検定は、変数は２つしかとれないので、「ファイル」「新しいスクリプト」を開いて、そこで 使うスクリプトを書いて実行するといいでしょう。　

いったん休憩・・・

>objects()

とすると、いままで作ったオブジェクトの一覧がでます。 それらのオブジェクトをすべてを消去するには、

>rm(list=ls())

というおまじないをしましょう。 一つだけなら

>rm(オブジェクトの名前) rm は remove の略ですね。

クラスター分析

ビッグファイブの５因子性格の得点で、クラスター分析（ウォード法）やってみましょう。

> hc <- hclust(dist(bf.df), "ward")

> plot(hc)

hclust は、クラスター分析の関数。 dist は、類似性の計算のための関数です。

大雑把に見て、４クラスターに分かれそうですね。memberというオブジェクトに、cutree関数を使って結果を入れます。

> member<-cutree(hc, k=4)

データの形式を確かめてみると、> class(member)[1] "integer"

ですが、これは番号の順番に特に意味はないので、名義尺度と考えたほうが良いでしょう。そこで、factor形式に変えます。（重要！）

> member<-factor(member) [1] 1 2 3 3 1 3 2 2 1 3 2 2 1 1 1 2 2 1 1 2 3 2 1 2 4 3 1 3 3 2 2 3 3 2 3 3 1[38] 3 3 2 2 1 1 2 2 2 2 2 2 1 3 3 2 3 1 1 3 2 2 3 1 1 3 3 1 1 3 1 4 2 4 2 3 3

16

186195 220193 229₁₆₄₁₇₆ 169206 203₂₂₅ 184221 167191 214₂₀₂₁₇₂ 223₂₅ 192₇₁ 137₁₁₁₇₆ 166168 179170 197₁₉₈₁₆₂ 217₂₂₆ 189228₁₉₄ 163224₁₇₃₉₈ 178185 161222 216₉₆ 199211₁₇₇₁₈₇ 196180 208₁₈₂₆₉ 174₁₂₉₅₆ 1241₈₅ 200175 201₇₇₆₅ 154₁₇₁₁₉₀18 212₂₇ 110₄₂₅₅ 135₄₃95 109₆₂₁₃ 119₁₅₁₈₆ 1171569₁₂₅ 152157₈₃ 116118₁₃₈₅₀

10568₈₄ 147₁₄₆₁ 143₁₅ 181₁₃₄₆₆ 1395₂₁₈₃₇87 2319 115₄₄₉₄ 104₃₀₄₁ 106₁₃₂₁₁ 127₁₄₂₂₃₀72 5358 145₅₉₈₈ 10222₂₀₄₀₇

318 159₁₆ 149₁₄₀₇₅ 113₁₄₁82₄₅ 12648 148153₁₂ 128₂₄₄₇

7846 122204 213₁₁₄₂₀₅ 215219₁₀₃₂₀₉

8979 227₁₆₅₁₇ 146₁₈₈₂₃₁

3497 112₇₀ 160₂₀₇₄₉ 183₂ 210₃₅₉₁

6774 1006₁₂₁₉₀

2839 123144₃₃₂₆ 150₁₃₃₇₃ 10821₆₀₃₈ 101₁₃₁₃₅₁

8093 107155 1201364₃₆₉₂ 52₉₉₈₁ 5764₃₂ 13054₁₀₂₉ 15863

0 500 1000 1500

ClusterDendrogram

hclust (*, "ward") dist(bf.df)

Height

[75] 2 4 1 2 2 3 3 2 1 1 1 1 1 2 2 3 3 3 3 2 1 4 2 4 3 3 3 2 2 2 1 2 3 3 1 1 4 [112] 2 2 2 1 1 1 1 1 3 3 2 3 1 1 2 2 2 1 3 3 2 3 1 1 3 4 1 1 2 2 2 1 3 2 2 1 2 [149] 2 3 1 1 2 1 3 1 1 3 2 2 4 4 4 4 2 4 4 4 4 4 1 4 4 4 1 4 4 4 4 4 1 4 2 4 4 [186] 4 4 2 4 1 4 4 4 4 4 4 4 4 4 1 1 4 4 2 2 4 2 4 2 2 4 1 2 4 2 4 4 1 2 4 4 4 [223] 4 4 4 4 2 4 4 2 2

Levels: 1 2 3 4

bf.df のフレームデータに、この member の変数を加えましょう。

> bfm.df<-cbind(bf.df,member)

で、同じ行数のデータを横につなげることができます。colum の c に bind ですね。

> bfm.df

で確認しましょう。

平均を示す図を描こう

ここで、「パッケージ」から「パッケージのインストール」をクリック。(インターネットにつないで いる状態でしてください。）

Japan(Aizu)あるいは(Tuskuba)を選択したうえで、gplots　を選択して、新しいパッケージをダウン ロードしましょう。

> library(gplots)

とすると、取り込みが完了します。

クラスターごとの平均の差を図で見てみましょう。

> plotmeans(bfopen~member, data=bfm.df)

> plotmeans(bfext~member, data=bfm.df)

> plotmeans(bfneur~member, data=bfm.df)

> plotmeans(bfagree~member, data=bfm.df)

> plotmeans(bfcon~member, data=bfm.df)

>

もし箱ひげ図でみたければ、以下のようにします。

> boxplot(bfopen~member,data=bfm.df)

一要因の分散分析

一要因の分散分析をしてみましょう。（注意：member を因子形式に変えていないと正しくいきませ ん！）

> summary(aov(bfopen~member,data=bfm.df))

Df Sum Sq Mean Sq F value Pr(>F) member 3 36090 12030 45.039 < 2.2e-16 *** Residuals 227 60631 267 ---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

（念のために、1e-2 は、1 x 10 のマイナス２乗ということで、0.01 1e2 は、1 x 10 の２乗ということで、１００ ということです。）

～～

別なやり方はこれ。var=T をしていしないか、あるいは var=F とすると、等分散でない仮定の元での 検定になります。

> oneway.test(bfopen~member,data=bfm.df,var=T) One-way analysis of means

data: bfopen and member

F = 45.0393, num df = 3, denom df = 227, p-value < 2.2e-16

～～～

by 関数を利用しながら、要約統計量を群ごとに出しておきましょう。

> by(bfm.df$bfopen, bfm.df$member, summary) INDICES: 1

Min. 1st Qu. Median Mean 3rd Qu. Max. 86.0 111.3 120.5 122.1 133.8 161.0

--- INDICES: 2

Min. 1st Qu. Median Mean 3rd Qu. Max. 96.0 114.3 126.0 124.3 133.0 156.0

--- INDICES: 3

Min. 1st Qu. Median Mean 3rd Qu. Max. 73.00 89.25 103.00 102.80 116.50 149.00

--- INDICES: 4

Min. 1st Qu. Median Mean 3rd Qu. Max. 88. 132.0 141.0 140.3 152.0 173.0

89.

標準偏差も出しておきましょう。

> by(bfm.df$bfopen, bfm.df$member, sd) INDICES: 1

[1] 16.01693

--- INDICES: 2

[1] 13.74355

--- INDICES: 3

[1] 18.02367

--- INDICES: 4

[1] 18.10518

>

それぞれ出すのが面倒ならば、データフレーム bfm.df だけで行いましょう。

> by(bfm.df, bfm.df$member, summary)

> by(bfm.df, bfm.df$member, sd) ちなみに

> by(bf.df, bfm.df$member, mean)

とやると、クラスターごとの平均値は出てきます。

～～

useful.r をいれてあると、

>by(bf.df, bfm.df$member, describe) が可能です。

～～

多重比較！

　群の間に有意な差がありましたので、どこに差があるかを検定したいですね。 一度、分散分析の結果をオブジェクトにしましょう。

>aov.open <- aov(bfopen~member,data=bfm.df) 分散分析の結果の要約は

>summary(aov.open)

http://phi.ypu.jp/statlib/stat.pdf　にある

中澤　港：「R による統計解析の基礎」ピアソン・エデュケーション，2003 年 10 月 15 日初版第１刷， A5判，183pp.，ISBN4-89471-757-3，1,800 円 の草稿の 111 ページから１１２ページを参照すると、

「現在では，かなり広い用途をもち，ノンパラメトリックな分析にも適応可能なホルム(Holm) の方法

（ボンフェローニの方法を改良して開発された方法）が第一に考慮されるべきである。」 とありますので、それでおこないましょう。

ホルムの多重比較

>pairwise.t.test(bf.df$bfopen, member, data=bfm.df)

Pairwise comparisons using t tests with pooled SD data: bf.df$bfopen and member

1 2 3 2 0.45 - - 3 3.0e-08 1.9e-10 - 4 3.0e-08 2.1e-07 < 2e-16

P value adjustment method: holm １群と３群、４群

２群と３群、４群

３群と４群　の間に差がありましたが、

クラスター１と２の間には開放性に差がありませんでしたね。

～～～～～～～

メモ：古典的な TukeyHSD による多重比較は、

> TukeyHSD(aov.open) 結果です。

Tukey multiple comparisons of means 95% family-wise confidence level

Fit: aov(formula = bfopen ~ member, data = bfm.df)

$member

diff lwr upr p adj 2-1 2.176355 -5.334131 9.68684 0.8766023 3-1 -19.377061 -27.728279 -11.02584 0.0000000 4-1 18.195402 10.306373 26.08443 0.0000001 3-2 -21.553416 -29.581781 -13.52505 0.0000000 4-2 16.019048 8.472619 23.56548 0.0000006 4-3 37.572464 29.188907 45.95602 0.0000000

同じく、クラスター１と２の間には開放性に差がありませんでしたね。

～～～～～～～～

さて、同じように、他の４つの性格得点ではどうなるかやってみましょう。＞宿題！ ４つのクラスターはそれぞれ、どんな人たちの集まりといえるでしょうか？

～～

～～ちょっと高度に図を描いてみましょう～～

> cl1.mean<-mean(bf.df[member==1,])

> cl2.mean<-mean(bf.df[member==2,])

> cl3.mean<-mean(bf.df[member==3,])

> cl4.mean<-mean(bf.df[member==4,])

> cl.mean<-cbind(cl1.mean,cl2.mean,cl3.mean,cl4.mean)

> matplot(cl.mean,type="l")

> matplot(cl.mean, main="クラスター別の性格因子得点",pch="1234", type="b")

～～～～

分散分析は、以下の式でもいいようですね。 anova(lm(bfopen~member, data=bfm.df) 同じ結果になります。

～～

クラスターに分けたわけですが、これらの群に、特性不安得点の差はあるのでしょうか？

> bfa.df<-cbind(bfm.df, traitanx)

> bfa.df

bfa.df という新しいデータフレームをつくりました。 平均値の図は、次で出ます。

>plotmeans(traitanx~member, connect=F, data=bfa.df) connect=F を追加したので、折れ線の図にはなっていません。

>avo.bfa<-aov(traitanx~member, data=bfa.df)

>summary(avo.bfa)

>pairwise.t.test(bfa.df$traitanx, member, data=bfa.df)

1

1

1

1

1

1 2 3 4 5

80100120140

クラスター別の性格因子得点

cl.mean

2

2

2 2

2

3

3

3

3

3 4

4

4

4

4

一要因の分散分析は、これでいいですね。

－ーーーーー２元配置は別な機会にします。

・回帰分析

　次は、linear model を使って回帰分析をしてみましょう。

　回帰分析と重回帰分析の説明は、http://www1.doshisha.ac.jp/~mjin/R/

の「

R と回帰分析 」http://www1.doshisha.ac.jp/~mjin/R/13.html と「

　ベックの抑うつ尺度得点(bdi)を予測していきます。最初は、EPI の神経症尺度得点(epiNeur)から、 次に、特性不安得点(traitanx)を追加してみます。#のあとはコメントで、計算のときは無視されま す。

Model1 <- lm(bdi~epiNeur) #simple regression of beck depression on Neuroticism summary(model1) # basic statistical summary

ここで、グラフの出力設定を変えておきます。

#pass parameters to the graphics device op <- par(mfrow = c(2, 2), # 2 x 2 pictures on one plot pty = "s") # square plotting region, # independent of device size

plot(model1) #diagnostic plots in the graphics window こんな結果になりますね。

Call:

lm(formula = bdi ~ epiNeur) Residuals:

Min 1Q Median 3Q Max -11.9548 -3.1577 -0.7707 2.0452 16.4092 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) -0.32129 0.73070 -0.44 0.66 epiNeur 0.68200 0.06353 10.74 <2e-16 *** ---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.721 on 229 degrees of freedom

Multiple R-Squared: 0.3348, Adjusted R-squared: 0.3319 F-statistic: 115.3 on 1 and 229 DF, p-value: < 2.2e-16

以上で、抑うつへ神経症傾向の回帰です。

Model2 <- lm(bdi~epiNeur+traitanx) #add in trait anxiety summary(model2) #basic output

plot(model2)

特性不安を加えた結果です。 どちらが良いでしょうか？

anova(model1,model2) #compare the difference between the two models Analysis of Variance Table

Model 1: bdi ~ epiNeur

Model 2: bdi ~ epiNeur + traitanx

Res.Df RSS Df Sum of Sq F Pr(>F) 1 229 5103.3 2 228 4214.3 1 889.1 48.101 4.133e-11 *** ---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 え～っと、どうみるのかな（謎？？？

## At end of plotting, reset to previous settings: par(op)

標準偏回帰係数について

標準偏回帰係数 β は、別途に計算しましょう。（下を参考に）

あるいは、使う変数の点数をすべて標準化してから重回帰をおこないましょう。 scale 関数で、z得点に変換できます。（mean=0,SD=1)

>zepiNeur <- scale(epiNeur)

>ztraitanx <- scale(traitanx)

> zbdi <- scale(bdi)

> model2=lm(zbdi~zepiNeur+ztraitanx)

> summary(model2)

Estimate Std. Error t value Pr(>|t|) (Intercept) 9.507e-17 4.898e-02 1.94e-15 1.00000 zepiNeur 2.164e-01 7.167e-02 3.019 0.00282 ** ztraitanx 4.971e-01 7.167e-02 6.935 4.13e-11 *** ---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7444 on 228 degrees of freedom Multiple R-Squared: 0.4507, Adjusted R-squared: 0.4459 F-statistic: 93.53 on 2 and 228 DF, p-value: < 2.2e-16

～～～

参考：

http://phi.ypu.jp/swtips/R.html#MISC 標準化偏回帰係数を得る

●ベクトルとして計算させると楽である。例えば，従属変数を y，独立変数を x1, x2, x3 とすれ ば，

res <- lm(y~x1+x2+x3)で線型回帰を行った後，

sdd <- c(0,sd(x1),sd(x2),sd(x3))として各独立変数の不偏標準偏差ベクトルを作り（0 は切片 用），

stb <- coef(res)*sdd/sd(y)として，偏回帰係数ベクトルに不偏標準偏差ベクトルを掛けて，従 属変数の不偏標準偏差で割ってやれば，stb に標準化偏回帰係数のベクトルが得られる。

～～～～～～ 参考２：

多重共線性のチェック（トレランスと VIF）

これは重回帰分析をするときにチェックしたほうがよいポイントです。 http://aoki2.si.gunma-u.ac.jp/R/ の以下のページを見てください。

# 多重共線性のチェック（トレランス）

# 多重共線性のチェック（従属性）

http://aoki2.si.gunma-u.ac.jp/lecture/Regression/mreg/mreg6.html に説明があり、 http://aoki2.si.gunma-u.ac.jp/R/tolerance.html に、使い方の説明があります。

VIF (Variance Inflation Factor 変動インフレーション因子)が 10 以上なら多重共線性が起こってお り、2 以下なら起こっていないと判断するようです。

> library(DAAG) # DAAG ライブラリを使う

> result <- lm(bdi~bfext+bfneur+bfagree+bfcon+bfopen)

> vif(result) # これが本来ほしかった結果 VIF bfext bfneur bfagree bfcon bfopen 1.4651 1.1356 1.5914 1.2858 1.5284

～～～～～～

では、ビッグファイブの５性格特性から、ベックの抑うつ得点への予測を検討しましょう。

> bfd.df<-cbind(bf.df,person.data$bdi)

ここで、http://aoki2.si.gunma-u.ac.jp/R/にある 総当たり法による重回帰分析

http://aoki2.si.gunma-u.ac.jp/R/All_possible_subset_selection.html の関数をとりいれてやってみましょう。

以下のソースを全部、コピーしてコマンドラインへ。

（詳しい使い方は、web ページを参照） ソース

All.possible.subset.selection <- function(df, adjusted=TRUE, limit=10) {

df <- df[complete.cases(df),] if ((nv <- ncol(df)-1) > limit) {

stop("Too many independent variables.") }

name <- names(df) n <- 2^nv-1

str <- character(nv)

result1 <- result2 <- numeric(n) result3 <- character(n)

for (i in 1:n) { k <- i m <- 0

for (j in 1:nv) { if (k%%2) { m <- m+1

str[m] <- name[j] }

k <- k %/% 2 }

form <- reformulate(str[1:m], name[nv+1])

result1[i] <- (result <- summary(lm(form, df)))$r.square result2[i] <- result$adj.r.square

result3[i] <- paste(c(name[nv+1], "~", paste(str[1:m], collapse=" + ")), collapse="

") }

o <- rev(order(if (adjusted) result2 else result1))

cbind("R square"=result1[o], "Adjusted R square"=result2[o], "formula"=result3[o]) }

その上で、

> result <- All.possible.subset.selection(bfd.df)

> print(result, quote=FALSE)# 余分な引用符が出力されないように quote=FALSE を指定する R square Adjusted R square

[1,] 0.286239056249838 0.276806092235519 [2,] 0.287793730787998 0.27518831009398 [3,] 0.286878849528703 0.274257236246025 [4,] 0.289413570946956 0.273622761412443 [5,] 0.270598367396201 0.260958698242846 [6,] 0.267105320574921 0.26067641987821 [7,] 0.268866117911653 0.259203555593305 [8,] 0.267907426802641 0.258232194557742 [9,] 0.270663118512867 0.257754501141414 [10,] 0.269019195330606 0.256081481973625 [11,] 0.257291445678644 0.250776458360035 [12,] 0.25844168789243 0.248641357776471 [13,] 0.241903313049384 0.235253342111221

[14,] 0.244579038432932 0.234595501495922 [15,] 0.231897369837557 0.225159627467711 [16,] 0.217311007194546 0.213893151330767 [17,] 0.0410872064910799 0.0326756907585455 [18,] 0.0420771396288114 0.0294173661437295 [19,] 0.0414000494102277 0.0287313275962658 [20,] 0.0366963270992104 0.0282462948807823 [21,] 0.032043619774868 0.0278167360184265 [22,] 0.0426033641525542 0.0256582909517143 [23,] 0.0325715476597082 0.024085333165495 [24,] 0.0367025405798106 0.0239717371513498 [25,] 0.026679944127585 0.0181420489006340 [26,] 0.0200088896687153 0.0157294525056967 [27,] 0.0194330682429096 0.0151511165758481 [28,] 0.0267054440343939 0.0138425203872713 [29,] 0.0205189590039890 0.0119270200478836 [30,] 0.0196328127715813 0.0110331006029109 [31,] 0.00585385855166141 0.00151260902568628

formula [1,] person.data$bdi ~ bfcon + bfneur + bfopen [2,] person.data$bdi ~ bfcon + bfext + bfneur + bfopen [3,] person.data$bdi ~ bfagree + bfcon + bfneur + bfopen [4,] person.data$bdi ~ bfagree + bfcon + bfext + bfneur + bfopen [5,] person.data$bdi ~ bfext + bfneur + bfopen [6,] person.data$bdi ~ bfneur + bfopen [7,] person.data$bdi ~ bfcon + bfext + bfneur [8,] person.data$bdi ~ bfagree + bfneur + bfopen [9,] person.data$bdi ~ bfagree + bfext + bfneur + bfopen [10,] person.data$bdi ~ bfagree + bfcon + bfext + bfneur [11,] person.data$bdi ~ bfcon + bfneur [12,] person.data$bdi ~ bfagree + bfcon + bfneur [13,] person.data$bdi ~ bfext + bfneur [14,] person.data$bdi ~ bfagree + bfext + bfneur [15,] person.data$bdi ~ bfagree + bfneur [16,] person.data$bdi ~ bfneur [17,] person.data$bdi ~ bfcon + bfext [18,] person.data$bdi ~ bfagree + bfcon + bfext [19,] person.data$bdi ~ bfcon + bfext + bfopen [20,] person.data$bdi ~ bfagree + bfcon [21,] person.data$bdi ~ bfcon [22,] person.data$bdi ~ bfagree + bfcon + bfext + bfopen [23,] person.data$bdi ~ bfcon + bfopen [24,] person.data$bdi ~ bfagree + bfcon + bfopen [25,] person.data$bdi ~ bfagree + bfext [26,] person.data$bdi ~ bfagree [27,] person.data$bdi ~ bfext [28,] person.data$bdi ~ bfagree + bfext + bfopen [29,] person.data$bdi ~ bfagree + bfopen [30,] person.data$bdi ~ bfext + bfopen [31,] person.data$bdi ~ bfopen

>

上のようなことも力技で可能ですが、意味がある変数を取り入れましょう！

～～～～ 終わり。

～～参考：統計処理の羅針盤～

http://phi.ypu.jp/swtips/compass.html

第２部：データの入力から、尺度の構成までのミ

ニミニヒント

・質問紙をつくる際には、０点をつけるようにするのは極力避けましょう。あとで欠損値との区別が

紛らわしくなるからです。

・回答忘れがないように調査の時には、できるだけしっかりとお願いしましょう。

・既存の尺度を使った場合も、因子分析と信頼性係数は確かめてから、どのように得点化するかを

考えましょう。

・分析方法まで考えないと質問紙をつくったり、実験を行ったりはできません。しっかり分析方法の

見通しまで立てましょう。

集めるデータの数について

　もし因子分析を考えているのならば、項目数の１０倍あるのが望ましいです。無理ならば、５～６

倍の人数は欲しいところです。

調査項目の作成については、以下を参照のこと。（必ず見てくださいね）

http://www.ec.kagawa-u.ac.jp/~hori/chosa05/chosa05index.html#item

http://www.okayama-

u.ac.jp/user/le/psycho/member/hase/education/2000/_0seminar/index.html#quest

http://www.human.tsukuba.ac.jp/~miyamoto/Questionnaire.htm

http://staff.miyakyo-u.ac.jp/~m-taira/Lecture/questionnaire2001.html

http://www.nurs.or.jp/~suzukirt/psychology/scale_develop.html

・質問紙を回収したら～

・・まずは、ナンバリングのスタンプを押して、回収した質問紙にすべて通し番号をつけていきま

しょう。それから、データの入力です。

通し番号(ID)から順番にデータを入力していきます。エクセルや OpenOffice の Calc に入力しま

す。

その際に、横の列には変数を、縦の行には各調査協力がはいっていきます。

ID GENDER v1 v2 v3

1 M 1 2 2 一人目のデータ

2 F 2 2 4 二人目のデータ

3 F 2 4 5 三人目のデータ

という具合にです。

欠損値は、空白にするか、-999 のような明白な数字をいれます。（とりあえず空白でいいでしょう）

欠損値があまりにも多いデータは利用をあきらめます。たとえば３分の１以上記入がない場合。

また、回答があっても、すべて１とか５とかのように、明らかに回答がいい加減だと判断される場合

は、ケースを削除します。

欠損値の扱い方にはいろいろあります。

ある分析をするときに、使われる変数の中に一つでも欠損値があったならば、その人のデータは

利用しないというリスト消去。欠損値が入るものにペアになるものを使わずに計算する方法。欠損

値を推定して、その数値（平均値や中央値や最頻値）を代入していく方法などです。

とりあえずのお勧めは、最頻値を代入です。あとあとの計算がしやすくなります。

R では、最頻値をもとめる関数はありませんので、代わりに stem()関数を使って調べてください。

あるいは、エクセル(Calc)上で、=mode(範囲)　で求めてください。

あとは、表計算上で、欠損値のセルに、求めた最頻値をそれぞれ代入していきましょう。

生のデータを見なおすためにも、手作業でやっていきましょう。

オリジナルのデータは保存しておき、欠損値を代入したデータは新しい名前で保存しておきましょ

う。

～～～

～～補足：エクセルのアドインを使っての因子分析

一番簡単なのは、以下のアドインをエクセルで使うことです。欠損値を処理できますし、因子分

析、アルファ係数も出せます。

Excel アドイン http://homepage3.nifty.com/hideakim/

因子分析について詳しいのは、以下のページを参考にしてください。探索的因子分析リンク集（日本語中心） http://www.ec.kagawa-u.ac.jp/~hori/spss/factorlink.html

因子分析についての一般的注意は以下のページを参考。

http://www4.ocn.ne.jp/~murakou/factor.htm

ほんとうの初心者は以下の本がお勧め。

誰も教えてくれなかった因子分析―数式が絶対に出てこない因子分析入門

尺度の構成へ

せっかくですから R を使って因子分析

～～大事な話～～

　因子分析は、あくまで利用した質問項目の中での構造を探索的に試行錯誤して調べようとするものです。と いうことは、利用する質問項目に制限された枠組みのなかでの因子探索です。いいかえれば、最初にどのよう な質問項目を設定するかが決定的だということです。質問項目の設定は十分の吟味をする必要がありますし、 それが心理学研究では最重要なことかもしれません。

　また後でも述べますが、因子分析という統計的手法をしたから因子が決定されるということはありません。あく まで研究者がおこなう解釈によって因子が決定されます。その決定された因子が、理論的な文脈において、十 分に他の人たちに説得性があるかが重要となります。

～～～

分析をしたいデータセットを x にいれます。

小塩のページにあるデータを、エクセルか、calc にコピーをしましょう。

http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/08_folder/da08_02.html

番号 外向性 社交性 積極性 知性 信頼性 素直さ

1 3 4 4 5 4 4

2 6 6 7 8 7 7

3 6 5 7 5 5 6

4 6 7 5 4 6 5

5 5 7 6 5 5 5

6 4 5 5 5 6 6

7 6 6 7 6 4 4

8 5 5 4 5 5 6

9 6 6 6 7 7 6

10 6 5 6 6 5 5

11 5 4 4 5 5 5

12 5 5 6 5 4 5

13 6 6 5 5 6 5

14 5 5 4 4 5 3

15 5 6 4 5 6 6

16 6 6 6 4 4 5

17 4 4 3 6 5 6

18 6 6 7 4 5 5

19 5 3 4 3 5 4

20 4 6 6 3 5 4

OpenOffice でダブルクリックすると表計算になります。

それから、番号は余分ですので、それ以外のデータをコピーしてから

>x <- read.delim(“clipboard”)

>x

外向性 社交性 積極性 知性 信頼性 素直さ

1 3 4 4 5 4 4

2 6 6 7 8 7 7

3 6 5 7 5 5 6

4 6 7 5 4 6 5

5 5 7 6 5 5 5

6 4 5 5 5 6 6

7 6 6 7 6 4 4

8 5 5 4 5 5 6

9 6 6 6 7 7 6

10 6 5 6 6 5 5

11 5 4 4 5 5 5

12 5 5 6 5 4 5

13 6 6 5 5 6 5

14 5 5 4 4 5 3

15 5 6 4 5 6 6

16 6 6 6 4 4 5

17 4 4 3 6 5 6

18 6 6 7 4 5 5

19 5 3 4 3 5 4

20 4 6 6 3 5 4

その上で、x の相関係数を求めて、それから固有値を計算して、y にいれます。

> y <-eigen(cor(x))

> names(y)

[1] "values" "vectors"

> y$values

出力例：

[1] 2.6911757 1.5214560 0.7145326 0.4821828 0.3340264 0.2566265

これで固有値がでます。

図を描いてみましょう。

> plot(y$values, type="b")

因子数を 決めるには、固有値が１以上の基準（カイザーガットマン基準)や、固有値が落ちるところ

で 決める基準(スクリープロット基準)などを参考にしていきます。

！　もっとも大事なことは、「解釈可能性」です！　

この場合、２としましょう。

以下のソースをコピーして、R にはりつけると、pfa 関数がつくられ、主因子法による因子分析が可

能になります。出典は、青木先生の url です。

http://aoki2.si.gunma-u.ac.jp/R/pfa.html

ソース

pfa <- function(dat, method=c("Varimax", "Biquartimax", "Quartimax", "Equimax",

"None"), eps1=1e-5, eps2=1e-5, max1=999, max2=999, factors=0) {

method <- match.arg(method) dat <- dat[complete.cases(dat),] nr <- nrow(dat)

nc <- ncol(dat)

r0 <- r <- if (nr == nc) dat else cor(dat) communality0 <- sqrt(1-1/diag(solve(r))) diag(r) <- communality0

result <- eigen(r) eval <- result$values evec <- result$vectors if (factors == 0) {

factors <- sum(eval >= 1) }

converged <- FALSE for (i in 1:max1) {

eval <- eval[1:factors] evec <- evec[,1:factors] evec <- t(sqrt(eval)*t(evec)) r <- r0

communality <- apply(evec^2, 1, sum)

if (all(abs(communality-communality0) < eps1)) { converged <- TRUE

break }

communality0 <- communality diag(r) <- communality result <- eigen(r) eval <- result$values evec <- result$vectors }

if (converged == FALSE) { warning("Not converged.") }

else if (any(communality >= 1)) { warning("Communality >= 1.")

}

if (factors == 1 || method == "None") { fl <- evec*sqrt(communality)

w <- solve(r0)%*%evec

scores <- (scale(dat)*sqrt(nr/(nr-1)))%*%w

invisible(list(rotation=method, correlation.matrix=r0, communality=communality,

before.rotation.fl=evec, before.rotation.eval=eval,

after.rotation.fl=NA, after.rotation.eval=NA, scores=scores)) }

else {

fl <- evec/sqrt(communality) eig <- rep(0, factors)

ov <- 0

wt <- switch (method, "Varimax" = 1,

"biQuartimax" = 0.5, "Quartimax" = 0, "Equimax" = nc/2) fnp <- nc

for (loop in 1:max2) {

for (k in 1:(factors-1)) { for (k1 in (k+1):factors) { x <- fl[,k]

y <- fl[,k1] xy <- x^2-y^2 a <- sum(xy) b <- 2*sum(x*y)

c <- sum(xy^2-4*x^2*y^2) d <- 4*sum(x*y*xy)

dd <- d-2*wt*a*b/fnp

theta <- atan(dd/(c-wt*(a^2-b^2)/fnp)) if(sin(theta)*dd < 0) {

if (theta > 0) { theta <- theta-pi }

else {

theta <- theta+pi }

}

theta <- theta/4 cs <- cos(theta) si <- sin(theta) fljk <- fl[,k] fljk1 <- fl[,k1]

fl[,k] <- fljk*cs+fljk1*si

fl[,k1] <- fljk1*cs-fljk*si }

}

v <- sum((t(fl)^2-apply(fl^2, 2, sum)*wt/fnp)^2)

if (abs(v-ov) < eps2) { break

}

ov <- v }

fl <- fl*sqrt(communality) w <- solve(r0)%*%fl

scores <- (scale(dat)*sqrt(nr/(nr-1)))%*%w

invisible(list(rotation=method, correlation.matrix=r0, communality=communality,

before.rotation.fl=evec, before.rotation.eval=eval,

after.rotation.fl=fl, after.rotation.eval=apply(fl^2, 2, sum), scores=scores))

} }

# 相関係数行列の整形出力のための関数 cor.matrix.out <- function(result) {

printf <- function(fmt, ...) cat(sprintf(fmt, ...)) r <- result$correlation.matrix

nv <- nrow(r) printf(" ") for (i in 1:nv) {

printf(" Var%03i", i) }

printf("\n") for (i in 1:nv) {

printf(" Var%03i", i) for (j in 1:nv) {

printf("%7.3f", r[i, j]) }

printf("\n") }

}

# 因子負荷量などの整形出力のための関数

fl.matrix.out <- function(result, before=FALSE) {

printf <- function(fmt, ...) cat(sprintf(fmt, ...)) out2 <- function(str, fmt, vec, n) {

printf("%7s", str) for (j in 1:n) { printf(fmt, vec[j]) }

printf("\n") }

communality <- result$communality

if (before || is.na(result$after.rotation.fl)) { fl <- result$before.rotation.fl

eval <- result$before.rotation.eval label="E-value"

if (result$rotation == "None") { printf("Result without rotation\n") }

else {

printf("Before %s rotation\n", result$rotation) }

} else {

fl <- result$after.rotation.fl eval <- result$after.rotation.eval label="Sq.Sum"

printf("After %s rotation\n", result$rotation) }

nv <- nrow(fl) nc <- ncol(fl) printf(" ") for (i in 1:nc) {

printf(" Fac%03i", i) }

printf(" Communality\n") for (i in 1:nv) {

printf(" Var%03i", i) for (j in 1:nc) {

printf("%7.3f", fl[i, j]) }

printf("%7.3f\n", communality[i]) }

out2(label, "%7.3f", eval, nc)

out2("Cont.", "%7.1f", eval/nv*100, nc)

out2("Cum.", "%7.1f", cumsum(eval/nv*100), nc) }

さっきのデータをつかってみましょう。

> result<-pfa(x)

> result

$rotation

[1] "Varimax"

$correlation.matrix　　←相関係数

外向性 社交性 積極性 知性 信頼性 素直さ

外向性 1.0000000 0.4865985 0.59742680 0.2423646 0.27631579 0.2188622

社交性 0.4865985 1.0000000 0.55796302 0.1250652 0.31685482 0.1725433

積極性 0.5974268 0.5579630 1.00000000 0.2407219 0.03733918 0.1466444

知性 0.2423646 0.1250652 0.24072187 1.0000000 0.43625629 0.6271032

信頼性 0.2763158 0.3168548 0.03733918 0.4362563 1.00000000 0.5836325

素直さ 0.2188622 0.1725433 0.14664439 0.6271032 0.58363246 1.0000000

$communality　←　共通性

[1] 0.5362094 0.4561217 0.6941481 0.4579444 0.4414802 0.8202977

$before.rotation.fl　←　回転前の因子負荷量

[,1] [,2]

[1,] -0.6321337 -0.3696165

[2,] -0.5654597 -0.3692926

[3,] -0.6095860 -0.5679375

[4,] -0.5883707 0.3343117

[5,] -0.5730835 0.3362371

[6,] -0.7098819 0.5624636

$before.rotation.eval

[1] 2.269470 1.136731

$after.rotation.fl　←　回転後の因子負荷量

[,1] [,2]

[1,] -0.20372323 -0.70335358

[2,] -0.15560769 -0.65719704

[3,] -0.05078255 -0.83160644

[4,] -0.65682534 -0.16286445

[5,] -0.64706833 -0.15093947

[6,] -0.90206049 -0.08114526

$after.rotation.eval

[1] 1.732126 1.674076

$scores　　　←　因子得点

[,1] [,2]

1 1.020792566 1.5183719

2 -2.023495927 -1.0206819

3 -0.501863422 -0.8340721

4 -0.019117969 -0.5395410

5 0.156678493 -0.5965202

6 -0.798849597 0.6438710

7 1.044371316 -1.2134759

8 -0.787240394 0.8812035

9 -1.225857589 -0.7051758

10 0.025510276 -0.4888099

11 0.006457166 0.9469725

12 0.362812188 -0.1080777

13 -0.131817938 -0.3218104

14 1.661687552 0.4550749

15 -0.970941738 0.5956135

16 0.434611557 -0.6386810

17 -0.967798954 1.8792744

18 0.382101755 -1.1888862

19 1.070419366 0.9829852

20 1.261541293 -0.2476348

>

相関係数、共通性、回転前と回転後の因子負荷量（マイナスになっていますが、全部ひっくりかえ

した意味として取って構いません）、また因子得点がでてきますね。（メモ：共通性＋独自性は、１

です。）

小塩の SPSS の結果と比べてください。

さて、プロマックス回転をしたいならば、

http://aoki2.si.gunma-u.ac.jp/R/factanal2.html

これは、最尤法（さいゆうほう)により因子を求めるやり方で、因子軸の回転はプロマックス回転，バリマックス回転か ら選ぶことが可能となってます。

メモ：最尤法が今は一番良いと言われてるそうです。うまく計算がいくならお勧めです。

（コピーして R にはりこみましょう)

ソース

factanal2 <- function(dat, factors=0, rotation=c("promax", "varimax", "none"), scores=c("none", "regression", "Bartlett"), from.uniquenesses=FALSE,

verbose=TRUE) {

dat <- dat[complete.cases(dat),] rotation <- match.arg(rotation) scores <- match.arg(scores)

p <- ncol(dat) if (factors == 0) {

factors <- max(1, floor((2*p+1-sqrt((2*p+1)^2-4*(p^2-p)))/2)) }

result<-factanal(dat, factors=factors, rotation=rotation, scores=scores) SS.loadings <- apply(result$loadings^2, 2, sum)

SS.loadings <- c(SS.loadings, sum(SS.loadings)) Proportion <- SS.loadings/p*100

Cum.Prop. <- cumsum(Proportion)