myscratchpad

a=gl(2, 10, labels=c("a1", "a2"))
b=rep(gl(2, 5, labels=c("b1", "b2")), 2)
value=c(1:5, 2:6, 6:10, 4:8)

aovres <- aov(value~a*b)
summary(aovres)
model.tables(aovres, "means")

# 1要因のプロット
library(gplots)
plotmeans(value~a)

# 2要因のプロット
interaction.plot(x.factor=a, trace.factor=b, response=value)

John FoxのIntroduction to Statistical Computing in R
http://socserv.socsci.mcmaster.ca/jfox/Courses/R-course/index.html

William RevelleのAn introduction to psychometric theory with applications in R 。"書きかけ"のe-book
http://www.personalitytheory.com/r/book/

Psych Web
http://www.psychwww.com/


UCLAのセミナー
http://www.ats.ucla.edu/stat/R/seminars/Repeated_Measures/repeated_measures.htm

J. J. McDowelの授業。Excelファイル
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/psych560.htm

分散分析の授業
http://psych.wisc.edu/moore/R_Analysis_of_Variance_Handouts_html.html
http://psych.wisc.edu/moore/610HandoutIndex.html

Karl Wuenschのsasコード
http://core.ecu.edu/psyc/wuenschk/SAS/SAS-Programs.htm
いろいろリンク
http://core.ecu.edu/psyc/wuenschk/links.htm
anova
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc

Howell, D. C. (2010). Statistical methods for psychology (7th ed.). Belmont, CA: Cengage Wadsworth.
amazon
Publisher

Tabachnick, B. G., & Fidell, L. S. (2007). Using multivariate statistics (5th ed.). Boston: Allyn & Bacon.
amazon
Publisher

# 心理統計学の基礎から
dat <- read.delim("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1258223909")

library(psych)
dcrp <- describe.by(dat$envy, dat$dominance.e)
print(dcrp, digits=3)

value <- dat$envy
dmn <- dat$dominance.e

gmean <- mean(value)
fmean <- tapply(value, list(dmn), mean)
a <- nlevels(dmn)
lvs <- levels(dmn)
na1 <- sum(dmn==lvs[1]); na2<- sum(dmn==lvs[2]); na3<- sum(dmn==lvs[3])
n <- length(value)

sst <- sum((value-gmean)^2) # 全体平方和

ssa <- (na1*(fmean[1]-gmean)^2)+(na2*(fmean[2]-gmean)^2)+(na3*(fmean[3]-gmean)^2) # 群間平方和。各水準の個々のデータの予測値を各水準の平均値としpy1-gmean, py2-gmean...の総和を水準ごとに出して足す。sum(15*(fmean-gmean)^2) でも同じ。級間平方和ともいう

sse <- var(value[dmn==lvs[1]])*(14/15)*15+var(value[dmn==lvs[2]])*(14/15)*15+var(value[dmn==lvs[3]])*(14/15)*15 # 群内平方和。各水準の"標本"分散から平均の情報を外し (*15=平方和にし) 、総和する。級内平方和ともいう

ssa+sse # sst
sst

summary(aov(value~dmn))
sst
ssa # namesが残っているが気にしない
sse

# 相関比
sqrt(ssa/sst)

ggplotはlatticeやRの組み込み作図機能よりよほど優れているらしい
 (元のpdf)

開発者: http://had.co.nz/ggplot2/
  ggplot: http://had.co.nz/ggplot/
  introduction: http://had.co.nz/ggplot/ggplot-introduction.pdf 
  One hour ggplot2 workshop: http://had.co.nz/ggplot2/resources/2007-vanderbilt.pdf
CRANのpdf: http://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf
yeroon.net: http://yeroon.net/ggplot2/
  UCLAの人がRの一部機能をWebアプリケーション化したらしい。これならマウスでインタラクティブにggplot2を使えてすごく便利…に見える
  youtubeのデモ1: http://www.youtube.com/watch?v=7XGN6OSCq6E&hd=1
  youtubeのデモ2: http://www.youtube.com/watch?v=_haIgb4nFFY&hd=1
Taglibro de H: http://ito-hi.blog.so-net.ne.jp/archive/200909-1
もうカツ丼でいいよな: http://d.hatena.ne.jp/Rion778/searchdiary?word=%2A%5Bggplot2%5D
r-bloggersで検索: http://www.r-bloggers.com/?s=ggplot2
SHO FUJIHARA'S PAGE: http://sites.google.com/site/shofujihara/r/rde-o-ekaki

  とりあえずintroductionとOne hour ggplot2 workshopからか…

Van Selst, M., & Jolicoeur, P. (1994). A solution to the effect of sample size on outlier elimination. The quarterly journal of experimental psychology, 47A, 631-650. を参考

source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1270317393")

ssz <- c(4,5,6,7,8,9,10,11,12,13,14,15,20,25,30,35,50,100)
nos <- 10000
mrv <- c(8, 6.20, 5.30, 4.80, 4.475,4.25,4.11,4.00,3.92,3.85,3.80,3.75,3.64,3.595,3.55,3.54,3.51,3.50)
nrv <- c(1.458,1.68,1.841,1.961,2.050,2.12,2.173,2.22,2.246,2.274,2.31,2.326,2.391,2.410,2.4305,2.45,2.48,2.50)

cop <- 2.5
system.time(rawres <- rawfun(ssz, nos, cop)) # 生データ平均値
system.time(medres <- medfun(ssz, nos, cop)) # 生データ中央値

system.time(lores <- lofun(ssz, nos, cop)) # 全体対数変換のみ
system.time(rores <- rofun(ssz, nos, cop)) # 全体対数変換、全体カットオフ
system.time(nrres.lo <- nrcs.lo(ssz, nos, cop)) # 全体対数変換、個別カットオフ
system.time(nrres <- nrcs(ssz, nos, cop)) # 非再帰法
cop <- 3.5
system.time(mrres <- mrcs(ssz, nos, cop)) # 修正再帰法
hybres <- (nrres+mrres)/2 # ハイブリッド法
cop <- nrv
system.time(nrreswmc <- nrcs(ssz, nos, cop)) # 非再帰基準変動
cop <- mrv
system.time(mrreswmc <- mrcs(ssz, nos, cop)) # 修正再帰基準変動

win.graph()
ylmt <- c(500, 650)
plot(ssz, rawres, type="l", xaxt="n", xlim=c(5,100), ylim=ylmt, xlab="Sample Size", yaxt="n", pch=1, col=1, lty=1) # 生データ平均値
axis(1, at=ssz, labels=as.character(ssz))
axis(2, ylim=ylmt)
#axis(2, at=seq(from=ylmt[1], to=ylmt[2], 5), labels=as.character(seq(from=ylmt[1], to=ylmt[2], 5)))
points(ssz, medres, type="l", col=1, pch=1, lty=3) # 生データ中央値
points(ssz, nrres, type="b", col=2, pch=2) # 非再帰
points(ssz, mrres, type="b", col=3, pch=3) # 修正再帰
points(ssz, hybres, type="b", col=4, pch=4) # ハイブリッド
points(ssz, nrreswmc, type="b", col=5, pch=5) # 非再帰基準変動
points(ssz, mrreswmc, type="b", col=6, pch=6) # 修正再帰基準変動
legend("topleft", lty=c(1,3,2), col=c(1,1,2), pch=c(NA, NA, 2), c("rawmean", "rawmedian", "non-recursive"), box.lty=0)
legend("top", lty=c(3:6), col=c(3:6), pch=c(3:6), c("modified recursive", "hybrid", "n.r.moving criterion", "m.r.moving criterion"), box.lty=0)

ylmt <- log(ylmt)
par(usr=c(par("usr")[1:2], ylmt[1], ylmt[2]))
axis(4, ylim=ylmt)
points(ssz, lores, type="b", col=1, ylim=ylmt, lty=2, pch=1)
points(ssz, rores, type="b", col=2, ylim=ylmt, lty=2, pch=2) # 全体カットオフ対数変換
points(ssz, nrres.lo, type="b", col=3, ylim=ylmt, lty=2, pch=3) # 修正再帰法対数変換
legend("topright", lty=c(2,2,2), col=1:3, pch=1:3, c("all.log", "log grand", "log individual"), box.lty=0)
#対数ってこう書いていいんかなー

## 出力pdf

## 表にする
tb1 <- rbind(rawres, medres, nrres, mrres, hybres, nrreswmc, mrreswmc)
tb2 <- rbind(rawres, explo=exp(lores), exprores=exp(rores), expnreslo=exp(nrres.lo), lores, rores, nrres.lo)
colnames(tb1)  <- ssz
colnames(tb2) <- ssz
options(digits=4)
tb1
tb2