myscratchpad

multilevelパッケージのsobel関数
http://davidakenny.net/dtt/mediate.htmのデータ。感謝

library(foreign)
dat <- data.frame(read.spss("http://davidakenny.net/dtt/morse_et_al.sav"))
names(dat) <- abbreviate(names(dat))
dmtr <- as.integer(dat$trtm=="Experimentals")

Treatment <- dmtr
DayHoused <- dat$hsng
HousingContacts <- dat$sr_h

library(psych)
print(corr.test(data.frame(Treatment, DayHoused, HousingContacts)), digits=3)

dat2 <- data.frame(Treatment, DayHoused, HousingContacts)
sapply(dat2, mean)
sd(dat2)

library(multilevel)
(sbl <- sobel(Treatment, HousingContacts, DayHoused))
(res <- rbind(c=sbl$Model.1[2,], a=sbl$Model.3[2,], b=sbl$Model.2[3,], "c'"=sbl$Model.2[2,]))

# 覚書
# 1. 予測変数が目的変数を予測する(6.558)
# 2. 予測変数が媒介変数を予測する(0.289)
# 3. 予測変数と媒介変数を投入して目的変数を予測する。このとき、媒介変数が目的変数を予測し(10.710) 、 予測変数、媒介変数を投入したとき、1. と比べて予測変数の係数が低い(3.468)


# sobel 関数の中身は、重回帰分析を3回やってるだけ
  summary(lm(DayHoused~Treatment)) # 1．　予測変数で目的変数を予測
  summary(lm(HousingContacts~Treatment)) # 2．　予測変数で媒介変数を予測
  summary(lm(DayHoused~Treatment+HousingContacts)) # 3．　予測変数と媒介変数を投入して目的変数を予測。予測変数 (Treatment) が有意じゃなくなって、媒介変数 (HousingContacts) が有意になればOK


#      Housing Contacts
#             /\
#            /  \
#           /    \
#          /      \
# 0.289*  /        \  10.710*
#        /          \
#       /            \
#      /              \
#     /                \
#Treatment-----------> Days Housed
#        3.468 (6.558*)


# Aroian, Goodmanのz値を出してみる。
# 参考: http://people.ku.edu/~preacher/sobel/sobel.htm

a <- sbl$Model.3[2,1] # HousingContacts~Treatment　モデルの、HousingContacts (媒介変数) に対するTreatment (予測変数) の回帰係数
b <- sbl$Model.2[3,1] #  DayHoused～Treatment+HousingContacts モデルの、HousingContacts (媒介変数) の回帰係数
a*b # いわゆる間接効果。これを以下の標準誤差で標準化し、z得点化して検定するのがsobel test
sga <- sbl$Model.3[2,2]
sgb <- sbl$Model.2[3,2]

sobeldn <- sqrt(a^2*sgb^2+b^2*sga^2)
aroiandn <- sqrt(a^2*sgb^2+b^2*sga^2+sga^2*sgb^2)
Goodmandn <- sqrt(a^2*sgb^2+b^2*sga^2-sga^2*sgb^2)

# Sobel
a*b/sobeldn
# Aroian
a*b/aroiandn
# Goodman
a*b/Goodmandn