×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
ここがわかりやすかったので、覚書。というかコピペ。官舎
dv <- c(8,6,7,7,6,4,5,5,6,3,5,3,3,6,2,3,4,2,2,3,6,7,5,4,6)
grp <- gl(5,5)
aovres <- aov(dv~grp)
summary(aovres)
model.tables(aovres)
# デフォルトの対比
options()$contrasts
contrasts(grp)
# ヘルマート対比で検定
coef.hel <- contr.helmert(5)
coef.hel
contrasts(grp) <- coef.hel
aovres.h <- aov(dv~grp)
summary(aovres.h, split=list(grp=c(1,2,3,4)))
# grp: C1 1 12.100 12.100 9.0299 0.0069971 ** grp1 vs. grp2
# grp: C2 1 12.033 12.033 8.9801 0.0071289 ** grp1&2 vs. grp3
# grp: C3 1 19.267 19.267 14.3781 0.0011441 ** grp1&2&3 vs. grp4
# grp: C4 1 4.840 4.840 3.6119 0.0718780 . grp1&2&3&4 vs. grp5
# 多項式対比で検定
coef.ply <- contr.poly(5)
coef.ply
contrasts(grp) <- coef.ply
aovres.p <- aov(dv~grp)
summary(aovres.p, split=list(grp=c(1,2,3,4)))
# 自分で対比を設定する
coef.s <- cbind(c(1,-1,1,0,-1), c(-1,-1,-1,4,-1))
coef.s
sum(coef.s)
contrasts(grp) <- coef.s
aovres.s <- aov(dv~grp)
summary(aovres.s, intercept=T, split=list(grp=c(1,2)))
# grp: C1 1 0.20 0.20 0.1493 0.7033289 # grp1&3 vs. grp2&5
# grp: C2 1 23.04 23.04 17.1940 0.0004994 *** # grp4 vs. grp1&2&3&5
dv <- c(8,6,7,7,6,4,5,5,6,3,5,3,3,6,2,3,4,2,2,3,6,7,5,4,6)
grp <- gl(5,5)
aovres <- aov(dv~grp)
summary(aovres)
model.tables(aovres)
# デフォルトの対比
options()$contrasts
contrasts(grp)
# ヘルマート対比で検定
coef.hel <- contr.helmert(5)
coef.hel
contrasts(grp) <- coef.hel
aovres.h <- aov(dv~grp)
summary(aovres.h, split=list(grp=c(1,2,3,4)))
# grp: C1 1 12.100 12.100 9.0299 0.0069971 ** grp1 vs. grp2
# grp: C2 1 12.033 12.033 8.9801 0.0071289 ** grp1&2 vs. grp3
# grp: C3 1 19.267 19.267 14.3781 0.0011441 ** grp1&2&3 vs. grp4
# grp: C4 1 4.840 4.840 3.6119 0.0718780 . grp1&2&3&4 vs. grp5
# 多項式対比で検定
coef.ply <- contr.poly(5)
coef.ply
contrasts(grp) <- coef.ply
aovres.p <- aov(dv~grp)
summary(aovres.p, split=list(grp=c(1,2,3,4)))
# 自分で対比を設定する
coef.s <- cbind(c(1,-1,1,0,-1), c(-1,-1,-1,4,-1))
coef.s
sum(coef.s)
contrasts(grp) <- coef.s
aovres.s <- aov(dv~grp)
summary(aovres.s, intercept=T, split=list(grp=c(1,2)))
# grp: C1 1 0.20 0.20 0.1493 0.7033289 # grp1&3 vs. grp2&5
# grp: C2 1 23.04 23.04 17.1940 0.0004994 *** # grp4 vs. grp1&2&3&5
PR
dat <- data.frame(
ps=paste("p", 1:27, sep=""),
trt=gl(3, 9, labels=c("a1", "a2", "a3")), # "notreat", "spaced", "massed"
grp=rep(gl(3, 3, labels=c("b1", "b2", "b3")),3), #"washers", "checkers", "seekers"
value=c(5,4,5,4,5,3,3,5,4,5,6,3,4,6,3,12,10,13,10,12,16,11,10,12,13,12,10)
)
agdat <- aggregate(dat[4], list(dat[,2], dat[,3]), mean)
lmres <- lm(value~trt*grp, dat)
library(car)
Anovares <- Anova(lmres)
Anovares
smry <- summary(lmres)
mse <- smry$sigma^2 # 残差の平均平方和
by(dat, dat$trt, function(x) Anova(lm(value~grp, data=x), error=lmres))
by(dat, dat$grp, function(x) Anova(lm(value~trt, data=x), error=lmres))
vnm <- paste(agdat[,1], agdat[,2], sep="")
for (i in 1:9) assign(vnm[i], agdat[,3][i])
cv <- c(1,-1,1,-1) # 対比係数
n <- 3 # 各群の人数
nm <- (sum(cv*c(a1b1, a2b1, a1b3, a2b3)))^2 # (a1b1-a2b1+a1b3-a2b3) b1水準のa1, a2の差とb3水準のa1, a2の差の差をとって二乗しただけ
dn <- (sum((cv^2)*mse))/n
(fv <- nm/dn)
(tv <- sqrt(fv))
1-pf(fv, df1=1, df2=18) # 分子の自由度は1。分母の自由度は分散分析の誤差自由度
# F(1, 18) = 18.03, p < .001. b1群とb3群はa1 -> a2への反応に違いがあるということ (原文: The significant F indicates that Washers (b1) and Seekers (b3) respond differently at the two levels of Factor A (a1 and a2) .
(1-pt(tv, 18))*2
参考:
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc
http://www.sprweb.org/articles/Keselman98.pdf
ps=paste("p", 1:27, sep=""),
trt=gl(3, 9, labels=c("a1", "a2", "a3")), # "notreat", "spaced", "massed"
grp=rep(gl(3, 3, labels=c("b1", "b2", "b3")),3), #"washers", "checkers", "seekers"
value=c(5,4,5,4,5,3,3,5,4,5,6,3,4,6,3,12,10,13,10,12,16,11,10,12,13,12,10)
)
agdat <- aggregate(dat[4], list(dat[,2], dat[,3]), mean)
lmres <- lm(value~trt*grp, dat)
library(car)
Anovares <- Anova(lmres)
Anovares
smry <- summary(lmres)
mse <- smry$sigma^2 # 残差の平均平方和
by(dat, dat$trt, function(x) Anova(lm(value~grp, data=x), error=lmres))
by(dat, dat$grp, function(x) Anova(lm(value~trt, data=x), error=lmres))
vnm <- paste(agdat[,1], agdat[,2], sep="")
for (i in 1:9) assign(vnm[i], agdat[,3][i])
cv <- c(1,-1,1,-1) # 対比係数
n <- 3 # 各群の人数
nm <- (sum(cv*c(a1b1, a2b1, a1b3, a2b3)))^2 # (a1b1-a2b1+a1b3-a2b3) b1水準のa1, a2の差とb3水準のa1, a2の差の差をとって二乗しただけ
dn <- (sum((cv^2)*mse))/n
(fv <- nm/dn)
(tv <- sqrt(fv))
1-pf(fv, df1=1, df2=18) # 分子の自由度は1。分母の自由度は分散分析の誤差自由度
# F(1, 18) = 18.03, p < .001. b1群とb3群はa1 -> a2への反応に違いがあるということ (原文: The significant F indicates that Washers (b1) and Seekers (b3) respond differently at the two levels of Factor A (a1 and a2) .
(1-pt(tv, 18))*2
参考:
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc
http://www.sprweb.org/articles/Keselman98.pdf
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
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
この記事のとは違う
対数変換してから非再帰法+基準値置換でやる
Ratcliff, R. (1993). Methods for dealing with reaction time outliers. Psychological Bulletin, 114, 510-532.
っぽいかんじ
nn <- 50
rt <- round(c((rnorm(nn, mean=300, sd=20)+(300*rexp(nn)))))
## 関数定義。反応時間のベクトルrt, 基準cop (デフォルト2.5SD) を与える。ifm=Tにすると他の情報も出す。
elrt <- function(rt, cop=2.5, ifm=FALSE) {
rt <- na.omit(as.vector(rt))
rtl <- log(rt)
cop <- cop
meanv <- mean(rtl)
sdv <- sd(rtl)
cv.l <- meanv-(cop*sdv)
cv.u <- meanv+(cop*sdv)
rtl2 <- rtl
tf.l <- rtl<cv.l
tf.u <- rtl>cv.u
rtl2[tf.l] <- cv.l # 外れ値をカットオフ値で置換
rtl2[tf.u] <- cv.u
nart <- rt
nart[tf.l] <- NA
nart[tf.u] <- NA
rplrt <- rt
rplrt[tf.l] <- exp(cv.l) # 生rtの外れ値に対しカットオフ値を指数変換して置換
rplrt[tf.u] <- exp(cv.u)
logrt <- rtl2
reslist <- list()
reslist$logrt <- logrt
reslist$nart <- nart
reslist$rplrt <- rplrt
reslist$outliers.l <- rt[tf.l]
reslist$outliers.u <- rt[tf.u]
reslist$cop <- paste(cop, "SD", sep="")
reslist$cv <- c(cv.l, sum(tf.l), cv.u, sum(tf.u))
names(reslist$cv) <- c("lower", "nofl", "upper", "nofu")
ifelse(ifm==FALSE, return(rtl2), return(reslist))
}
(logrt <- elrt(rt)) # 対数変換・カットオフ置換のみ
(elres <- elrt(rt, ifm=T)) # logrtに加え、外れ値をNAにしたもの、指数変換で置き換えたもの、カットオフ値とその個数
対数変換してから非再帰法+基準値置換でやる
Ratcliff, R. (1993). Methods for dealing with reaction time outliers. Psychological Bulletin, 114, 510-532.
っぽいかんじ
nn <- 50
rt <- round(c((rnorm(nn, mean=300, sd=20)+(300*rexp(nn)))))
## 関数定義。反応時間のベクトルrt, 基準cop (デフォルト2.5SD) を与える。ifm=Tにすると他の情報も出す。
elrt <- function(rt, cop=2.5, ifm=FALSE) {
rt <- na.omit(as.vector(rt))
rtl <- log(rt)
cop <- cop
meanv <- mean(rtl)
sdv <- sd(rtl)
cv.l <- meanv-(cop*sdv)
cv.u <- meanv+(cop*sdv)
rtl2 <- rtl
tf.l <- rtl<cv.l
tf.u <- rtl>cv.u
rtl2[tf.l] <- cv.l # 外れ値をカットオフ値で置換
rtl2[tf.u] <- cv.u
nart <- rt
nart[tf.l] <- NA
nart[tf.u] <- NA
rplrt <- rt
rplrt[tf.l] <- exp(cv.l) # 生rtの外れ値に対しカットオフ値を指数変換して置換
rplrt[tf.u] <- exp(cv.u)
logrt <- rtl2
reslist <- list()
reslist$logrt <- logrt
reslist$nart <- nart
reslist$rplrt <- rplrt
reslist$outliers.l <- rt[tf.l]
reslist$outliers.u <- rt[tf.u]
reslist$cop <- paste(cop, "SD", sep="")
reslist$cv <- c(cv.l, sum(tf.l), cv.u, sum(tf.u))
names(reslist$cv) <- c("lower", "nofl", "upper", "nofu")
ifelse(ifm==FALSE, return(rtl2), return(reslist))
}
(logrt <- elrt(rt)) # 対数変換・カットオフ置換のみ
(elres <- elrt(rt, ifm=T)) # logrtに加え、外れ値をNAにしたもの、指数変換で置き換えたもの、カットオフ値とその個数
## 基本的に、解析等が終了しその結果をhtmlファイルにまとめておくためのもの
dat <- data.frame(s=paste("s", 1:10), val=rnorm(10))
library(psych)
library(R2HTML)
HTMLStart(outdir=getwd(), file="report",extension="html", echo=FALSE, HTMLframe=TRUE)
HTML.title("My Report", HR=1)
HTML.title("Description of my data", HR=3)
describe(dat)
HTMLhr()
HTML.title("X Y Scatter Plot", HR=2)
plot(dat$val)
HTMLplot()
HTMLStop()
Quick-R: http://www.statmethods.net/interface/output.html
CRAN: http://cran.r-project.org/web/packages/R2HTML/index.html
開発者 (?) : http://www.stat.ucl.ac.be/R2HTML/
dat <- data.frame(s=paste("s", 1:10), val=rnorm(10))
library(psych)
library(R2HTML)
HTMLStart(outdir=getwd(), file="report",extension="html", echo=FALSE, HTMLframe=TRUE)
HTML.title("My Report", HR=1)
HTML.title("Description of my data", HR=3)
describe(dat)
HTMLhr()
HTML.title("X Y Scatter Plot", HR=2)
plot(dat$val)
HTMLplot()
HTMLStop()
Quick-R: http://www.statmethods.net/interface/output.html
CRAN: http://cran.r-project.org/web/packages/R2HTML/index.html
開発者 (?) : http://www.stat.ucl.ac.be/R2HTML/
