2

R の Winbugs を使用してベイジアン分析を行っています。2 つの Winbugs スクリプトを 1 つに結合する必要がありますが、エラー メッセージ ( Variable x2 is not defined in model or in data set) が表示されます。winbugs のコードは次のとおりです。

model{
# Model’s likelihood
 for (i in 1:n) {
    tto[i] ~ dnorm( mu[i], tau ) # stochastic componenent
    b[i] ~ dnorm(0.0, tau2)
    # link and linear predictor
    mu[i] <- 1 - (beta.concern2*concern2[i] + beta.concern3*concern3[i] + b[i])
 }

 for (i in 1:1002) {
    # Linear regression on logit
    logit(p[i]) <- beta.concern2*x2[i,1] + beta.concern2*x2[i,2]

    # Likelihood function for each data point
    y2[i] ~ dbern(p[i])
 }

  s2<-1/tau
  s <-sqrt(s2)

  a2<-1/tau2
  a <-sqrt(a2)

  }

ここx2で、 は 1002*2 行列で、yはベクトルです

これは、データを定義する R コードです。

 combined.data <- list(n=n,tto=tto,concern2=concern2,
                       concern3=concern3,y2=y2, x2=x2)

誰が何が悪いのか知っていますか？

4

1 に答える 1

1

ここでかなりの仮定を立てるつもりです...

おそらく、変数間の関係を示す図を追加できます。これは決定論的と確率論的です。BUGS でモデルを作成するときに、これが役に立ちます。また、すべてのデータの次元、その意味、nおよびモデル化しているものと関心のあるノードに関するコンテキストまたは詳細を把握しておくと役立ちます。

これはy長さ 1002 のバイナリ (0,1) ベクトルであり、x2[,1]およびx2[,2](ここx1では , x2) およびconcern2, concern3(ここc2では , c3) およびttoie

nrow(x2) == 1002

nrow==10を使用して作業するサンプル データを次に示します。

y <- sample(x=c(0,1), size=10, replace=TRUE, prob=c(0.5,0.5))
x2 <- matrix(rnorm(20), nrow=10, ncol=2)
c2 <- rnorm(10)
c3 <- rnorm(10)
tto <- rnorm(10)

ロジットの の両方の値についてbeta.concern2(ここでは )の値を決定しようとしているようです。2 つの異なる予測子に対して同じパラメーターを使用する理由がわかりません。これがタイプミスの場合は、代わりにパラメーターとして指定します。これをお客様のニーズに合わせて調整できることを願っています。b2x2b2b3

b2これらの, b3(確率的) とc2, (与えられた)の値の積はc3、変数 を生成するために使用されmu、これには誤差項もあります。b[i]( (ここで ) は正規分布の誤差項であると推測しb1[i]ています。) 次にtto、 は の値に依存する正規分布の変数でありmu、それ自体に誤差項があります。誤差項の精度はどちらの場合も等しくなるように設定しました。

したがって、そのようなモデルの場合:

require(rjags)

### The data
dataList <- list(
    x1 = x2[,1],
    x2 = x2[,2],
    y = y,
    c2 = c2,
    c3 = c3,
    tto = tto,
    nRowX = nrow(x2)
    )

### make sure logistic model can be fitted
f1 <- stats::glm(dataList$y ~ dataList$x1 + dataList$x2 -1, family=binomial(logit))
show(f1)

### set some approximate initial values
b1Init <- 0.1 # arbitrary
b2Init <- f1$coef[2]
b3Init <- f1$coef[3]

initsList <-  list(
 b1 = b1Init,
 b2 = b2Init,
 b3 = b3Init)

### Model: varying parameters (b2, b3) per observation; 2x error terms
modelstring <- "
    model {
        for(i in 1:nRowX){
            tto[i] ~ dnorm(mu[i], prec)
            mu[i] <- 1 - (b1 + b2*c2[i] + b3*c3[i])
            y[i] ~ dbern(L[i]) # L for logit
            L[i] <- 1/(1+exp(- ( b2*x1[i] + b3*x2[i]) ))
        }
        b1 ~ dnorm(0, prec) # precision
        prec <- 1/sqrt(SD) # convert to Std Deviation
        SD <- 0.5
        b2 ~ dnorm(0, 1.4) # arbitrary
        b3 ~ dnorm(0, 1.4)
    }
"

writeLines(modelstring,con="model.txt")

parameters <- c("b1","b2","b3") # to monitor
adaptSteps <- 1e4 # "tune in" samplers
burnInSteps <- 2e4 # "burn in" samplers
nChains <- 3
numSavedSteps <-2e3
thinSteps <- 1 # Steps to "thin" (1=keep every step).
nPerChain <- ceiling(( numSavedSteps * thinSteps ) / nChains) # Steps per chain

rm(jagsModel) # in case already present

jagsModel <- rjags::jags.model(
 "model.txt", data=dataList,
 inits=initsList, n.chains=nChains,
 n.adapt=adaptSteps)

stats::update(jagsModel, n.iter=burnInSteps)

### MCMC chain
MCMC1 <- as.matrix(rjags::coda.samples(
 jagsModel, variable.names=parameters,
 n.iter=nPerChain, thin=thinSteps))

### Extract chain values
b2Sample <- as.vector(MCMC1[,grep("b2",colnames(MCMC1))])
于 2013-01-25T22:42:22.990 に答える