First time posting here so I hope I have been able to tick all the boxes for requirements. I am running a network meta-analysis using code provided by the NICE technical support documents via R and R2WinBUGS. Model works well using mean difference, but when I go to convert to hedges G I suddenly get an undefined real error that seems to be the result of a non positive definite covariance matrix. I'm sure I am missing something obvious, but would appreciate any help possible. Please see code below.
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)
After speaking with the author of the original code, I learned that the issue was in how I calculated V for the SMD. Covariance of the md is just variance in the control arm, but for smd it is ~ the inverse of the sample size in the control arm. Once this was calculated correctly my matrixes behaved.
Source: http://methods.cochrane.org/sites/methods.cochrane.org.cmi/files/public/uploads/S8-L%20Problems%20introduced%20by%20multi-arm%20trials%20-%20full%20network%20meta-analysis.pdf