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

2007/1/11

目次

第1部:前書き 2

Rの初歩的操作 3

・データの入力 3

<別のデータのやりとり方法> 4

person.dataを利用して 5

平均値の算出 6

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

要約統計量 8

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

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

散布図を描こう 9

相関係数 9

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

by関数を用いる 10

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

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

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

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

t検定 13

ピアソンの相関係数 13

相関係数の差の検定 14

散布図 15

相関係数の検定 15

クラスター分析 15

平均を示す図を描こう 17

一要因の分散分析 17

多重比較! 19

ホルムの多重比較 19

・回帰分析 22

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

参考: 24

〜〜参考:統計処理の羅針盤〜 27

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

・質問紙を回収したら〜 27

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

尺度の構成へ 29

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

尺度の信頼性の検討 41

尺度得点の出し方 42

Rcommanderを使ってみよう 45

第3部:クロス集計表の分析 46

お勧め本! 46

比率の検定 46

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

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

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

<おまけです> 51

分散分析について 51

リンク 51

一要因の分散分析 52

〜〜ちょっと細かい話〜〜 55

2要因の分散分析 60

多重比較 65

2要因の分散分析(混合計画) 72

〜〜〜付記:重要な事柄〜〜 76

第4部:Rでt検定 77

対応のあるt検定 81

参考リンク 84


第1部:前書き

おかしなところがあったらぜひ緒賀( 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研究事例で学ぶSPSSAmosによる心理・調査データ解析 東京図書
田中 敏 ()  1996 実践心理データ解析―問題の発想・データ処理・論文の作成 新曜社

Rについては、以下の本とwebをお勧めします。

データ解析環境「R」―定番フリーソフトの基本操作からグラフィックス、統計解析まで    IO 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(関数名)を入れましょう。


・データの入力

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


 変数名も含めて、データ部分をコピーすると、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は、5因子性格の調和性得点。

bfconは、5因子性格の生真面目さ得点。

bfextは、5因子性格の外向性得点。

bfneurは、5因子性格の神経症傾向得点。

bfopenは、5因子性格の開放性得点。


<メモ>

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.dataepiE得点の平均値が返ってきます。


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


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

[1] 13.3


小数点2桁にしたいならば、上を1を2にするだけです。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)


です。もちろん、( )内に変数を指定できます。


小数点2桁にしたいならば

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


箱ひげ図を描いてみよう

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

>boxplot(person.data)


figureは、「メタファイル」形式でコピーをして、OpenOfficedrawにいったん貼り付けて、「切り離し」をすると、自由に加工ができます!
・コマンドラインで、↑のキイを押すと、前に打ち込んだものが使えて重宝します。


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


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

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


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


person.dataの中から、epiEからepiNeurまでの変数をとりだして、epi.dfというオブジェクトにいれました。epi.dfdata.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(データ名,条件式,関数) です。
以下は、虚偽尺度が3点より小さい人と、そうでない人を分けて、要約統計量を出す例です。
(条件式は、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

これでずいぶん便利になりましたね。




他にグラフィックで以下の関数がつかえるようになりました。下の2つを試してください。

multi.hist(person.data) 

	#draw the histograms for all variables
  
pairs.panels(person.data[1:5])
	#draws scatterplots, histograms, and shows correlations 

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

〜〜

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

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

>person.nolie.df<-subset(person.data,epilie<3)

とりあえず作っただけです。あとで使いません。
>rm(person.nolie.df)
で消去しましょう。

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

もし性別の変数があれば、同じようにして男だけのデータ、女だけのデータを取り出すことが可能です。
たとえば、originaldataというなかの変数sexMFがはいっていたとしたら次のようになります。
(data.frameでは、文字列は自動的に因子形式として扱われます。)

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

あるいは、女性だけを取り出すには、次のやりかたでもOKです。(注意! , を忘れないように。)
>femaledata.df<-originaldata[originaldata$sex==”F”,]

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

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

> x<-c(1:9) #xに1から9までの数字を入れたベクトルにする。combinationc です。
> 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を使うと、1 あるいは 0 という数字に変わります。
さいごに、1を足して、1、2、3 に変わりました。

epilieepiNeurの関係は、
>cor(epilie,epiNeur)
[1] -0.2508926

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

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


次に、epiEbfextの相関係数を確認

> cor(epiE,bfext)

[1] 0.5434974


そして、epiEbfagreeの相関係数を確認

> 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が低くても、相関が強いというわけではありません。単に無相関でないということしか言えません!)


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


いったん休憩・・・

>objects()

とすると、いままで作ったオブジェクトの一覧がでます。

それらのオブジェクトをすべてを消去するには、

>rm(list=ls())

というおまじないをしましょう。

一つだけなら

>rm(オブジェクトの名前)

rmremoveの略ですね。



クラスター分析

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


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

> plot(hc)



hclustは、クラスター分析の関数。

distは、類似性の計算のための関数です。












大雑把に見て、4クラスターに分かれそうですね。

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

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

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


> 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のマイナス2乗ということで、0.01

1e2は、1 x 10の2乗ということで、100 ということです。)

〜〜

別なやり方はこれ。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.

  1. 132.0 141.0 140.3 152.0 173.0


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


> 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による統計解析の基礎」ピアソン・エデュケーション,20031015日初版第1刷,A5判,183pp.ISBN4-89471-757-31,800円 の草稿の111ページから112ページを参照すると、「現在では,かなり広い用途をもち,ノンパラメトリックな分析にも適応可能なホルム(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


1群と3群、4群

2群と3群、4群

3群と4群 の間に差がありましたが、


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


〜〜〜〜〜〜〜

メモ:古典的な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


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

〜〜〜〜〜〜〜〜




さて、同じように、他の4つの性格得点ではどうなるかやってみましょう。>宿題!

4つのクラスターはそれぞれ、どんな人たちの集まりといえるでしょうか?



〜〜

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

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



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


−ーーーーー2元配置は別な機会にします。


・回帰分析

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

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

の「Rと回帰分析 」http://www1.doshisha.ac.jp/~mjin/R/13.html と「Rと重回帰分析」http://www1.doshisha.ac.jp/~mjin/R/14.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に標準化偏回帰係数のベクトルが得られる。
〜〜〜〜〜〜
参考2:
多重共線性のチェック(トレランスと 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 

〜〜〜〜〜〜

では、ビッグファイブの5性格特性から、ベックの抑うつ得点への予測を検討しましょう。
> 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


第2部:データの入力から、尺度の構成までのミニミニヒント


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

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


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


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


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

 もし因子分析を考えているのならば、項目数の10倍あるのが望ましいです。無理ならば、5〜6倍の人数は欲しいところです。


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

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)から順番にデータを入力していきます。エクセルやOpenOfficeCalcに入力します。

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

ID GENDER v1 v2 v3

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

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

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

という具合にです。


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


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

また、回答があっても、すべて1とか5とかのように、明らかに回答がいい加減だと判断される場合は、ケースを削除します。


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

ある分析をするときに、使われる変数の中に一つでも欠損値があったならば、その人のデータは利用しないというリスト消去。欠損値が入るものにペアになるものを使わずに計算する方法。欠損値を推定して、その数値(平均値や中央値や最頻値)を代入していく方法などです。



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

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





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


因子数を決めるには、固有値が1以上の基準(カイザーガットマン基準)や、固有値が落ちるところで決める基準(スクリープロット基準)などを参考にしていきます。

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


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





以下のソースをコピーして、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


>


相関係数、共通性、回転前と回転後の因子負荷量(マイナスになっていますが、全部ひっくりかえした意味として取って構いません)、また因子得点がでてきますね。(メモ:共通性+独自性は、1 です。)

小塩の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)
    Cum.Prop.[factors+1] <- NA
    Communality <- if (from.uniquenesses == TRUE) 1-result$uniquenesses else apply(result$loadings^2, 1, sum)
    result$cosmetic <- rbind(cbind(result$loadings, Communality), SS.loadings, Proportion, Cum.Prop.)    
    if (verbose == TRUE) {
        test <- c(result$STATISTIC, result$dof, result$PVAL)
        names(test) <- c("Chi sq.", "d.f.", "P value")
        print(paste("H0:", factors, "factors are sufficient."), quote=FALSE)
        print(test)
        print(paste("Factor loadings(rotation:", rotation, ")", sep=""), quote=FALSE)
        print(result$cosmetic)
        if (scores != "none") {
            print(paste("Factor scores(", scores, ")", sep=""), quote=FALSE)
            print(result$scores)
        }
    }
    invisible(result)
}

小塩のデータを分析してみましょう。
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/08_folder/da08_03.html

から、エクセルのデータをとりいれてください。

その上で、

> factanal2(x,factors=2)
[1] H0: 2 factors are sufficient.
   Chi sq.       d.f.    P value 
23.5811573 26.0000000  0.5999016 
[1] Factor loadings(rotation:promax)
                                                     		 Factor1     Factor2 Communality
F1_悩みを話し合えるような友人ができた            	0.6744282764  0.02671390  0.45556713
F2_たくさんの友人と一緒に遊ぶようになった        	0.7239945308 -0.21193606  0.56908498
F3_一生つきあっていけるような友人ができた        	0.4731980142  0.17938080  0.25609383
F4_グループで色々なことをするようになった       	-0.1482557366  0.40162922  0.18328579
F5_言いたいことを何でも言い合える友だちができた  0.1300917236  0.11609465  0.03040182
F6_みんなで一緒にいることが多くなった            	0.1274229194  0.35695006  0.14364994
F7_お互いに信頼できる友人ができた                	0.3831518386 -0.04508295  0.14883780
F8_たくさんの人と知り合いになった                	0.0002106384  0.59171507  0.35012677
F9_友達と心から理解し合えるようになった          	0.6071022245  0.17073966  0.39772514
F10_友達グループの一員になった                   	0.0104258902  0.81375437  0.66230488
SS.loadings                        	1.7735653477  1.42351275  3.19707810
Proportion                          17.7356534765 14.23512748 31.97078096
Cum.Prop.                           17.7356534765 31.97078096          NA
> 

注意:
promaxの因子間相関はRでは直接算出されていません。(JMPの斜交回転も同様ですので、必要ならSPSSあるいはエクセルのアドインを使いましょう。)

最近はpromax回転が勧められていることが多いのですが、個人的見解として、varimax回転で十分ではないかと思っています。promax回転のほうが解釈しやすい結果を得られやすいということなのですが、それは項目作成の段階で項目の吟味と洗練を十分にさせていたらクリアーできると思います。また、因子得点を使うよりも、α係数を出して、下位尺度得点を使う場合のほうが多い(独断か?)ので、下位尺得点間の相関を検討したほうが話が分かりやすいと思っています。でも結局のところ、データをうまく説明できることが大切ですから、いろいろと方法も回転も試行錯誤していくことが大事でしょう。たとえば、まずはバリマックスでやってみて、因子の意味がわかりにくかったら斜交回転をためしていくというようにです。


実際には、分析を進めていくには・・・・下の小塩先生の
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/da09_01.html
のように、試行錯誤していきます。

天井効果とフロアー効果のチェックには、以下でできます。
>mean(データ・フレーム)+sd(データ・フレーム)
mean(データ・フレーム)+sd(データ・フレーム)


> result <- pfa(x, factors=3)
> names(result)
[1] "rotation"             "correlation.matrix"   "communality"          "before.rotation.fl"   "before.rotation.eval" "after.rotation.fl"   
[7] "after.rotation.eval"  "scores"              
> result$after.rotation.fl
              [,1]        [,2]        [,3]
 [1,] -0.375220538 -0.07243387 -0.19864645
 [2,] -0.237844594 -0.04551322 -0.46257321
 [3,] -0.083470323 -0.52587951  0.07371994
 [4,]  0.077251990 -0.72368670 -0.37545857
 [5,] -0.331907365 -0.31998542 -0.40501972
 [6,]  0.105466360  0.49435425 -0.10585969
 [7,] -0.008339662 -0.23784657 -0.45809823
 [8,]  0.081420902 -0.07099637 -0.08756984
 [9,]  0.020539507  0.19300210 -0.55512605
[10,] -0.371322695 -0.11907617 -0.26214357
[11,] -0.348968366  0.01050690 -0.37001365
[12,] -0.248946799 -0.18468131 -0.40649953
[13,] -0.082494655 -0.07089483 -0.50444964
[14,] -0.093291323 -0.55662368 -0.35205662
[15,] -0.623913408 -0.12702767 -0.12159470
[16,] -0.522076782  0.20307303 -0.06739247
[17,] -0.431529797 -0.10856063 -0.36481338
[18,] -0.512858186 -0.05937604 -0.04970596
[19,] -0.244938506 -0.19254474 -0.31543140
[20,] -0.464779083 -0.05469920 -0.47180889
[21,]  0.099675561 -0.03666271  0.40268051
[22,]  0.538542643  0.30376842 -0.10686309
[23,] -0.337927424 -0.33415528 -0.26733222
> round(result$after.rotation.fl,3)
        [,1]   [,2]   [,3]
 [1,] -0.375 -0.072 -0.199
 [2,] -0.238 -0.046 -0.463
 [3,] -0.083 -0.526  0.074
 [4,]  0.077 -0.724 -0.375
 [5,] -0.332 -0.320 -0.405
 [6,]  0.105  0.494 -0.106
 [7,] -0.008 -0.238 -0.458
 [8,]  0.081 -0.071 -0.088
 [9,]  0.021  0.193 -0.555
[10,] -0.371 -0.119 -0.262
[11,] -0.349  0.011 -0.370
[12,] -0.249 -0.185 -0.406
[13,] -0.082 -0.071 -0.504
[14,] -0.093 -0.557 -0.352
[15,] -0.624 -0.127 -0.122
[16,] -0.522  0.203 -0.067
[17,] -0.432 -0.109 -0.365
[18,] -0.513 -0.059 -0.050
[19,] -0.245 -0.193 -0.315
[20,] -0.465 -0.055 -0.472
[21,]  0.100 -0.037  0.403
[22,]  0.539  0.304 -0.107
[23,] -0.338 -0.334 -0.267
>

いくつかの項目を落として、新しいデータフレームをつくる。
8番、19番、23番を落としてみます。

> x2 <- transform(x[,c(1:7,9:18,20:22)])

あるいは

> x2 <- transform(x[,-c(8,19,23)])

上のいずれかのやりかたでOK

> result <- pfa(x2,factors=3)
> round(result$after.rotation.fl,3)
        [,1]   [,2]   [,3]
 [1,] -0.369 -0.063 -0.204
 [2,] -0.227  0.003 -0.454
 [3,] -0.066 -0.498  0.071
 [4,]  0.067 -0.676 -0.374
 [5,] -0.339 -0.310 -0.401
 [6,]  0.087  0.497 -0.096
 [7,] -0.014 -0.196 -0.425
 [8,]  0.016  0.198 -0.575
 [9,] -0.393 -0.127 -0.265
[10,] -0.345  0.017 -0.367
[11,] -0.270 -0.218 -0.414
[12,] -0.101 -0.089 -0.509
[13,] -0.111 -0.627 -0.381
[14,] -0.584 -0.127 -0.124
[15,] -0.567  0.200 -0.020
[16,] -0.462 -0.119 -0.327
[17,] -0.470 -0.036 -0.040
[18,] -0.501 -0.054 -0.441
[19,]  0.107 -0.031  0.419
[20,]  0.552  0.344 -0.114

通し番号が変わったので注意してください。

以下のほうが分かりやすいでしょう。

> > factanal2(x2, factors=3)
[1] H0: 3 factors are sufficient.
    Chi sq.        d.f.     P value 
136.1603418 133.0000000   0.4078537 
[1] Factor loadings(rotation:promax)
                                                                              		Factor1     Factor2      Factor3 Communality
A01_人の話を黙って聞いていられない                                       		0.3338075025  0.13928700 -0.012273603   0.1309790
A02_自分はすぐ調子に乗るタイプだ                                         		0.1012982485  0.45519291  0.005729681   0.2174948
A03_一人で行動することが好きではない                                     		0.0653527173 -0.17006647  0.471365950   0.2553794
A04_友達がやるからという理由で何をするか決めることが多い       	 	-0.2260173028  0.31007797  0.788105487   0.7683424
A05_気分によって考えがコロコロ変わる                                     		0.2507305397  0.32813969  0.300839330   0.2610458
A06_分からないことがあったとき.まずは自分で考えたい\0      		-0.1153280008  0.20863513 -0.437849224   0.2485411
A07_買い物に行ったとき.予定していなかった物まで買ってしまうことがある\0 	-0.1207874244  0.42066340  0.281775929   0.2709450
A09_夏休みの宿題は終わりがけにあわててやっていた                        	-0.1640962420  0.63075164 -0.078364838   0.4309163
A10_自分の失敗を周りのせいにすることが多い                               	0.3472582730  0.18547936  0.089065419   0.1629235
A11_約束の時間によく遅刻をしてしまう                                     		0.2907224402  0.31770518 -0.066459201   0.1898729
A12_同じ失敗を繰り返してしまうことが多い                                 		0.1894009383  0.34413708  0.208360922   0.1977173
A13_勉強をしているときに他のことに気をとられやすい                      	-0.0212535386  0.50708194  0.130680567   0.2746612
A14_自分一人で何かを決めることが苦手だ                                  	 	0.0002202487  0.28142865  0.665039408   0.5214796
A15_思っていることや感情が.ついつい表情に出てしまう\0                	0.5952679414  0.00858751 -0.015873974   0.3546697
A16_人に指図されるのは嫌いだ                                             		0.6499970910 -0.03610451 -0.370567521   0.5611200
A17_自分のやりたくないことはやらないことが多い                           		0.4427186867  0.25198751  0.025415099   0.2601435
A18_授業中そわそわして落ち着かない                                       		0.4778510376 -0.02556823 -0.077617529   0.2350198
A20_自分に都合が良くない考えは受け入れないことが多い                    	 0.4442337411  0.37584960 -0.007347853   0.3386605
A21_お金は先のことを考えて使っている                                    		-0.0226219414 -0.41322823  0.007505726   0.1713257
A22_カッとなることがあっても冷静になって考え直すことができる            	-0.6559797286  0.28787864 -0.156818127   0.5377754
SS.loadings                                                 2.3655057980  2.12820903  1.895298059   6.3890129
Proportion                                                 11.8275289899 10.64104513  9.476490297  31.9450644
Cum.Prop.                                                  11.8275289899 22.46857412 31.945064417          NA


〜〜
ところで、もともとRで実装されているのは、factanal()です。

> result <-factanal(x2,factors=3,rotation=”promax”)
> names(result)
 [1] "converged"    "loadings"     "uniquenesses" "correlation"  "criteria"     "factors"      "dof"          "method"      
 [9] "STATISTIC"    "PVAL"         "n.obs"        "call"        

result$何々 で、細かい結果が示されます。

> print(result$loading, cutoff=0.001, sort=T)

 cutoffのオプションで、表示しない絶対値を指定できます。
  sortのオプションで、並び替えをできますが・・・いまひとつ変な感じです。

出力は以下のとおり。

いったん、witerに貼り付けて、範囲選択をしたうえで次のことをしてください。

編集→検索と置換→詳細オプションで、正規表現にチェックを入れる→検索のボックスに、 +  (半角スペースにプラス)、置換ボックスに、\t(円記号にt 意味はタブです)→すべて置換→ コピーして、Calcに形式を指定して貼り付け(書式設定のないテキストを選んでください) 必要な修正を加えた上で、コピーして、writerに貼り付けます。整形したのは以下のとおりです。
               Factor1 Factor2 Factor3
SS loadings       2.38    2.31    1.81
Proportion Var    0.12    0.12    0.09
Cumulative Var    0.12    0.23    0.32
> 
こんなこともできます。

>  plot(result$loading)                                 
>  identify(result$loading,labels=names(x2))           #

グラフのポイントをクリックしてみてください。
終わるときは、グラフ画面を閉じましょう。

尺度の信頼性の検討

何のことなのかは、以下のページを見てください。
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/da09_02.html

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

特に、因子得点と尺度得点は違いますので、注意してください。
因子得点を使う場合は、α係数を求める必要性はありません。
〜〜
余談:
因子分析においてアルファ係数に代わるものとして、オメガがあります。
http://personality-project.org/r/r.omega.html

〜〜

尺度得点の出し方

〜〜メモ
 エクセルやCalcで出すときには、欠損値の扱いに注意しましょう。
たとえば以下のようにします。
=IF(COUNTBLANK(A1:C1)>0; ""; SUM(A1:C1))
これは、欠損値(ブランク)があれば、ブランク、そうでなければ合計を求めるという式です。

逆転項目については、以下の式で、新たに変数をつくります。(5件法の場合)
=IF(D1<>"";6-D1;"")
〜〜

Rには、α係数を求める関数は実装されていませんので、以下のものを使います。
http://personality-project.org/r/alpha.html

ソース
#define a function to calculate coefficient alpha
alpha.scale=function (x,y)   #find coefficient alpha given a scale and a data.frame of the items in the scale
        {
                n=length(y)          #number of variables
                Vi=sum(diag(var(y,na.rm=TRUE)))     #sum of item variance
                Vt=var(x,na.rm=TRUE)                #total test variance
                ((Vt-Vi)/Vt)*(n/(n-1))}              #alpha
alpha.scaleという関数を得ました。

小塩の以下のデータを用いてみます。
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/09items_folder/data_youji.xls

表計算上で、コピーして。

> oshio.df <- read.delim("clipboard")

> names(oshio.df)
 [1] "A01_人の話を黙って聞いていられない"                                     
 [2] "A02_自分はすぐ調子に乗るタイプだ"                                       
 [3] "A03_一人で行動することが好きではない"                                   
 [4] "A04_友達がやるからという理由で何をするか決めることが多い"               
 [5] "A05_気分によって考えがコロコロ変わる"                                   
 [6] "A06_分からないことがあったとき.まずは自分で考えたい\0"                  
 [7] "A07_買い物に行ったとき.予定していなかった物まで買ってしまうことがある\0"
 [8] "A08_他人のペースに合わせることができる"                                 
 [9] "A09_夏休みの宿題は終わりがけにあわててやっていた"                       
[10] "A10_自分の失敗を周りのせいにすることが多い"                             
[11] "A11_約束の時間によく遅刻をしてしまう"                                   
[12] "A12_同じ失敗を繰り返してしまうことが多い"                               
[13] "A13_勉強をしているときに他のことに気をとられやすい"                     
[14] "A14_自分一人で何かを決めることが苦手だ"                                 
[15] "A15_思っていることや感情が.ついつい表情に出てしまう\0"                  
[16] "A16_人に指図されるのは嫌いだ"                                           
[17] "A17_自分のやりたくないことはやらないことが多い"                         
[18] "A18_授業中そわそわして落ち着かない"                                     
[19] "A19_欲しいと思ったものは手に入らないと気が済まない"                     
[20] "A20_自分に都合が良くない考えは受け入れないことが多い"                   
[21] "A21_お金は先のことを考えて使っている"                                   
[22] "A22_カッとなることがあっても冷静になって考え直すことができる"           
[23] "A23_楽しいときにはどこでもはしゃぐ"                                     
> 


http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/da09_02.html

にあるように、

「第3因子に高い負荷量を示した項目は

 A04A14A06(逆転),A03の4項目。

 A06は「負の負荷量」を示しているので,「逆転項目」とする」尺度をつくって、α係数を出してみましょう。



変数名が長くて扱いにくいので、表計算で変数名を変えてから、もう一度読み込みましょう。

A1からA23という変数名にしましょう。

(したの表をOpenOfficeを使っていたら、ダブルクリックしましょう。)




> oshio.df <- read.delim("clipboard")

> names(oshio.df)

[1] "A1" "A2" "A3" "A4" "A5" "A6" "A7" "A8" "A9" "A10" "A11" "A12"

[13] "A13" "A14" "A15" "A16" "A17" "A18" "A19" "A20" "A21" "A22" "A23"

>



仮にS3という尺度名にします。

まず

>attach(oshio.df)

として、直接に変数名を扱えるようにします。

S3の得点を求めましょう。A&は逆転項目ですので、6から引き算をした数を加えます。

> S3 <- A4+A14-A6+A3+6
> S3
  [1] 12  7  8 15 11  8 10 15  9  7 13  5 11  9 13  6 13 10 12 13 18 14 12 11  6
 [26] 13 14  7 11 11  8  6 14  9 10 12 16  8 14 15  7 12 14 10 19 13 13  8 15 12
 [51]  4  8  8 13  4 10  7 14 12 10 10 12 14  8 14  9 11 10  9  8  7 15 11 15  8
 [76] 10 10 11  5  7  8 11 13 16 12 12  9 15 11  8 10 10 10 12 11 13 13 12 14 19
[101] 16  5  9  9 13 10 14 11 11 14 11
> 

S3α係数を求めるためには、以下のデータフレームを作りましょう。

> S3.df <- data.frame(A4,A14,A6,A3)
> summary(S3.df)
       A4             A14              A6              A3       
 Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
 1st Qu.:2.000   1st Qu.:2.000   1st Qu.:3.000   1st Qu.:2.000  
 Median :3.000   Median :3.000   Median :4.000   Median :2.000  
 Mean   :3.072   Mean   :3.054   Mean   :3.658   Mean   :2.477  
 3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:3.000  
 Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.000  
> 

S3を構成する項目間相関をみてみましょう。

> round(cor(S3.df,use="pair"),2)    
       A4   A14    A6    A3
A4   1.00  0.57 -0.22  0.35
A14  0.57  1.00 -0.31  0.23
A6  -0.22 -0.31  1.00 -0.29
A3   0.35  0.23 -0.29  1.00
> 

A&が逆転項目になるのはわかりますね。

S3尺度得点と項目の相関を見てみましょう。

> round(cor(S3.df,S3,use="pair"),2)  
     [,1]
A4   0.76
A14  0.75
A6  -0.62
A3   0.68
> 

さて、α係数を求めると以下のようになります。

> alpha.scale(S3,S3.df) 
[1] 0.6605016
> 

このalpha.scale関数では、スケールを構成するデータフレームの中に逆転項目をそのまま残しても良いです。
大事なのは、尺度得点を求めるときに、ちゃんと逆転項目を変換していることです。

Useful.rを取り入れていたらdescribeを使ってみましょう。

> describe(S3)
  V   n  mean   sd median min max range skew  se
1 1 111 10.95 3.12     11   4  19    15 0.06 0.3
> 

Rcommanderを使ってみよう

SPSS風のGUI画面で使いたいという人には、Rcommanderというパッケージがあります。

Rcmdrを選んでインストールしてください。

コンソール画面で

> library(Rcmdr)

と打つと、GUI画面が立ち上がります。


使い方は、

http://oscar.gsid.nagoya-u.ac.jp/tech/wiki/wiki.cgi?page=R%A4%CE%BB%C8%A4%A4%CA%FD が参考になります。


それをインストールすると、reliability関数が使えます。

> help(reliability)

> data(DavisThin)

> reliability(cov(DavisThin))

Alpha reliability = 0.8919

Standardized alpha = 0.8888


Reliability deleting each item in turn:

Alpha Std.Alpha r(item, total)

DT1 0.8997 0.8991 0.4512

DT2 0.8720 0.8681 0.7267

DT3 0.8598 0.8566 0.8139

DT4 0.8902 0.8861 0.5698

DT5 0.8660 0.8630 0.7716

DT6 0.8601 0.8576 0.8117

DT7 0.8758 0.8728 0.6934

>

第3部:クロス集計表の分析


お勧め本!

クイック・データアナリシス―10秒でできる実践データ解析法
田中 敏 (), 中野 博幸 () 出版社: 新曜社 ; ISBN: 4788509199 ; (October 2004)



この本いいですね〜。絶対おすすめします!


比率の検定

1ページの「今回の授業は今までの授業よりよかったですか」の結果を分析してみましょう。

生徒数が20名で、ハイが15名、イイエ5名です。



二項検定 binom.test


これは、正確な比率の検定です。


> binom.test(15,20)


Exact binomial test


data: 15 and 20

number of successes = 15, number of trials = 20, p-value = 0.04139

alternative hypothesis: true probability of success is not equal to 0.5

95 percent confidence interval:

0.5089541 0.9134285

sample estimates:

probability of success

0.75


両側検定(デフォールト)で、5%水準で有意に「ハイ」が多いですね。

もし片側だと、その半分の水準になります。



> binom.test(15, 20, alternative="greater")


Exact binomial test


data: 15 and 20

number of successes = 15, number of trials = 20, p-value = 0.02069

alternative hypothesis: true probability of success is greater than 0.5

95 percent confidence interval:

0.5444176 1.0000000

sample estimates:

probability of success

0.75


さて6人と14人ではどうなるでしょうか?>宿題

なお、オプションで母比率は変えることもできます。


P34の例を分析してみましょう。


女性

男性

40代管理職の人数

12人

10人

英語教師の母比率

7割

3割


> binom.test(12,22, p=0.7,alternative="less")


Exact binomial test


data: 12 and 22

number of successes = 12, number of trials = 22, p-value = 0.0916

alternative hypothesis: true probability of success is less than 0.7

95 percent confidence interval:

0.0000000 0.7286868

sample estimates:

probability of success

0.5454545


片側で、9.16%ですね。


この分析の実践は、上記の本にいろいろ参考になること(符号検定)が載っています。

加えて、説明のある自動集計検定は、とても便利そうです。

以下のリンクのJAVA-SCRIPT STARを使ってみてください。

http://www.kisnet.or.jp/nappa/software/star/


~~おまけ

普通に行われているΧ二乗検定は、

chisq.test()

です。以下のリンクをみてください。

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

でも、近似値計算ではない、フィッシャーの検定をできるだけ行いましょう。


あと比率の検定には、

prop.test() というのもあります。

http://www.okada.jp.org/RWiki/index.php?cmd=read&page=Ķ%CC%F5%A1%A7%C8%E6Ψ%A4κ%B9%A4θ%A1%C4%EA&word=prop.test

〜〜


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

次は、クロス集計表の分析です。

フィッシャーの直接確率検定 でやっていきましょう。



>fisher.test() は、引数にmatrix(行列)形式のデータをとります。

そこで、まずはmatrixを作成する方法を学びましょう。


1,2,3,4,5,6の数字を2行3列の行列にすることで学びます。

まず、ひとつのベクトルにします。これはc(1,2,3,4,5,6)です。それにmatrix関数を使います。


>matrix(c(1,2,3,4,5,6), nrow=2, ncol=3, byrow=T)


結果は以下のとおり。

[,1] [,2] [,3]

[1,] 1 2 3

[2,] 4 5 6


nrowは、number of row ということで行数を指定し、

ncolは、number of colum ということで列数を指定し、

byrow=T は、行ごとに、言い換えると、横ならびで入れていきます。省略すれば、縦ならびになります。(1から6までの数字の並びを、1:6 とも指定できます。)


> matrix(c(1,2,3,4,5,6), nrow=2, ncol=3)

[,1] [,2] [,3]

[1,] 1 3 5

[2,] 2 4 6


あるいは、

> matrix(1:6, nrow=2, ncol=3)

[,1] [,2] [,3]

[1,] 1 3 5

[2,] 2 4 6


でも同じです。見やすいのは、byrow=Tのオプションをつけたほうでしょう。

Nrowncolはどちらかを省略できます。


> matrix(1:6, nrow=2, byrow=T)

[,1] [,2] [,3]

[1,] 1 2 3

[2,] 4 5 6

> matrix(1:6, ncol=3, byrow=T)

[,1] [,2] [,3]

[1,] 1 2 3

[2,] 4 5 6

>

どちらも同じ行列が作られました。


さて、フィッシャーの直接確立検定をします。

次は、私の指導生が卒業論文でとったデータです。対象の内的想起の4タイプと、携帯電話依存の4タイプの行列です。

> x <- matrix(c(11,13,6,8,12,7,5,6,7,8,4,1,4,4,2,2),ncol=4,byrow=T)

> x

[,1] [,2] [,3] [,4]

[1,] 11 13 6 8

[2,] 12 7 5 6

[3,] 7 8 4 1

[4,] 4 4 2 2

列に内的対象想起のタイプ、行に携帯電話依存のタイプとなっています。

> fisher.test(x)

以下にエラーfisher.test(x) : FEXACT error 6.

LDKEY is too small for this problem.

Try increasing the size of the workspace.

>

エラーが出ました。計算するのに必要なメモリーを確保していなかったようです。

workspace=3000000 というように、メモリー量を指定して行ってみましょう。


> fisher.test(x,workspace=3000000)


Fisher's Exact Test for Count Data


data: x

p-value = 0.8859

alternative hypothesis: two.sided


今度は、計算されました・・・・が、関係は認められなかったですね。


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


「テレビゲームの暴力表現は青少年に悪影響を及ぼすか?」を100名の男子高校生に質問をしたあと、ゲームで覚えた技を使いたくて男性を殺害した事件が発生した。その事件の前と後で、同じ100人に、それぞれ質問を繰り返した結果が以下の表となったそうです。同じ人に質問を繰り返したので、4つの回答パターンにわかれます。その人数です。

例題から分かる心理統計学

この本もお勧めです。



影響あり(事件前)

影響なし(事件前)

影響あり(事件後)

49

16

65

影響なし(事件後)

7

28

35

56

44

100

例題からわかる心理統計学 p46より



> x <- matrix(c(49, 16, 7, 28), ncol=2, byrow=T)

> x

[,1] [,2]

[1,] 49 16

[2,] 7 28

> mcnemar.test(x, correct=F)


McNemar's Chi-squared test


data: x

McNemar's chi-squared = 3.5217, df = 1, p-value = 0.06057


>

6%ですので、殺人事件の報道をきいて意見が変わったとは言えないという結果ですね。


correctのオプションの意味は自分で調べてください。)


カテゴリデータの検定を参考にしてください。


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


モザイク図を描きたいときには、下のページが参考になります。


http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%A5%B0%A5%E9%A5ե%A3%A5å%AF%A5%B9%BB%B2%B9ͼ%C2%CE%A1%A7%A5⥶%A5%A4%A5%AF%A5ץ%ED%A5å%C8&word=%A5⥶%A5%A4%A5%AF



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

κ (カッパ) 係数を計算する際は,パッケージ vcd を呼び出した後で関数 Kappa() を実行.

信頼性係数 ICC (Intraclass correlation coefficient) を求める場合は,パッケージ irr 中の関数 icc() を用い

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





<おまけです>


心理統計をする上での重要なポイント


http://www.okayama-u.ac.jp/user/le/psycho/member/hase/articles/1994/9407Hasegawa/9407Hasegawa.html

心理学研究法再考

(1)基礎的統計解析の誤用をなくすための30のチェック項目



分散分析について

リンク

基本の説明があります。

http://www.shiga-med.ac.jp/~koyama/stat/com-anova.html


2要因以上になるとややこしくなる分散分析ですが、以下のリンクを参考にしてください。

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

上は、分散分析をするときの注意点です。


「心理学のためのデータ解析テクニカルブック」の例題をRで書いてみたというものが下にあります。

http://mat.isc.chubu.ac.jp/R/tech.html


雑誌に載ったテキストの分散分析の部分です。

http://www1.doshisha.ac.jp/~mjin/R/12.pdf


心理学者用のページです。

繰り返しのある場合の分散分析の例も載っています。

http://personality-project.org/r/r.guide.html#anova

特にそのなかの「R and Analysis of Variance」は以下のリンク。

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


#Statistical Computing Seminars Repeated Measures Analysis using SAS

#http://www.ats.ucla.edu/stat/sas/seminars/sas_repeatedmeasures/default.htm

それをRで書いたコードが以下にあるます。

http://cat.zero.ad.jp/~zak52549/R/anova/repeated.txt


Rでラクラク分散分析」という研修テキストです。

http://cse.niaes.affrc.go.jp/minaka/R/Ranova-nlbc04.pdf


交互作用が出たときの理解のためのヒントです。

http://www.juen.ac.jp/psych/nakayama/anova.html


Practical Statistics with R ・・・参考になるテキストです。

http://www.brandonu.ca/zoology/rutherford/RBookWebSite/RBook_pdf.html




以上のリンク先を参考にしましょう。

〜〜〜〜〜〜


さて、基本からはじめましょう。


一要因の分散分析


http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/04_folder/da04_02.html

にあるデータを使います。(被験者間データです。)







というデータです。ダブルクリックして、Calcで開いて(OpenOfficeを使ってる場合)、クリップボードにコピーしましょう。


> x.df <- read.delim("clipboard")

> x.df

番号 条件 結果

1 1 1 4

2 2 1 1

3 3 1 3

4 4 1 2

5 5 1 2

6 6 1 4

7 7 1 3

8 8 2 6

9 9 2 8

10 10 2 5

11 11 2 9

12 12 2 8

13 13 2 7

14 14 2 7

15 15 3 4

16 16 3 3

17 17 3 4

18 18 3 6

19 19 3 5

20 20 3 5

21 21 3 5

>



まず最初にやること、条件の変数を、因子タイプに変えること!

> x.df$条件 <- factor(x.df$条件)

> class(x.df$条件)

[1] "factor"


分散分析は英語では、analysis of variance ですね。頭文字をとると、aov。その名前のオブジェクトを使います。


> result <- aov(結果~条件, data=x.df)


分析結果を、resultにいれておきました。どんな種類の結果がresultの中にあるかをnames()で調べます。

Aov( )の中の式は、~ チルダを使います。



> names(result)

[1] "coefficients" "residuals" "effects" "rank"

[5] "fitted.values" "assign" "qr" "df.residual"

[9] "xlevels" "call" "terms" "model"


通常、必要な結果は、summary( )で、以下の出てきます。

> summary(result)

Df Sum Sq Mean Sq F value Pr(>F)

条件 2 69.238 34.619 25.964 4.961e-06 ***

Residuals 18 24.000 1.333

---

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

>


繰り返し強調しますが、要因となるのは、factor( )で、因子タイプに変えてから分析してください。でないと、おかしなことになります。


要因ごとの、基本統計量を describebyを使って確認しておきましょう。


> by(x.df$結果, x.df$条件, describe)

INDICES: 1

V n mean sd median min max range skew se

1 1 7 2.71 1.11 3 1 4 3 -0.15 0.42

------------------------------------------------------------

INDICES: 2

V n mean sd median min max range skew se

1 1 7 7.14 1.35 7 5 9 4 -0.22 0.51

------------------------------------------------------------

INDICES: 3

V n mean sd median min max range skew se

1 1 7 4.57 0.98 5 3 6 3 -0.17 0.37

>


有意な結果になっていますので、Tukeyの下位検定をしてみます。分散分析の結果を入れたresultTukeyHSD( )の中にいれます。


> TukeyHSD(result)

Tukey multiple comparisons of means

95% family-wise confidence level


Fit: aov(formula = 結果 ~ 条件, data = x.df)


$条件

diff lwr upr p adj

2-1 4.428571 2.8533421 6.0038007 0.0000032

3-1 1.857143 0.2819136 3.4323722 0.0196279

3-2 -2.571429 -4.1466579 -0.9961993 0.0015998


>



どの条件間においても差が認められますね。


フィッシャーのLSDを使いたいならば、

> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”none”)


Bonferroniの手法ならば、

> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”bonferroni”)


Holmの手法ならば

> pairwise.t.test(x.df$結果, x.df$条件, p.adjust.method=”holm”)



グラフにしたいのならば、まず条件ごとの平均値を、y にいれましょう。

> y <- by(x.df$結果, x.df$条件, mean)

> y

INDICES: 1

[1] 2.714286

------------------------------------------------------------

INDICES: 2

[1] 7.142857

------------------------------------------------------------

INDICES: 3

[1] 4.571429


次に、plot( )を使います。


> plot(y)


折れ線にしたければ、line の”l” typeで指定します。


> plot(y , type="l")




〜〜ちょっと細かい話〜〜


分散分析は、各群の分散が等質性を前提としています。その前提をチェックするには、bartlett.test() を用いることができます。

> bartlett.test(x.df$結果, x.df$条件)


Bartlett test of homogeneity of variances


data: x.df$結果 and x.df$条件

Bartlett's K-squared = 0.5877, df = 2, p-value = 0.7454


上の結果は等分散ということですが、もし等分散でなかったならば、oneway( )で検定できます。オプションで、var=F と明示するか、書かずに分析してみましょう。


なお、正規性の検定は、shapiro.test() Shapiro-Wilk 検定 を, ks.test() Kolmogorov-Smirnov 検定 を行えます。分散分析は、正規性については頑強だといわれてます。



〜〜〜〜


さて、次は被験者内要因の場合をみていきましょう。

同じく小塩のデータを使います。

「要因の分散分析(被験者内計画)

という例です。


> x.df <- read.delim("clipboard")

> x.df

subject r30 r60 r90 r120

1 1 42 39 36 34

2 2 22 17 15 8

3 3 35 32 25 25

4 4 34 30 22 20

5 5 40 33 28 23

6 6 34 30 23 28

>


この形でRは分析しませんので、データの形を変えます。


> subject <- c(rep(1,4), rep(2,4), rep(3,4), rep(4,4), rep(5,4), rep(6,4))

> subject

[1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6

>



> condition <- rep(c(1:4), 6)

> condition

[1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

>


> sakusi <- c(42,39,36,34,

+ 22,17,15, 8,

+ 35,32,25,25,

+ 34,30,22,20,

+ 40,33,28,23,

+ 34,30,23,28)

> sakusi

[1] 42 39 36 34 22 17 15 8 35 32 25 25 34 30 22 20 40 33 28 23 34 30 23 28

>


> subject <- factor(subject)

> condition <- factor(condition)


要因は、factorタイプに変えます。


> x <- data.frame(subject, condition, sakusi)

> x

subject condition sakusi

1 1 1 42

2 1 2 39

3 1 3 36

4 1 4 34

5 2 1 22

6 2 2 17

7 2 3 15

8 2 4 8

9 3 1 35

10 3 2 32

11 3 3 25

12 3 4 25

13 4 1 34

14 4 2 30

15 4 3 22

16 4 4 20

17 5 1 40

18 5 2 33

19 5 3 28

20 5 4 23

21 6 1 34

22 6 2 30

23 6 3 23

24 6 4 28

>


被験者内要因が入るときには、+Error( / )を付け加えます。


〜〜〜ちょっと詳しく〜〜

Mauchlyの球面性検定のやり方は謎??なので、仮定を満たしたものとしてやってみます。

   mauchly.test(stats) Mauchly's Test of Sphericity があり、SSD(stats)を使って前処理したあとで使うようですが・・ちょっと使い方を分かりかねています。

  また、http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#SECTION000710000000000000000 6.10 Sphericity にthe Greenhouse-Geisser and the Huynh-Feldt corrections イプシロン算出の方法の記載があります。Greenhous-Geisser と Huygh-Feldt の修正は、自由度にイプシロンをかけるということです。それでF値が決まり、pf( )関数を使って、p値を求めます。


http://www.psych.upenn.edu/~baron/rpsych/rpsych.html Errorの使い方を含めて分散分析のやり方が詳しくあります。

〜〜〜


> result <- aov(sakusi~condition+Error(subject/condition))

> summary(result)


Error: subject

Df Sum Sq Mean Sq F value Pr(>F)

Residuals 5 1058.37 211.67


Error: subject:condition

Df Sum Sq Mean Sq F value Pr(>F)

condition 3 491.46 163.82 32.855 7.766e-07 ***

Residuals 15 74.79 4.99

---

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

>


角度の主効果は(3,15)=32.86で,0.1%水準で有意でした。


平均値を出してグラフにしてます。


> y <- by(x$sakusi, condition, mean)

> y

INDICES: 1

[1] 34.5

------------------------------------------------------------

INDICES: 2

[1] 30.16667

------------------------------------------------------------

INDICES: 3

[1] 24.83333

------------------------------------------------------------

INDICES: 4

[1] 23



> plot(y, type="o")


>library(gplots)

で、plotmeansが使えるならば以下のもOK


> plotmeans(sakusi~condition)


次の下位検定ですが、Bonferroniの修正(補正)をほどこしましょう。

修正方法は、6回の対応のあるt検定ですので、p X 6 で修正します。

> x.df <- read.delim("clipboard")

> x.df

subject r30 r60 r90 r120

1 1 42 39 36 34

2 2 22 17 15 8

3 3 35 32 25 25

4 4 34 30 22 20

5 5 40 33 28 23

6 6 34 30 23 28

>

としたあとで、以下のスクリプトを一気にコンソールに流し込んでみましょう。


result.t <- t.test(r30,r60, paired=T)

result.t$p.value*6

result.t <- t.test(r30,r90, paired=T)

result.t$p.value*6

result.t <- t.test(r30,r120, paired=T)

result.t$p.value*6

result.t <- t.test(r60,r90, paired=T)

result.t$p.value*6

result.t <- t.test(r60,r120, paired=T)

result.t$p.value*6

result.t <- t.test(r90,r120, paired=T)

result.t$p.value*6


結果は以下のとおり


> result.t <- t.test(r30,r60, paired=T)

> result.t$p.value*6

[1] 0.005322308

> result.t <- t.test(r30,r90, paired=T)

> result.t$p.value*6

[1] 0.001551205

> result.t <- t.test(r30,r120, paired=T)

> result.t$p.value*6

[1] 0.006571619

> result.t <- t.test(r60,r90, paired=T)

> result.t$p.value*6

[1] 0.01774532

> result.t <- t.test(r60,r120, paired=T)

> result.t$p.value*6

[1] 0.01621808

> result.t <- t.test(r90,r120, paired=T)

> result.t$p.value*6

[1] 1.982593

>


確率は1を超えることはありえませんので、r90r120p値は、1.000と理解することにしましょう。

〜〜〜

ちなみに、JavaScript-STAR で分析した結果は以下の通り。LSDによる多重比較結果がでています。

== Analysis of Variance ==


S.V SS df MS F

------------------------------------------------

Sub 1058.3750 5 211.6750

------------------------------------------------

A 491.4583 3 163.8194 32.86 **

SxA 74.7917 15 4.9861

------------------------------------------------

Total 1624.6250 23 +p<.10 *p<.05 **p<.01


== Multiple Comparisons by LSD ==


(MSe= 4.9861, * p<.05)

(LSD= 2.7479)

------------------------------------------------

A1 > A2 *

A1 > A3 *

A1 > A4 *

A2 > A3 *

A2 > A4 *

A3 = A4 n.s.

------------------------------------------------


_/_/_/ Analyzed by JavaScript-STAR _/_/_/

〜〜〜


2要因の分散分析

同じく小塩のページから・・


2要因の分散分析(どちらも被験者間要因)







> x.df <- read.delim("clipboard")

> x.df

番号 失敗経験 完全主義 抑うつ

1 1 1 1 10

2 2 1 1 13

3 3 1 1 21

4 4 1 1 16

5 5 1 2 16

6 6 1 2 19

7 7 1 2 13

8 8 1 2 8

9 9 2 1 15

10 10 2 1 16

11 11 2 1 12

12 12 2 1 15

13 13 2 2 21

14 14 2 2 23

15 15 2 2 16

16 16 2 2 19

17 17 3 1 21

18 18 3 1 14

19 19 3 1 24

20 20 3 1 20

21 21 3 2 31

22 22 3 2 36

23 23 3 2 24

24 24 3 2 34

>


> attach(x.df)


変数名を直接扱えるようにしました。


> 失敗経験 <- factor(失敗経験)

> 完全主義 <- factor(完全主義)

> class(失敗経験)

[1] "factor"

> class(完全主義)

[1] "factor"

>


必ず、factorタイプに要因は変更!


交互作用を考えて分析してみます。(交互作用を入れなくても分析はできます。)


> y <- aov(抑うつ~失敗経験+完全主義+失敗経験*完全主義)


注意: 抑うつ~失敗主義*完全主義  という形も同じことをあらわします。



> summary(y)

Df Sum Sq Mean Sq F value Pr(>F)

失敗経験 2 528.08 264.04 15.6727 0.0001143 ***

完全主義 1 165.37 165.37 9.8162 0.0057477 **

失敗経験:完全主義 2 156.25 78.13 4.6373 0.0237486 *

Residuals 18 303.25 16.85

---

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

>


失敗経験の主効果、完全主義の主効果の両方とも認められますが、交互作用があるみたいです。

グラフを見てみましょう。

> interaction.plot(失敗経験, 完全主義, 抑うつ, fixed=TRUE)

>





それぞれのセル/ブロックの平均値をだして、グラフを再確認しましょう。

model.tables( )を使います。引数の”means”を忘れないように。


> model.tables(y, "means")

Tables of means

Grand mean

19.04167


失敗経験

失敗経験

1 2 3

14.500 17.125 25.500


完全主義

完全主義

1 2

16.417 21.667


失敗経験:完全主義

完全主義

失敗経験 1 2

1 15.00 14.00

2 14.50 19.75

3 19.75 31.25


上の前半部分は、以下でも出ます。


> by(抑うつ,失敗経験,describe)

INDICES: 1

V n mean sd median min max range skew se

1 1 8 14.5 4.38 14.5 8 21 13 0 1.55

-------------------------------------------------------------------

INDICES: 2

V n mean sd median min max range skew se

1 1 8 17.12 3.6 16 12 23 11 0.3 1.27

-------------------------------------------------------------------

INDICES: 3

V n mean sd median min max range skew se

1 1 8 25.5 7.56 24 14 36 22 0.04 2.67

> by(抑うつ, 完全主義, describe)

INDICES: 1

V n mean sd median min max range skew se

1 1 12 16.42 4.21 15.5 10 24 14 0.3 1.22

-------------------------------------------------------------------

INDICES: 2

V n mean sd median min max range skew se

1 1 12 21.67 8.49 20 8 36 28 0.28 2.45

>

〜〜〜

JavaScript-Starで分析してみると

ABs

Sippai

3

Kanzen

2

4 4 4 4 4 4

10 13 21 16

16 19 13 8

15 16 12 15

21 23 16 19

21 14 24 20

31 36 24 34


以上、入力画面

以下、結果。

[ ABs-Type Design ]


== Means & SDs ( of samples') ==


A= Sippai

B= Kanzen

-----------------------------------------------------

A B N Mean S.D.

-----------------------------------------------------

1 1 4 15.0000 4.0620

1 2 4 14.0000 4.0620


2 1 4 14.5000 1.5000

2 2 4 19.7500 2.5860


3 1 4 19.7500 3.6315

3 2 4 31.2500 4.5484


-----------------------------------------------------


== Analysis of Variance ==


A(3) = Sippai

B(2) = Kanzen

-----------------------------------------------------

S.V SS df MS F

-----------------------------------------------------

A 528.0833 2 264.0417 15.67 **

(A at B1 67.1667 2 33.5833 1.99 ns)

(A at B2 617.1667 2 308.5833 18.32 **)

-----------------------------------------------------

B 165.3750 1 165.3750 9.82 **

(B at A1 2.0000 1 2.0000 0.11 ns)

(B at A2 55.1250 1 55.1250 3.27 +)

(B at A3 264.5000 1 264.5000 15.70 **)

-----------------------------------------------------

AxB 156.2500 2 78.1250 4.64 *

Sub 303.2500 18 16.8472

-----------------------------------------------------

Total 1152.9583 23 +p<.10 *p<.05 **p<.01




== Multiple Comparisons by LSD ==


(MSe= 16.8472, * p<.05)

(LSD= 6.0975)

-----------------------------------------------------

A at B1 Level

-----------------------------------------------------

A at B2 Level

-----------------------------------------------------

A1 = A2 n.s.

A1 < A3 *

A2 < A3 *

-----------------------------------------------------

-----------------------------------------------------


_/_/_/ Analyzed by JavaScript-STAR _/_/_/


〜〜〜


多重比較

事後の多重比較は、Practical Statistics with R

http://www.brandonu.ca/zoology/rutherford/RBookWebSite/RBook_pdf.html

のなかの、Multiple Comparisons を参照しながらおこなってみます。


心理学では、単純主効果の分析が多いようですが、簡便で分かりやすいと思われる、すべての群間での比較をしてみます。


まずはTukeyHSD() 分散分析の結果を( )の中に入れます。


> TukeyHSD(y)

Tukey multiple comparisons of means

95% family-wise confidence level


Fit: aov(formula = 抑うつ ~ 失敗経験 * 完全主義)


$失敗経験

diff lwr upr p adj

2-1 2.625 -2.612724 7.862724 0.4245182

3-1 11.000 5.762276 16.237724 0.0001214

3-2 8.375 3.137276 13.612724 0.0019288


$完全主義

diff lwr upr p adj

2-1 5.25 1.729548 8.770452 0.0057477


$"失敗経験:完全主義"

diff lwr upr p adj

2:1-1:1 -0.50 -9.723756 8.723756 0.9999747

3:1-1:1 4.75 -4.473756 13.973756 0.5867657

1:2-1:1 -1.00 -10.223756 8.223756 0.9992404

2:2-1:1 4.75 -4.473756 13.973756 0.5867657

3:2-1:1 16.25 7.026244 25.473756 0.0003203

3:1-2:1 5.25 -3.973756 14.473756 0.4847027

1:2-2:1 -0.50 -9.723756 8.723756 0.9999747

2:2-2:1 5.25 -3.973756 14.473756 0.4847027

3:2-2:1 16.75 7.526244 25.973756 0.0002252

1:2-3:1 -5.75 -14.973756 3.473756 0.3895912

2:2-3:1 0.00 -9.223756 9.223756 1.0000000

3:2-3:1 11.50 2.276244 20.723756 0.0099886

2:2-1:2 5.75 -3.473756 14.973756 0.3895912

3:2-1:2 17.25 8.026244 26.473756 0.0001588

3:2-2:2 11.50 2.276244 20.723756 0.0099886


>


失敗経験が多くて、完全主義が強い人(3:2)の人は、失敗経験が少なくても、多くても、完全主義が弱い人(1:1, 2:1, 3:1)の人たちよりも抑うつ的になっており、また完全主義が強いけれども、失敗経験が少ない人(1:2, 2:2)よりも抑うつ的だという結果ですね。


Holmを使った多重比較

さて、別な方法で群間の比較をしてみます。

最初にすべきことは、失敗経験と完全主義の組み合わせのグループをつくることです。ifelse( ) を使って群分けをしましょう。( )の最初は、条件式で、もし真ならば、第2引数にして、そうでないならば、第3引数にするということです。忘れないように、factor形式に直しておきます。それから、横に並べてくっつけましょう。



group <- ifelse(失敗経験=="1" & 完全主義=="1", 1,

ifelse(失敗経験=="2" & 完全主義=="1", 2,

ifelse(失敗経験=="3" & 完全主義=="1", 3,

ifelse(失敗経験=="1" & 完全主義=="2", 4,

ifelse(失敗経験=="2" & 完全主義=="2", 5,

ifelse(失敗経験=="3" & 完全主義=="2", 6, NA))))))

group <- factor(group)

x.df <- cbind(x.df, group)


以下のようなデータフレームになりますね。


> x.df

番号 失敗経験 完全主義 抑うつ group

1 1 1 1 10 1

2 2 1 1 13 1

3 3 1 1 21 1

4 4 1 1 16 1

5 5 1 2 16 4

6 6 1 2 19 4

7 7 1 2 13 4

8 8 1 2 8 4

9 9 2 1 15 2

10 10 2 1 16 2

11 11 2 1 12 2

12 12 2 1 15 2

13 13 2 2 21 5

14 14 2 2 23 5

15 15 2 2 16 5

16 16 2 2 19 5

17 17 3 1 21 3

18 18 3 1 14 3

19 19 3 1 24 3

20 20 3 1 20 3

21 21 3 2 31 6

22 22 3 2 36 6

23 23 3 2 24 6

24 24 3 2 34 6

>


一要因の時と同じように、pairwise.t.testを使います。

noneを指定すると、FisherLSDでの検定。有意な差がでやすいが、タイプIのエラーが問題となります。


> pairwise.t.test(抑うつ, group, p.adjust.method="none")


Pairwise comparisons using t tests with pooled SD


data: 抑うつ and group


1 2 3 4 5

2 0.86514 - - - -

3 0.11908 0.08720 - - -

4 0.73443 0.86514 0.06306 - -

5 0.11908 0.08720 1.00000 0.06306 -

6 2.6e-05 1.8e-05 0.00091 1.3e-05 0.00091


P value adjustment method: none


Bonferroniの補正で計算すると、以下のようになります。

補正については、以下の説明が分かりやすいです。

http://chiryo.phar.nagoya-cu.ac.jp/javastat/BONFERRONI-P.PDF


> pairwise.t.test(抑うつ, group, p.adjust.method="bonferroni")


Pairwise comparisons using t tests with pooled SD


data: 抑うつ and group


1 2 3 4 5

2 1.00000 - - - -

3 1.00000 1.00000 - - -

4 1.00000 1.00000 0.94592 - -

5 1.00000 1.00000 1.00000 0.94592 -

6 0.00039 0.00027 0.01370 0.00019 0.01370


P value adjustment method: bonferroni



Bonferroniは、逆に検出力が弱くなるので、発展させたHolmで分析してみます。

pの値が、LSD>Holm>Bonferroniの順になってることが分かりますね。


> pairwise.t.test(抑うつ, group, p.adjust.method="holm")


Pairwise comparisons using t tests with pooled SD


data: 抑うつ and group


1 2 3 4 5

2 1.00000 - - - -

3 0.71446 0.69763 - - -

4 1.00000 1.00000 0.63061 - -

5 0.71446 0.69763 1.00000 0.63061 -

6 0.00034 0.00025 0.01096 0.00019 0.01096


P value adjustment method: holm

>


あるいは、一元配置の分散分析の結果をTukeyHSDに入れて、下位検定をおこなうことも可能でしょう。


> TukeyHSD(aov(抑うつ~group))

Tukey multiple comparisons of means

95% family-wise confidence level


Fit: aov(formula = 抑うつ ~ group)


$group

diff lwr upr p adj

2-1 -0.50 -9.723756 8.723756 0.9999747

3-1 4.75 -4.473756 13.973756 0.5867657

4-1 -1.00 -10.223756 8.223756 0.9992404

5-1 4.75 -4.473756 13.973756 0.5867657

6-1 16.25 7.026244 25.473756 0.0003203

3-2 5.25 -3.973756 14.473756 0.4847027

4-2 -0.50 -9.723756 8.723756 0.9999747

5-2 5.25 -3.973756 14.473756 0.4847027

6-2 16.75 7.526244 25.973756 0.0002252

4-3 -5.75 -14.973756 3.473756 0.3895912

5-3 0.00 -9.223756 9.223756 1.0000000

6-3 11.50 2.276244 20.723756 0.0099886

5-4 5.75 -3.473756 14.973756 0.3895912

6-4 17.25 8.026244 26.473756 0.0001588

6-5 11.50 2.276244 20.723756 0.0099886


>

「交互作用の分析(単純主効果の検定)


ということなので、x.dfを分割しながら検定を行ってみましょう。


以下の最初のやり方は、小塩のSPSSやり方と違い、単純に輪切りをしておこなっています。

(このやり方ですると、誤差項は全体ではありません。「水準別誤差項」を用いた分析です。)



完全主義が低群”1”のケースのみと取り出すには、以下のとおりです。


> x.df[x.df$完全主義=="1",]

番号 失敗経験 完全主義 抑うつ

1 1 1 1 10

2 2 1 1 13

3 3 1 1 21

4 4 1 1 16

9 9 2 1 15

10 10 2 1 16

11 11 2 1 12

12 12 2 1 15

17 17 3 1 21

18 18 3 1 14

19 19 3 1 24

20 20 3 1 20

>


そこで、


> x1.df <- x.df[x.df$完全主義=="1",]

> summary(aov(x1.df$抑うつ ~ factor(x1.df$失敗経験)))

Df Sum Sq Mean Sq F value Pr(>F)

factor(x1.df$失敗経験) 2 67.167 33.583 2.3659 0.1494

Residuals 9 127.750 14.194


完全主義が小さい人では、失敗経験は抑うつとは関係がなさそうですね。(誤差の自由度が小塩SPSS結果と違うことに注意)


次に完全主義が強い人のなかで検定してみましょう。

> x2.df <- x.df[x.df$完全主義=="2",]

> summary(aov(x2.df$抑うつ ~ factor(x2.df$失敗経験)))

Df Sum Sq Mean Sq F value Pr(>F)

factor(x2.df$失敗経験) 2 617.17 308.58 15.825 0.001131 **

Residuals 9 175.50 19.50

---

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


有意だったので、多重比較をしてみます。


> pairwise.t.test(x2.df$抑うつ, factor(x2.df$失敗経験), p.adjust.method="bonferroni" )


Pairwise comparisons using t tests with pooled SD


data: x2.df$抑うつ and factor(x2.df$失敗経験)


1 2

2 0.2960 -

3 0.0011 0.0152


P value adjustment method: bonferroni


あるいは・・・

> pairwise.t.test(x2.df$抑うつ, factor(x2.df$失敗経験), p.adjust.method="holm" )


Pairwise comparisons using t tests with pooled SD


data: x2.df$抑うつ and factor(x2.df$失敗経験)


1 2

2 0.0987 -

3 0.0011 0.0101


P value adjustment method: holm

>


完全主義が強い人においては、失敗経験が多い人は、失敗経験の少ない人や中程度の人と比べると、抑うつが高いようです。


以下は2群の検定になりますので、t検定してみます。

> t.test(x3.df$抑うつ ~ factor(x3.df$完全主義), var.equal=T)


Two Sample t-test


data: x3.df$抑うつ by factor(x3.df$完全主義)

t = 0.3015, df = 6, p-value = 0.7732

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-7.115489 9.115489

sample estimates:

mean in group 1 mean in group 2

15 14



有意な差はないようです。





> x4.df <- x.df[x.df$失敗経験=="2",]

> t.test(x4.df$抑うつ ~ factor(x4.df$完全主義), var.equal=T)


Two Sample t-test


data: x4.df$抑うつ by factor(x4.df$完全主義)

t = -3.0417, df = 6, p-value = 0.02275

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-9.473434 -1.026566

sample estimates:

mean in group 1 mean in group 2

14.50 19.75


有意差が認められました。


> x5.df <- x.df[x.df$失敗経験=="3",]

> t.test(x5.df$抑うつ ~ factor(x5.df$完全主義), var.equal=T)


Two Sample t-test


data: x5.df$抑うつ by factor(x5.df$完全主義)

t = -3.4223, df = 6, p-value = 0.01410

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-19.722376 -3.277624

sample estimates:

mean in group 1 mean in group 2

19.75 31.25



これも有意差が認められました。


〜〜〜

小塩のSPSSでのやり方(「プールされた誤差項 」を用いる方法は以下のようになります。まずSPSSでされた抑うつの単純主効果の検定結果は以下のようです。違いは誤差項の平方和と自由度、平均平方は、最初に分散分析をおこなったもの(全体の誤差項)を使っています。下は、その数値をコピーしたものです。


Df Sum Sq Mean Sq

Residuals 18 303.25 16.85


 完全主義が低群と高群の自由度2についてのS.S. と M.S.は、最初に群わけをした上でおこなった分散分析に利用したものを使います。それぞれの結果を引用してくると以下のとおりです。




小塩のページからSPSSの結果 


     Df Sum Sq Mean Sq

factor(x1.df$失敗経験) 2 67.167 33.583

actor(x2.df$失敗経験) 2 617.17 308.58


F値の求め方は、それぞれの平均平方を、誤差の平均平方で割った値ですね。

低群に場合

> 33.583/16.85

[1] 1.993056


高群の場合

> 308.58/16.85

[1] 18.31335


F値が求められました。有意確率は、pf(F値, 自由度1, 自由度2)関数で求められます。


> 1-pf((33.583/16.85), 2, 18)

[1] 0.1652405


> 1-pf(( 308.58/16.85), 2, 18)

[1] 4.579388e-05


同じように、失敗経験の水準ごとの単純主効果も求められます。やってみてください。


〜〜〜


2要因の分散分析(混合計画)



次は、2要因の分散分析(混合計画)

http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/05_folder/da05_03.html


従属変数は、血圧で、独立変数は,偽薬(0)・Aを投与(1)・Bを投与(2)の3群(被験者間要因)と投与前後(被験者内要因)という例です。





これを、以下のような形式に変えます。





> SPFp.q <- read.delim("clipboard")

> attach(SPFp.q)

> Subj <- factor(Subj)

> Drug <- factor(Drug)

> Time <- factor(Time)

> print(summary(aov(Blood ~ Drug * Time + Error(Subj:Drug + Subj:Drug:Time), SPFp.q)))


Error: Subj:Drug

Df Sum Sq Mean Sq F value Pr(>F)

Drug 2 312.47 156.23 0.651 0.539

Residuals 12 2880.00 240.00


Error: Subj:Drug:Time

Df Sum Sq Mean Sq F value Pr(>F)

Time 1 202.8 202.8 8.0476 0.014985 *

Drug:Time 2 351.8 175.9 6.9802 0.009755 **

Residuals 12 302.4 25.2

---

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

Warning message:

エラーモデルは特異です in: aov(Blood ~ Drug * Time + Error(Subj:Drug + Subj:Drug:Time),

>


下の小塩と同じ結果が得られています。

投与と薬品の交互作用が5%水準で有意:(2,12)=6.98, p<.05

投与の主効果は5%水準で有意:(1,12)=8.05, p<.05

薬品の主効果はみられない:(2,12)=.651, n.s.


あるいは、aovのなかでfactor指定も可。

summary(aov(Blood ~ factor(Drug) * factor(Time) + Error(Subj:factor(Drug) + Subj:factor(Drug):factor(Time)), SPFp.q))


グラフは、attach していたら以下。

interaction.plot(Time, Drug, Blood)


group <- factor(ifelse(Time==1 & Drug==1,1,

ifelse(Time==1 & Drug==2,2,

ifelse(Time==1 & Drug==3,3,

ifelse(Time==2 & Drug==1,4,

ifelse(Time==2 & Drug==2,5,

ifelse(Time==2 & Drug==3,6,NA)))))))


上のように、グループ分けをして、LSD法で6つの群間の差を見てみましょう。(BonferroniHolmでは有意差はでませんでした。)


> pairwise.t.test(Blood, group, p.adjust.method="none")


Pairwise comparisons using t tests with pooled SD


data: Blood and group


1 2 3 4 5

2 0.913 - - - -

3 0.935 0.978 - - -

4 0.745 0.828 0.807 - -

5 0.684 0.607 0.626 0.466 -

6 0.074 0.060 0.063 0.038 0.158


P value adjustment method: none

>


どうにか薬品Bがプラセボよりも効果があった感じですか。


小塩のページでは、単純主効果を調べることを推奨しています。可能な人はトライしてみてください。


〜〜


分散分析の下位検定 http://www5e.biglobe.ne.jp/~tbs-i/psy/semi/anova/ のページには、交互作用の下位検定が4種類紹介されています。


http://home.hiroshima-u.ac.jp/nittono/ANOVA_JSPP2004.pdf の「心理生理学データの分散分析 」にも分かりやすい説明があります。



〜〜〜付記:重要な事柄〜〜


http://mat.isc.chubu.ac.jp/R/tech.html#RB 言語Rによる分散分析のページには、「心理学のためのデータ解析テクニカルブック」例題が書かれおり、とても参考になります。


http://www.okada.jp.org/RWiki/index.php?cmd=read&page=R%A4%CE%C5%FD%B7%D7%B2%F2%C0%CF%B4%D8%BF%F4Tips&word=%CA%BF%CA%FD%CF%C2#content_1_1 に、aov, lm, glm における「線形モデル当てはめ関数におけるモデル公式」が有益な情報です。



 ところで、平方和の計算方法にはいくつか種類があります。

 説明は、http://androids.happy.nu/ja/anova.html が参考になりますが、セル(ブロック)内の数が、アンバランスの場合は、type=II type=III が奨められてます。http://www.aichi-gakuin.ac.jp/~chino/anova/chapter1/sec1-3-7.html では、「非釣り合い型デザインでも、交互作用項を含まない主効果 のみのデザインであれば、Type II 平方和とType III 平方和は一致するが、 交互作用項を含んだデザインでは、両者は異なる(竹内監修、高橋ら、1990)。 このような場合、高橋らはType II 平方和を用いることを薦めている。 」とあります。

 また http://phi.med.gunma-u.ac.jp/swtips/Rnewsarchives.html#SAS-SS ここには、「Type IIIは繰り返し数が不揃いのときにデータ数の少ないセルを他のセルと同等とみなす目的で使うもの。 」「高橋・大橋・芳賀「SASによる実験データの解析」(東大出版会)によると,数量化一類をするときや,乱塊法の場合や,MANOVAの場合や欠損値がある 場合はType IIの使用が薦められるとあるので,とりあえずType IType IIだけ出せれば充分ではないかと思う。 」「結論としては,Rで,因子が直交していなくてセルごとの繰り返し数が不揃いの二元配置分散分析をしたいときは,library(car)としてから, Anova(lm(Y~C1+C2+C1:C2))を使えば偏平方和が計算されるので,そうすればいいということのようである。 」とあります。

 一方で、http://mat.isc.chubu.ac.jp/fpr/fpr1999/0051.html には、「要は「Type II の平方和では、セルサイズによる重みが影響してくるのに対し、Type III の平方和ではセルサイズによる重みをつけない効果を推定してそれを捨象している」というのが2つの平方和の違いなのです。セルサイズによる重みを与えた方がいいのか、無視した方がいいのか、議論は決着していません。私は無視した(つまり Type III )方がいいと思っているのですが、Type II派の人は「セルサイズが極端に違う場合、たとえばあるセルのサイズがゼロであっても等しい重みをかけるというのか。そんな方法に統計的妥当性は認めない!」と主張しています。」ともあります。


  備忘録:TypeIIIの分散分析をRでするには、library(car)のあとに、Anova(lm(Y~C1+C2+C1:C2), type=”III”) になります。

 ちなみにSPSSdefaulttypeIIIだそうです。


 私の考えでは、http://www4.ocn.ne.jp/~murakou/anova.htm にかかれてある「いくつかの種類の平方和で同じ結果が得られたなら,それはセルの人数比が不揃いであっても,自信を持ってその結果を報告してもよいということです. 」というところが納得できるところだと思います。


第4部:Rでt検定


http://kogolab.jp/elearn/hamburger/chap4/sec1.html


にある例題をRでみていきます。

最初は、対応の無い二標本 t 検定の例です。


ワクワクバーガーを食べた女子高生8名の点数をaに、モグモグバーガーを食べた女子高生8名(最初の8名とは別人)の点数をbにいれます。


> a <- c(70,75,70,85,90,70,80,75)

> b <- c(85,80,95,70,80,75,80,90)


まずは分布をみてみましょう。


> stem(a)


The decimal point is 1 digit(s) to the right of the |


7 | 000

7 | 55

8 | 0

8 | 5

9 | 0


> stem(b)


The decimal point is 1 digit(s) to the right of the |


7 | 05

8 | 0005

9 | 05


ヒストグラムを見ましょう。


> hist(a)

> hist(b)


箱ひげ図は。


> boxplot(a)

> boxplot(b)


ひとつにまとめて表示しましょう。


> boxplot(a, b, names=c("ワクワク","モグモグ"))



それぞれの平均は・・・


> mean(a)

[1] 76.875

> mean(b)

[1] 81.875


小数点2桁にして、見やすくしましょう。


> round(mean(a),2)

[1] 76.88

> round(mean(b),2)

[1] 81.88


useful.rをいれていたらdescribeを使って基礎統計量を確認できます。

> describe(a)

V n mean sd median min max range skew se

1 1 8 76.88 7.53 75 70 90 20 0.54 2.66

> describe(b)

V n mean sd median min max range skew se

1 1 8 81.88 7.99 80 70 95 25 0.2 2.82

>


t検定をする前に、それぞれの分散が等しいがどうかを調べます。

帰無仮説は、等分散である。(分散に差がない)

対立仮説は、分散に差がある。(等分散ではない)


> var.test(a,b)


F test to compare two variances


data: a and b

F = 0.8881, num df = 7, denom df = 7, p-value = 0.8796

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

0.1778034 4.4360383

sample estimates:

ratio of variances

0.8881119


等しいという帰無仮説は棄却されませんでしたので、等分散のようですね。


t検定は、t.test( )なのですが、注意をしないと以下のようになります。


> t.test(a,b)


Welch Two Sample t-test


data: a and b

t = -1.2881, df = 13.951, p-value = 0.2187

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-13.327988 3.327988

sample estimates:

mean of x mean of y

76.875 81.875


 上のようにオプションをつけないと、等分散でない時の検定がおこなわれます。(Welch(ウェルチ)の検定)


 先に調べたように、等分散であるならば、オプション var.equal=T を引数(ひきすう)にくわえ、t検定をします。

 通常、心理学研究の場合、等分散であるという仮定をしていく場合が多い気がします。


帰無仮説:味に差がない。

対立仮説:味に差がある。


> t.test(a,b, var.equal=T)


Two Sample t-test


data: a and b

t = -1.2881, df = 14, p-value = 0.2186

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-13.325244 3.325244

sample estimates:

mean of x mean of y

76.875 81.875



味に差がないという帰無仮説は棄却されません。


次の方法もぜひ覚えておいてください。普通に質問紙調査をすると、dataframe形式で扱います。



> group <- c(rep("A",8),rep("B",8))


repは、repeatですね。 最初の引数を、次の引数の数だけ繰り返します。

Aが、ワクワクバーガー、Bがモグモグバーガーということです。

> group

[1] "A" "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B" "B"


> aji <- c(a,b)

最初に、ワクワクバーガーの得点を並べて、次にモグモグバーガーの得点を並べます。


> aji

[1] 70 75 70 85 90 70 80 75 85 80 95 70 80 75 80 90


> buger.df <- data.frame(group,aji)

> buger.df

group aji

1 A 70

2 A 75

3 A 70

4 A 85

5 A 90

6 A 70

7 A 80

8 A 75

9 B 85

10 B 80

11 B 95

12 B 70

13 B 80

14 B 75

15 B 80

16 B 90


上のようなデータフレームになります。

~は、チルダです。


> t.test(aji ~ group, data = buger.df, var.equal=T)


Two Sample t-test


data: aji by group

t = -1.2881, df = 14, p-value = 0.2186

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-13.325244 3.325244

sample estimates:

mean in group A mean in group B

76.875 81.875


この形式の場合の箱ひげ図は以下のようにします。

> boxplot(aji~group, data=buger.df)



> library(gplots)

要求されたパッケージ gtools をロード中です

要求されたパッケージ gdata をロード中です


次のパッケージを付け加えます: 'gplots'



The following object(s) are masked from package:stats :


lowess


> plotmeans(aji~group, data=buger.df)


useful.rを取り入れていて、describeを使えるならば・・・


> by(buger.df$aji, buger.df$group, describe)

INDICES: A

V n mean sd median min max range skew se

1 1 8 76.88 7.53 75 70 90 20 0.54 2.66

---------------------------------------------------------------------------------------------------

INDICES: B

V n mean sd median min max range skew se

1 1 8 81.88 7.99 80 70 95 25 0.2 2.82

>



対応のあるt検定

次は、

http://kogolab.jp/elearn/hamburger/chap5/sec1.html


対応のあるt検定です。

「女子高生8人を集めて、2種類のハンバーガーを食べて比較してもらい、それぞれのハンバーガーに点数をつけてもらいました。」という例題です。


> a <- c(90,75,75,75,80,65,75,80)

> b <- c(95,80,80,80,75,75,80,85)

> describe(a)

V n mean sd median min max range skew se

1 1 8 76.88 7.04 75 65 90 25 0.22 2.49

> describe(b)

V n mean sd median min max range skew se

1 1 8 81.25 6.41 80 75 95 20 1.02 2.27

>


http://cse.naro.affrc.go.jp/takezawa/r-tips/r/64.html を参考にします。


オプションに、paired = T をつけます。


> t.test(a, b, paired = TRUE)


Paired t-test


data: a and b

t = -2.9656, df = 7, p-value = 0.02094

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-7.8633933 -0.8866067

sample estimates:

mean of the differences

-4.375


>

ワクワクとモグモグの評価点の平均には差があるという結果になりました。


もし、どちらかがおいしいという仮説が前もってあるならば、オプションをつけることで、片側検定ができます。


> t.test(a, b, paired = TRUE, alternative = "less")


Paired t-test


data: a and b

t = -2.9656, df = 7, p-value = 0.01047

alternative hypothesis: true difference in means is less than 0

95 percent confidence interval:

-Inf -1.580038

sample estimates:

mean of the differences

-4.375


あるいは以下のとおり。

> t.test(a - b, alternative = "less")


One Sample t-test


data: a - b

t = -2.9656, df = 7, p-value = 0.01047

alternative hypothesis: true mean is less than 0

95 percent confidence interval:

-Inf -1.580038

sample estimates:

mean of x

-4.375


データフレームにしてやってみましょう。


まず、協力者を8名番号をつけます。

>subject <- c(1:8)

>pair.t.df <- data.frame(subject, a, b)


> subject <- c(1:8)

> pair.t.df <- data.frame(subject, a, b)

> pair.t.df

subject a b

1 1 90 95

2 2 75 80

3 3 75 80

4 4 75 80

5 5 80 75

6 6 65 75

7 7 75 80

8 8 80 85

>

これで、データフレームができました。


> t.test(a, b, paired=T, data=pair.t.df)


Paired t-test


data: a and b

t = -2.9656, df = 7, p-value = 0.02094

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-7.8633933 -0.8866067

sample estimates:

mean of the differences

-4.375



http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/03_folder/da03_03.html#%91%CE%89%9E%82%CC%82%A0%82%E9%82%94%8C%9F%92%E8

にある、10名に対してある授業を行った前と後にテストを行い,成績を算出した。授業前後で成績が伸びているといえるかという例です。


> before <- c(5,4,5,4,3,4,3,6,5,6)

> after <- c(7,6,8,7,6,6,5,9,9,8)

> t.test(before, after, paired=T)


Paired t-test


data: before and after

t = -11.7589, df = 9, p-value = 9.151e-07

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-3.100182 -2.099818

sample estimates:

mean of the differences

-2.6


>

参考リンク

以下のページも参考になりそうですので、お勧めしておきます。


http://www.sky.sannet.ne.jp/mikeneko/psychology/main.htm

●心理統計の基礎

├心理学と統計

├平均値・標準偏差

├偏差値

├相関関数

├統計的検定とは?

├統計的検定の考え方

t検定とは?

t検定の手順

├分散分析と多重比較

├独立変数と従属変数

├主効果と交互作用

└分散分析の手順



1