myscratchpad

# ちょっと復習
## 心理統計学の基礎、p. 9 変化量データ。
y1 <- c(4, 3, -3, 4, 1, 4, 0, 6, 2, -6, 2, 6, 6, -2, 1, 1, 6, -2, 7, -1) # 男子
y2 <- c(9, 0, 6, 5, 5, 5, 2, 11, 5, 3, 4, 7, 4, 10, -2, 6, 1, 2, 4, 4) # 女子
n1 <- 20
n2 <- 20
m1 <- mean(y1)
m2 <- mean(y2)
(tres <- t.test(y1, y2, var.equal=T))

# 標本分散の場合
s1 <- var(y1)*((n1-1)/(n1))
s2 <- var(y2)*((n2-1)/(n2))
sps <- sqrt(((s1*n1+s2*n2)/(n1+n2-2)))
(m1-m2)/sps

# 不偏分散の場合
v1 <- var(y1)
v2 <- var(y2)
spe <- sqrt(((v1*(n1-1)+v2*(n2-1))/(n1+n2-2)))
(m1-m2)/spe

# つまるところ、プールされた標準偏差の分子は2群それぞれの平方和を足したもので、それを(n1+n2-2) = 自由度　で割る
# Peasonのrは対応のないt検定だけ

## t.testの結果オブジェクト、もしくはt値とサンプルサイズを与えて効果量を算出する
tef <- function(tres=NULL, tv, n1, n2) {
  if (is.null(tres)) {
      tv <- tv
      df <- n1+n2-2
      }
      else {
       tv <- tres$statistic
       df <- tres$parameter
      }
   #dv1 <- (2*tv)/sqrt(df)
   dv2 <- (2*tv)/sqrt(df+2)
   r1 <- sqrt(tv^2/(tv^2+df))
   #print("# G*Powerと一致するのはdv2")
   #c(dv1=dv1, dv2=dv2, r=r1)
   res <- c(dv2, r1)
   names(res) <- c("Cohen's d", "Pearson's r")
   res
 }

tef(tres)
tef(tv=2.20, n1=50, n2=60)

# 芝・南風原 (1990). 行動科学における統計解析法　p251　より
n1v <- c(50,30,29,25,20,20,15,15)
n2v <- c(60,35,30,30,26,24,18,15)
tvv <- c(2.20,0.91,0.76,1.15,-0.06,0.93,1.64,0.38)

dvs <- list()
for (i in 1:8) {
  n1i <- n1v[i]
  n2i <- n2v[i]
  tvi <- tvv[i]
  x <- tef(tv=tvi, n1=n1i, n2=n2i)
  dvs[[i]] <- x
  }
dvs


# 効果量dを他の色んな効果量に変換するExcelシート
http://www.psychsystems.net/Manuals/StatsCalculators/Effect_Size_Calculator%2017.xls

# psychometricパッケージのCI.Rsqlm関数を用いる
# CI.Rsqでもいいが、lmオブジェクトを受け取れるCI.Rsqlmの方が楽そう

# データ準備
x1 <- rnorm(100)
x2 <- rnorm(100)
y <- x1+x2^2+rnorm(100)
lmres <- lm(y~x1+x2)
summary(lmres)

# 信頼区間
library(psychometric)
CI.Rsqlm(lmres)
## 回帰係数の信頼区間は組み込みのconfintで
confint(lmres)

psychometricパッケージの関数一覧
http://rss.acs.unt.edu/Rdoc/library/psychometric/html/00Index.html

開発者のThomas D. Fletcher
http://www.umsl.edu/~fletchert/quant/

# 相関係数の信頼区間はcor.testで出る
x <- rnorm(100)
y <- rnorm(100)
cor.test(x,y)
## あと、psychometoricパッケージのCIrでも出る。これはrを直接与える

参考は偏イータ二乗と同じ

source("http://psychology3.anu.edu.au/people/smithson/details/CIstuff/Nonct.ssc")
tv <- 2.141
df.t <- 58
n1 <- 30
n2 <- 30

# 効果量Cohen's d
(d <- tv/(sqrt((n1*n2)/(n1+n2))))
# 関連度r
(rv <- sqrt((tv^2)/(tv^2+df.t)))

# t値の信頼区間
lot(tv, df.t, 0.95)
hit(tv, df.t, 0.95)

# dの信頼区間
ncplow <- lot(tv, df.t, 0.95)[1]
ncpupp <- hit(tv, df.t, 0.95)[1]
(dlow <- ncplow/sqrt(n1*n2/(n1+n2)))
(dupp <- ncpupp/sqrt(n1*n2/(n1+n2)))
# 関連度rの信頼区間。単なる点双列相関なのでcor.testでもいいが
library(psychometric)
(CIr(rv, df.t+2, 0.95))


# Cohenのdについてメモ
(d <- (m1-m2)/sqrt(sp))
## 普通は分散とNからspを求める
var1 <- 2.001^2 # 普通はローデータを使おう
var2 <- 3.073^2
## sp (プールされた分散) について。とりあえずこれ使おう
sqrt((var1+var2)/2) # よく見る効果量の分母
sqrt((((n1-1)*var1)+((n2-1)*var2))/(n1+n2)) # 各群の人数が違うとき
### 心理統計学の基礎 p. 162では違う数式をあげている。色々回ってみるとこれはレアケースのようだ
sqrt((var1*n1+var2*n2)/(n1+n2-2)) # 心理統計学の基礎 p. 162
sqrt((var1*n1+var2*n2)/(n1+n2)) # 心理統計学の基礎 p. 162で、分母の-2をとった (自由度じゃなくてNの総計にした) ら#よく見る…と一致する
sqrt((((n1-1)*var1)+((n2-1)*var2))/(n1+n2-2)) # 各群の人数が違うとき。これは-2すると他と一致する

参考
1) Scripts and Software for Noncentral Confidence Interval and Power Calculations
2) Exploring Effect Size and Measures of Association

# 参考1)より、関数lof, hifを利用する
source("http://psychology3.anu.edu.au/people/smithson/details/CIstuff/Fscript.SSC")

# F値、df1 (分子自由度) 、df2 (分母自由度) を設定しておく
df1 <- 3
df2 <- 461
Fv <- 92.901

# とりあえず偏イータ二乗をだす
(p.eta <- (df1*Fv) / (df2 + (df1*Fv)))

# lof, hif関数利用
lof(Fv, df1, df2, 0.95)
hif(Fv, df1, df2, 0.95)

# 偏イータ二乗の信頼区間にする
ncplow <- lof(Fv, df1, df2, 0.95)[1]
ncpupp <- hif(Fv, df1, df2, 0.95)[1]
(rsq <- df1 * Fv / (df2 + df1 * Fv)) # これは偏イータ二乗と同じ
(rsqlow <- ncplow / (ncplow + df1 + df2 + 1))
(rsqupp <- ncpupp / (ncpupp + df1 + df2 + 1))
# lof, hifだけでは偏イータ二乗の信頼区間は出ないので、sasの方のコードを援用

偏イータ二乗とイータ二乗の信頼区間の求め方は同じなんだろうか…

# MBESSというパッケージがあるらしい
http://www.nd.edu/~kkelley/site/MBESS.html

リンクだけ。そのうち調べよう

Loftus, G. R., & Masson, M. E. J. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin & Review, 1, 476-490.


Cousineau, D. (2007). Confidence intervals in within-subjects designs: A simpler solution to Loftus and Masson's method. Tutorials in Quantitative Methods for Psychology, 1, 42-45.
Masson, M. E. J. & Loftus, G. R. (2003). Using confidence intervals for graphically based data interpretation. Canadian Journal of Experimental Psychology, 57, 203-220.

S-PLUSで計算する関数があった。Rのことも書いてある。親切な人がいるものである
Wright, D. B. (2007). Graphing within-subjects confidence intervals using SPSS and S-Plus. Behavior Research Methods, 39, 82-85.
解説HP: http://www.fiu.edu/~dwright/WSCI/wscihome.htm