#======================================================
#	6 一般化線形モデル
#======================================================
#======================================================
#	6.3 ロジスティック回帰
#======================================================
#------------<二項分布>------------
y <- 0:20
#図6.1
plot(y, dbinom(y, size=20, prob=0.2), type="h", main=expression(mu==0.2))
#図6.2
plot(y, dbinom(y, size=20, prob=0.5), type="h", main=expression(mu==0.5))

#------------<適用例１>------------

dose<-c(1.6907,1.7242,1.7552,1.7842,
	1.8113,1.8369,1.8610,1.8839) #気体二硫化炭素
n<-c(59,60,62,56,63,59,62,60)       #観測数
y<-c(6,13,18,28,52,53,61,60)        #死亡数

beetle.glm <- glm(cbind(y, n-y) ~ dose, family=binomial(link=logit))
summary(beetle.glm)

resid <- resid(beetle.glm)
quantile(resid)

#説明された逸脱度
(beetle.glm$null.deviance-beetle.glm$deviance)/beetle.glm$null.deviance

#関数化
dev.explained <- function(model) {
    null.dev <- model$null.deviance
    res.dev <- model$deviance
    ans <- (null.dev - res.dev)/null.dev
    ans
}
dev.explained(beetle.glm)
 
#図6.3（通常の診断プロット）
op <- par(mfrow=c(2, 2))
plot(beetle.glm)
par(op)

図6.4（応答変数の尺度）
x <- seq(0, 2, length=200)
pred <- predict(beetle.glm, newdata=data.frame(dose=x), type="response")
plot(dose, y/n, xlim=c(1.6, 2), ylim=c(0, 1), pch=16, col="steelblue")
lines(x, pred)

#検算
y.n <- y/n #確率
mu <- beetle.glm$fitted　#応答の当てはめ値
w <- n　                 #観測数（重み）
devy <- y.n　
nz <- y.n != 0　         #確率０を除く
devy[nz] <- y.n[nz] * log(y.n[nz])
nz <- (1 - y.n) != 0　   #確率１を除く
devy[nz] <- devy[nz] + (1 - y.n[nz])*log(1 - y.n[nz]) #飽和モデルの対数尤度
devmu <- y.n * log(mu) +(1 - y.n)*log(1 - mu)　       #推定モデルの対数尤度
devi <- 2 * (devy -devmu)　                           #残逸脱度
sign(y.n - mu) * sqrt(abs(devi) * w)　                #各逸脱度残差
sum(w * devi)　          #残逸脱度の重み付き合計を計算する。

#------------<適用例２>------------

#初期モデル
names(肝臓病)
肝臓病.glm0 <- glm(肝臓病 ~ .^2, data = 肝臓病, family=binomial)

#モデル選択
library(MASS)
肝臓病.glm1 <- stepAIC(肝臓病.glm0)

#要約
summary(肝臓病.glm1)

#説明された逸脱度
dev.explained(肝臓病.glm1)

#通常の診断プロット
op <- par(mfrow=c(2, 2))
plot(肝臓病.glm1)
par(op)

#応答変数の尺度
#図6.5
肝臓病.fit <- predict(肝臓病.glm1)
肝臓病.pred <- predict(肝臓病.glm1, type="response")
ord <- order(肝臓病.fit)
plot(肝臓病.fit[ord], 肝臓病.pred[ord], type="l", 
	xlab="線形予測子", ylab="確率", xlim=c(-2, 10))
rug(肝臓病.fit[肝臓病$肝臓病!=1], side=1, col="springgreen")
rug(肝臓病.fit[肝臓病$肝臓病==1], side=3, col="steelblue")
grid()

#判別結果
ans <- table(肝臓病$肝臓病, (肝臓病.pred>0.5))
dimnames(ans)[[1]] <- c("肝臓病以外（真）", "肝臓病（真）")
dimnames(ans)[[2]] <- c("肝臓病以外（予測）", "肝臓病（予測）")
ans
(error <- (sum(ans)-sum(diag(ans)))/sum(ans)) #誤判別率
sum(肝臓病$肝臓病==1)/length(肝臓病$肝臓病)   #肝臓病比率（実績）
sum(肝臓病.pred>0.5)/length(肝臓病.pred)      #肝臓病比率（予測）

#======================================================
#	6.4 ネストされたロジスティック回帰
#======================================================

library(mlbench)
library(car)
data(DNA)
names(DNA)
table(DNA$Class)       #応答変数の確認
temp <- DNA[ ,61:181]  #データの整理

#応答変数の作成

mRNA <- recode(temp$Class, "'n'='no'; else ='yes'")
table(mRNA, temp$Class)
exon <- recode(temp$Class, "'ie'='no'; 'ei'='yes'; else =NA")
table(exon, temp$Class)

#因子に戻す（少し時間がかかる。）

n <- nrow(temp)                 　#観測数
dna <- matrix(0, nrow=n, ncol=40) #初期化
for(i in 1:n) {
	for(j in 1:40) {
		if(temp[i, (j-1)*3+1]==1)
			dna[i, j] <- "A"
		else 
		if(temp[i, (j-1)*3+2]==1)
			dna[i, j] <- "C"
		else
		if(temp[i, (j-1)*3+3]==1)
			dna[i, j] <- "G"
		else
			dna[i, j] <- "T"
	}
}
dna <- data.frame(dna, Class=temp$Class)

#確認
temp[1:5, 1:15]
dna[1:5, 1:5]

#サンプリング（例）
set.seed(456)
sample(1:10, 5)

#サンプリング（実際）
set.seed(456)
n <- nrow(dna)
train <- sample(1:n, ceiling(n/2))

#モデルの当てはめ
mRNA.glm0 <- glm(mRNA[train] ~ ., data=dna[train, ], family=binomial)
mRNA.glm1 <- stepAIC(mRNA.glm0)
summary(mRNA.glm1)
dev.explained(mRNA.glm1)

exon.glm0 <- glm(exon[train] ~ ., data=dna[train, ], family=binomial, 
    maxit=50)
exon.glm1 <- stepAIC(exon.glm0)
summary(exon.glm1)
dev.explained(exon.glm1)

#応答変数の尺度における予測
p.mRNA <- predict(mRNA.glm1, newdata=dna[-train, ], type="response")
p.exon <- predict(exon.glm1, newdata=dna[-train, ], type="response")
p.ei <- p.mRNA * p.exon        #エクソン
p.ie <- p.mRNA * (1 - p.exon)  #イントロン
p.n <- 1 - p.mRNA              #それ以外

p.all <- cbind(p.ei, p.ie, p.n)　#まとめる
head(p.all)

p.max <- apply(p.all, 1, which.max) #どの確率が最大か　
p.dna <- character(length(p.max))　 #文字ベクトルを初期化
p.dna[p.max==1] <- "ei"　　#エクソンの確率が最大である観測値
p.dna[p.max==2] <- "ie"　　#イントロンの確率が最大である観測値
p.dna[p.max==3] <- "n"　　 #それ以外の確率が最大である観測値
ans <- table(DNA$Class[-train], p.dna)　　#分類結果のクロス集計
dimnames(ans)[[1]] <- c("エクソン（真）", "イントロン（真）", 
    "それ以外（真）")
dimnames(ans)[[2]] <- c("エクソン（予測）", "イントロン（予測）", 
    "それ以外（予測）")
ans
(error <- (sum(ans)-sum(diag(ans)))/sum(ans)) #誤判別率

#======================================================
#	6.5 その他の二項分布モデル
#======================================================

dose<-c(1.6907,1.7242,1.7552,1.7842,
	1.8113,1.8369,1.8610,1.8839) #気体二硫化炭素
n<-c(59,60,62,56,63,59,62,60) #観測数
y<-c(6,13,18,28,52,53,61,60)  #死亡数
beetle.glm <- glm(cbind(y, n-y) ~ dose, family=binomial(link=logit))

#---------------<プロビット回帰>---------------

#モデルの当てはめ
beetle.probit <- glm(cbind(y, n-y) ~ dose, family=binomial(link=probit))
summary(beetle.probit)

#説明された逸脱度
dev.explained(beetle.probit)

#図6.6（通常の診断プロット）
op <- par(mfrow=c(2, 2))
plot(beetle.probit)
par(op)

図6.7（応答変数の尺度）
x <- seq(0, 2, length=200)
pred.probit <- predict(beetle.probit, newdata=data.frame(dose=x), type="response")
pred.logit <- predict(beetle.glm, newdata=data.frame(dose=x), type="response")
plot(dose, y/n, xlim=c(1.6, 2), ylim=c(0, 1), pch=16, col="steelblue")
lines(x, pred.probit)
lines(x, pred.logit, lty=2)
legend("topleft", legend=c("probit", "logit"), lty=1:2)

#---------------<補対数対数回帰>---------------

#モデルの当てはめ
beetle.cloglog <- glm(cbind(y, n-y) ~ dose, family=binomial(link=cloglog))
summary(beetle.cloglog)

#説明された逸脱度
dev.explained(beetle.cloglog)

#図6.8（通常の診断プロット）
op <- par(mfrow=c(2, 2))
plot(beetle.cloglog)
par(op)

図6.9（応答変数の尺度）
x <- seq(0, 2, length=200)
pred.cloglog <- predict(beetle.cloglog, newdata=data.frame(dose=x), type="response")
pred.logit <- predict(beetle.glm, newdata=data.frame(dose=x), type="response")
plot(dose, y/n, xlim=c(1.6, 2), ylim=c(0, 1), pch=16, col="steelblue")
lines(x, pred.cloglog)
lines(x, pred.logit, lty=2)
legend("topleft", legend=c("cloglog", "logit"), lty=1:2)

#======================================================
#	6.4 ポアソン回帰
#======================================================
#------------<ポアソン分布>------------
y <- 0:15
#図6.10
plot(y, dpois(y, lambda=0.8), type="h", main=expression(lambda==0.8))
#図6.11
plot(y, dpois(y, lambda=4), type="h", main=expression(lambda==4))

#------------<適用例：電子回路の故障発生件数>------------
data(solder, package="rpart") #データのロード

summary(solder)　#要約
solder$Panel <- as.factor(solder$Panel) #因子へ変換
summary(solder)　#要約

#図6.12
hist(solder$skips, col="springgreen", main="電子回路の不良箇所",
	xlab="不良箇所の数", ylab="度数")

#図6.13
plot(skips ~ Mask, data=solder, col="skyblue")

#------------<初期モデルの当てはめ>------------
glm.pois0 <- glm(skips ~ .^2, data=solder, family=poisson)

#------------<最終モデルの当てはめ>------------
library(MASS) #stepAICのため
glm.pois1 <- stepAIC(glm.pois0)
summary(glm.pois1)

#説明された逸脱度
dev.explained(glm.pois1)

anova(glm.pois1, test="Chi")

#------------<診断プロット>------------
#図6.14
par(mfrow=c(2, 2))
plot(glm.pois1)
par(op)

#図6.15
lp <- predict(glm.pois1)　                 #線形予測子の推定結果
res <- predict(glm.pois1, type="response") #応答変数の尺度
ord <- order(lp)
plot(lp, solder$skips, col="orange", xlab="線形予測子", ylab="応答変数")
lines(lp[ord], res[ord], col="steelblue", lwd=3)

#------------<IRLS>------------
irls.pois <- function(y, X, maxit=50) {
maxit <- maxit    #最大反復回数
eps <- 1e-7     #収束判断値

mu <- y         #muの初期化
eta <- log(mu)  #etaの初期化=poisson()$link=log

for(i in 1:maxit) {
    z <- (y-mu)/mu + eta  #擬似データ：g'(mu)=log(mu)'=1/muから
    w <- as.numeric(mu)   #重み：V(mu)=mu, g'(mu)^2=1/mu^2から
    wz <- sqrt(w)*z       #重み付き擬似データ=応答変数
    WX <- sqrt(w)*X       #重み付き予測変数
    fit <- lm(wz ~ WX-1)  #重み付き線形回帰
    beta <- coef(fit)     #回帰係数

    eta.old <- eta
    eta <- X %*% beta     #新しい線形予測子eta
    mu <- exp(eta)        #新しい推定値mu
    cat(paste(i,"回目：デヴィアンス=", deviance(fit),"\n", sep=""))
    if(max(abs(eta-eta.old)) < eps*max(abs(eta))) break
    }
fit
}

y <- solder$skips             #応答変数
X <- model.matrix(glm.pois1)  #モデル行列
fit.irls <- irls.pois(y, X)   #IRLS当てはめ
coef(fit.irls)[1:5]           #IRLS当てはめ係数の最初のいくつか
coef(glm.pois1)[1:5]          #関数glmによる当てはめ係数の最初のいくつか

#======================================================
#	6.5 ガンマ回帰
#======================================================
#------------<ガンマ分布>------------
#図6.16
y <- seq(0, 40, length=200)
op <- par(mfrow=c(2, 2))
plot(y, dgamma(y, shape=1, rate=1), type="l", 
	main=expression(paste(nu==1,",",nu/mu==1)))
plot(y, dgamma(y, shape=10, rate=1), type="l", 
	main=expression(paste(nu==10,",",nu/mu==1)))
plot(y, dgamma(y, shape=1, rate=0.5), type="l", 
	main=expression(paste(nu==1,",",nu/mu==0.5)))
plot(y, dgamma(y, shape=10, rate=0.5), type="l", 
	main=expression(paste(nu==10,",",nu/mu==0.5)))
par(op)

#------------<適用例：動物のレム睡眠時間>------------
data(msleep, package="ggplot2") #データセットmsleepのロード
summary(msleep)

#図6.17
hist(msleep$sleep_rem, col="springgreen", 
	main="レム睡眠の時間の分布", xlab="レム睡眠時間", ylab="度数")

sleep.lm <- lm(sleep_rem ~ vore + awake, data=msleep)
summary(sleep.lm)

#図6.18
op <- par(mfrow=c(2, 2))
plot(sleep.lm)
par(op)

sleep.Gam <- glm(sleep_rem ~ vore + awake, na.action=na.exclude,
	data=msleep, family=Gamma(link=log))
summary(sleep.Gam)
dev.explained(sleep.Gam)    #説明された逸脱度
anova(sleep.Gam, test="Chi")
#図6.19
op <- par(mfrow=c(2, 2))
plot(sleep.Gam)
par(op)

#応答変数の尺度
#図6.20
lp <- predict(sleep.Gam)　#線形予測子の推定結果
res <- predict(sleep.Gam, type="response") #応答変数の尺度
ord <- order(lp)
plot(lp, msleep$sleep_rem, col="orange", pch=16,
	xlab="線形予測子", ylab="応答変数")
lines(lp[ord], res[ord], col="steelblue", lwd=3)
