Функция ранга в зазубринах - PullRequest
Новая
       36

Функция ранга в зазубринах

1 голос
/ 11 апреля 2019

У меня проблема с тем, чтобы заставить работать функцию 'rank' в Jags. Ниже приведен код, модель и данные, адаптированные из Winbugs. Точка, которая не работает: 

  myorder[i] <-  rank(aux.u[],i)

, который выдает следующую ошибку: 

  RUNTIME ERROR:
  Incorrect number of parameters in function rank

Мне известно, что у Джагса ранг функции отличается. Итак, когда я пытаюсь перекодировать его, используя что-то вроде: 

  myorderx <-  rank(aux.u[]);myorder[i]<-myorderx[i];

Я получаю ошибку: 

  Attempt to redefine node myorderx[1:4]

Любой совет приветствуется. Большое спасибо, Лисо

Модель 

rm(list = ls())
library(coda)
library(rjags)
library(pracma);
library(nmathresh)

setwd("/Users/test")

nummodel <- function(){

  for(i in 1:NS)  {
    w[i,1] <- 0
    delta[i,t[i,1]] <- 0
    mu[i] ~ dnorm(0,.0001)                        # Fixed study effects
    for (k in 1:na[i])   {
      r[i,k] ~ dbin(p[i,t[i,k]],n[i,k])                # Binomial likelihood for data
      logit(p[i,t[i,k]])<-mu[i] + delta[i,t[i,k]] }    # Model for log odds parameters
    for (k in 2:na[i])   {
      delta[i,t[i,k]] ~ dnorm(md[i,t[i,k]],precd[i,t[i,k]])   # Random treatment effects
      md[i,t[i,k]] <-  d[t[i,k]] - d[t[i,1]]  + sw[i,k]                   # Mean of LOR distributions
      precd[i,t[i,k]] <- prec[t[i,1],t[i,k]]*2*(k-1)/k     # Precision of LOR distributions
      w[i,k] <- (delta[i,t[i,k]]  - d[t[i,k]] + d[t[i,1]])           # Adjustment, multi-arm RCTs
      sw[i,k] <-sum(w[i,1:(k-1)])/(k-1) }        # Cumulative adjustment, multi-arm RCTs
  }

  d[1]<-0
  for (k in 2:NT)  {d[k] ~ dnorm(0,.0001) }    # Vague priors for basic parameters

  for(j in 1:3){
    prec[j,j]<-1
    for(k in (j+1):4){
      prec[j,k]<-1/tausq[j,k]
      prec[k,j]<-prec[j,k] }}

  for(k in 1:4){
    v.a[k]~dlnorm(-3.31,0.346)    # Informative prior chosen to imply (approximately) chosen data-based prior for contrast heterogeneity variance
    sd.a[k]<-sqrt(v.a[k])
  } 

  pi.half<-1.5708

  for(i in 1:3) {for(j in (i+1):4){
    g[j,i]<-0
    tausq[i,j]<-v.a[i]+v.a[j]-2*rho.star[i,j]*sd.a[i]*sd.a[j]
    tau[i,j] <- sqrt(tausq[i,j]) }}

  # Implementing random permutation
  for(i in 1:4) { for(j in 1:4) {
    rho[i,j]<-inprod(g[, i],g[, j])   
    rho.star[i,j]<- rho[ myorder[i], myorder[j] ]    }}

  for(i in 1:4) {  
    aux.u[i] ~ dunif(0, 1) 
    myorder[i] <-  rank(aux.u[],i)
    #myorderx <-  rank(aux.u[]);myorder[i]<-myorderx[i];
  }

  # Constructing entries of upper-triangular matrix for Cholesky decomposition
  g[1,1]<-1
  g[2,2]<-sin.a[1,2]
  g[3,3]<-sin.a[1,3]*sin.a[2,3]
  g[4,4]<-sin.a[1,4]*sin.a[2,4]*sin.a[3,4]
  g[1,2]<-cos.a[1,2]
  g[1,3]<-cos.a[1,3]
  g[1,4]<-cos.a[1,4]
  g[2,3]<-sin.a[1,3]*cos.a[2,3]
  g[2,4]<-sin.a[1,4]*cos.a[2,4]
  g[3,4]<-sin.a[1,4]*sin.a[2,4]*cos.a[3,4]

  # Beta prior for cos(a[i,j])
  for (i in 1:3) {
    for (j in (i+1):4)  {
      cos.a[i,j] ~ dbeta(0.93,1.07)       #  From Table 3
      sin2.a[i,j] <- 1 - pow(cos.a[i,j] , 2)
      sin.a[i,j]  <- pow(sin2.a[i,j] , 1/2)    }  }

  # Defining pairwise ORs
  for (c in 1:(NT-1))
  {  for (k in (c+1):NT)  
  {  lor[c,k] <- d[k] - d[c]
  log(or[c,k]) <- lor[c,k]   }    
  }
}
write.model(nummodel, "Model.txt")
model.file1 = paste(getwd(),"Model.txt", sep="/")

Определение некоторых параметров MCMC для JAGS 

nchains  <- 3; # How Many Chains?
nadapt<-200
nburnin  <- 200; # How Many Burn-in Samples?
nsamples <- 300;  # How Many Recorded Samples?
nthin <- 20;

Данные 

r1<-c(9,11,75,2,58,0,3,1,6,79,18,64,5,20,0,8,95,15,78,69,20,7,12,9)
n1<-c(140,78,731,106,549,33,100,31,39,702,671,642,62,234,20,116,1107,187,584,1177,49,66,76,55)
r2<-c(23,12,363,9,237,9,31,26,17,77,21,107,8,34,9,19,143,36,73,54,16,32,20,3)
n2<-c(140,85,714,205,1561,48,98,95,77,694,535,761,90,237,20,149,1031,504,675,888,43,127,74,26)
r3<-c(10,29,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
n3<-c(138,170,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
t1<-c(1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,3,3)
t2<-c(3,3,3,3,3,3,3,3,3,2,2,3,3,3,4,2,3,3,3,3,3,4,4,4)
t3<-c(4,4,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
na<-c(3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2)

datasrc <- data.frame(r1,r2,r3,n1,n2,n3,t1,t2,t3,na)

datastruct<-list(NS=length(datasrc$t1),NT=18,na=c(datasrc$na),
                 t=structure(.Data=c(datasrc$t1,datasrc$t2,datasrc$t3), .Dim=c(length(datasrc$t1),3)),
                 r=structure(.Data=c(datasrc$r1,datasrc$r2,datasrc$r3), .Dim=c(length(datasrc$t1),3)),
                 n=structure(.Data=c(datasrc$n1,datasrc$n2,datasrc$n3), .Dim=c(length(datasrc$t1),3)))

Запустить модель 

parameters = c("lor");
model.file1 = "/Users/test/Model.txt"
mod1 <- jags.model(file =model.file1, data=datastruct, n.chains=nchains, n.adapt=nadapt);

1 Ответ

1 голос
/ 12 апреля 2019

Похоже, проблема в том, что myorderx вложено в цикл for.В результате myorderx записывается один раз, а затем записывается более трех раз, поскольку i находится в диапазоне от 1 до 4. Все, что вам нужно сделать, это вывести его из цикла for: 

for(i in 1:4) {  
 aux.u[i] ~ dunif(0, 1)
}
myorder <- rank(aux.u[])

Помните, что функция ранга получает и возвращает вектор, вы в настоящее время задаете ему вектор и пытаетесь вернуть результат в скаляр.В приведенном здесь коде myorder имеет те же размеры, что и aux.u, что означает, что оба являются векторами.

...