Упорядоченный пробит в пробках Использование масштабированного инверсного желаемого

Question

Я пытаюсь использовать следующий код (адаптированный из кода, приведенного в книге Гельмана и Хилла), для оценки переменной модели/перехвата модели пробит-кода в Jags. Тем не менее, он дает мне «Наблюдаемый узел несовместим с незаметными родителями при инициализации. Попробуйте установить соответствующие начальные значения». Где я иду не так? Может кто-нибудь, пожалуйста, помогите мне? Заранее спасибо !!

rm(list=ls(all=TRUE)); 
options(warn=-1) 
library(mvtnorm) 
library(arm) 
library(foreign) 
library("R2jags") 
library(MCMCpack) 

set.seed(1) 

standardizeCols = function(dataMat) { 
    zDataMat = dataMat 
    for (colIdx in 1:NCOL(dataMat)) { 
     mCol = mean(dataMat[,colIdx]) 
     sdCol = sd(dataMat[,colIdx]) 
     zDataMat[,colIdx] = (dataMat[,colIdx] - mCol)/sdCol 
    } 
    return(zDataMat) 
} 

keep<-1 
nobs = 150; 
nis<-sample(1:40,nobs,replace=T) # number obs per subject 
id<-rep(1:nobs,nis) 
N<-length(id) 
corr_beta = 0.6; 
Sigma_beta = matrix(c(1, corr_beta, corr_beta, corr_beta, 
       corr_beta, 1, corr_beta, corr_beta, 
       corr_beta, corr_beta, 1, corr_beta, 
       corr_beta, corr_beta, corr_beta, 1), ncol=4); 
betas <- rmvnorm(n=N, mean=c(-1.45, 0.90, 0.25, -2.3), sigma=Sigma_beta); 

#Generate the data 
x3 = matrix(0, nrow=N,ncol=3); 
y3 = matrix(0, nrow=N,ncol=1); 

for (i in 1:N) { 
    error_v = rnorm(1,0,1); 
    x3[i,1] = rnorm(1,0,1); 
    x3[i,2] = rnorm(1,0,1); 
    x3[i,3] = rnorm(1,0,1); 
    y3[i,1] = betas[id[i], 1] + betas[id[i], 2]*x3[i,1] + betas[id[i], 3]*x3[i,2] + betas[id[i], 4]*x3[i,3] + error_v; 
} 
cutoff=c(-100, 0, 1.5, 2.4, 100) 
k=length(cutoff)-1; 
Y3<-cut(y3, br = cutoff, right=TRUE, include.lowest = TRUE, labels = FALSE) 
Y3=Y3 
X3=x3 
m1=max(Y3) 

y = as.vector(Y3) 
n = length(y) 
J<-length(unique(id)) 
X = cbind(1, standardizeCols(X3)) 
nPred = NCOL(X) 
subjects<-as.vector(as.numeric(id)) 
K=nPred 
W <- diag (K) 

# MCMC settings 
ni <- 5000; nb <- 2500; nt <- 6; nc <- 3 

tau1u=c(0,1,2) 
jags_data <- list ("n", "J", "K", "y", "subjects", "X", "W", "m1") 

inits <- function(){ 
    list (B.raw=array(rnorm(J*K),c(J,K)), mu.raw=rnorm(K), sigma.y=runif(1), Tau.B.raw=rwish(K+1,diag(K)), xi=runif(K)) 
} 

params <- c ("B", "mu", "sigma.B", "rho.B", "tau1u") 

cat("model { 

    for (i in 1:n){ 
     y.hat[i] <- inprod(B[subjects[i],],X[i,]) 
     y[i] ~ dcat(p[i,]) 
     estar[i]~dnorm (y.hat[i], tau.y); 


     for (j in 1:(m1-1)) { 
      Q1[i,j]<-pnorm(tau1[j]-estar[i],0,1) 
     } 

     p[i,1] <- Q1[i,1] 
     for(j in 2:(m1-1)) { 
      p[i,j] <- Q1[i,j] - Q1[i,j-1] 
     } 
     p[i,m1] <- 1 - Q1[i,m1-1] 
    } 

    tau.y <- pow(sigma.y, -2) 
    sigma.y ~ dunif (0, 100) 

    # thresholds (unordered priors) 
    for(j in 1:(m1-1)){ 
     tau1u[j] ~ dnorm(0,.01) 
    } 

    # ordered thresholds 
    tau1 <- sort(tau1u) 

    for (j in 1:J){ 
     for (k in 1:K){ 
      B[j,k] <- xi[k]*B.raw[j,k] 
     } 
     B.raw[j,1:K] ~ dmnorm (mu.raw[], Tau.B.raw[,]) 
    } 
    for (k in 1:K){ 
     mu[k] <- xi[k]*mu.raw[k] 
     mu.raw[k] ~ dnorm (0, .0001) 
     xi[k] ~ dunif (0, 100) 
    } 

    Tau.B.raw[1:K,1:K] ~ dwish (W[,], df) 
    df <- K+1 
    Sigma.B.raw[1:K,1:K] <- inverse(Tau.B.raw[,]) 

    for (k in 1:K){ 
     for (k.prime in 1:K){ 
      rho.B[k,k.prime] <- Sigma.B.raw[k,k.prime]/sqrt(Sigma.B.raw[k,k]*Sigma.B.raw[k.prime,k.prime]) 
     } 
     sigma.B[k] <- abs(xi[k])*sqrt(Sigma.B.raw[k,k]) 
    } 
}", fill=TRUE, file="wishart2.txt") 

# Start Gibbs sampler 
outj <- jags(jags_data, inits=inits, parameters.to.save=params, model.file="wishart2.txt", n.thin=nt, n.chains=nc, n.burnin=nb, n.iter=ni) 

Answer 1

Ваша начальная функция возвращает значения случайных чисел от нормальных и равномерных распределений, которые он появляется не достаточно близко чувственных значений, чтобы не-0 значение кзади рассчитывается. Я думаю, вам нужно более тщательно выбирать начальные значения и, возможно, основываться на значениях, генерируемых в данных, чтобы обеспечить компиляцию модели. Получают ли Gelman и Hill исходные значения для своей модели, с которой вы могли бы начать?

Обновление: вы также можете попробовать удалить аргумент 'inits = inits', чтобы позволить JAGS выбирать свои собственные начальные значения, которые работают для большинства (хотя и не для всех) моделей. Я не использую R2JAGS, хотя я не уверен, что это разрешено для функции jags (но это для rjags и runjags).