myscratchpad

# 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

#　対応のない1要因
# データ生成。テクニカルブックp. 87
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9, 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))

# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)

# lmで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a, dat)
# carのAnovaで分析
library(car)
Anovares <- Anova(lmres, type=3)
Anovares

# aovと比較
summary(aov(result~a, dat))


#　対応のある1要因
# データ生成。テクニカルブックp. 92
dat <- data.frame(s=factor(1:8), a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))

# 平均と合計
sapply(dat[2:5], mean)
sapply(dat[2:5], sum)

# 分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1, a2, a3, a4)~1, dat)
# carのAnovaで分析
library(car)
afact <- factor(c("a1", "a2", "a3", "a4"))
idat <- data.frame(afact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact, type=3)
summary(Anovares)
(AnovaTable <- Anovastat(Anovares))

# aovと比較
library(reshape)
dat2 <- melt(dat, id="s")
names(dat2) <- c("s", "a", "result")
summary(aov(result~a+s+Error(a:s), dat2))

carパッケージ、Anova関数の分散分析表とかはprintで結果を出力しているだけ。
データフレームとか行列になってると、後からF値や自由度を取り出して効果量なんかを出しやすいので関数 "Anovastat" を作ってみた。
少し改良して"Astat"にした。
といっても、car:::summary.Anova.mlmやcar:::print.linear.hypothesis.mlmを適当に書き換えただけ。Rjpwikiの記事に感謝。

## 関数の読み込み
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1264349365")
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1273264153")
# 我ながらきたないコードだなぁ…

## 使い方

# 分析サンプル
library(car)
hsb2 <- read.table('http://www.ats.ucla.edu/stat/R/faq/hsb2.csv', header=T, sep=",")
result1 <- Anova(lm(cbind(write,read)~female+math+science+socst, hsb2))


## Anovastatで結果の取り出し
Tableres <- Anovastat(result1)
Tableres <- Astat(result1)
Tableres
# manovaはWilksをデフォルトで出力する。単変量Anovaには偏イータ二乗、f二乗も出力する。
Tableres$Anova # 単変量分散分析の結果。混合要因のときのみ出力される (上の分析例では出力されない) 。通常はUAnova (単変量分散分析表) 、Mauchly, GG(Greenhouse-Geisser) 、HF (Huynh-Feldt) がmatrixで出力される。
# Tableres$Anova$UAnovaとかで取り出す
Tableres$Manova # 多変量分散分析の結果
Tableres$Manova[[2]] # とかやれば個々の分散分析表がデータフレームで取り出せる

## 被験者間要因のみのときはそもそもデータフレームなので上の関数は必要ない
result2 <- Anova(lm(write~female*math, hsb2))
class(result2)