########################################### #Kline, R. B. (2004). #Effect size estimation in one-way designs #In Beyond significance testing: Reforming data analysis methods in behavioral research (pp.163-202). #Washington, DC: American Psychological Association. ########################################### source("http://blue.zero.jp/yokumura/R/kline/onewayes.r") #--------------------------------# #Step1: データセットの作成 (Table 6.3 p. 173) #--------------------------------# ##(a) データ tab6.3 <- data.frame(matrix( c(9,12,13,15,16, 8,12,11,10,14, 10,11,13,11,15), ncol=3 )) names(tab6.3) <- c("a","b","c") ##(b) 平均・分散・標本サイズ m <- apply(tab6.3, 2, mean) s <- apply(tab6.3, 2, var) n <- apply(tab6.3, 2, length) sig.level <- .05 ##(c) 対比の設定 contr1 <- c(1, 0, -1) contr2 <- c(1/2, -1, 1/2) psi1 <- sum(contr1 * m) #6.2式 psi2 <- sum(contr2 * m) #--------------------------------# #Step2: 対応のない1元配置分散分析(table 6.4 p.173) #--------------------------------# ##(a) 群の数,標本サイズ,grand meanの設定 a <- length(m) N <- sum(n) Mt <- mean(m) ##(b) between effect(df, MS, SS) dfa <- a-1 MSa <- sum(n*(m-Mt)^2)/(a-1) #2.29 SSa <- dfa * MSa ##(c) within effect(df, MS, SS) dfw <- N-a #2.28 MSw <- sum((n-1)*s)/sum(n-1) #2.30 SSw <- dfw * MSw ##(d) f値・p値 f.value <- MSa/MSw #2.27 p.value <- pf(f.value, dfa, dfw, lower.tail=F) ##(e) total df/total SS dft <- dfa + dfw SSt <- SSa + SSw ##(f) contrasts effect SSpsi1 <- psi1^2/sum(contr1^2/n) #6.8 SSpsi2 <- psi2^2/sum((contr2^2)/n) f.psi1 <- SSpsi1/MSw f.psi2 <- SSpsi2/MSw p.psi1 <- pf(f.psi1, 1, dfw, lower.tail=F) p.psi2 <- pf(f.psi2, 1, dfw, lower.tail=F) g.psi1 <- psi1/sqrt(MSw) g.psi2 <- psi2/sqrt(MSw) ##(g) descriptive measure eta.sq <- SSa/SSt #6.15 eta.psi1 <- SSpsi1/SSt eta.psi2 <- SSpsi2/SSt part.eta.psi1 <- SSpsi1/(SSpsi1+SSw) #6.18 part.eta.psi2 <- SSpsi2/(SSpsi2+SSw) ##(h) standard error of contrasts and confidence interval std.psi1 <- sqrt(MSw * sum(contr1^2/n)) #6.7 std.psi2 <- sqrt(MSw * sum(contr2^2/n)) upper.psi1 <- psi1 + std.psi1 * qt(p=sig.level/2, df=dfw,lower.tail=F) #6.12 lower.psi1 <- psi1 + std.psi1 * qt(p=sig.level/2, df=dfw,lower.tail=T) #6.12 upper.psi2 <- psi2 + std.psi2 * qt(p=sig.level/2, df=dfw,lower.tail=F) #6.12 lower.psi2 <- psi2 + std.psi2 * qt(p=sig.level/2, df=dfw,lower.tail=T) #6.12 upper.g.psi1 <- upper.psi1/sqrt(MSw) lower.g.psi1 <- lower.psi1/sqrt(MSw) upper.g.psi2 <- upper.psi2/sqrt(MSw) lower.g.psi2 <- lower.psi2/sqrt(MSw) t.psi1 <- sqrt(f.psi1) t.psi2 <- sqrt(f.psi2) ncp.lower1 <- TNONCT(t.psi1, dfw, prob=1-sig.level/2) ncp.upper1 <- TNONCT(t.psi1, dfw, prob=sig.level/2) ncp.lower2 <- TNONCT(t.psi2, dfw, prob=1-sig.level/2) ncp.upper2 <- TNONCT(t.psi2, dfw, prob=sig.level/2) upper.delta.psi1 <- ncp.upper1 * sqrt(sum(contr1^2/n)) #6.14 lower.delta.psi1 <- ncp.lower1 * sqrt(sum(contr1^2/n)) upper.delta.psi2 <- ncp.upper2 * sqrt(sum(contr1^2/n)) lower.delta.psi2 <- ncp.lower2 * sqrt(sum(contr1^2/n)) ##(i) table6.4 output <- data.frame(matrix(NA,ncol=8, nrow=5)) rownames(output) <- c("Between", "psi1", "psi2", "Within", "Total") colnames(output) <- c("SS", "df", "MS", "F", "p.value","g.psi", "eta.sq", "part.eta.sq") output$SS <- c(SSa,SSpsi1,SSpsi2,SSw,SSt) output$df <- c(dfa, 1, 1, dfw, dft) output$MS <- c(MSa,SSpsi1,SSpsi2,MSw,NA) output$F <- c(f.value, f.psi1, f.psi2, NA,NA) output$p.value <- c(p.value, p.psi1, p.psi2,NA,NA) output$g.psi <- c(NA,g.psi1, g.psi2,NA,NA) output$eta.sq <- c(eta.sq, eta.psi1, eta.psi2, NA,NA) output$part.eta.sq <- c(NA,part.eta.psi1, part.eta.psi2, NA,NA) class(output) <- c("anova", "data.frame") output #--------------------------------# #Step3: 対応のある1元配置分散分析(table 6.5 p.173) #--------------------------------# ##(a) データの設定 m <- apply(tab6.3, 2, mean) s <- apply(tab6.3, 2, var) n <- apply(tab6.3, 2, length) Mcov <- mean(cov(tab6.3)[lower.tri(cov(tab6.3))]) ##(b) 対比の設定 contr1 <- c(1, 0, -1) contr2 <- c(1/2, -1, 1/2) psi1 <- sum(contr1 * m) #6.2式 psi2 <- sum(contr2 * m) ##(c) 群の数,標本サイズ,grand meanの設定 a <- length(m) Mt <- mean(m) N <- sum(n) ##(d) between effect(df, MS, SS) dfa <- a-1 MSa <- sum(n*(m-Mt)^2)/(a-1) #2.29 SSa <- dfa * MSa ##(e) within effect(df, MS, SS) dfw <- N-a #2.28 MSw <- sum((n-1)*s)/sum(n-1) #2.30 SSw <- dfw * MSw ##(f) A*S error MSas <- MSw-Mcov dfas <- (a-1)*(n[[1]]-1) SSas <- dfas * MSas ##(g) Subject effect dfs <- n[[1]]-1 SSs <- SSw - SSas MSs <- SSs/dfs ##(h) f値・p値 f.value <- MSa/MSas #2.32 p.value <- pf(f.value, dfa, dfas, lower.tail=F) ##(i) total df/total SS dft <- dfa + dfw SSt <- SSa + SSw ##(j) standard error of contrasts var.d.psi1 <- as.vector(var(as.matrix(tab6.3) %*% as.matrix(contr1))) #6.9 分子 var.d.psi2 <- as.vector(var(as.matrix(tab6.3) %*% as.matrix(contr2))) ##(k) contrasts effect SSpsi1 <- psi1^2/sum(contr1^2/n) #6.8 SSpsi2 <- psi2^2/sum((contr2^2)/n) f.psi1 <- SSpsi1/MSas f.psi2 <- SSpsi2/MSas p.psi1 <- pf(f.psi1, 1, dfas, lower.tail=F) p.psi2 <- pf(f.psi2, 1, dfas, lower.tail=F) g.psi1 <- psi1/sqrt(MSw) g.psi2 <- psi2/sqrt(MSw) ratio.psi1 <- psi1/sqrt(var.d.psi1) ratio.psi2 <- psi2/sqrt(var.d.psi2) ##(l) descriptive measure eta.sq <- SSa/SSt #6.15 eta.psi1 <- SSpsi1/SSt eta.psi2 <- SSpsi2/SSt part.eta <- SSa /(SSa + SSas) part.eta.psi1 <- SSpsi1/(SSpsi1+SSas) #6.18 part.eta.psi2 <- SSpsi2/(SSpsi2+SSas) ##(m) confidence interval std.psi1 <- sqrt(var.d.psi1/n[[1]]) #6.9 std.psi2 <- sqrt(var.d.psi2/n[[1]]) upper.psi1 <- psi1 + qt(sig.level/2, dfs, lower.tail=F) * std.psi1 #6.12 lower.psi1 <- psi1 + qt(sig.level/2, dfs, lower.tail=T) * std.psi1 upper.psi2 <- psi2 + qt(sig.level/2, dfs, lower.tail=F) * std.psi2 lower.psi2 <- psi2 + qt(sig.level/2, dfs, lower.tail=T) * std.psi2 upper.gpsi1 <- upper.psi1/sqrt(MSw) lower.gpsi1 <- lower.psi1/sqrt(MSw) upper.gpsi2 <- upper.psi2/sqrt(MSw) lower.gpsi2 <- lower.psi2/sqrt(MSw) ##(n) table6.5 output <- data.frame(matrix(NA,ncol=9, nrow=7)) rownames(output) <- c("Between", "psi1", "psi2", "Within", "Subj", "A*S(error)","Total") colnames(output) <- c("SS", "df", "MS", "F", "p.value","g.psi", "ratio","eta.sq", "part.eta.sq") output$SS <- c(SSa,SSpsi1,SSpsi2,SSw, SSs, SSas,SSt) output$df <- c(dfa, 1, 1, dfw, dfs, dfas, dft) output$MS <- c(MSa,SSpsi1,SSpsi2,MSw, MSs, MSas,NA) output$F <- c(f.value, f.psi1, f.psi2, NA,NA,NA,NA) output$p.value <- c(p.value, p.psi1, p.psi2,NA,NA,NA,NA) output$g.psi <- c(NA,g.psi1, g.psi2,NA,NA,NA,NA) output$ratio <- c(NA,ratio.psi1, ratio.psi2, NA,NA,NA,NA) output$eta.sq <- c(eta.sq, eta.psi1, eta.psi2, NA,NA,NA,NA) output$part.eta.sq <- c(part.eta,part.eta.psi1, part.eta.psi2, NA,NA,NA,NA) class(output) <- c("anova", "data.frame") output #--------------------------------# #Step4: example of omega and intera class correlation (p.187-) #--------------------------------# a <- 3; n <-30 SSa <- 60; SSw <- 478.5; SSt <- 538.5 dfa <- a-1 MSa <- SSa/(a-1) #2.29 MSw <- SSw/((n*a)-a) #2.28 f.value <- MSa/MSw #2.27 eta.sq <- SSa/SSt #6.15 varcomp.a <- ((a-1)/(a*n))*(MSa-MSw) #table 6.7 omega.sq <- (dfa*(MSa-MSw))/(SSt+MSw) #6.31 rho.inter <- (MSa-MSw)/(MSa + (n-1)*MSw) #6.32 #--------------------------------# #Step5: 平均・標準偏差・標本サイズから1元配置の分散分析関数 #--------------------------------# #--------------------------------# #oneway.anova.second(m, sd, n, sig.level=.05, unbiased=T, graph=F, contr=NULL, digits=3) #--------------------------------# #library(gtools)がインストール済みであること #m = 平均値; sd = 標準偏差; n = 標本サイズ; sig.level =有意水準(デフォルトは,α=.05) #unbiased # unbiased=Tの場合は,不偏分散を元にした標準偏差を使っていることを指定する(デフォルト) # unbiased=Fの場合は,標本分散を元にした標準偏差を使っていることを指定する #graph # graph=Fの場合は,検出力曲線を描画しない(デフォルト)母集団効果量は,Cohen の規準を採用 # graph=Tの場合は,検出力曲線を描画する #contr= 対比(指定しなければ,全ての組み合わせ) #digits=表示桁数(デフォルトは3) #--------------------------------# ##(a) table 6.4 m <- c(13, 11, 12) sd <- sqrt(c(7.5, 5, 4)) n <- 5 contr <- matrix(c(1,0,-1,1/2,-1,1/2), ncol=3, nrow=2, byrow=T) oneway.anova.second(m, sd, n, contr=contr) ##(b) Table 6.12 oneway.anova.second( m=c(218.9, 221.1, 218.7), sd=c(28.2, 26.3, 27.5), n =28, sig.level=.05) ##attention selective oneway.anova.second( m=c(532, 484.4, 478.6), sd=c(65.4, 57.9, 48.4), n =28, sig.level=.05) ##learning and abstract thinking verbal oneway.anova.second( m=c(4.46, 3.71, 3.29), sd=c(.79, 1.15, 1.12), n =28, sig.level=.05) ##learning and abstract thinking visual oneway.anova.second( m=c(4.61, 4, 4.11), sd=c(.96, 1.41, 1.13), n =28, sig.level=.05) ##learning and abstract thinking abstract thinking oneway.anova.second( m=c(25.96,29.46,29.5), sd=c(4.1,4.19,3.64), n =28, sig.level=.05) ##(c) Table 6.13 n <- c(129,211,171,78,38,40) m <- c(12.02,10.95,10.18,9.58,9.89,9.05) sd <- c(2.95,3.1,2.91,2.84,3.34,3.82) contr <- c(1/3,1/3,1/3,-1/3,-1/3,-1/3) oneway.anova.second(m, sd, n, contr=contr)