myscratchpad

##　心理統計学の基礎p226 より。感謝

x1 <- c(12,12,7,17,14,9,10,13,15,12,12,15,11,14,17,17,16,15,15,10,12,9,12,12,19,11,14,15,15,15,16,15,12,10,11,12,15,13,15,12,12,12,13,17,13,11,14,16,12,12)
x2 <- c(2,2,2,3,2,2,3,3,3,1,3,3,2,2,4,2,4,3,4,2,2,1,2,2,4,2,3,2,3,3,2,3,2,2,3,1,2,3,2,2,2,3,3,3,2,3,2,4,2,2)
y <- c(6,11,11,13,13,10,10,15,11,11,16,14,10,13,12,15,16,14,14,8,13,12,12,11,16,9,12,13,13,14,12,15,8,12,11,6,12,15,9,13,9,11,14,12,13,9,11,14,16,8)


dat <- data.frame(x1, x2, y)
library(psych)
describe(dat)[2:4]
cor(dat)

r12 <- cor(x1,x2)
ry1 <- cor(x1,y)
ry2 <- cor(x2,y)

# 部分相関 (semipartial correlation, part correlation)
## x2からx1の影響を除き、そのx2とyとの相関を求める
ry21 <- (ry2 -ry1*r12)/sqrt(1-r12^2)
ry21
## x1からx2の影響を除き、そのx1とyとの相関を求める
ry12 <- (ry1 -ry2*r12)/sqrt(1-r12^2)
ry12
## 関数にしておこう
spcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、そのxとyの相関\n")
  return((ryx -rye*rxe)/sqrt(1-rxe^2))
}
spcor(x2, y, x1)
ry21

# 偏相関 (partial correlation)
## x2からx1の影響を除き、yからもx1の影響を除き、そのうえでx2とyとの相関をとる
pry21 <- (ry2-ry1*r12)/((sqrt(1-ry1^2))*(sqrt(1-r12^2)))
pry21
## 部分相関から偏相関を求める
ry21
pry21*(sqrt(1-ry1^2))
## 偏相関から部分相関を求める
pry21
ry21/sqrt(1-ry1^2)
## 関数にしておこう
pcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、yからelの影響を除き、そのxとyの相関\n")
  return((ryx-rye*rxe)/((sqrt(1-rye^2))*(sqrt(1-rxe^2))))
}
pcor(x2, y, x1)
pry21

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

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 ~ cGPA+college, dat) # 交互作用なしモデル
res.int <- lm(salary ~ cGPA * college, dat) # 交互作用ありモデル
summary(res.idv)
summary(res.int)


## 標準化偏回帰係数を求める。ダミー変数をつくり、交互作用項を計算し、全部標準化して同じ分析をするだけ
# ダミー変数の作成。青木先生のコードを参考にさせていただいた。感謝。
lvs <- levels(dat$college) # 水準名を確認し、オブジェクトに格納する
la.dm <- as.integer(dat$college==lvs[1])
e.dm <- as.integer(dat$college==lvs[2])
bus.dm <- as.integer(dat$college==lvs[3])

dat2 <- data.frame(la.dm, e.dm, bus.dm, dat[3:5])
dat3 <- data.frame(scale(dat2))　# 変数の標準化

## 交互作用なしモデルの標準化偏回帰係数
res.idv2 <- lm(scale(salary)~scale(GPA)+e.dm+bus.dm+la.dm, dat2)
res.idv3 <- lm(salary~GPA+e.dm+bus.dm+la.dm, dat3)
summary(res.idv2) # ダミー変数以外を標準化
summary(res.idv3) # 全て標準化
## 交互作用ありモデルの標準化偏回帰係数
# ダミー変数以外を標準化する。ただし、この場合ダミー変数以外の係数が1を超えることがある
summary(lm(scale(salary)~scale(cGPA)*(e.dm+bus.dm), dat2))
summary(lm(salary~GPA*(e.dm+bus.dm), dat3))# dat3は全て標準化したもの


これは間違いらしい
# 交互作用項を先につくっておく。
gpa.la <- la.dm*dat$cGPA # ここは"cGPA"にしておく。GPAだと間違い(?)
gpa.e <- e.dm*dat$cGPA
gpa.bus <- bus.dm*dat$cGPA
dat3.1 <- data.frame(dat2, gpa.la, gpa.e, gpa.bus)
dat4 <- data.frame(scale(dat3.1))
res.int2 <- lm(salary~GPA+e.dm+bus.dm+gpa.e+gpa.bus, dat4)
summary(res.int2)

datf <- read.table("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1258223909", header = T)
## 心理統計学の基礎、p312より
datf$dominance <- factor(datf$dominance)

options(contrasts = c("contr.sum", "contr.sum"))
covres.ind <- lm(envy ~ dominance + cov, data = datf)
covres.int <- lm(envy ~ dominance * cov, data = datf)
anova(covres.ind, covres.int)

library(car)
anova(covres.int)
Anova(covres.int, type = 3) # 教科書どおりはこれ
anova(covres.ind)
Anova(covres.ind, type = 3)

# プロットの準備
lmres.app  <- lm(envy ~ cov, data = datf, subset = (dominance == 1)) # 容姿
lmres.acbg <- lm(envy ~ cov, data = datf, subset = (dominance == 2)) # 学歴
lmres.afl  <- lm(envy ~ cov, data = datf, subset = (dominance == 3)) # 豊かさ

# プロット
par(family = "Japan1GothicBBB")
plot(1, xlim = c(0, 16), ylim = c(0, 10), xlab = "cov.", ylab = "envy")
points(envy ~ cov, pch = 1, data = datf, subset = (dominance == 1))
abline(lmres.app, lty = 1)
points(envy ~ cov, pch = 17, data = datf, subset = (dominance == 2))
abline(lmres.acbg, lty = 4)
points(envy ~ cov, pch = 4, data = datf, subset = (dominance == 3))
abline(lmres.afl, lty = 2)
legend("topleft", c("容姿", "学歴", "豊かさ"), pch = c(1, 17, 4), lty = c(1, 4, 2), box.lty = 0)

dat <- read.table("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1256404631",header=T)

dat1 <- dat

mod1sd <- sd(dat1$mod1)
mod1ab <- dat1$mod1 - mod1sd
mod1bl <- dat1$mod1 - (-1*mod1sd)
mod2sd <- sd(dat1$mod2)
mod2ab <- dat1$mod2 - mod2sd
mod2bl <- dat1$mod2 - (-1*mod2sd)
dat2 <- cbind(dat1, mod1ab, mod1bl, mod2ab, mod2bl)

lmres.whl <- lm(dv ~ idv*mod1*mod2, data = dat2)
summary(lmres.whl)

lmres.1bl.2bl <- lm(dv ~ idv*mod1bl*mod2bl, data = dat2)
lmres.1bl.2ab <- lm(dv ~ idv*mod1bl*mod2ab, data = dat2)
lmres.1ab.2bl <- lm(dv ~ idv*mod1ab*mod2bl, data = dat2)
lmres.1ab.2ab <- lm(dv ~ idv*mod1ab*mod2ab, data = dat2)
summary(lmres.1bl.2bl) #LL
summary(lmres.1bl.2ab) #LH
summary(lmres.1ab.2bl) #HL
summary(lmres.1ab.2ab) #HH

# 標準化値を調べる
dat3 <- dat1
dat3$idv <- scale(dat1$idv)
dat3$mod1 <- scale(dat1$mod1)
dat3$mod2 <- scale(dat1$mod2)
dat3$dv <- scale(dat1$dv)
smod1sd <- sd(dat3$mod1)
smod1ab <- dat3$mod1 - smod1sd
smod1bl <- dat3$mod1 - (-1*smod1sd)
smod2sd <- sd(dat3$mod2)
smod2ab <- dat3$mod2 - smod2sd
smod2bl <- dat3$mod2 - (-1*smod2sd)
dat4 <- cbind(dat3, smod1ab, smod1bl, smod2ab, smod2bl)

slmres.whl <- lm(dv ~ idv*mod1*mod2, data = dat4)
summary(slmres.whl)
slmres.1bl.2bl <- lm(dv ~ idv*smod1bl*smod2bl, data = dat4)
slmres.1bl.2ab <- lm(dv ~ idv*smod1bl*smod2ab, data = dat4)
slmres.1ab.2bl <- lm(dv ~ idv*smod1ab*smod2bl, data = dat4)
slmres.1ab.2ab <- lm(dv ~ idv*smod1ab*smod2ab, data = dat4)
summary(slmres.1bl.2bl) # LL
summary(slmres.1bl.2ab) # LH
summary(slmres.1ab.2bl) # HL
summary(slmres.1ab.2ab) # HH