ベイズ統計 第6回演習

  • ここでは、慶應大学SFCでの開講科目「ベイズ統計」での演習で用いている分析例やRコード・JAGSコードなどを、授業履修者のために示しています。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rパッケージやデータの更新などにより、本書執筆時とはデータやRコードを予告無しに変更する場合があります。 本書記載のR分析例以外は、R最新版等で動作確認できていないRコードも存在しますが、参考のために記載しておきます(旧バージョンで動作することがあります)。

第6回R演習の内容

  • 線形回帰モデル及びマルチレベル・モデルのベイズ推定について理解を深める
  • データ

データの作成と表示

# パッケージの読み込み
library(nlme) # マルチレベルモデル
# library(lme4)  # マルチレベルモデル
library(mgcv) # 一般化加法モデル
library(MASS)
# データの読み込み
lph <- read.table("lph.csv", sep=",", header=T)
summary(lph)

回帰モデル

# 線形回帰モデル
# 通常最小二乗法
# GISデータの属性からモデル推定
lph.lm <- lm(LPH~POPD+EMP3D,data=lph)
summary(lph.lm)
# 誤差の算出
lph.lm.resid <- resid(lph.lm)
kanto$lm.resid <- lph.lm.resid
#ベイズ線形回帰モデル
# MCMCpackのMCMCregressを使う方法
# パッケージの呼び出し
library(MCMCpack)
lph.mcmc <- MCMCregress(LPH~POPD+EMP3D,data=lph, b0=0, B0=1e-6, c0=1e-2, d0=1e-2, mcmc=100000, burnin=10000)
summary(lph.mcmc)
# 図10.3
plot(lph.mcmc)
# Jagsを使う方法
# パッケージの読み込み
library(R2jags)
# JAGS用データ作成 
# 説明変数・被説明変数・誤差項
# 変数の指定
# 被説明変数
y <- as.vector(lph$LPH)
# 地域数
n <- length(y)
# 説明変数
x <- as.matrix(cbind(
lph$POPD,
lph$EMP3D
))
# 線形回帰モデルの誤差項(最小二乗法)
model.lm <- lm(y~x[,1]+x[,2])
summary(model.lm)
#
####
# JAGSコード(lm.txt)
model{
for(i in 1:n){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- b0+b1*x[i,1]+b2*x[i,2] 
}
b0 ~ dnorm(0,1.0E-6)
b1 ~ dnorm(0,1.0E-6)
b2 ~ dnorm(0,1.0E-6)
tau ~ dgamma(0.01, 0.01)
sigma <- 1/sqrt(tau)
}
#
####
# Rコード
# JAGS変数設定
# データ
data <- list("n", "y", "x")
# MCMC初期値(事前情報)
in1 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
b2=model.lm$coefficients[3], tau=1)
in2 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
b2=model.lm$coefficients[3], tau=1)
in3 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
b2=model.lm$coefficients[3], tau=1)
inits <- list(in1,in2,in3)
# パラメータ
parameters <- c("b0", "b1", "b2", "tau", "sigma")
# モデルファイル
model.file <- system.file(package="R2jags", "model", "lm.txt")
# MCMC
lm.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
n.burnin=1000, n.chains=3, model.file=model.file)
print(lm.jags, digits=3)
# jpeg(filename="fig1.jpg")
plot(lm.jags)
# dev.off()
traceplot(lm.jags)
# jpeg(filename="fig2.jpg")
# traceplot(lm.jags, varname="b0")
# dev.off()
lm.fit <- update(lm.jags)
print(lm.fit, digits=3)

マルチレベルモデル

# lmer()関数を使ったマルチレベルモデルの推定
# パッケージ
library(lme4)
# 固定効果:傾き、ランダム効果:切片
lph.lme1 <- lmer(LPH~POPD+EMP3D+(1|AREA), data=lph)
summary(lph.lme1)
ranef(lph.lme1)
# random.effects(lph.lme1)
# 固定効果:切片、ランダム効果:傾き
lph.lme2 <- lmer(LPH~1+(0+POPD+EMP3D|AREA), data=lph)
summary(lph.lme2)
# random.effects(lph.lme2)
ranef(lph.lme2)
# mcmcsamp()関数を使ったマルチレベルモデルのベイズ推定
lph.mcmc1 <- mcmcsamp(lph.lme1, n=1000)
print(lph.mcmc1)
densityplot(lph.mcmc1)
#
lph.mcmc2 <- mcmcsamp(lph.lme2, n=1000)
print(lph.mcmc2)
densityplot(lph.mcmc2)
# ベイズ推定
# bayesmパッケージのrhierLinearModel()関数を使う推定方法
# 市区町村単位でパラメータ推定しているため計算に時間を要します
library(bayesm)
# mcmcの設定
R=10000
keep=1
# 事前分布を設定
reg=levels(factor(lph$JCODE))
nreg=length(reg)
nvar=3 	#説明変数2つ+定数項
# 変数を設定
regdata=NULL
for (j in 1:nreg) {
       y=lph$LPH[lph$JCODE==reg[j]]
       iota=c(rep(1,length(y)))
       X=cbind(iota, lph$POPD[lph$JCODE==reg[j]],
               lph$EMP3D[lph$JCODE==reg[j]])
       regdata[[j]]=list(y=y,X=X)}
Z=matrix(c(rep(1,nreg)),ncol=1)
Data1=list(regdata=regdata,Z=Z)
Mcmc1=list(R=R,keep=1)
set.seed(66)
# モデル推定
out=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)
# 推定結果の表示
# summary(out$Deltadraw, burnin=1000)
# summary(out$Vbetadraw, burnin=1000)
summary(t(out$betadraw[1,,]), burnin=1000)
plot(out$betadraw, burnin=1000)
#
out_mat <- matrix(0,nrow=nreg,ncol=nvar+3)
colnames(out_mat) <- c("JCODE","POPD","EMP3D", "const","b.POPD","b.EMP3D")
out_mat <- as.data.frame(out_mat)
out_mat$JCODE <- lph$JCODE
out_mat$POPD  <- lph$POPD
out_mat$EMP3D <- lph$EMP3D
for(i in 1:nreg){ 
	for(j in 1:nvar){
		out_mat[i,j+3] <- mean(out$betadraw[i,j,1001:10000])
}}

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-07-12 (土) 02:27:17 (1930d)