myscratchpad

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)