################################################ #南風原（2002）第6章 #Last update: 2006/01/09 13:21:20 ################################################ #Step1：t分布を用いた正確な検定 p.161-64 ################################################ ##(a)データを読み込み，検定する2つの変数をy1,y2と命名 setwd("C:/data/haebara/chap01") data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T) 逸脱男 <- subset(data1,性別=="男",select=c(変化量)) #男性だけのデータを抽出 逸脱女 <- subset(data1,性別=="女",select=c(変化量)) #女性だけのデータを抽出 y1 <- 逸脱男$変化量 y2 <- 逸脱女$変化量 ##(b)標本分散の定義 2.12式（svarと命名） svar <- function(x){ #BeginFnc 1/length(x)*sum((x-mean(x))^2) } #EndFnc ##(c)(6.7)式 分子 mean(y1)-mean(y2) ##(d)母集団標準偏差σの推定量 (6.4)式 s.asta <- sqrt( (length(y1) * svar(y1) + length(y2) * svar(y2))/ (length(y1) +length(y2)-2) ) ##(e)平均値差の標準誤差の推定量(6.6)式 mean.std <- s.asta * sqrt(1/length(y1) + 1/length(y2)) ##(f)t値(6.7)式 t.value <- (mean(y1)-mean(y2))/mean.std ##(g)自由度38のt分布の上側確率.025に対応する値（棄却域） #qt(p=有意確率,df=自由度,lower.tail=F 上側確率) qt(p=0.025,df=38,lower.tail=F) ##(h)t.test(x,y, var.equal=T) 関数を使うと一発で求まる #---------------------------------------------------- #var.equal=T→2群の分散が等質である（t検定） #var.equal=F→2群の分散が等質でない（ウェルチの検定） #---------------------------------------------------- t.test(y1,y2,var.equal=T) ################################################ #Step2：効果量を用いた独立な2群の平均値差の検定（両側検定）における # 棄却の限界値 付表6 を再現 p.381 ################################################ ##(a)付表4 t分布の上側確率に対応する値を再確認 qt(p=0.025,df=18,lower.tail=F) #qt(p=有意確率,df=自由度,lower.tail=F 上側確率) ##(b)両群のサンプルサイズが等しいと考えると 6.8式は #---------------------------------------------------- # t= d*sqrt( (n1^2)/(n1*2) ) # = d * sqrt(n1/2) と展開できる #従って， d = t * sqrt(2/n1) となる． #以上の結果から，t分布の限界値を利用して，効果量の場合に変換すると #---------------------------------------------------- n <- c( seq(2,20), seq(22,40,by=2),seq(50,100,by=10),seq(200,400,by=100) ) #各群のサンプルサイズを指定 x1 <- qt(p=0.05, n*2-2,lower.tail=F) * sqrt(2/n) x2 <- qt(p=0.025, n*2-2,lower.tail=F) * sqrt(2/n) x3 <- qt(p=0.005, n*2-2,lower.tail=F) * sqrt(2/n) data.frame(x1,x2,x3, row.names=n) #data.frame()で変数をまとめる ################################################ #Step3：独立な2群の平均値差の検定の検定力 # （有意水準.05の両側検定）を再現 付図3 p.388 ################################################ ##(a)検定力を求める関数を定義（関数の中身は覚えなくても大丈夫） power.t <- function(delta,n1,n2,sig.level){ pt(qt(p=sig.level/2,df=n1+n2-2,lower.tail=T), df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=T)+ pt(qt(p=sig.level/2,df=n1+n2-2,lower.tail=F), df=n1+n2-2,ncp=delta*sqrt((n1*n2)/(n1+n2)),lower.tail=F) } ##(b)power.t(delta=母集団効果量, # n1=第1群の標本サイズ,n2=第2群の標本サイズ,sig.level=有意水準(両側検定)) power.t(delta=1,n1=20,n2=20,sig.level=0.05) #母集団効果量=1,各群の標本サイズ=20,有意水準=0.05（両側検定） ##(c)描画 delta <- seq(-2,2,by=0.01) #x軸の指定 power.t(delta=delta,n1=10,n2=10,sig.level=0.05) #テスト plot(power.t(delta=delta,n1=5,n2=5,sig.level=0.05)~ delta, ylim=c(0,1),type="l", xlab="母集団効果量", ylab="検定力") #n=5の場合 lines(power.t(delta=delta,n1=20,n2=20,sig.level=0.05)~ delta,col=2) #n=20の場合 lines(power.t(delta=delta,n1=200,n2=200,sig.level=0.05)~ delta,col=3) #n=200の場合 graphics.off() #グラフを閉じる ################################################ #Step4：効果量を用いた検定，検定力とサンプルサイズ p.164-65 ################################################ ##(a)効果量を用いた検定 d <- (mean(y1)-mean(y2))/s.asta #(6.9)式 d qt(p=0.025, 20*2-2,lower.tail=F) * sqrt(2/20) #棄却の限界値 #(b)検定力とサンプルサイズ power.t(delta=1,n1=20,n2=20,sig.level=0.05) #母集団効果量=1,各群の標本サイズ=20,有意水準=0.05（両側検定） power.t(delta=0.5,n1=20,n2=20,sig.level=0.05) #母集団効果量=0.5,各群の標本サイズ=20,有意水準=0.05（両側検定） power.t(delta=1,n1=25,n2=15,sig.level=0.05) #母集団効果量=1,各群の標本サイズ=25と15,有意水準=0.05（両側検定） power.t(delta=0.5,n1=25,n2=15,sig.level=0.05) #母集団効果量=0.5,各群の標本サイズ=25と15,有意水準=0.05（両側検定） ########################################## #Step5：平均値差および効果量の区間推定 p.167-71 ########################################## ##(a)平均値差の区間推定 mean.std #平均値差の標準誤差の推定量(6.6)式 mean(y1)-mean(y2) #平均値差 qt(0.025,df=38,lower.tail=F) delta.L <- (mean(y1)-mean(y2)) - qt(0.025,df=38,lower.tail=F) * mean.std #(6.13)式 delta.U <- (mean(y1)-mean(y2)) + qt(0.025,df=38,lower.tail=F) * mean.std #(6.13)式 delta.L delta.U t.test(y1,y2,var.equal=T) #t.test()でも95 percent confidence intervalが算出される ##(b)効果量の区間推定 n1 <- length(y1); n2 <- length(y2) #各群の標本サイズを求める d.sem <- sqrt( ((n1+n2)/(n1*n2)) + (d^2)/(2*(n1+n2-2)) ) #dの標準誤差の推定値 6.17式 d - 1.96 * d.sem d + 1.96 * d.sem ##(c)6.13/14式を関数として定義する d.lower <- function(n1,n2,d){ d.sem <- sqrt( ((n1+n2)/(n1*n2)) + (d^2)/(2*(n1+n2-2)) )#dの標準誤差の推定値 6.17式 d - 1.96 * d.sem } d.upper <- function(n1,n2,d){ d.sem <- sqrt( ((n1+n2)/(n1*n2)) + (d^2)/(2*(n1+n2-2)) ) #dの標準誤差の推定値 6.17式 d+1.96*d.sem } d.lower(n1=20,n2=20,d=-.77) #テスト d.upper(n1=20,n2=20,d=-.77) #テスト ##(d)効果量の標本分布の上側および下側確率.025に対応する値 付図4 を再現 p.388 d <- seq(-4,4,by=.01) #x軸を定義 plot(d.lower(n1=5,n2=5,d=d)~ d, type="l", xlab="標本効果量", ylab="母集団効果量", xlim=c(-3,3), ylim=c(-2,2), col=1 ) lines(d.upper(n1=5, n2=5, d=d)~ d, col=1) #標本サイズが5の場合 lines(d.lower(n1=10, n2=10, d=d)~ d, col=2) #標本サイズが10の場合 lines(d.upper(n1=10, n2=10, d=d)~ d, col=2) #標本サイズが10の場合 lines(d.lower(n1=50, n2=50, d=d)~ d, col=3) #標本サイズが50の場合 lines(d.upper(n1=50, n2=50, d=d)~ d, col=3) #標本サイズが50の場合 lines(d.lower(n1=100,n2=100,d=d)~ d, col=4) #標本サイズが100の場合 lines(d.upper(n1=100,n2=100,d=d)~ d, col=4) #標本サイズが100の場合 lines(d.lower(n1=200,n2=200,d=d)~ d, col=5) #標本サイズが200の場合 lines(d.upper(n1=200,n2=200,d=d)~ d, col=5) #標本サイズが200の場合 graphics.off() #図を閉じる ##(e)d.sample()で，信頼区間を利用したサンプルサイズの決定 ## （関数の中身は覚えなくても大丈夫） d.sample <- function(d,w,conf.level,interval=c(4,10e10)){ temp <- function(N){ (d + qnorm( (1-conf.level)/2, lower.tail=F) * sqrt(2/N + d^2/(4*(N-1))))- (d - qnorm( (1-conf.level)/2, lower.tail=F) * sqrt(2/N + d^2/(4*(N-1))))-w } ceiling(uniroot(temp,interval)[[1]]) } d.sample(d=.8,w=0.5,conf.level=.95) #d=予想される効果量，w=望まれる信頼区間の幅，conf.level=信頼係数 ########################################## #Step6：対応のある2群の平均値差の検定と推定 ########################################## ##(a)表6-1 男子20人と小6と中2のときの逸脱行動得点と変化量のデータを作成 逸脱男 <- subset(data1,性別=="男",select=c(小6,中2,変化量)) 逸脱男 mean(逸脱男) sd(逸脱男) ##(b)「変化量」の標本平均の標準誤差の推定値 6.25式 attach(逸脱男) #attach()で変数名を使えるようにする Sv <- sd(変化量)/sqrt(length(変化量)) ##(c) 検定統計量tの値 6.29式 mean(変化量)/Sv #教科書との値が僅かに違うが，丸め誤差の影響である． ##(d) 自由度19のt分布の上側確率.025に対応する値 qt(p=0.025,df=19,lower.tail=F) ##(e)平均値差の区間推定 6.30/31式 delta.lower <- mean(変化量) - qt(p=0.025,df=19,lower.tail=F) * Sv delta.upper <- mean(変化量) + qt(p=0.025,df=19,lower.tail=F) * Sv ##(f)t.test(x,y, var.equal=T, paired=T)で，(b)〜(e)の作業を1行で計算する t.test(中2,小6,var.equal=T, paired=T) detach(逸脱男) ##作業が終了したらdetach()を使う ########################################## #Step7:独立な2群の比率差の検定 ########################################## ##(a)データの読み込みと修正 table6.2 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/prop.csv",header=T) table6.2 #データの確認 attach(table6.2) #attach()で変数名を使えるようにする table(性別,経験の有無) #table(x,y)でクロス表を作る #男子と女子の順番が逆になってしまっている 性別 #levelsをみると，「女子，男子」という順番になっている #factor(data, levels=c(任意の順番))で，levelsを修正する factor(性別,levels=c("男子","女子")) 性別 <- factor(性別,levels=c("男子","女子")) #修正結果がOKならば，変数名を保存する 性別 #OK table(性別,経験の有無) #完成 ##(b) vcdライブラリにある #------------------------------------------- # table2d_summary(table(x,y), percentages=F, conditionals=c("none","row","column")) # percentages=T → 比率が求まる # conditionals="none"（デフォルト）→ 合計（この場合，40名）を中心にした比率を求める # conditionals="row"　→ 行の合計を中心とした比率を求める（男女共に20名） # conditionals="column" → 列の合計を中心とした比率を求める（ある = 16名，ない=24名） #------------------------------------------- # を使う library(vcd) table2d_summary(table(性別,経験の有無)) #比率なし（表6-2と同一） table2d_summary(table(性別,経験の有無) ,percentages=T) #各セルの度数/セルの度数の合計が%に表示 table2d_summary(table(性別,経験の有無) ,percentages=T, conditionals="row") #row%=度数/一番右側の合計 table2d_summary(table(性別,経験の有無) ,percentages=T, conditionals="column") #col%=度数/一番下側の合計 ##(c)tableの要素を取り出して，比率の検定をする table(性別,経験の有無)[1,1] #「男子」で「ある」と回答した人数 # table()[行,列]で指定する p1 <- table(性別,経験の有無)[1,1]/20 #p.179上（男子） p2 <- table(性別,経験の有無)[2,1]/20 #p.179上（女子） p <- 16/40 n1 <- 20 n2 <- 20 z <- (p1-p2)/sqrt(p*(1-p)*(1/n1+1/n2)) #6.37式 z qnorm(.025,lower.tail=F) #標準正規分布の上側確率.025に対応する値 ##(e) prop.test, chisq.test, table2d_summary (2群の比率差の検定に使える関数) library(epitools) prop.test(table(性別, 経験の有無), correct = F) #強み: 信頼区間が出る table2d_summary(table(性別,経験の有無)) #強み: クロス集計表と同時に見れる chisq.test(table(性別,経験の有無), correct =F) #強み: p値のシミュレーションができる epitab(table(性別,経験の有無), pvalue = "chi2") #library(epitools) z^2 #z^2 = chisq detach(table6.2) #作業が終了したらdetach()を使う ########################################## #Step8:対応のある2群の比率差の検定 p.181-83 ########################################## ##(a)データを読み込み table6.3 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/prop2.csv",header=T) attach(table6.3) ##(b)6.38式を計算 n12 <- table(逸脱行動A,逸脱行動B)[1,2] m <- table(逸脱行動A,逸脱行動B)[1,2] + table(逸脱行動A,逸脱行動B)[2,1] z <- (n12-m/2) / (sqrt(m)/2) z qnorm(.025,lower.tail=F) #標準正規分布の上側確率.025に対応する値 ##(c) mcnemar.test z^2 mcnemar.test(table(table6.3), correct=F) ########################################## #Step9:連関の検定 p.184-90 ########################################## ##(a) epitoolsライブラリにある #------------------------------------------- # expected(table(x,y)) #------------------------------------------- #を使うと推定期待度数が求まる expected(table(逸脱行動A,逸脱行動B)) #表6.5（かっこ）内 table6.5 <- table(逸脱行動A,逸脱行動B) table6.5 detach(table6.3) #作業が終了したらdetach()を使う ##(b) 6.39式，6.40式を確認 e11 <- (sum(table6.5[1,]) * sum(table6.5[,1]))/40 e12 <- (sum(table6.5[1,]) * sum(table6.5[,2]))/40 e21 <- (sum(table6.5[2,]) * sum(table6.5[,1]))/40 e22 <- (sum(table6.5[2,]) * sum(table6.5[,2]))/40 chisq <- (table6.5[1,1]-e11)^2/e11 + (table6.5[1,2]-e12)^2/e12 + (table6.5[2,1]-e21)^2/e21 + (table6.5[2,2]-e22)^2/e22 ##(c) クラメルの連関係数とファイ係数 V <- sqrt(chisq/40) #6.42式 V #6.44式 phi <- ( table6.5[1,1] * table6.5[2,2] - table6.5[1,2] * table6.5[2,1] )/ sqrt( sum(table6.5[1,]) * sum(table6.5[2,]) * sum(table6.5[,1]) * sum(table6.5[,2]) ) phi ##(d) assocstats assocstats(table6.5) #library(vcd) ##(e) 付表8を確認 p.383 #------------------------------------------- #qchisq(p=確率, df=自由度, lower.tail = F（上側確率）) #------------------------------------------- qchisq(p=.05,df=1,lower.tail=F) df <- seq(1,30) x1 <- qchisq(p=.1,df=df,lower.tail=F) x2 <- qchisq(p=.05,df=df,lower.tail=F) x3 <- qchisq(p=.01,df=df,lower.tail=F) data.frame(x1,x2,x3,row.names=df) #data.frame()で変数をまとめる