myscratchpad

# 交互作用重回帰の別の例。以下参考。感謝
# http://wilhelmhofmann.info/wordpress/?page_id=120
library(foreign)
dat <- read.spss("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1286213208", to.data.frame=T, max.value.labels=99)
head(dat)

# 記述統計
library(psych)
describe(dat)
x <- dat[c(6,2,3,4)]

attach(dat)
heightc <- scale(height, scale=F)
agec <- scale(age, scale=F)
summary(lm(lifesat~heightc*agec))

# simple slope
agech <- agec-sd(agec)
hmod <- lm(lifesat~heightc*agech)
summary(hmod)
agecl <- agec+sd(agec)
lmod <- lm(lifesat~heightc*agecl)
summary(lmod)

# プロット
x <- agec
y <- lifesat
plot(x,y, col=0)
abline(hmod, lty=1)
abline(lmod, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)


# 正しい標準化値。交互作用項を作成してから標準化するのは誤り
hsc <- scale(height)
asc <- scale(age)
lsc <- scale(lifesat)
summary(lm(lsc~hsc*asc))

# simple slopeの標準化値
asch <- asc-1
hmods <- lm(lsc~hsc*asch)
summary(hmods) # 高群
agecl <- asc+1
lmods <- lm(lsc~hsc*agecl)
summary(lmods) # 低群

# 標準化値のプロット
win.graph()
x <- asc
y <- lsc
plot(x,y, col=0)
abline(hmods, lty=1)
abline(lmods, lty=2)
legend("topright", legend=c("ageabove", "agebelow"), lty=1:2)

# そのうちやろう
  # 統制変数


# 連続 x カテゴリカル
summary(lm(lifesat~heightc*gender))
# simple slope
 mm <- factor(gender, levels=c("male", "female")) # 男性をリファレンスカテゴリにする。ダミーコードではmale=0, female=1になる。
 ff <- factor(gender, levels=c("female", "male"))
 summary(lm(lifesat~heightc*mm)) # 男性
 summary(lm(lifesat~heightc*ff)) # 女性

# 標準化値
# The unstandardized solution from the output is the correct solution (Friedrich, 1982)!
summary(lm(lsc~hsc*gender))
 # .507 = estimated difference in regression weights between groups
# simple slope. カテゴリカル変数は標準化しないこと
 summary(lm(lsc~hsc*mm)) # 男性
 summary(lm(lsc~hsc*ff)) # 女性

# プロット
plot(heightc, lifesat)
abline(lm(lifesat~heightc*mm), lty=1)
abline(lm(lifesat~heightc*ff), lty=2)
legend("bottomright", legend=c("male", "female"), lty=1:2)

# x軸とライン分けを変えたプロット
 lmres <- lm(lifesat~heightc*gender) # 全体のモデル
 coef(lmres)
 inc <- coef(lmres)[1] # 全体の切片
 cfidv <- coef(lmres)[2] # 独立変数の回帰係数
 cfmod <- coef(lmres)[3] # 調整変数の回帰係数
 cfint <- coef(lmres)[4] # 独立変数の回帰係数
 sdidv <- sd(heightc)
 # 全体の重回帰分析でfemaleをリファレンスカテゴリにしている場合
 fh <- inc+sdidv*cfidv # 切片＋（平均＋1SD)*回帰係数
 fl <- inc-sdidv*cfidv # 切片＋（平均-1SD)*回帰係数
 mh <- inc+sdidv*cfidv + cfmod+((sdidv)*cfint) # 切片＋（平均＋1SD)*回帰係数 + 調整変数の回帰係数＋((平均＋1SD)*交互作用項の回帰係数)
 ml <- inc-sdidv*cfidv + cfmod-((sdidv)*cfint) # 切片＋（平均-1SD)*回帰係数 + 調整変数の回帰係数＋((平均-1SD)*交互作用項の回帰係数)

 mat <- matrix(c(fl, ml, fh, mh), nr=2, dimnames=list(c("female", "male"), c("heightlow", "heighthigh")))
 mat
 matplot(mat, type="l", xaxt="n")
 axis(1,
      at = 1:2,
      labels=rownames(mat),
     )
 legend("bottomleft", legend=colnames(mat), lty=1:2, col=1:2)
# 棒グラフにしてみる
win.graph(); barplot(t(mat), beside=T, legend=c("heightlow", "heighthigh"), ylab="lifesat")