1
初めての投稿ですので、すべてのボックスにチェックを入れてください。私はRICEとR2WinBUGSを介してNICEのテクニカルサポートの文書で提供されているコードを使って、ネットワークのメタ解析を行っています。モデルは平均差を使ってうまく機能しますが、ヘッジGに変換すると、私は突然、非正定値共分散行列の結果であると思われる未定義の実際の誤差を得ます。私は明らかな何かを見逃していると確信していますが、可能な助けを感謝します。以下のコードをご覧ください。WinBUGSでMDからSMDに切り替えると、未定義の実際のトラップエラーが発生する
library(R2WinBUGS)
re_normal_gaus <- function() # this code for this model was adapted from WinBUGS code from the multi-parameter Evidence Synthesis Research Group at the University of Bristol: Website: www.bris.ac.uk/cobm/research/mpes
{
for(i in 1:ns2) { # LOOP THROUGH 2-ARM STUDIES
y[i,2] ~ dnorm(delta[i,2],prec[i,2]) # normal likelihood for 2-arm trials
resdev[i] <- (y[i,2]-delta[i,2])*(y[i,2]-delta[i,2])*prec[i,2] #Deviance contribution for trial i
}
for(i in (ns2+1):(ns2+ns3)) { # LOOP THROUGH 3-ARM STUDIES
for (k in 1:(na[i]-1)) { # set variance-covariance matrix
for (j in 1:(na[i]-1)) { Sigma[i,j,k] <- V[i]*(1-equals(j,k)) + var[i,k+1]*equals(j,k) }
}
Omega[i,1:(na[i]-1),1:(na[i]-1)] <- inverse(Sigma[i,,]) #Precision matrix
# multivariate normal likelihood for 3-arm trials
y[i,2:na[i]] ~ dmnorm(delta[i,2:na[i]],Omega[i,1:(na[i]-1),1:(na[i]-1)])
#Deviance contribution for trial i
for (k in 1:(na[i]-1)){ # multiply vector & matrix
ydiff[i,k]<- y[i,(k+1)] - delta[i,(k+1)]
z[i,k]<- inprod2(Omega[i,k,1:(na[i]-1)], ydiff[i,1:(na[i]-1)])
}
resdev[i]<- inprod2(ydiff[i,1:(na[i]-1)], z[i,1:(na[i]-1)])
}
for(i in (ns2+ns3+1):(ns2+ns3+ns4)) { # LOOP THROUGH 4-ARM STUDIES
for (k in 1:(na[i]-1)) { # set variance-covariance matrix
for (j in 1:(na[i]-1)) { Sigma2[i,j,k] <- V[i]*(1-equals(j,k)) + var[i,k+1]*equals(j,k) }
}
Omega2[i,1:(na[i]-1),1:(na[i]-1)] <- inverse(Sigma2[i,,]) #Precision matrix
# multivariate normal likelihood for 4-arm trials
y[i,2:na[i]] ~ dmnorm(delta[i,2:na[i]],Omega2[i,1:(na[i]-1),1:(na[i]-1)])
#Deviance contribution for trial i
for (k in 1:(na[i]-1)){ # multiply vector & matrix
ydiff[i,k]<- y[i,(k+1)] - delta[i,(k+1)]
z[i,k]<- inprod2(Omega2[i,k,1:(na[i]-1)], ydiff[i,1:(na[i]-1)])
}
resdev[i]<- inprod2(ydiff[i,1:(na[i]-1)], z[i,1:(na[i]-1)])
}
for(i in 1:(ns2+ns3+ns4)){ # LOOP THROUGH ALL STUDIES
w[i,1] <- 0 # adjustment for multi-arm trials is zero for control arm
delta[i,1] <- 0 # treatment effect is zero for control arm
for (k in 2:na[i]) { # LOOP THROUGH ARMS
var[i,k] <- pow(se[i,k],2) # calculate variances
prec[i,k] <- 1/var[i,k] # set precisions
dev[i,k] <- (y[i,k]-delta[i,k])*(y[i,k]-delta[i,k])*prec[i,k]
}
#resdev[i] <- sum(dev[i,2:na[i]]) # summed residual deviance contribution for this trial
for (k in 2:na[i]) { # LOOP THROUGH ARMS
delta[i,k] ~ dnorm(md[i,k],taud[i,k]) # trial-specific treat effects distributions
md[i,k] <- d[t[i,k]] - d[t[i,1]] + sw[i,k] # mean of treat effects distributions (with multi-arm trial correction)
taud[i,k] <- tau *2*(k-1)/k # precision of treat effects distributions (with multi-arm trial correction)
w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]]) # adjustment for multi-arm RCTs
sw[i,k] <- sum(w[i,1:k-1])/(k-1) # cumulative adjustment for multi-arm trials
}
}
totresdev <- sum(resdev[]) #Total Residual Deviance
d[1]<-0 # treatment effect is zero for reference treatment
for (k in 2:nt){ d[k] ~ dnorm(0,.0001) } # vague priors for treatment effects
sd ~ dunif(0,5) # vague prior for between-trial SD
tau <- pow(sd,-2) # between-trial precision = (1/between-trial variance)
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
# Extra code for all mean differences, rankings
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
# pairwise mean differences for all possible pair-wise comparisons, if nt>2
for (c in 1:(nt-1)) {
for (k in (c+1):nt) {
meandif[c,k] <- (d[k] - d[c])
# pairwise comparison between all treatments
better[c,k] <- 1 - step(d[k] - d[c])
}
}
# ranking calculations
for (k in 1:nt) {
# assumes differences<0 favor the comparator ===> number of elements in d[] that are less than or equal to d[k]
rk[k] <- rank(d[],k)
#calculate probability that treat k is best ===> k is the best treatment if rk[k] = 1
best[k] <- equals(rk[k],1)
for(h in 1:nt) {
prob[k,h]<- equals(rk[k],h)
}
}
for(k in 1:nt) {
for(h in 1:nt) {
cumeffectiveness[k,h]<- sum(prob[k,1:h])
}
}
#SUCRAS#
for(i in 1:nt) {
SUCRA[i]<- sum(cumeffectiveness[i,1:(nt-1)]) /(nt-1)
}
} # END Program
write.model(re_normal_gaus, "re-normal-gaus.txt")
MODELFILE.re <- c("re-normal-gaus.txt")
md = read.table(textConnection("t_1 t_2 t_3 t_4 y_2 y_3 y_4 se_2 se_3 se_4 V na
1 5 NA NA 1.6 NA NA 1.7801767 NA NA 1.381371651 2
1 8 NA NA -0.2 NA NA 0.365808405 NA NA 0.075789474 2
2 3 NA NA 2.7 NA NA 0.48894018 NA NA 0.1378125 2
1 4 NA NA -2.2 NA NA 1.255423019 NA NA 0.695652174 2
2 3 NA NA 0 NA NA 1.060660172 NA NA 0.5625 2
2 3 NA NA 2.7 NA NA 0.555150008 NA NA 0.154285714 2
1 7 NA NA -2 NA NA 0.441189176 NA NA 0.039192632 2
1 2 NA NA -0.5 NA NA 0.970562384 NA NA 0.457142857 2
4 9 NA NA 0.1 NA NA 0.807757815 NA NA 0.188072678 2
1 5 NA NA 2.5 NA NA 1.011075035 NA NA 0.465454545 2
2 3 NA NA 2.6 NA NA 0.921954446 NA NA 0.49 2
1 10 NA NA -1.3 NA NA 0.850881895 NA NA 0.4 2
2 3 NA NA 2.6125 NA NA 0.439282085 NA NA 0.14546875 2
1 4 NA NA 2 NA NA 1.306118643 NA NA 0.989117513 2
1 6 NA NA -2.6 NA NA 0.550381686 NA NA 0.15842 2
1 4 NA NA -3.2 NA NA 0.194924242 NA NA 0.0162 2
1 2 NA NA -3.6 NA NA 0.380131556 NA NA 0.07225 2
5 11 12 NA -2.4 -2.2 NA 0.79304967 0.848808689 NA 0.347142857 3
1 4 7 NA -4.7 -0.8 NA 0.420823389 0.45845681 NA 0.065641026 3
2 3 6 NA 0.3 0.48 NA 0.584636639 0.703153611 NA 0.20402 3
1 2 3 4 -3.2 -3 -1 1.219830171 1.113092025 0.785493475 0.361 4
"),header = T)
smd = read.table(textConnection("t_1 t_2 t_3 t_4 y_2 y_3 y_4 se_2 se_3 se_4 V na
1 5 NA NA 0.284314186 NA NA 0.387759426 NA NA NA 2
1 8 NA NA -0.088247814 NA NA 0.424250037 NA NA NA 2
2 3 NA NA 1.363769324 NA NA 0.621067031 NA NA NA 2
1 4 NA NA -0.507895157 NA NA 0.404635828 NA NA NA 2
2 3 NA NA 0 NA NA 0.702664963 NA NA NA 2
2 3 NA NA 1.5093214 NA NA 0.91184626 NA NA NA 2
1 7 NA NA -0.843577542 NA NA 0.461709287 NA NA NA 2
1 2 NA NA -0.123574144 NA NA 0.383877128 NA NA NA 2
4 9 NA NA 0.038564349 NA NA 0.495157679 NA NA NA 2
1 5 NA NA 0.732129077 NA NA 0.451400734 NA NA NA 2
2 3 NA NA 1.001922271 NA NA 0.553706463 NA NA NA 2
1 10 NA NA -0.654394486 NA NA 1.502905599 NA NA NA 2
2 3 NA NA 1.843306589 NA NA 1.078204274 NA NA NA 2
1 4 NA NA 0.553077868 NA NA 0.428210865 NA NA NA 2
1 6 NA NA -1.464178893 NA NA 1.155770916 NA NA NA 2
1 4 NA NA -3.346020348 NA NA 1.497631686 NA NA NA 2
1 2 NA NA -2.097219595 NA NA 0.613586425 NA NA NA 2
5 11 12 NA -0.906267636 -0.784775956 NA 0.535679681 0.518303517 NA 0.347142857 3
1 4 7 NA -2.488768168 -0.386547304 NA 0.734950178 0.595461098 NA 0.065641026 3
2 3 6 NA 0.15904499 0.211580576 NA 1.013288527 0.858087687 NA 0.20402 3
1 2 3 4 -1.100360888 -1.182907814 -0.545284283 1.059726917 1.270240674 1.588768054 0.361 4
"),header = T)
treatments = read.table(textConnection("description t
A 1
B 2
C 3
D 4
E 5
F 6
G 7
H 8
I 9
J 10
K 11
L 12
"),header = T)
t = as.matrix(md[1:4])
# number of treatments
nt = length(treatments[[1]])
y = as.matrix(cbind(rep(NA, length(t[,1])), md[5:7]))
se = as.matrix(cbind(rep(NA, length(t[,1])), md[8:10]))
na = as.vector(md[,12])
V = as.vector(md[,11])
ns2 = length(subset(na, na==2))
ns3 = length(subset(na, na==3))
ns4 = length(subset(na, na==4))
data_md = list(nt=nt, ns2=ns2, ns3=ns3, ns4=ns4,t=t, y=y, se=se, na=na, V=V)
y_smd = as.matrix(cbind(rep(NA, length(t[,1])), smd[5:7]))
se_smd = as.matrix(cbind(rep(NA, length(t[,1])), smd[8:10]))
data_smd = list(nt=nt, ns2=ns2, ns3=ns3, ns4=ns4,t=t, y=y_smd, se=se_smd, na=na, V=V)
params = c("meandif", 'SUCRA', 'best', 'totresdev', 'rk', 'dev', 'resdev', 'prob', "better","sd")
bugs(data_md, NULL, params, model.file= MODELFILE.re,
n.chains = 3, n.iter = 40000, n.burnin = 20000, n.thin=1,
bugs.directory = bugsdir, debug=F)
bugs(data_smd, NULL, params, model.file= MODELFILE.re,
n.chains = 3, n.iter = 40000, n.burnin = 20000, n.thin=1,
bugs.directory = bugsdir, debug=F)