分层狄利克雷回归中的随机截距锯齿

问题描述

我有以下数据结构：

• y：3 列，表示多年来观察到的死亡比例。
• x1：GDP - 与每年相关的连续变量
• x2：与死亡相关的年龄

这里是模型规格：

Model

这里是模拟数据：

for line in iter(command.stderr.readline,b''):
    line = line.rstrip().decode('utf8')
    for x in re.split("(?<=.)(?=copying)",line):
        logger.debug(x)

Jags 模型

library(tidyverse)

Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)

# X2 AGES
n.ages <- length(Ages)

# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)

# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)


Year.tot <- rep(Year,length(unique(Ages)))

DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),C1=C_Porp[,1],C2=C_Porp[,2],C3=C_Porp[,3]) 

# X1 GDP OVER YEARS
GDP <- rnorm(n.year,15)
GDP <- as.data.frame(cbind(GDP,Year))

# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)

我收到此错误：

library(runjags)
dirlichet.model = 
  
  
  "model {
  #setup priors for each species
  for(j in 1:N.c){
    m0[j] ~ dnorm(0,1.0E-3) #intercept prior
    m1[j] ~ dnorm(0,1.0E-3) #      mat prior
  }
  
  #implement dirlichet
  for(i in 1:N){
  
    y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 
     
    for(j in 1:N.c){
    
      a0[i,j] ~ dnorm(mu[i,j],0.001) 
      
      log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])
    
    }
    
     # Effect of site:
    for (s in 1:S){
      
      alpha[s] ~ dnorm(0,0.001)
      
    }}
    
  }
"

jags.data <- list(y = D[,c(3,4,5)],mat = D[,5],Ages=D[,N = nrow(D[,5)]),N.c = ncol(D[,S=(length(unique(D[,1]))))

jags.out <- run.jags(dirlichet.model,data=jags.data,adapt = 500,burnin = 5000,sample = 10000,n.chains=5,monitor=c('m0','m1',"mu"))
out <- summary(jags.out)

对模型规范有什么建议吗？

解决方法

还有其他几个小错误。不过，这段代码现在应该可以工作了。问题是：

1. Age 被定义为从 0 到 50 乘以 5 的序列，但随后您使用 Age 来索引 alpha。我认为您真正想要的是 Age 的每个值都有一个不同的 alpha。我通过在数据中执行以下操作解决了这个问题：Ages=match(D[,1],unique(D[,1])),
2. alpha 被重新定义，因为 N 上的循环同时包含了 N.c 和 S 上的循环。通过在 N 上的循环开始之前关闭 S 上的循环，可以解决问题。
3. GDP 在数据中定义为 mat，因此我将 mat = D[,5] 替换为 GDP = D[,5]。

进行这些更改后，模型将运行。

library(tidyverse)

Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)

# X2 AGES
n.ages <- length(Ages)

# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)

# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)


Year.tot <- rep(Year,length(unique(Ages)))

DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),C1=C_Porp[,C2=C_Porp[,2],C3=C_Porp[,3]) 

# X1 GDP OVER YEARS
GDP <- rnorm(n.year,15)
GDP <- as.data.frame(cbind(GDP,Year))

# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)



library(runjags)
dirlichet.model = 
  
  
  "model {
  #setup priors for each species
  for(j in 1:N.c){
    m0[j] ~ dnorm(0,1.0E-3) #intercept prior
    m1[j] ~ dnorm(0,1.0E-3) #      mat prior
  }
  
  #implement dirlichet
  for(i in 1:N){
  
    y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 
     
    for(j in 1:N.c){
    
      a0[i,j] ~ dnorm(mu[i,j],0.001) 
      
      log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])
    
    }
  }
     # Effect of site:
    for (s in 1:S){
      
      alpha[s] ~ dnorm(0,0.001)
      
    }
    
  }
"

jags.data <- list(y = D[,c(2,3,4)],GDP = D[,5],Ages=match(D[,N = nrow(D[,c(3,4,5)]),N.c = ncol(D[,S=(length(unique(D[,1]))))

jags.out <- run.jags(dirlichet.model,data=jags.data,adapt = 500,burnin = 5000,sample = 10000,n.chains=5,monitor=c('m0','m1',"mu"))
out <- summary(jags.out)