|
使用以下程序（附件为该程序）做统计时，OpenBUGS系统报错“array index is greater than array bound for y”
求大神指点啊，新手第一次提问~
#Random effects model for multi-arm trials
# data are y=a table of the arm-means, sd=a table of the arm sd, N=a table of the arm sample size, t=a table with the names (numbers) of treatments, na=a vector with the number of arms in each study, ref=a number specifying which is the reference treatment
model{
for(i in 1:ns){
w[i,1] <-0
delta[i,t[i,1]]<-0
u[i] ~ dnorm(0,.0001)
for (k in 1:na[i]) {
se[i,t[i,k]]<-sd[i,t[i,k]]/sqrt(n[i,t[i,k]])
var[i,t[i,k]]<-se[i,t[i,k]]*se[i,t[i,k]]
prec[i,t[i,k]]<-1/var[i,t[i,k]]
#normal likelihood
y[i,t[i,k]]~dnorm(phi[i,t[i,k]],prec[i,t[i,k]])
phi[i,t[i,k]]<-(u[i]+delta[i,t[i,k]])*pooled.sd[i]
#calculate the pooled SD
nom1[i,k]<-n[i,t[i,k]]*sd[i,t[i,k]]*sd[i,t[i,k]] #nominator for the pooled sd
}
# model
ss[i]<-sum(n[i,1:nt])-nt+na[i] #total sample size in a study
nom[i]<-sum(nom1[i,1:na[i]])#nominator for the pooled sd
pooled.sd[i]<-sqrt(nom[i]/(ss[i]-na[i]))# pooled sd
for (k in 2:na[i]) {
delta[i,t[i,k]] ~ dnorm(md[i,t[i,k]],taud[i,t[i,k]]) # trial-specific SMD distributions
md[i,t[i,k]] <- d[t[i,k]] - d[t[i,1]] + sw[i,k] # mean of SMD distributions
taud[i,t[i,k]] <- tau *2*(k-1)/k #precision of SMD 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 for multi-arm trials
}
d[ref]<-0
for (k in 1:(ref-1)){d[k] ~ dnorm(0,.0001) }
for (k in (ref+1):nt){d[k] ~ dnorm(0,.0001) }
SD~dunif(0,5) # vague prior for random effects standard deviation
tau<-1/pow(SD,2)
# Collection of results###########
# pairwise SMDs
# for all comparisons
for (c in 1:(nt-1)) { for (k in (c+1):nt) { SMD[c,k] <- d[c] - d[k] } }
#compared to baseline
for (c in 1:nt) { SMD.ref[c] <-d[c] - d[ref] }
#predictions
for (c in 1:(ref-1)) { X[c]<-d[c] - d[ref]
predSMD.ref[c] ~dnorm( X[c],tau)}
for (c in (ref+1):nt) { X[c]<-d[c] - d[ref]
predSMD.ref[c] ~dnorm( X[c],tau)}
for (c in 1:(nt-1)) { for (k in (c+1):nt) { predSMD[c,k] ~ dnorm( SMD[c,k],tau) } }
for(k in 1:nt) {
order[k]<- rank(d[],k)
# this is when the outcome is positive - omit 'nt+1-' when the outcome is negative
most.effective[k]<-equals(order[k],1)
for(j in 1:nt) {
effectiveness[k,j]<- equals(order[k],j)
}
}
for(k in 1:nt) {
for(j in 1:nt) {
cumeffectiveness[k,j]<- sum(effectiveness[k,1:j])
}
}
#SUCRAS#
for(k in 1:nt) {
SUCRA[k]<- sum(cumeffectiveness[k,1:(nt-1)]) /(nt-1)
}
#Fit of the Model#
for(i in 1:ns) {
for(k in 1:na[i]) {
Darm[i,k]<-(y[i,t[i,k]]-phi[i,t[i,k]])*(y[i,t[i,k]]-phi[i,t[i,k]])/var[i,t[i,k]]
}
D[i]<- sum(Darm[i,1:na[i]])
}
D.bar<- sum(D[])
}
list(ns=36, nt=15, ref=1)
t[,1] t[,2] t[,3] y[,1] y[,2] y[,3] sd [,1] sd [,2] sd [,3] n [,1] n [,2] n [,3] na[]
1 2 NA -2.5 -5.9 NA 2.38 3.34 NA 10 10 NA 2
1 2 NA -4.4 -4 NA 4.31 4.71 NA 13 18 NA 2
1 3 NA -16.5 -21.7 NA 14.75 15.09 NA 85 89 NA 2
1 3 NA -11.91 -12.09 NA 15.93 16.36 NA 91 96 NA 2
4 13 NA -13.3 -17.8 NA 9.94 8.40 NA 58 63 NA 2
1 5 NA -6.72 -10.61 NA 4.53 3.21 NA 18 18 NA 2
1 5 NA -10.35 -9.95 NA 4.32 4.61 NA 25 17 NA 2
1 6 8 -21.6 -24.50 -22.6 13.74 13.86 13.44 117 219 112 3
1 6 8 -24.3 -24.3 -23.7 11.27 11.59 11.27 103 113 113 3
1 7 NA -20.2 -21.9 NA 13.50 13.50 NA 132 129 NA 2
1 7 NA -18.8 -22.1 NA 15.91 15.14 NA 157 154 NA 2
8 12 NA -10.95 -2.6 NA 2.61 0.8 NA 20 20 NA 2
1 8 NA -19.41 -22.64 NA 6.93 7.08 NA 112 109 NA 2
1 8 NA -8.14 -9.14 NA 4.54 4.02 NA 9 7 NA 2
1 8 NA -12.42 -13.91 NA 7.24 7.84 NA 19 21 NA 2
1 8 NA -10.5 -20.1 NA 14.84 13.19 NA 48 48 NA 2
1 8 NA -14.9 -22 NA 12.74 12.11 NA 105 109 NA 2
1 8 NA -22.63 -18.40 NA 12.48 12.47 NA 16 18 NA 2
8 15 NA -15 -15 NA 6 6 NA 30 29 NA 2
1 9 13 -9.09 -8.91 -10.74 6.71 6.80 6.66 87 94 90 3
1 9 NA -6.4 -11.5 NA 10.4 14.7 NA 7 9 NA 2
1 9 NA -1.1 -1.2 NA 0.78 0.60 NA 22 16 NA 2
10 13 NA -23.1 -23.3 NA 5.37 5.89 NA 26 26 NA 2
1 10 NA -14.69 -15.85 NA 12.55 12.54 NA 44 82 NA 2
1 10 NA -8.81 -13.48 NA 12.22 12.24 NA 41 83 NA 2
1 11 NA -22.5 -26.5 NA 14.46 14.46 NA 96 99 NA 2
1 11 NA -21.6 -21.74 NA 13.68 14.04 NA 93 180 NA 2
1 12 NA -17.6 -17 NA 8.49 9.99 NA 24 26 NA 2
1 12 NA -13.6 -16.6 NA 7.93 6.77 NA 19 12 NA 2
1 13 NA -12.80 -13.60 NA 10.30 10.91 NA 91 177 NA 2
1 13 NA -23.38 -22.58 NA 16 14.77 NA 100 101 NA 2
1 13 NA -11.9 -16.5 NA 13.20 13.19 NA 27 29 NA 2
1 14 NA -22.1 -25.9 NA 14.71 14.71 NA 88 93 NA 2
1 14 NA -25.6 -28.8 NA 15.71 15.71 NA 91 92 NA 2
1 15 NA -16.1 -18.1 NA 12.90 12.86 NA 73 68 NA 2
1 15 NA -22.6 -24.3 NA 13.52 12.66 NA 92 101 NA 2
END
#####初始值####
#####第一条链###
list(SD = 1,d = c(NA,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
u = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
####第二条链####
list(SD = 0.5,d = c(NA,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
u = c(-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3))
|
jg-xs1
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
发帖
雷达卡
