myscratchpad

dat <- read.table("http://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)

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

dat <- read.table("http://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)

データフレームdfの1列目がidv, 2列目がmod (以上独立変数)
4列目がdv (従属変数) として
以下は全て同じ結果を返す。列番号指定のときはas.matrixをつけないとダメ
ただし、個々の独立変数のF値等を見たいときは一番上

lm(dv ~ idv + mod, data = df)
lm(df[, 4] ~ as.matrix(df[c(1, 2)]))
lm(df[, 4] ~ as.matrix(df[, 1:2]))
lm(df[, 4] ~ as.matrix(df[1:2]))

結論：面倒でも独立変数は全部変数名と"+"でつなごう

大本。idv, mod, dv、全て連続変量
result <- lm(dv ~ idv*mod, df)
# これでidv, modの主効果、idv*modの交互作用項を検定している

#　結果
##　残差、偏回帰係数、その標準誤差、t値、p値、決定係数、調整済み決定係数、回帰式のF値
summary(result)
##　QuantPsycパッケージのlm.beta関数で標準化偏回帰係数を出してくれるようだ
lm.beta(result)


##　各独立変数の平方和とF値
anova(result)

##　VIFとトレランス
x <- df$idv*df$mod # 交互作用項をつくっておく
ndf <- cbind("idv"=c(df$idv), "mod"=c(df$mod), x) # 独立変数だけのデータフレームをつくっておく
diag(solve(cor(ndf))) # VIF計算
1/diag(solve(cor(ndf))) # トレランス

##　モデルの比較
result2 <- lm(dv ~ idv+mod, df) # 交互作用なしモデル
anova(result2, result)


simple slope (単純傾斜)
#　QuantPsycパッケージを利用する
## この関数の結果はSPSSによるものと一致しない。さらにanovaの分散分析表も一致しない。たぶん平方和の問題じゃないかと思う。別に一致する必要はないが...

library(QuantPsyc)
mresult <- moderate.lm(idv, mod, dv, df, mc = TRUE)
##　mcはデータを中心化しているかどうか
##　mresultは上記のresultと一致する
sresult <- sim.slopes(mresult, df$mod, zsd=1, mcz=FALSE)
##　zsdはsimple slopeをつくるとき、SDをいくつとるか。普通は1SD (だったと思う)
sresult # 結果の出力。各群の切片、偏回帰係数、標準誤差、信頼区間が出力される
graph.mod(sresult, idv, dv, df, title = "Sslope", xlab = "Centered X", ylab = "Raw Y", ylimit = 1.5)