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)