myscratchpad

# 交互作用重回帰の別の例。以下参考。感謝
# http://wilhelmhofmann.info/wordpress/?page_id=120
library(foreign)
dat <- read.spss("http://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")
 

回帰直線と信頼区間のプロット

適当にぐぐったら便利なものがあった。
http://www.iiap.res.in/astrostat/School07/R/confidence.band.r
http://www.iiap.res.in/astrostat/School07/R/reg.html

## confidence.bandという関数を読み込む
source("http://www.iiap.res.in/astrostat/School07/R/confidence.band.r")

x <- rnorm(30)
y <- rnorm(30)
lmres <- lm(y~x)
confidence.band(lmres)

## 単回帰にしか使えない。そのうちまじめに検証しよう

2章を忘れてた
# Kreft, I., & de Leeuw, J. (1998). Introducing multilevel modeling. Newbury Park, CA: Sage. (クレフト, I., デレウ, J. 小野寺孝義 (編訳) (2006). 基礎から学ぶマルチレベルモデル―入り組んだ文脈から新たな理論を創出するための統計手法　ナカニシヤ出版
# 第2章
以下のページを参考。感謝
http://blue.zero.jp/yokumura/R/multilevel/chap4.txt
http://www.ats.ucla.edu/stat/examples/imm/default.htm



library(foreign)
dat0 <- data.frame(read.spss("http://www.ats.ucla.edu/stat/spss/examples/mlm_imm/imm10.sav"))
names(dat0) <- tolower(names(dat0))
dat <- dat0[c(1,2,11,5)]

## 表2.1 学校内平均 p21
size <- as.vector(table(dat[1]))
dat.agmean <- aggregate(dat[c(3,4)], list(dat[,1]), mean)
tbl2.1 <- data.frame(school=1:10, size, dat.agmean[2:3])
tbl2.1

## 表2.2 散布度 共分散行列と相関係数 p21
by(dat[3:4], dat[1], function(x) cov(x)*((nrow(x)-1)/nrow(x)))
by(dat[3:4], dat[1], cor)

## 表2.3 10校の全体回帰 p24
summary(lm(math~1, dat))
summary(lm(math~homework, dat))

## 表2.4 10校の集計回帰 p24
summary(lm(math~1, dat.agmean, weights=size))
summary(lm(math~homework, dat.agmean, weights=size))

## 表2.5 10校に対する文脈モデル p25
dat$schid <- factor(dat$schid)
cx <- rep(dat.agmean[,3], size)
dat$cx <- cx
summary(lm(math~1, dat))
summary(lm(math~homework+cx, dat))

## 表2.6 10校のクロンバックモデル p26
dat$bw <- dat$homework-dat$cx
dat$bb <- dat$cx - mean(dat$homework)
summary(lm(math~1, dat))
summary(lm(math~bw+bb, dat))

## 表2.7 10校のANCOVAモデル p28 # 推定値に関しては教科書どおりの結果にならないが、基本的に切片からの差分が回帰係数とされているので問題なし、ということにしよう (predict(オブジェクト) でもいい) 。決定係数とかは同じ
lmres.null <- lm(math~schid, dat)
lmres.acv <- lm(math~schid+homework, dat)
summary(lmres.null)
summary(lmres.acv)

こっちの続き。参考は同じ
# Kreft, I., & de Leeuw, J. (1998). Introducing multilevel modeling. Newbury Park, CA: Sage. (クレフト, I., デレウ, J. 小野寺孝義 (編訳) (2006). 基礎から学ぶマルチレベルモデル―入り組んだ文脈から新たな理論を創出するための統計手法　ナカニシヤ出版
# 第4章
以下のページを参考。感謝
http://blue.zero.jp/yokumura/R/multilevel/chap4.txt
http://www.ats.ucla.edu/stat/examples/imm/default.htm

library(foreign)
dat0 <- data.frame(read.spss("http://www.ats.ucla.edu/stat/spss/examples/mlm_imm/imm23.sav"))
names(dat0) <- tolower(names(dat0))
data.frame(names(dat0))


library(psych)
describe(dat0)

dat.corr <- dat0[c(3,5,6,8)]
corr.test(dat.corr)

## 表4.2　ヌルモデル　p57
library(nlme)
library(multilevel)
model.null <- lme(math~1, random=~1|schid,  data=dat0, method="ML")
summary(model.null)
VarCorr(model.null)
getVarCov(model.null) # 分散
logLik(model.null)*-2 # 逸脱度
intervals(model.null) # 信頼区間. Randome effect が水準2、Within...が水準1
library(multilevel)
GmeanRel(model.null) # 級内相関。psychometricパッケージにも似たようなのがある


## 表4.3 p58
model1 <- lme(math~homework, random=~1|schid,  data=dat0, method="ML")
summary(model1)
VarCorr(model1)
getVarCov(model1) # 分散
logLik(model1)*-2 # 逸脱度
intervals(model1) # 信頼区間
GmeanRel(model1) # 級内相関

## 表4.4 p59
model2 <- lme(math~homework, random=~homework|schid,  data=dat0, method="ML")
summary(model2)
VarCorr(model2)
getVarCov(model2) # 分散
logLik(model2)*-2 # 逸脱度
intervals(model2) # 信頼区間
GmeanRel(model2) # 級内相関

## 表4.5 p61
model3 <- lme(math~homework+parented, random=~homework|schid,  data=dat0, method="ML")
summary(model3)
VarCorr(model3)
getVarCov(model3) # 分散
logLik(model3)*-2 # 逸脱度
intervals(model3) # 信頼区間
GmeanRel(model3) # 級内相関

## 表4.6 p62
model4 <- lm(math~homework+parented, data=dat0)
summary(model4)
logLik(model4)*-2 # 逸脱度

## 表4.7　p65
model2s2 <- lme(math~homework + scsize, random=~homework|schid,  data=dat0, method="ML")
summary(model2s2)
VarCorr(model2s2)
getVarCov(model2s2) # 分散
logLik(model2s2)*-2 # 逸脱度
intervals(model2s2) # 信頼区間
GmeanRel(model2s2) # 級内相関

## 表4.8　p66
model3s2 <- lme(math~homework + public, random=~homework|schid,  data=dat0, method="ML")
summary(model3s2)
VarCorr(model3s2)
getVarCov(model3s2) # 分散
logLik(model3s2)*-2 # 逸脱度
intervals(model3s2) # 信頼区間
GmeanRel(model3s2) # 級内相関

## 表4.9 p67
dat0$hp <- dat0$homework*dat0$public
dat0$hm <- dat0$homework*dat0$meanses
dat0$hr <- dat0$homework*dat0$ratio
cmat <- cor(dat0[,c("public", "meanses", "ratio", "homework", "hp", "hm", "hr")])
cmat[c(5:7,4),]


## 表4.10　p68
model4s2 <- lme(math~homework+public+homework:public, random=~homework|schid,  data=dat0, method="ML")
summary(model4s2)
VarCorr(model4s2)
getVarCov(model4s2) # 分散
logLik(model4s2)*-2 # 逸脱度
intervals(model4s2) # 信頼区間
GmeanRel(model4s2) # 級内相関

###########
# 表4.11 データなし
##########


## 表4.12　p72
model5s2 <- lme(math~homework+white+public, random=~homework|schid,  data=dat0, method="ML")
summary(model5s2)
VarCorr(model5s2)
getVarCov(model5s2) # 分散
logLik(model5s2)*-2 # 逸脱度
intervals(model5s2) # 信頼区間
GmeanRel(model5s2) # 級内相関

## 表4.13　p73
model6s2 <- lme(math~homework + public + white, random=~homework + white|schid,  data=dat0, method="ML")
summary(model6s2)
VarCorr(model6s2)
getVarCov(model6s2) # 分散
logLik(model6s2)*-2 # 逸脱度
intervals(model6s2) # 信頼区間
GmeanRel(model6s2) # 級内相関

## 表4.14　p74
model7s2 <- lme(math~homework+white+public+meanses, random=~homework|schid,  data=dat0, method="ML")
summary(model7s2)
VarCorr(model7s2)
getVarCov(model7s2) # 分散
logLik(model7s2)*-2 # 逸脱度
intervals(model7s2) # 信頼区間
GmeanRel(model7s2) # 級内相関

## 表4.15　p75
model8s2 <- lme(math~homework+white+meanses, random=~homework|schid,  data=dat0, method="ML")
summary(model8s2)
VarCorr(model8s2)
getVarCov(model8s2) # 分散
logLik(model8s2)*-2 # 逸脱度
intervals(model8s2) # 信頼区間
GmeanRel(model8s2) # 級内相関

## 表4.16　p76
model9s2 <- lme(math~homework+white+meanses+homework:meanses, random=~homework|schid,  data=dat0, method="ML")
summary(model9s2)
VarCorr(model9s2)
getVarCov(model9s2) # 分散
logLik(model9s2)*-2 # 逸脱度
intervals(model9s2) # 信頼区間
GmeanRel(model9s2) # 級内相関

## 表4.17　p77
model10s2 <- lme(math~homework+white+meanses+ses, random=~homework|schid,  data=dat0, method="ML")
summary(model10s2)
VarCorr(model10s2)
getVarCov(model10s2) # 分散
logLik(model10s2)*-2 # 逸脱度
intervals(model10s2) # 信頼区間
GmeanRel(model10s2) # 級内相関

## 表4.18　データなし

## 表4.19　p79
model1s3 <- lme(math~ses, random=~1|schid,  data=dat0, method="ML")
summary(model1s3)
VarCorr(model1s3)
getVarCov(model1s3) # 分散
logLik(model1s3)*-2 # 逸脱度
intervals(model1s3) # 信頼区間
GmeanRel(model1s3) # 級内相関

## 表4.20　p80
model2s3 <- lme(math~ses, random=~ses|schid,  data=dat0, method="ML")
summary(model2s3)
VarCorr(model2s3)
getVarCov(model2s3) # 分散
logLik(model2s3)*-2 # 逸脱度
intervals(model2s3) # 信頼区間
GmeanRel(model2s3) # 級内相関

## 表4.21　p81
model3s3 <- lme(math~ses+percmin, random=~1|schid,  data=dat0, method="ML")
summary(model3s3)
VarCorr(model3s3)
getVarCov(model3s3) # 分散
logLik(model3s3)*-2 # 逸脱度
intervals(model3s3) # 信頼区間
GmeanRel(model3s3) # 級内相関

## 表4.22　p82
model4s3 <- lme(math~ses+percmin+meanses, random=~1|schid,  data=dat0, method="ML")
summary(model4s3)
VarCorr(model4s3)
getVarCov(model4s3) # 分散
logLik(model4s3)*-2 # 逸脱度
intervals(model4s3) # 信頼区間
GmeanRel(model4s3) # 級内相関

## 表4.25　p86
model1s4 <- lme(math~homework+ratio, random=~homework|schid,  data=dat0, method="ML")
summary(model1s4)
VarCorr(model1s4)
getVarCov(model1s4) # 分散
logLik(model1s4)*-2 # 逸脱度
intervals(model1s4) # 信頼区間
GmeanRel(model1s4) # 級内相関

## 表4.26　p87
model2s4 <- lme(math~homework+homework:ratio, random=~homework|schid,  data=dat0, method="ML")
summary(model2s4)
VarCorr(model2s4)
getVarCov(model2s4) # 分散
logLik(model2s4)*-2 # 逸脱度
intervals(model2s4) # 信頼区間
GmeanRel(model2s4) # 級内相関

「基礎から学ぶマルチレベルモデル」をRで実行
# Kreft, I., & de Leeuw, J. (1998). Introducing multilevel modeling. Newbury Park, CA: Sage. (クレフト, I., デレウ, J. 小野寺孝義 (編訳) (2006). 基礎から学ぶマルチレベルモデル―入り組んだ文脈から新たな理論を創出するための統計手法　ナカニシヤ出版
# 第3章

uclaや奥村先生のサイトを参考にした。というかほとんどコピペ。感謝。
上記サイトではlme4パッケージのlmer関数が主だが、lme (nlmeパッケージ) を使えばいいじゃないか、的な記事もあったのでとりあえずlmeで書き直してみた。
RjpWikiのオブジェクト一覧も有用。感謝

library(foreign)
dat0 <- data.frame(read.spss("http://www.ats.ucla.edu/stat/spss/examples/mlm_imm/imm10.sav"))
names(dat0) <- tolower(names(dat0))
data.frame(names(dat0))
dat <- dat0[c(1,2,11,5,8,19)]
head(dat)
dat.agmean <- aggregate(dat[3:5], list(dat[,1]), mean)
dat.agmean

# 傾きを結果変数とするアプローチ
## 表3.2　10校の宿題に対する数学の成績のOLS回帰　p41
by(dat[3:4], dat[1], cor) # 相関係数
by(dat, dat$schid, function(x) coef(summary(lm(math~homework, x)))) # 学校ごとにデータを分割し、math~homeworkの回帰分析を行ったもの
summary(lm(math~homework, dat)) #p24 表2.3と同じ

## マクロ回帰の結果
bylmres <- by(dat, dat$schid, function(x) lm(math~homework, x))
a <- sapply(bylmres, coef)
intc <- a[1,]
slp <- a[2,]
pbl <- factor(dat.agmean[,4])
summary(lm(intc~pbl))
summary(lm(slp~pbl))

## 表3.3　ランダム係数モデル p42
# レベル1 Y=math, X = homework, i=個人, j=学校, ε=誤差
## β_0j, β_1j はそれぞれa_j, b_jとも表記される (教科書中) 。学校の切片よび傾きを表す。
## 最初の添え字 (i) はミクロ水準、後の添え字は (j) はマクロ水準を示す。つまりiはネストされる変数、jはネストする変数。この例では個人がi, 学校がj。また、例えば日ごとの反復測定ならiが日、jが個人となる
## 表3.3、水準1の分散はε_ijであり、重回帰分析でいうところの誤差項にあたる。VarCorr関数のResidualsにあたる。
# Y_ij = β_0j + β_1jX_ij + ε_ij

# レベル2 # γ₀₀=切片の全体の平均, γ₁₀=全体の傾き、u=誤差
## VarCorr関数で出力される分散成分、および教科書表3.3の「水準2」の「推定値」とはここのu_0j, u_1jの誤差分散を示す。
# β_0j = γ₀₀ + u_0j
# β_1j = γ₁₀ + u_1j
library(nlme)
lme.res <- lme(math ~ homework, random=~homework|schid, data=dat, method="ML")
summary(lme.res)
VarCorr(lme.res)
getVarCov(lme.res) # 分散
logLik(lme.res)*-2 # 逸脱度
intervals(lme.res) # 信頼区間

## 表3.4 p44
lme.res2 <- lme(math ~ homework+public, random=~homework|schid, data=dat, method="ML")
summary(lme.res2)
VarCorr(lme.res2)
getVarCov(lme.res2) # 分散
logLik(lme.res2)*-2 # 逸脱度
intervals(lme.res2) # 信頼区間


## 表3.5 p45
lme.res3 <- lme(math~homework+public+homework:public, random=~homework|schid,  data=dat, method="ML")
summary(lme.res3)
VarCorr(lme.res3)
getVarCov(lme.res3) # 分散
logLik(lme.res3)*-2 # 逸脱度
intervals(lme.res3) # 信頼区間

## 表3.6 p46
lme.res4 <- lme(math~homework, random=~homework|schid, data=dat, method="ML")
summary(lme.res4)
VarCorr(lme.res4)
getVarCov(lme.res4) # 分散
logLik(lme.res4)*-2 # 逸脱度
intervals(lme.res4) # 信頼区間
cf.v <- coef(lme.res4)
cf.v


## 図3.8
n<-nrow(cf.v)
x<-cbind(0, 7)
y<-cbind(seq(1:n), seq(1:n))
y[,1]<-cf.v[,1]
y[,2]<-cf.v[,1]+cf.v[,2]*7
plot(x, y[1,], type="l", ylab="predicte", ylim=c(20,100))
for (i in 2:n){points(x, y[i,], type="l", lty=i)}