一般化線形モデル

一般化線形モデルの解析はglm( ) により行う。書式は以下のとおり。

glm(formula, family=familytype(link=linkfunction), data=)

Family Default Link Function
binomial：二項分布 (link = "logit")
gaussian：正規分布 (link = "identity")
Gamma：ガンマ分布 (link = "inverse")
inverse.gaussian：逆正規分布 (link = "1/mu^2")
poisson：ポアソン分布 (link = "log")
quasi：擬似尤度 (link = "identity", variance = "constant")
quasibinomial：擬似二項分布 (link = "logit")
quasipoisson：擬似ポアソン分布 (link = "log")

オプションについて詳しくは help(glm) や help(family) を見よう。
ここでは一般線形モデルのうちロジスティック回帰、ポアソン回帰、生存時間分析について紹介する。.


ロジスティック回帰

ロジスティック回帰は連続量である予測変数から2値変数を予測するときに用いる。制約が少ないので判別分析よりも好まれる。予測変数にカテゴリカル変数を含めることもできる。

# データとspss結果はuclaのサイトから。感謝。
# 以下も参考: R Data Analysis Examples Logit Regression

# データ準備
library(foreign)
dat <- read.spss("http://www.ats.ucla.edu/stat/data/hsb2.sav", to.data.frame=T)
names(dat) <- tolower(names(dat))
honcomp <- as.integer(dat$write>=60)
ses2 <- factor(dat$ses, levels=c("high", "middle", "low")) # リファレンスカテゴリを変える
dat2 <- data.frame(dat, honcomp, ses2)

fit <- glm(honcomp~read+science+ses2, data=dat2, family=binomial("logit")) # ses2はカテゴリカル変数

# 係数 (B) 、係数の標準誤差 (S.E.) 、有意性 (Sig.) 。 () はSPSSでの表示。
summary(fit)

# 信頼区間
confint(fit)

# 回帰係数の指数。オッズ比 (Exp(B))
exp(coef(fit))

# オッズ比の信頼区間
exp(confint(fit))

# 予測値
predict(fit, type="response")

# 尤離度
residuals(fit, type="deviance")

# Wald統計量
library(aod)
# readのwald検定
wald.test(b=coef(fit), Sigma=vcov(fit), Terms=2)
# scienceのwald検定
wald.test(b=coef(fit), Sigma=vcov(fit), Terms=3)

# Nagelkerke R2 (中澤先生のサイトより。感謝。)
NagelkerkeR2 <- function(rr,n) {(1-exp((rr$dev-rr$null)/n))/(1-exp(-rr$null/n)) }
print(NagelkerkeR2(fit,nrow(dat)))

# 対数尤度
logLik(fit) # summary()のResidual deviance
-2*logLik(fit)

# プロット
x1 <- fit$linear.predictor # 予測変数
y <- dat2$honcomp # 目的変数
x2 <- seq(from=min(x1), to=max(x1), length=200) # 曲線描画用
y2 <- exp(x2)/(1+exp(x2)) # 曲線描画用

par(mfrow=c(2,1))
plot(x1, y, pch = ifelse(y=="1", 16, 1))
lines(x2, y2)
cdplot(factor(y)~x1)

anova(fit1,fit2, test="Chisq") とするとモデルの適合度を比較できる。
cdplot(F~x, data=mydata) は、2値変数Fを連続変量xで予測したときの条件つきの密度プロットを表示できる。

conditional density plot click to view


ポアソン回帰

連続変量から計数データを予測するときに用いる。

# 例によってUCLAのサイトから。感謝。結果はSPSSに一致。

# データ準備
library(foreign)
dat <- read.spss("http://www.ats.ucla.edu/stat/spss/dae/poisson_sim.sav", to.data.frame=T)
prog2 <- factor(dat$prog, levels=c("vocation", "academic", "general"))
dat2 <- data.frame(dat, prog2)

# 記述統計
library(psych)
describe.by(dat$num_awards, list(dat2$prog2))

fit <- glm(num_awards~prog2+math, data=dat2, family=poisson())
summary(fit)

# 係数の信頼区間、オッズ比、オッズ比の信頼区間
confint(fit)
exp(coef(fit))
exp(confint(fit))

過分散 (overdispersion) が生じたとき、すなわち残差尤離度が自由度を超えるような場合は、誤差分布にquasipoisson() を使用する。


生存時間分析

生存時間分析 (survival analysis, event history analysis or reliability analysisともいう) はイベント (e.g., 死亡、故障) とそこまでの時間、すなわち生存時間とほかの変数との関係をモデリングすることである。データは研究終了までイベントが起きない right censored であったり、イベントは起きていないが参加者がドロップアウトしたような不完全情報であったりすることもある。

生存時間分析にはglm関数ではなく survival パッケージの関数を使う。survival パッケージはone and two sample problems, parametric accelerated failure modelやコックス比例ハザードモデルを扱える。

データにはstart time, stop time, and status (1=イベント発生, 0=イベントなし) の情報がある。あるいはイベントまでの時間 (time to event) と status という場合もある。status=0とは先述のright cencoredを意味する。データはまずSurv関数をもちいてSurvオブジェクトにしておく必要がある。

survfit( ) 関数でひとつまたは複数の群の生存分布を推定できる。
survdiff( ) 関数により2つ以上の群の生存分布の差を検定できる。
coxph( ) 関数で予測変数に対するハザード関数のあてはまりをモデリングできる。

以下参考

# Mayo Clinic Lung Cancer Data
library(survival)

# データセットについて
help(lung)

# Surv オブジェクトの生成
survobj <- with(lung, Surv(time,status))

# 全サンプルの生存分布をプロット
# Kaplan-Meier 推定
fit0 <- survfit(survobj~1, data=lung)
summary(fit0)
plot(fit0, xlab="Survival Time in Days",
   ylab="% Surviving", yscale=100,
   main="Survival Distribution (Overall)")

# 性別で生存分布を比較
fit1 <- survfit(survobj~sex,data=lung)

# 性別ごとにプロット
plot(fit1, xlab="Survival Time in Days",
  ylab="% Surviving", yscale=100, col=c("red","blue"),
  main="Survival Distributions by Gender")
  legend("topright", title="Gender", c("Male", "Female"),
  fill=c("red", "blue"))

# 性差の検定
# 生存曲線 (ログランク検定)
survdiff(survobj~sex, data=lung)

# 年齢と医療得点から男性の生存を予測
MaleMod <- coxph(survobj~age+ph.ecog+ph.karno+pat.karno,
  data=lung, subset=sex==1)

# 結果の表示
MaleMod

# 比例ハザード仮定の評価
cox.zph(MaleMod)

survival distribution for total group survival distribution by gender click to view

Thomas Lumley's R newsの記事、Mai Zhou のUse R Software to do Survival Analysis and Simulation 、M. J. CrawleyのSurvival Analysis も参照。


マルチレベルモデル

線形混合モデル、一般化線形混合モデルとかいったりする。

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

以下参考


第2章

## 表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)


第3章

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 はそれぞれaj, bjとも表記される (教科書中) 。学校の切片よび傾きを表す。
## 最初の添え字 (i) はミクロ水準、後の添え字は (j) はマクロ水準を示す。つまりiはネストされる変数、jはネストする変数。この例では個人がi, 学校がj。また、例えば日ごとの反復測定ならiが日、jが個人となる
## 表3.3、水準1の分散はεijであり、重回帰分析でいうところの誤差項にあたる。VarCorr関数のResidualsにあたる。
# Yij = β0j + β1jXij + εij

# レベル2 # γ00=切片の全体の平均, γ10=全体の傾き、u=誤差
## VarCorr関数で出力される分散成分、および教科書表3.3の「水準2」の「推定値」とはここのu0j, u1jの誤差分散を示す。
# β0j = γ00 + u0j
# β1j = γ10 + u1j
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)}


第4章

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) # 級内相関