################################################ #南風原(2002)第9章 実験デザインと分散分析 ################################################ #Step1:デザインの特徴 p.270-72 ################################################ ##(a) 妬み感情の測定データの読み込み setwd("C:/data/haebara/") dat1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap9.csv", header=T) dat1 ##(b) 名義変数を適切な順序に並べ替える dat1$優越属性 #表示される,Levelsは「学歴 豊かさ 容姿」の順になっているが #「容姿 学歴 豊かさ」の順に変更したい #------------------------------------------------------# #factor(名義変数, levels=c(名義変数の順番を指定する)) #------------------------------------------------------# factor(dat1$優越属性,levels=c("容姿","学歴","豊かさ")) #「容姿 学歴 豊かさ」の順になっていることを確認 dat1$優越属性 <- factor(dat1$優越属性,levels=c("容姿","学歴","豊かさ")) #保存する dat1$相手の態度 #表示される,Levelsは「ふつう 敵対的 友好的」の順になっているが #「友好的 敵対的 ふつう」の順に変更したい factor(dat1$相手の態度,levels=c("友好的","敵対的","ふつう")) #「友好的 敵対的 ふつう」の順になっていることを確認 dat1$相手の態度 <- factor(dat1$相手の態度,levels=c("友好的","敵対的","ふつう")) #保存する ##(c)表9-2 優越属性の各条件における妬み感情の測定値の平均と標準偏差 #------------------------------------------------------# # 標準偏差関数の定義(ssd)と命名 ssd <- function(x){ sqrt(1/length(x)*sum((x-mean(x))^2)) } #------------------------------------------------------# attach(dat1) #attach()で使用するデータの変数名で指定できるようにする #------------------------------------------------------# # by(数値変数, 質的変数, 関数) # tapply(数値変数, 質的変数, 関数) #で各条件ごとの平均・標準偏差を算出できる #------------------------------------------------------# by(妬み感情,優越属性,mean) #優越属性ごとの平均 tapply(妬み感情,優越属性,mean) #優越属性ごとの平均 by(妬み感情,優越属性,ssd) #優越属性ごとの標準偏差 tapply(妬み感情,優越属性,ssd) #優越属性ごとの標準偏差 mean(妬み感情) #全体の平均 ssd(妬み感情) #全体の標準偏差 #------------------------------------------------------# # boxplot(従属変数~独立変数)で図を描く #------------------------------------------------------# boxplot(妬み感情~優越属性) #本来は,箱型図で条件ごとの要約統計量を確認すると良い graphics.off() #図を閉じる detach(dat1) #attach()したものをdetach()する ################################################ #Step2:平方和の分解 p.272-74 ################################################ ##(a) 全体の平方和の式を確認する(9.1式) attach(dat1) #attach()で使用するデータの変数名で指定できるようにする a <- 優越属性 #変数名が長いので,見やすくする為に「a」と書き換える y <- 妬み感情 #変数名が長いので,見やすくする為に「y」と書き換える detach(dat1) #attach()したものをdetach()する #------------------------------------------------------# # 分散関数の定義(svarと命名) svar <- function(x){ 1/length(x)*sum((x-mean(x))^2) } #------------------------------------------------------# mean(y) #全体の平均 svar(y) #全体の分散 svar(y) * length(y) #全体の平方和 sum((y-mean(y))^2) #全体の平方和(9.1式使用) SSt <- sum((y-mean(y))^2) #SStという名前で,全体の平方和の値を保存 ##(b) 群間平方和の式を確認する(9.2式) tapply(y,a,mean) #各条件の平均 tapply(y,a,mean)-mean(y) #各条件の平均-全体の平均 (tapply(y,a,mean)-mean(y))^2 #(各条件の平均-全体の平均)^2 tapply(y,a,length) #各条件の標本サイズ tapply(y,a,length)* (tapply(y,a,mean)-mean(y))^2 #各条件の標本サイズ * (各条件の平均-全体の平均)^2 sum(tapply(y,a,length)* (tapply(y,a,mean)-mean(y))^2) #Σ(各条件の標本サイズ * (各条件の平均-全体の平均)^2) SSa <- sum(tapply(y,a,length)* (tapply(y,a,mean)-mean(y))^2) #SSaという名前で,群間平方和の値を保存 ##(c) 群内平方和の式を確認する(9.3式) tapply(y,a,svar) #各条件の分散 tapply(y,a,length) #各条件の標本サイズ tapply(y,a,length) * tapply(y,a,svar) #各条件の標本サイズ*各条件の分散 sum(tapply(y,a,length) * tapply(y,a,svar)) #Σ(各条件の標本サイズ*各条件の分散) SSe <- sum(tapply(y,a,length) * tapply(y,a,svar)) #SSeという名前で,群内平方和の値を保存 ################################################ #Step3:相関比/自由度/仮定と検定統計量/分散分析表 p.274-79 ################################################ ##(a) 相関比を確認(9.4式) SSa/SSt #分散説明率 sqrt(SSa/SSt) #相関比η ##(b) 自由度 (9.5/9.6 式) nlevels(a) #水準数 nlevels(a) - 1 #要因の自由度 dfa <- nlevels(a) - 1 #dfa という名前で,要因の自由度の値を保存 length(y) #全体の標本サイズ length(y) - nlevels(a) #残差の自由度 dfe <- length(y) - nlevels(a) #dfeという名前で,残差の自由度の値を保存 length(y) - 1 #全体の自由度 dfa + dfe #2つの自由度の和は,全体の自由度になる ##(c) F値を確認(9.9式) ( SSa/dfa ) / ( SSe/dfe ) F.value <- ( SSa/dfa ) / ( SSe/dfe ) #F.valueという名前で,F値の値を保存 ##(d) 付表9aを確認 p.384 #------------------------------------------------------# # qf(p=有意水準, df1= 分子の自由度, df2= 分母の自由度, lower.tail=F(上側確率)) #------------------------------------------------------# qf(p=.05,df1=2,df2=40,lower.tail=F) qf(p=.05,df1=2,df2=42,lower.tail=F) qf(p=.05,df1=2,df2=50,lower.tail=F) ##(e) 群内平均平方 (9.10式)を確認 SSe/dfe MSe <- SSe/dfe #MSeという名前で,群内平均平方の値を保存 sqrt(MSe) #各水準共通と仮定される母標準偏差の推定値 ##(f) 分散分析表 #------------------------------------------------------# #1. aov(従属変数~独立変数,data=データセット) #2. summary(aov()) # で分散分析を実行する #------------------------------------------------------# aov(妬み感情~優越属性,data=dat1) table9.3 <- aov(妬み感情~優越属性,data=dat1) #table9.3という名前で,分散分析の結果を保存する summary(table9.3) #summary(保存したオブジェクト)で分散分析表を出す ################################################ #Step4:多重比較の考え方 p.279-84 ################################################ ##(a) 統計量q を算出する(9.11式) tapply(y,a,mean) #各群の標本平均 max(tapply(y,a,mean)) #各群の標本平均のうちの最大のもの min(tapply(y,a,mean)) #各群の標本平均のうちの最小のもの sqrt(MSe/15) #標本平均の標準誤差の推定量(9.12式) (max(tapply(y,a,mean))-min(tapply(y,a,mean)))/sqrt(MSe/15) #9.11式 ##(b) 付表10を確認 p.386 #------------------------------------------------------# # qtukey(p=有意水準, nmeans=比較する平均の総数, df= 残差の自由度, lower.tail=F(上側確率)) #------------------------------------------------------# qtukey(p=.05,nmeans=3, df=40,lower.tail=F) qtukey(p=.05,nmeans=3, df=42,lower.tail=F) qtukey(p=.05,nmeans=3, df=60,lower.tail=F) qtukey(p=.05,nmeans=3, df=42,lower.tail=F) * sqrt(MSe/15) #平均値の差の絶対値が,この値を越えれば5%水準で有意と判断できる ##(c) 有意差がある対を確認する tapply(y,a,mean) #各条件の平均 tapply(y,a,mean)[[1]] #1番目の条件の平均だけを取り出す a1 <- tapply(y,a,mean)[[1]] #各条件の平均をa1-a3と名前をつける a2 <- tapply(y,a,mean)[[2]] a3 <- tapply(y,a,mean)[[3]] abs(a1-a2) #1.56を超える abs(a2-a3) #1.56を超えない abs(a3-a1) #1.56を超えない ##(d) TuckeyHSDという多重比較用の関数を使って求める #------------------------------------------------------# #1. aov(従属変数~独立変数,data=データセット) #2. TukeyHSD(aov(), ordered=F) # 「orderd」は,平均値差の大きい水準の順番で並べ替えるかどうかを指定する #3. plot(TukeyHSD(aov(), ordered=F)) # で多重比較を実行する #------------------------------------------------------# TukeyHSD(aov(妬み感情~優越属性,data=dat1), ordered=T) #aov関数を使った場合 TukeyHSD(table9.3,ordered=T) #aov関数の結果をtable9.3という名前で保存して使った場合 TukeyHSD(table9.3,ordered=F) #orderedをFにした場合 plot(TukeyHSD(table9.3,ordered=T)) #多重比較の結果を図示する #95%信頼区間が0を含まなければ, #有意であると判断できる graphics.off() #図を閉じる ################################################ #Step5:完全無作為2要因デザイン p.284-94 ################################################ ##(a) 優越属性と相手の態度の組み合わせごとの妬み感情の測定値の平均(表9-4) dat1 #データを表示 attach(dat1) #attach()で使用するデータの変数名で指定できるようにする tapply(妬み感情,list(優越属性,相手の態度),mean) tapply(妬み感情,list(優越属性,相手の態度),ssd) tapply(妬み感情,list(優越属性),mean) tapply(妬み感情,list(相手の態度),mean) mean(妬み感情) ##(b) セル平均のプロットを描く(図9-1) #------------------------------------------------------# #interaction.plot(x.factor, trace.factor, response, xtick = F) # x.factor =x軸に持ってくる名義変数 # trace.factorは =トレースする名義変数 # response =従属変数 # xtick =x軸にチェックマークを付けるか否か(デフォルトは付けない) #------------------------------------------------------# interaction.plot(優越属性,相手の態度,妬み感情,xtick=T,ylab="妬み感情",ylim=c(0,7)) graphics.off() #図を閉じる #(c) 平方和の分割 p.288 y <- 妬み感情 #変数名が長いので,みやすくするために「y」と書き換える A <- 優越属性 #変数名が長いので,みやすくするために「A」と書き換える B <- 相手の態度 #変数名が長いので,みやすくするために「B」と書き換える detach(dat1) #attach()したものをdetach()する SSt <- sum((y-mean(y))^2) #全体の平方和(9.15式) a <- nlevels(A) #要因Aの水準数 b <- nlevels(B) #要因Bの水準数 n <- length(y)/(a*b) #各条件の標本サイズ SSa <- n*b*sum((tapply(y,A,mean) - mean(y))^2) #要因Aの平方和(9.16式) SSb <- n*a*sum((tapply(y,B,mean) - mean(y))^2) #要因Bの平方和(9.17式) SSe <- n*sum(tapply(y,list(A,B),ssd)^2) #残差の平方和(9.20式) ##(d) 交互作用の平方和を求める(9.19式) ##j=1, k=1の場合 tapply(y,list(A,B),mean)[1,1] tapply(y,A,mean)[[1]] #j=1 tapply(y,B,mean)[[1]] #k=1 mean(y) ##j=1, k=2の場合 tapply(y,list(A,B),mean)[1,2] tapply(y,A,mean)[[1]] #j=1 tapply(y,B,mean)[[2]] #j=3 mean(y) ##j=3,k=3の場合 tapply(y,list(A,B),mean)[3,3] tapply(y,A,mean)[[3]] #j=3 tapply(y,B,mean)[[3]] #k=3 mean(y) #これを,(j,k)=(1,1),(1,2),(1,3),(2,1),(2,2),(2,3),(3,1),(3,2),(3,3) #と9通りの場合を計算したい SSab <- matrix(1,a,b) #3 * 3 のサイズの行列を作る for(j in 1:a){ for(k in 1:b){ SSab[j,k] <- (tapply(y,list(A,B),mean)[j,k]-tapply(y,A,mean)[[j]]- tapply(y,B,mean)[[k]]+ mean(y) )^2 } } SSab sum(SSab) n * sum(SSab) #交互作用の平方和 SSab <- n * sum(SSab) #交互作用の平方和をSSabという名前で保存する ##(e) 平方和の分割を確認(9.14式) SSt SSa + SSb + SSab + SSe ##(f) 自由度 df.ab <- (a-1) * (b-1) #交互作用の自由度(9.21式) df.e <- length(y)-a*b #残差の自由度(9.22式) ##(g) F値 Fa <- (SSa/(a-1)) / (SSe/df.e) #要因Aの主効果(9.25) Fb <- (SSb/(b-1)) / (SSe/df.e) #要因Bの主効果(9.26) Fab <- (SSab/df.ab) / (SSe/df.e) #要因AとBの交互作用(9.27) ##(h) 検定 Fab qf(p=0.05,df1=4,df2=36,lower.tail=F) ##(i) もちろん関数を使えば簡単に計算できる #------------------------------------------------------# #1. aov(従属変数~ A + B + A:B,data=データセット) #2. summary(aov()) # で分散分析を実行する # ・A:Bは,AとBの交互作用を意味する # ・A*Bは,A + B + A:Bを意味する #------------------------------------------------------# table9.5 <- aov(妬み感情~優越属性 * 相手の態度,data=dat1) summary(table9.5) #表9-5 ##(j) 1要因デザインとの比較 p.293 summary(table9.3) summary(table9.5) ################################################ #Step6:単純主効果 p.294-96 ################################################ ##(a)単純主効果の検定 subset(dat1,相手の態度=="友好的",select=c(優越属性,妬み感情)) boxplot(妬み感情~優越属性, data=dat1,subset=相手の態度=="友好的") summary(aov(妬み感情~優越属性, data=dat1,subset=相手の態度=="友好的")) ##優越属性のSum Sq(平方和に注目する) ##summary(table9.5)の Mean Sq 残差の平均平方に注目する (3.6/2)/1.29 #9.30式 qf(p=0.05,df1=2,df=36,lower.tail=F) graphics.off() #図を閉じる ################################################### #Step7:対応のある1要因デザイン p.299-306 ################################################### ##(a)データの読み込みと確認(表9-6) dat2 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap9_2.csv",header=T) dat2 #データの確認 subset(dat2,select=c(容姿:平均)) mean(subset(dat2,select=c(容姿:平均))) ##(b) 全体の平方和の計算 (9.39)式 dat2.1 <- subset(dat2,select=c(容姿:豊かさ)) a <- ncol(dat2.1) #水準数 n <- nrow(dat2.1) #ブロックの数 y.mean <- mean(mean(dat2.1)) #全体の平均 ##i=1, j=1の場合 (dat2.1[1,1]-y.mean)^2 ##i=1, j=2の場合 (dat2.1[1,2]-y.mean)^2 ##i=15, j=3の場合 (dat2.1[15,3]-y.mean)^2 SSt <- matrix(1,nrow=n,ncol=a) #15 * 3 のサイズの行列を作成 for(i in 1:n){ for(j in 1:a){ SSt[i,j] <- (dat2.1[i,j]-y.mean)^2 } } sum(SSt) SSt <- sum(SSt) #全体の平方和をSStという名前で保存 ##(c) ブロック×要因Aの平方和の計算 (9.40)式 ##i=1, j=1の場合 dat2.1[1,1] apply(dat2.1[1,],1,mean) mean(dat2.1[,1]) SSba <- matrix(1,nrow=n,ncol=a) #15 * 3 のサイズの行列を作成 for(i in 1:n){ for(j in 1:a){ SSba[i,j] <- (dat2.1[i,j] - apply(dat2.1[i,],1,mean) - mean(dat2.1[,j]) + y.mean)^2 } } SSba <- sum(SSba) #ブロック×要因Aの平方和をSSbaという名前で保存 ##(d) 要因Aの平方和の計算 (9.41)式 mean(dat2.1) mean(dat2.1) - y.mean (mean(dat2.1) - y.mean)^2 sum((mean(dat2.1) - y.mean)^2) SSa <- n * sum((mean(dat2.1) - y.mean)^2) #要因Aの平方和をSSaという名前で保存 ##(e) ブロックの平方和の計算 (9.42)式 apply(dat2.1,1,mean) apply(dat2.1,1,mean) - y.mean (apply(dat2.1,1,mean) - y.mean)^2 sum((apply(dat2.1,1,mean) - y.mean)^2) SSblock <- a * sum((apply(dat2.1,1,mean) - y.mean)^2) #ブロックの平方和をSSblockという名前で保存 ##(f) 9.43式を確認 SSa + SSblock + SSba SSt ##(g) 検定統計量 (9.47)式 f.value <- ( SSa/(a-1) ) / ( SSba/((n-1)*(a-1)) ) qf(p=0.05,df1=(a-1), df2=(n-1)*(a-1), lower.tail=F) ##(h) 分散分析表 #------------------------------------------------------# #1. stack(data) # でデータセットを整形する #2. summary(aov()) # で分散分析を実行する # Error()で9.45式の確率変数を指定する #------------------------------------------------------# subset(dat2,select=c(容姿:豊かさ)) #1要因にしたい stack( subset(dat2,select=c(容姿:豊かさ)) ) dat2.2 <- stack( subset(dat2,select=c(容姿:豊かさ)) ) dat2$ブロック #ブロックをdat2.2に対応させたい rep(dat2$ブロック,3) #水準数(3)分,ブロックを繰り返せばよい dat2.3 <- data.frame(rep(dat2$ブロック,3),dat2.2) #合成させる names(dat2.3) <- c("ブロック","得点","優越属性") #変数の名前をつける dat2.3$ブロック <- factor(dat2.3$ブロック) #ブロックを名義変数に変換する summary(aov(得点~優越属性 + Error(ブロック + 優越属性:ブロック), data=dat2.3)) #表9-7の分散分析表 ################################################### #Step8:共分散分析 p.309- ################################################### ##(a) データの読み込み(表9-8) dat3 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap9_3.csv",header=T) ##(b) 図9-5を描画 (パッケージ使用) #------------------------------------------------------# # library(car) #scatterplot(従属変数~独立変数 | 質的変数, smooth=F, data=データセット名) #------------------------------------------------------# library(car) scatterplot(妬み感情~共変数|優越属性,smooth=F,data=dat3) graphics.off() ##(c) 図9-5を描画 subset(dat3,優越属性=="容姿",select=c(共変数,妬み感情)) a <- subset(dat3,優越属性=="容姿",select=c(共変数,妬み感情)) b <- subset(dat3,優越属性=="豊かさ",select=c(共変数,妬み感情)) c <- subset(dat3,優越属性=="学歴",select=c(共変数,妬み感情)) plot(a$妬み感情~a$共変数, xlim=c(0,14),ylim=c(0,10),xlab="共変数",ylab="妬み感情") abline(lm(a$妬み感情~a$共変数)) lines(b$妬み感情~b$共変数,type="p", pch=4) #pch=4で,×を指定 abline(lm(b$妬み感情~b$共変数),lty=2) #lty=2で,点線を指定 lines(c$妬み感情~c$共変数,type="p", pch=17) #pch=17で,▲を指定 abline(lm(c$妬み感情~c$共変数),lty=4) #lty=4で,-・-・-・-・-を指定 graphics.off() ##(d) 回帰係数の等質性の検定 Fullmod <- lm(妬み感情~優越属性*共変数,data=dat3) anova(Fullmod) ##交互作用が有意である場合は,共分散分析をしてはいけない ##ただし,anova()は,typeI の平方和を返す library(car) Anova(Fullmod,type="II") Anova(Fullmod,type="III") ##Anova()のtype III の平方和は,SASと定義が異なる ##(e) 共分散分析の実行(交互作用項を除くだけ) Ancovamod <- lm(妬み感情~優越属性+共変数,data=dat3) anova(Ancovamod) #表9-10 anova(Ancovamod, Fullmod) #表9-9に登場する #属性×共変数のTypeIIIの平方和 ##(f) 別解(type IIIの平方和の求め方) lm(妬み感情~優越属性*共変数,data=dat3) #因子の対比を"contr.sum"とする model1 <- lm(妬み感情~優越属性*共変数,data=dat3,contrasts=list(優越属性="contr.sum")) Anova(model1,type="III")