ESTRELA連載 第82回

注)コードは現在作成中です

今回使うデータ

JAGSの準備

  • R2jagsを参照してください。

データの作成と表示

  • パッケージの読み込み
    library(spdep)
    library(R2jags)
    library(bayesm)
  • データの読み込み
    lph <- read.table("lph.csv", sep=",", header=T)
    summary(lph)

JAGS用データ作成

  • 説明変数・被説明変数・誤差項
    # 変数の指定
    # 被説明変数
    y <- as.vector(lph$LPH)
    # 地域数
    n <- length(y)
    # 地域セグメント数
    m <- max(lph$seg)
    # 説明変数
    x <- as.matrix(cbind(
    lph$POPD,
    lph$EMP3D,
    lph$seg,
    lph$ID
    ))
    # 線形回帰モデルの誤差項(最小二乗法)
    model.lm <- lm(y~x[,1]+x[,2])
    summary(model.lm)
    # e <- as.matrix(resid(model.lm))
  • 空間隣接行列の作成
    coords <- as.matrix(cbind(lph$Easting,lph$Northing))
  • 空間重み付け行列の作成
    nb <- tri2nb(coords)
    # nb.mat <- nb2mat(nb, style="B")
    nb.mat <- nb2mat(nb, style="W")
  • 空間ウェイト付き変数
    yL <- nb.mat %*% y
    xL <- nb.mat %*% x

空間自己回帰モデル(空間自己回帰ラグモデル:Spatial Lag Model)

  • JAGSコード(slag1.txt)
    # spatial lag model
    model{  
    for (i in 1 : n) {          
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- rho*yL[i,1]+b0+b1*x[i,1]+b2*x[i,2]
    }
    b0 ~ dnorm(0.0, 1.0E-6)
    b1 ~ dnorm(0.0, 1.0E-6)
    b2 ~ dnorm(0.0, 1.0E-6)
    tau ~ dgamma(0.001, 0.001)
    rho ~ dnorm(0.0, 1.0E-6)
    sigma <- 1/sqrt(tau)
    }
  • Rコード
    # JAGS変数設定
    # データ
    data <- list("n", "y", "x", "yL")
    # MCMC初期値(事前情報)
    in1 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], rho=0, tau=1)
    in2 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], rho=0.5, tau=1)
    in3 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], rho=1, tau=1)
    inits <- list(in1,in2,in3)
    # パラメータ
    parameters <- c("b0", "b1", "b2", "rho", "tau", "sigma")
    # モデルファイル
    model.file <- system.file(package="R2jags", "model", "slag1.txt")
    # MCMC
    slag.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
    print(slag.jags, digits=3)
    plot(slag.jags)
    traceplot(slag.jags)
    slag.fit <- update(slag.jags)
    print(slag.fit, digits=3)

Spatial Durbin Model

  • JAGSコード(sdm1.txt)
    # spatial durbin model
    model{  
    for (i in 1 : n) {          
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- rho*yL[i,1]+b0+b1*x[i,1]+b2*x[i,2]+g1*xL[i,1]+g2*xL[i,2]
    }
    b0 ~ dnorm(0.0, 1.0E-6)
    b1 ~ dnorm(0.0, 1.0E-6)
    b2 ~ dnorm(0.0, 1.0E-6)
    g1 ~ dnorm(0.0, 1.0E-6)
    g2 ~ dnorm(0.0, 1.0E-6)
    tau ~ dgamma(0.001, 0.001)
    rho ~ dnorm(0.0, 1.0E-6)
    sigma <- 1/sqrt(tau)
    }
  • Rコード
    # JAGS変数設定
    # データ
    data <- list("n", "y", "x", "yL", "xL")
    # MCMC初期値(事前情報)
    in1 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], g1=model.lm$coefficients[2],
    g2=model.lm$coefficients[3], rho=0, tau=1)
    in2 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], g1=model.lm$coefficients[2],
    g2=model.lm$coefficients[3], rho=0.5, tau=1)
    in3 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], g1=model.lm$coefficients[2],
    g2=model.lm$coefficients[3], rho=1, tau=1)
    inits <- list(in1,in2,in3)
    # パラメータ
    parameters <- c("b0", "b1", "b2", "g1", "g2", "rho", "tau", "sigma")
    # モデルファイル
    model.file <- system.file(package="R2jags", "model", "sdm1.txt")
    # MCMC
    sdm.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
    print(sdm.jags, digits=3)
    plot(sdm.jags)
    traceplot(sdm.jags)
    sdm.fit <- update(sdm.jags)
    print(sdm.fit, digits=3)

Spatial Error Model

  • JAGSコード(sem1.txt)
    # spatial error model
    model{  
    for (i in 1 : n) {          
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- lambda*yL[i,1]+b0*(1-lambda)+b1*x[i,1]+b2*x[i,2]
    -b1*lambda*xL[i,1]-b2*lambda*xL[i,2]
    }
    b0 ~ dnorm(0.0, 1.0E-6)
    b1 ~ dnorm(0.0, 1.0E-6)
    b2 ~ dnorm(0.0, 1.0E-6)
    tau ~ dgamma(0.001, 0.001)
    lambda ~ dunif(0,1)
    sigma <- 1/sqrt(tau)
    }
  • Rコード
    # JAGS変数設定
    # データ
    data <- list("n", "y", "x", "yL", "xL")
    # MCMC初期値(事前情報)
    in1 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], lambda=0, tau=1)
    in2 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], lambda=0.5, tau=1)
    in3 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], lambda=1, tau=1)
    inits <- list(in1,in2,in3)
    # パラメータ
    parameters <- c("b0", "b1", "b2", "lambda", "tau", "sigma")
    # モデルファイル
    model.file <- system.file(package="R2jags", "model", "sem1.txt")
    # MCMC
    sem.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
    print(sem.jags, digits=3)
    plot(sem.jags)
    traceplot(sem.jags)
    sem.fit <- update(sem.jags)
    print(sem.fit, digits=3)
  • エラーメッセージ
    Warning message:
    In readLines(file) :
      '/Library/Frameworks/R.framework/Resources/library/R2jags/model/sem1.txt' で 不完全な最終行が見つかりました

マルチレベルモデル

  • bayesmパッケージのrhierLinearModel?()関数を使う推定方法
  • 空間重み付け行列に対するcross validationは行わない
  • 市区町村単位でパラメータ推定しているため計算に時間を要します
    # 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)
    # モデル推定
    out1=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)
    # 推定結果の表示
    summary(out1$Deltadraw, burnin=1000)
    summary(out1$Vbetadraw, burnin=1000)
    summary(t(out1$betadraw[1,,]), burnin=1000)

Geographically Weighted Model

  • bayesmパッケージのrhierLinearModel?()関数を使う推定方法
  • 空間重み付け行列に対するcross validationは行わない
  • 市区町村単位でパラメータ推定しているため計算に時間を要します
    # library(bayesm)
    # mcmcの設定
    R=10000
    keep=1
    # 事前分布を設定
    reg=levels(factor(lph$JCODE))
    nreg=length(reg)
    nvar=3 	#説明変数2つ+定数項
    # 変数を設定
    lph$POPDw <- nb.mat %*% lph$POPD
    lph$EMP3Dw <- nb.mat %*% lph$EMP3D 
    regdata=NULL
    for (j in 1:nreg) {
           y=lph$LPH[lph$JCODE==reg[j]]
           iota=c(rep(1,length(y)))
           X=cbind(iota, lph$POPDw[lph$JCODE==reg[j]],
                   lph$EMP3Dw[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)
    # モデル推定
    out1=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)
    # 推定結果の表示
    summary(out1$Deltadraw, burnin=1000)
    summary(out1$Vbetadraw, burnin=1000)
    summary(t(out1$betadraw[1,,]), burnin=1000)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-22 (水) 17:53:34 (3318d)