myscratchpad

dat <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1256390904",header=T)

dat1 <- dat
dat1$idv <- dat1$idv - mean(dat1$idv)
dat1$mod <- dat1$mod - mean(dat1$mod)

modsd <- sd(dat1$mod)
modab <- dat1$mod - modsd
modbl <- dat1$mod - (-1*modsd)
dat2 <- cbind(dat1, modab, modbl)

lmres.ind <- lm(dv ~ idv + mod, data = dat2)
lmres.whl <- lm(dv ~ idv * mod, data = dat2)
anova(lmres.ind, lmres.whl)
summary(lmres.whl)
anova(lmres.whl)

# 単純傾斜
lmres.modab <- lm(dv ~ idv * modab, data = dat2)
lmres.modbl <- lm(dv ~ idv * modbl, data = dat2)

summary(lmres.modab) # 高群
summary(lmres.modbl) # 低群

# プロット
x <- dat2$idv
y <- dat2$dv
above.line <- summary(lmres.modab)$coef[1:2,1] # 高群の切片と傾き
below.line <- summary(lmres.modbl)$coef[1:2,1] # 低群の切片と傾き
plot(x, y)
abline(above.line, col = "red")
abline(below.line, col = "blue")
legend(locator(1), legend = c("mod.above", "mod.below"), lwd = 0, box.lty = 0, col = c("red", "blue"))


# 標準化偏回帰係数とか
co.whl <- summary(lmres.whl)$coefficients[,1]
co.ab <- summary(lmres.modab)$coefficients[,1]
co.bl <- summary(lmres.modbl)$coefficients[,1]

sd.idv <- sd(dat2$idv)
sd.dv <- sd(dat2$dv)

co.whl*(sd.idv/sd.dv)
co.ab*(sd.idv/sd.dv)
co.bl*(sd.idv/sd.dv)

# 回帰係数の標準誤差と信頼区間
cose.ab <- summary(lmres.modab)$coefficients[,2]
cose.bl <- summary(lmres.modbl)$coefficients[,2]

## 標準誤差
cose.ab
cose.bl

## 信頼区間
cose.df <- summary(lmres.modbl)$fstatistic[3]
t95 <- abs(qt(0.975, cose.df))

co.ab - t95*cose.ab # 下側
co.ab + t95*cose.ab # 上側

co.bl - t95*cose.bl # 下側
co.bl + t95*cose.bl # 上側

# 独立変数を全て標準化してから実行すると回帰係数とかは全部標準化...になる
dat1$idv <- scale(dat1$idv)
dat1$mod <- scale(dat1$mod)
dat1$dv <- scale(dat1$dv)

dat <- data.frame(
s=factor(c("s1", "s2", "s3", "s4", "s5", "s6", "s7")),
a=factor(c(rep("a1", 4), rep("a2", 3))),
b1c1=c(2,6,5,7,1,2,5), b1c2=c(5,7,9,9,2,1,4), b1c3=c(9,10,13,14,1,5,3),
b2c1=c(6,6,8,10,5,3,4), b2c2=c(3,6,5,8,5,6,6), b2c3=c(6,7,5,6,5,5,9)
)

options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(b1c1, b1c2, b1c3, b2c1, b2c2, b2c3) ~ a, data = dat)
library(car)
bfact <- factor(rep(c("b1", "b2"), c(3, 3)))
cfact <- factor(rep(c("c1", "c2", "c3"), 2))
idata <- data.frame(bfact, cfact)
Anovares <- Anova(lmres, idata = idata , idesign = ~ bfact*cfact, type = 3)
summary(Anovares)

# aovと比較
library(reshape)
dat.aov <- melt(dat, id=c("s", "a"))
b <- c(rep(c("b1", "b2"), each=21))
c <- c(rep(c("c1", "c2", "c3"), each=7), rep(c("c1", "c2", "c3"), each=7))
dat.aov <- data.frame(dat.aov[c(1, 2, 4)], b, c)
names(dat.aov) <- c("s", "a", "result", "b", "c")
summary(aov(result ~ a * b * c + Error(s:a + b:s:a + c:s:a + b:c:s:a), dat.aov))

carパッケージのAnova関数で混合要因の分散分析
データはデータ解析テクニカルブック
SPFp.qデザイン　データ数が同じ(p.107)

dat <- data.frame(
s=factor(c("s1","s2","s3","s4","s5")),
a=factor(c(rep("a1",5), rep("a2",5))),
b1=c(3,3,1,3,5,3,5,2,4,6), b2=c(4,3,4,5,7,2,6,3,6,4), b3=c(6,6,6,4,8,3,2,3,6,5), b4=c(5,7,8,7,9,2,3,3,4,6)
)
dat
aggregate(dat[3:6], list(dat[,2]), sum)

# メインの分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(b1, b2, b3, b4) ~ a , contrasts = list(a = contr.sum), data = dat)

bfact <- as.factor(c("b1", "b2", "b3", "b4"))
idat <- data.frame(bfact)
library(car)
Anovares <- Anova(lmres, idata = idat, idesign = ~ bfact, type = 3)
summary(Anovares)

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

# 作図のため、平均値を出して行列にする
agdat <- aggregate(dat[3:6], list(dat[,2]), mean)
mdat <- as.matrix(agdat[2:5])
# 折れ線で作図
matplot(mdat, type="b", pch=1:4, xaxt="n")
axis(1, at=1:2, labels=c("a1",  "a2"))
legend("topright", legend=c("b1", "b2", "b3", "b4"), lty=1:4, pch=1:4, col=1:4)
# 棒グラフで作図
rownames(mdat) <- c("a1", "a2")
palette(grey.colors(nrow(t(mdat)))) # a要因を横軸にしたいのでtで転置
barplot(t(mdat), beside=T)
legend("topright", rownames(t(mdat)), fill = 1:4)


## 事後検定　テクニカルブックとは違って水準別
## b要因 (被験者内要因) の水準ごと
lmab1 <- lm(b1 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab1 <- Anova(lmab1, type = 3)
lmab2 <- lm(b2 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab2 <- Anova(lmab2, type = 3)
lmab3 <- lm(b3 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab3 <- Anova(lmab3, type = 3)
lmab4 <- lm(b4 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab4 <- Anova(lmab4, type = 3)

Anv.ab1
Anv.ab2
Anv.ab3
Anv.ab4

## 事後検定
## a要因 (被験者間要因) の水準ごと、aの１
levels(dat$a)
data1 <- subset(dat, a=="a1")
data1
lmba1 <- lm(cbind(b1, b2, b3, b4) ~ 1, contrasts=contr.sum, data = data1)
Anv.a1 <- Anova(lmba1, idata = idat, idesign = ~ bfact2, type = 3)
summary(Anv.a1)

## a要因 (被験者間要因) の水準ごと、aの2
data2 <- subset(dat, a=="a2")
data2
lmba2 <- lm(cbind(b1, b2, b3, b4) ~ 1, contrasts=contr.sum, data = data2)
Anv.a2 <- Anova(lmba2, idata = idat, idesign = ~ bfact2, type = 3)
summary(Anv.a2)

##多重比較 単純効果のあったa1で
#holm法
datpa <- stack(data1[3:6])
datpa$s <- rep(data1[,2], 4)
pttres <- pairwise.t.test(datpa$values, datpa$ind, p.adjust.method="holm", paired=TRUE)
pttres

#holm法2
p12 <- t.test(data1$b1, data1$b2, paired = T)
p13 <- t.test(data1$b1, data1$b3, paired = T)
p14 <- t.test(data1$b1, data1$b4, paired = T)
p23 <- t.test(data1$b2, data1$b3, paired = T)
p24 <- t.test(data1$b2, data1$b4, paired = T)
p34 <- t.test(data1$b3, data1$b4, paired = T)
p.adjust(c(p12$p.value,p13$p.value,p14$p.value,p23$p.value,p24$p.value,p34$p.value), method = "holm")

Rで基礎的な多変量解析を行う、というフリー文書へのリンク
http://www.plymouth.ac.uk/pages/dynamic.asp?page=staffdetails&id=phewson#
このなかの"Multivariate Statistics with R"
著者のPaul HewsonはCRAN Task View: Multivariate StatisticsのMaintainer

関連リンク: http://www.opentextbook.org/

Aiken & West (1991) のp. 130-131
連続変量とカテゴリカル変量の重回帰分析で、単純傾斜を検定する。
カテゴリカル変量のダミーコードの割り当て方を変えて分析するだけ

dat <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1255449130", header = T)
dat$college <- factor(dat$college,levels=c("LA","E","BUS"))
dat$cGPA <- dat$GPA - mean(dat$GPA) # 連続変量の中心化
res.idv <- lm(salary ~ college + cGPA, dat) # 交互作用なしモデル
res.int <- lm(salary ~ college * cGPA, dat) # 交互作用ありモデル
anova(res.idv, res.int) # 交互作用項投入に効果があるか検定

## 交互作用ありモデルの概要
summary(res.int)
anova(res.int)

## 事後検定
# LA群の単純傾斜
datla <- dat
datla$college <- factor(datla$college, levels = c("LA","E","BUS"))
res.la <- lm(salary ~ college * cGPA, datla)
summary(res.la)

# E群の単純傾斜
date <- dat
date$college <- factor(date$college, levels = c("E","LA","BUS"))
res.e <- lm(salary ~ college * cGPA, date)
summary(res.e)

# BUS群の単純傾斜
datbus <- dat
datbus$college <- factor(datbus$college, levels = c("BUS","LA","E"))
res.bus <- lm(salary ~ college * cGPA, datbus)
summary(res.bus)

# 各々の分析のcGPAの係数、検定等が単純傾斜検定に相当する
summary(res.la)$coef["cGPA", ]
summary(res.e)$coef["cGPA", ]
summary(res.bus)$coef["cGPA", ]

plot(dat$GPA, dat$salary, pch = 1:3, type = "p", xlab = "GPA", ylab = "Salary", xlim = c(1.5, 4.5), ylim = c(16000, 30000), col = 1:3)
abline(lm(dat$salary[dat$college == "LA"] ~ dat$GPA[dat$college == "LA"]), lty = 1, col = 1, lwd = 3)
abline(lm(dat$salary[dat$college == "E"] ~ dat$GPA[dat$college == "E"]), lty = 2, col = 2, lwd = 3)
abline(lm(dat$salary[dat$college == "BUS"] ~ dat$GPA[dat$college == "BUS"]), lty = 3, col = 3, lwd = 3)
legend("bottomright", c("Liberal Arts", "Engineering", "Business"), pch = 1:3, lty = 1:3, col = 1:3, lwd = 1.5, box.lty = 0)