myscratchpad

datf <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1258223909", header = T)
## 心理統計学の基礎、p312より
datf$dominance <- factor(datf$dominance)

options(contrasts = c("contr.sum", "contr.sum"))
covres.ind <- lm(envy ~ dominance + cov, data = datf)
covres.int <- lm(envy ~ dominance * cov, data = datf)
anova(covres.ind, covres.int)

library(car)
anova(covres.int)
Anova(covres.int, type = 3) # 教科書どおりはこれ
anova(covres.ind)
Anova(covres.ind, type = 3)

# プロットの準備
lmres.app  <- lm(envy ~ cov, data = datf, subset = (dominance == 1)) # 容姿
lmres.acbg <- lm(envy ~ cov, data = datf, subset = (dominance == 2)) # 学歴
lmres.afl  <- lm(envy ~ cov, data = datf, subset = (dominance == 3)) # 豊かさ

# プロット
par(family = "Japan1GothicBBB")
plot(1, xlim = c(0, 16), ylim = c(0, 10), xlab = "cov.", ylab = "envy")
points(envy ~ cov, pch = 1, data = datf, subset = (dominance == 1))
abline(lmres.app, lty = 1)
points(envy ~ cov, pch = 17, data = datf, subset = (dominance == 2))
abline(lmres.acbg, lty = 4)
points(envy ~ cov, pch = 4, data = datf, subset = (dominance == 3))
abline(lmres.afl, lty = 2)
legend("topleft", c("容姿", "学歴", "豊かさ"), pch = c(1, 17, 4), lty = c(1, 4, 2), box.lty = 0)