WinBUGS Weibull Netzwerk Meta-Analyse

Question

Ich arbeite derzeit an einer Meta-Analyse von Überlebensdaten über mehrere klinische Studien.

Um dies zu tun, habe ich Code aus einer veröffentlichten Analyse mit der gleichen Methode. Wenn ich diesen Code mit den Daten der veröffentlichten Analyse ausführe, kann ich die Ergebnisse nicht replizieren. Tatsächlich stimmen die Ergebnisse nicht mit irgendeiner vernünftigen Schätzung überein.

Der Code selbst (ohne die Daten) sollte korrekt sein, da er direkt von den Autoren stammt. Ich nehme an, das Problem hat w/initial Werte oder Parameter, wie die Probenahme ausgeführt wird, aber nach dem Spielen w/viele Anfangswerte, die Länge der Burn in, Ausdünnung, etc ... Ich habe keine aussagekräftigen Ergebnisse erhalten .

Ich würde mich freuen, jemand Vorschläge über das Ausführen (Anfangswerte, etc ...), damit es richtig ausgeführt wird. Wenn es im Code Probleme gibt oder wenn die Daten so eingerichtet sind, dass sie nicht mit dem Code übereinstimmen, wäre es hilfreich, dies zu wissen.

Als eine Randnotiz, mache ich die Analysen mit R2WinBUGs, obwohl ich die gleiche Art von Problemen mit WinBUGs allein bekommen haben.

ein bisschen mehr Hintergrund auf der Methode:

Der Weg, um dies durch Abschätzen des Unterschiedes in der Form und Größe zwischen Behandlungen über mehrere Studien unter Verwendung von Zufallseffekten Parametern einer reparametrisierten Weibull-Verteilung funktioniert, ist.

Die Weibull-Verteilung ist so umparametriert, dass der Logarithmus der Hazard Rate ein + b * log (t) ist, wobei a ein Skalierungsparameter und b ein Shape-Parameter ist. Daraus können Sie die Wahrscheinlichkeit Funktion einer bestimmten Anzahl von Fehlern aus einer bestimmten Anzahl von Patienten über ein Intervall berechnen.

Leider ist der Artikel öffentlich, aber wenn Sie den Zugriff hier ist der Link zu erreichen: http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=

Eine kurze Zusammenfassung der Variablen in das Modell eingegeben:

NT: Anzahl der separaten Behandlungen enthalten.

N: Anzahl der Zeilen im Hauptdatensatz. NS: Anzahl von Studien

s: Studie, dass die Datenzeile entspricht (dies ist nummeriert 1: 6)

r: Anzahl der Patienten in Intervall für diese Behandlung versagt/Studie

n : bei Risikopatienten zu Beginn der Intervallzahl der für diese Behandlung /Studie

t: Behandlung, die diese Datenzeile entspricht (nummeriert von 1: 3)

b: Gibt an, welche treatme nt ist die Baseline, mit der andere verglichen werden (für jede Zeile auf 1 gesetzt).

bs: Behandlung, die der Steuerarm dieser Studie ist

bt: Behandlung, die die Forschung Arm dieser Studie ist

WinBUGS Code (einschließlich Daten):

#Winbugs code for random effects networks meta-analysis model 
Model 
{ 
    for (i in 1:N) 
    { # N=number of data points in dataset 
    #likelihood 
    r[i]~ dbin(p[i],n[i]) 
    p[i]<-1-exp(-h[i]*dt[i]) # hazard h over interval [t,t+dt] # expressed as deaths per unit person-time (e.g. months) 
    #random effects model 
    log(h[i])<-nu[i]+log(time[i])*theta[i] 
    nu[i]<-mu[s[i],1]+delta[s[i],1]*(1-equals(t[i],b[i])) 
    theta[i]<-mu[s[i],2]+ delta[s[i],2]*(1-equals(t[i],b[i])) 
    } 
    for(k in 1 :NS) 
    { # NS=number of studies in dataset 
    delta[k,1:2]~dmnorm(md[k,1:2],omega[1:2,1:2]) 
    md[k,1]<-d[ts[k],1]-d[bs[k],1] 
    md[k,2]<-d[ts[k],2]-d[bs[k],2] 
    } 
    # priors 
    d[1,1]<-0 
    d[1,2]<-0 
    for(j in 2 :NT) 
    { # NT=number of treatments 
    d[j,1:2] ~ dmnorm(mean[1:2],prec2[,]) 
    } 
    for(k in 1 :NS) 
    { 
    mu[k,1:2] ~ dmnorm(mean[1:2],prec2[,]) 
    } 
    omega[1:2, 1:2] ~ dwish(R[1:2,1:2],2) 
} 
# Winbugs data set 
list(N=242, NS=6, NT=3, 
mean=c(0,0), 
prec2 = structure(.Data = c(
0.0001,0, 
0,0.0001), .Dim = c(2,2)), 
R = structure(.Data = c(
0.01,0, 
0,0.01), .Dim = c(2,2)) 
) 

s[] r[] n[] t[] b[] time[] dt[] 
1 15 152 3 1 3 3 
1 11 140 3 1 6 3 
1 8 129 3 1 9 3 
1 9 121 3 1 12 3 
1 9 112 3 1 15 3 
1 3 83 3 1 18 3 
1 4 80 3 1 21 3 
1 5 76 3 1 24 3 
1 2 71 3 1 27 3 
1 2 41 3 1 30 3 
1 1 39 3 1 33 3 
1 3 38 3 1 36 3 
1 2 35 3 1 39 3 
1 1 33 3 1 42 3 
1 3 32 3 1 45 3 
1 3 29 3 1 48 3 
1 2 26 3 1 51 3 
1 1 24 3 1 54 3 
1 1 23 3 1 57 3 
1 1 22 3 1 60 3 
1 10 149 1 1 3 3 
1 11 140 1 1 6 3 
1 9 128 1 1 9 3 
1 5 119 1 1 12 3 
1 6 114 1 1 15 3 
1 3 72 1 1 18 3 
1 5 70 1 1 21 3 
1 4 65 1 1 24 3 
1 7 61 1 1 27 3 
1 2 34 1 1 30 3 
1 2 32 1 1 33 3 
1 3 30 1 1 36 3 
1 2 27 1 1 39 3 
1 2 25 1 1 42 3 
1 1 23 1 1 45 3 
1 2 22 1 1 48 3 
1 1 19 1 1 51 3 
1 2 19 1 1 54 3 
1 1 17 1 1 57 3 
1 0 16 1 1 60 3 
2 4 125 2 1 3 3 
2 4 121 2 1 6 3 
2 2 117 2 1 9 3 
2 5 114 2 1 12 3 
2 2 109 2 1 15 3 
2 3 107 2 1 18 3 
2 2 104 2 1 21 3 
2 4 94 2 1 24 3 
2 4 90 2 1 27 3 
2 3 81 2 1 30 3 
2 4 78 2 1 33 3 
2 3 61 2 1 36 3 
2 5 58 2 1 39 3 
2 1 48 2 1 42 3 
2 2 47 2 1 45 3 
2 3 41 2 1 48 3 
2 0 38 2 1 51 3 
2 3 29 2 1 54 3 
2 3 26 2 1 57 3 
2 2 18 2 1 60 3 
2 0 16 2 1 63 3 
2 1 10 2 1 66 3 
2 0 9 2 1 69 3 
2 0 3 2 1 72 3 
2 0 3 2 1 75 3 
2 0 3 2 1 78 3 
2 15 196 1 1 3 3 
2 9 179 1 1 6 3 
2 10 170 1 1 9 3 
2 9 162 1 1 12 3 
2 9 153 1 1 15 3 
2 5 141 1 1 18 3 
2 5 136 1 1 21 3 
2 10 121 1 1 24 3 
2 5 111 1 1 27 3 
2 7 92 1 1 30 3 
2 7 85 1 1 33 3 
2 4 71 1 1 36 3 
2 6 67 1 1 39 3 
2 4 53 1 1 42 3 
2 5 49 1 1 45 3 
2 6 36 1 1 48 3 
2 3 30 1 1 51 3 
2 2 26 1 1 54 3 
2 2 24 1 1 57 3 
2 0 13 1 1 60 3 
2 1 13 1 1 63 3 
2 1 11 1 1 66 3 
2 1 10 1 1 69 3 
2 0 6 1 1 72 3 
2 0 6 1 1 75 3 
2 0 6 1 1 78 3 
3 6 113 2 1 3 3 
3 4 105 2 1 6 3 
3 3 101 2 1 9 3 
3 1 97 2 1 12 3 
3 9 96 2 1 15 3 
3 4 84 2 1 18 3 
3 2 80 2 1 21 3 
3 4 74 2 1 24 3 
3 3 70 2 1 27 3 
3 2 59 2 1 30 3 
3 0 57 2 1 33 3 
3 6 51 2 1 36 3 
3 2 45 2 1 39 3 
3 1 37 2 1 42 3 
3 3 36 2 1 45 3 
3 1 27 2 1 48 3 
3 1 26 2 1 51 3 
3 2 25 2 1 54 3 
3 7 116 1 1 3 3 
3 6 111 1 1 6 3 
3 4 105 1 1 9 3 
3 3 99 1 1 12 3 
3 9 96 1 1 15 3 
3 5 85 1 1 18 3 
3 5 80 1 1 21 3 
3 3 68 1 1 24 3 
3 7 65 1 1 27 3 
3 8 48 1 1 30 3 
3 4 40 1 1 33 3 
3 2 33 1 1 36 3 
3 0 31 1 1 39 3 
3 1 28 1 1 42 3 
3 2 27 1 1 45 3 
3 3 20 1 1 48 3 
3 1 17 1 1 51 3 
3 0 16 1 1 54 3 
4 10 167 2 1 3 3 
4 5 149 2 1 6 3 
4 6 145 2 1 9 3 
4 3 138 2 1 12 3 
4 4 135 2 1 15 3 
4 5 128 2 1 18 3 
4 2 122 2 1 21 3 
4 2 120 2 1 24 3 
4 7 104 2 1 27 3 
4 9 98 2 1 30 3 
4 5 89 2 1 33 3 
4 2 57 2 1 36 3 
4 2 55 2 1 39 3 
4 4 53 2 1 42 3 
4 2 49 2 1 45 3 
4 2 26 2 1 48 3 
4 1 24 2 1 51 3 
4 1 23 2 1 54 3 
4 1 11 2 1 57 3 
4 0 10 2 1 60 3 
4 0 10 2 1 63 3 
4 2 164 1 1 3 3 
4 5 153 1 1 6 3 
4 7 148 1 1 9 3 
4 6 141 1 1 12 3 
4 12 135 1 1 15 3 
4 6 119 1 1 18 3 
4 4 113 1 1 21 3 
4 3 109 1 1 24 3 
4 5 98 1 1 27 3 
4 2 94 1 1 30 3 
4 2 92 1 1 33 3 
4 4 55 1 1 36 3 
4 3 50 1 1 39 3 
4 1 48 1 1 42 3 
4 2 47 1 1 45 3 
4 1 22 1 1 48 3 
4 1 21 1 1 51 3 
4 0 20 1 1 54 3 
4 1 7 1 1 57 3 
4 0 6 1 1 60 3 
4 0 6 1 1 63 3 
5 12 152 2 1 3 3 
5 7 135 2 1 6 3 
5 9 128 2 1 9 3 
5 8 120 2 1 12 3 
5 7 112 2 1 15 3 
5 1 77 2 1 18 3 
5 3 76 2 1 21 3 
5 2 73 2 1 24 3 
5 4 71 2 1 27 3 
5 5 45 2 1 30 3 
5 3 40 2 1 33 3 
5 2 37 2 1 36 3 
5 3 35 2 1 39 3 
5 3 32 2 1 42 3 
5 3 32 2 1 45 3 
5 1 32 2 1 48 3 
5 9 149 1 1 3 3 
5 4 131 1 1 6 3 
5 5 127 1 1 9 3 
5 8 122 1 1 12 3 
5 11 114 1 1 15 3 
5 5 76 1 1 18 3 
5 5 71 1 1 21 3 
5 5 66 1 1 24 3 
5 6 61 1 1 27 3 
5 3 35 1 1 30 3 
5 4 32 1 1 33 3 
5 1 28 1 1 36 3 
5 1 27 1 1 39 3 
5 6 26 1 1 42 3 
5 5 26 1 1 45 3 
5 0 26 1 1 48 3 
6 22 179 2 1 3 3 
6 13 151 2 1 6 3 
6 3 138 2 1 9 3 
6 5 135 2 1 12 3 
6 1 130 2 1 15 3 
6 5 104 2 1 18 3 
6 7 99 2 1 21 3 
6 6 92 2 1 24 3 
6 6 66 2 1 27 3 
6 7 60 2 1 30 3 
6 4 53 2 1 33 3 
6 0 30 2 1 36 3 
6 2 29 2 1 39 3 
6 3 27 2 1 42 3 
6 1 24 2 1 45 3 
6 0 16 2 1 48 3 
6 1 15 2 1 51 3 
6 0 14 2 1 54 3 
6 1 14 2 1 57 3 
6 0 14 2 1 60 3 
6 13 178 1 1 3 3 
6 7 160 1 1 6 3 
6 7 153 1 1 9 3 
6 10 146 1 1 12 3 
6 10 136 1 1 15 3 
6 2 97 1 1 18 3 
6 5 95 1 1 21 3 
6 3 90 1 1 24 3 
6 5 57 1 1 27 3 
6 2 52 1 1 30 3 
6 6 50 1 1 33 3 
6 3 37 1 1 36 3 
6 1 34 1 1 39 3 
6 2 33 1 1 42 3 
6 4 31 1 1 45 3 
6 0 17 1 1 48 3 
6 0 17 1 1 51 3 
6 1 17 1 1 54 3 
6 0 16 1 1 57 3 
6 0 16 1 1 60 3 
END 


ts[] bs[] 
3 1 
2 1 
2 1 
2 1 
2 1 
2 1 
END 

Answer 1

Letztendlich konnte ich das Modell nicht ordnungsgemäß in WinBUGS ausführen. Allerdings war ich in der Lage, das Modell auf STAN zu portieren und eine sehr enge Übereinstimmung zu erzielen. Siehe unten für den STAN Code:

data { 


int<lower=1> N; 
    int<lower=1> NS; 
    int<lower=1> NT; 

    cov_matrix[2] prec2; 
    matrix[2,2] R; 
    vector[2] means; 

    int<lower=0> bs[NS]; 
    int<lower=0> ts[NS]; 

    int<lower=0> s[N]; 
    int<lower=0> t[N]; 
    int<lower=0> n[N]; 
    int<lower=0> r[N]; 
    real<lower=0> dt[N]; 
    real<lower=0> time[N]; 
} 
parameters { 
    vector[2] mu[NS]; 
    vector[2] delta[NS]; 
    vector[2] dj[NT-1]; 
    cov_matrix[2] omega; 
} 
transformed parameters{ 
    real<lower=0,upper=1> p[N]; 
    real<lower=0> h[N]; 
    real nu[N]; 
    real theta[N]; 
    vector[2] md[NS]; 
    vector[2] d[NT]; 

    d[1][1] <- 0; 
    d[1][2] <- 0; 
    for(j in 2:NT){ 
    d[j] <- dj[j-1]; 
    } 
    for(k in 1:NS){ 
    md[k] <- d[ts[k]] - d[bs[k]]; 
    } 
    for(i in 1:N){ 
    if(t[i] == 1){ 
     nu[i] <- mu[s[i]][1]; 
     theta[i] <- mu[s[i]][2]; 
    }else{ 
     nu[i] <- mu[s[i]][1] + delta[s[i]][1]; 
     theta[i] <- mu[s[i]][2] + delta[s[i]][2]; 
    } 
    h[i] <- exp(nu[i] + log(time[i]) * theta[i]); 
    p[i] <- 1 - exp(- h[i] * dt[i]); 
    } 
} 
model { 
    omega ~ inv_wishart(2,R); 
    for(j in 1:(NT-1)){ 
    dj[j] ~ multi_normal(means,prec2); 
    } 
    for(k in 1:NS){ 
    delta[k] ~ multi_normal(md[k],omega); 
    mu[k] ~ multi_normal(means,prec2); 
    } 
    for(i in 1:N){ 
    r[i] ~ binomial(n[i],p[i]); 
    } 
} 
generated quantities{ 
    vector[N] log_lik; 
    for (l in 1:N) { 
    log_lik[l] <- binomial_log(r[l], n[l], p[l]); 
    } 
} 

Beachten Sie, dass aufgrund der Unterschiede zwischen STAN/BUGS, Verwirrung über Präzision/Varianz machen kann es verwirrend, was für Daten einzugeben. Für R und Prec2, ich geladen:

dataList$R <- matrix(c(0.01,0,0,0.01),nrow=2,ncol=2,byrow=TRUE) 
dataList$prec2 <- matrix(c(10^4,0,0,10^4),nrow=2,ncol=2,byrow=TRUE) 

Answer 2

Einige Vorschläge für Sie, die hoffentlich hilfreich sein werden:

1. Ich denke, Ihre Frage könnte besser beantwortet werden in https://stats.stackexchange.com/; Wenn Sie die Frage dorthin verschieben möchten, sollten Sie sie jedoch neu schreiben, anstatt einen Code-Dump zu veröffentlichen.

2. WinBUGS ist mehrere Jahre alt, und es ist 8 Jahre seit seiner letzten Aktualisierung! Ich denke, du solltest es vergessen, denn es gibt mehrere bessere Alternativen.

3. Sie können praktisch den gleichen Code versuchen in Jags oder Stan, in dem sowohl in R über rJags und RStan verwendbar sind. Stan ist besonders wichtig, weil es HCMC verwendet, das in vielen Situationen konvergiert, die WinBUGS nicht hat.

4. Ich glaube, Sie das einfache Buch lesen sollten: Doing Bayesian Data Analysis 2e von John K. Kruschke der Lage sein, selbst Bayesian Datenanalyse zu verstehen und zu tun.