九州理学療法士・作業療法士合同学会2017in宮崎で紹介したスクリプト

モンテカルロ法によるｔ検定の検定力分析
モンテカルロ法については以下を参照ください

yoshida931.hatenablog.com

tval <- numeric(length = 10000)　　#ｔ値(空のベクトルを設定)
cs <- 0
for(i in 1:10000){
 x <- rnorm(n = 30, mean = 10, sd = 3)  　# xの乱数30個生成
  y <- rnorm(n = 30, mean = 12, sd = 3) 　 # yの乱数30個生成
   res <- t.test(x, y, var.equal = T)　　　　　　　#xとｙのｔ検定の結果
    tval[i] <- res[[1]]　　　　　　　　　　　　　　　#ｔ値のみを抜き取る
     cs <- cs + ifelse(res[[3]] <  0.05, 1,0)　
}

#cs：ｐ値が0.05より小さい場合には1を足し、大きい場合には0を足す
#ｔ検定を1000回繰り返した結果、0.05より小さいP値の個数となる

cs/10000

データの要約、散布図、エラーバー付き平均値の一括処理

x<-iris$Sepal.Length[51:75]
y<-iris$Petal.Length[101:125]


sumxy <- function(x, y, TEST01){
 
  library(pwr)

  
  cat1 <- c("Min.", "1st Qu", "Median", "Mean", "3rd Qu.", "Max.")   
  xs <- summary( x ) 　　　　　　　　　　　　　　　　　　　　　　　
  ys <- summary(y)
  
  cat(cat1,"\n","x",xs,"\n","y",ys,"\n\n",file = "TEST01");
  tx <- t.test(x);　　　　　　　　　　　　　
  cat2<- c("tvalue","pvalue","95%CI");　　
  txt <- tx$statistic;　　　　　　　　　
  txp <- tx$p.value;　　　　　　　　　　
  txc1 <- tx$conf.int[1];　　
  txc2 <- tx$conf.int[2];　　　　　　
  x1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txt,txp,txc1,txc2); 
  
  ty <- t.test(y); 　　　
  tyt <- ty$statistic;
  typ <- ty$p.value;
  tyc1 <- ty$conf.int[1];
  tyc2 <- ty$conf.int[2];
  y1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f", tyt, typ, tyc1, tyc2);
  
  txy <- t.test(x,y,var.equal = T);
  txyt <- txy$statistic;　　　　　　
  txyp <- txy$p.value;　　　　　　　　　
  txyc1 <- txy$conf.int[1];　　　　　　　
  txyc2 <- txy$conf.int[2];　　　　　　　
  xy1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txyt,txyp,txyc1,txyc2);
  
  d <- (mean(x) - mean(y)) / sqrt((var(x) + var(y)) / 2);
  p <- pwr.t.test(n = length(x), d = d, sig.level = 0.05);
  
  xy <- c(x,y);
  zu1 <- pdf(file = "散布図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  plot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu1
  
  zu2 <- pdf(file = "エラーバー.pdf");
  se2 <- c(sqrt(var(x)/(length(x))), sqrt(var(y)/(length(y))));
  m2 <- c(mean(x), mean(y));
  plot(m2, xlim=c(0.5,2.5), type="b", ylim=c(5,7), xaxt="n",xlab = "", ylab = "");
  arrows(1:2,m2-se2,
         1:2,m2+se2,
         lwd=1.0,angle=90,length=0.1,code=3);
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu2
  
  zu3 <- pdf(file = "箱ひげ図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  boxplot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu3
  
  cat(cat2, "\n", "x", x1, "\n", "y", y1, "\n", "xy", xy1, "\n",
      "x:1標本t検定, y:1標本t検定, xy:独立2群t検定", "\n", "\n",
      "独立2群t検定の検定力", p[[4]], "\n", file = "TEST01", append =T);
}

sumxy(x, y, "TEST01")
#　出力されたファイルはDocumentsにあると思います