0

我有以下数据结构：

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

这里的型号规格：

模型

这里模拟数据：

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,50,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[,1],N = nrow(D[,c(3,4,5)]), 
                  N.c = ncol(D[,c(3,4,5)]),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)

我收到了这个错误：

Errore: The following error occured when compiling and adapting the model using rjags:
 Error in rjags::jags.model(model, data = dataenv, n.chains = length(runjags.object$end.state),  : 
  RUNTIME ERROR:
Compilation error on line 12.
Index out of range taking subset of  alpha

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

4

1 回答 1

1

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

  1. Age被定义为从 0 到 50 由 5 组成的序列，但随后您使用Age索引alpha. 我认为您真正想要的是每个Age. 我通过在数据中执行以下操作解决了这个问题：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[,1],
                   C2=C_Porp[,2],
                   C3=C_Porp[,3]) 

# X1 GDP OVER YEARS
GDP <- rnorm(n.year,50,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[,1], unique(D[,1])),
                  N = nrow(D[,c(3,4,5)]), 
                  N.c = ncol(D[,c(3,4,5)]),
                  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)
于 2021-01-13T21:11:22.500 回答