WinBUGS Meta-analisi della rete Weibull

Question

Attualmente sto lavorando a una meta-analisi dei dati di sopravvivenza in diversi studi clinici.

Per fare ciò, ho il codice di un'analisi pubblicata che utilizza la stessa metodologia. Tuttavia, quando si esegue questo codice utilizzando i dati dell'analisi pubblicata, non sono in grado di replicare i risultati. In realtà, i risultati non riescono a convergere in alcun tipo di stima ragionevole.

Il codice stesso (esclusi i dati) deve essere corretto poiché proviene direttamente dagli autori. Presumo che il problema debba fare w/i valori iniziali oi parametri di come viene eseguito il campionamento, ma dopo aver giocato w/molti valori iniziali , la lunghezza di burn in, thinning, ecc ... Non ho ottenuto risultati significativi .

Gradirei tutti i suggerimenti su come eseguire questo (valori iniziali, ecc.) Per farlo funzionare correttamente. In alternativa, se ci sono problemi nel codice o se i dati sono impostati in un modo che non corrisponde al codice, sarebbe utile sapere.

Come nota a margine, sto facendo le analisi usando R2WinBUGs, anche se ho ottenuto lo stesso tipo di problemi usando WinBUG da solo.

Un po 'di fondo in più sul metodo:

Il modo in cui funziona è stimando la differenza di forma e dimensioni parametri di una distribuzione di Weibull riparametrizzato tra trattamenti attraverso molteplici studi che utilizzano effetti casuali.

La distribuzione di Weibull è riparametrizzata in modo che il registro della frequenza di rischio sia un + b * log (t) dove a è un parametro di scala eb è un parametro di forma . Da questo, è possibile calcolare la probabilità della funzione di un determinato numero di guasti su un determinato numero di pazienti su un intervallo.

Purtroppo, l'articolo è pubblico, ma se è possibile ottenere l'accesso qui è il link: http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=

un breve riassunto delle variabili inserite nel modello:

NT: Numero di trattamenti separati inclusi.

N: numero di righe nel set di dati principale. NS: Numero di studi

s: Studio che la linea di dati corrisponde a (questo è numerato 1: 6)

r: numero di pazienti che non in intervallo per questo trattamento/studio

n : numero di pazienti a rischio a inizio intervallo per questo trattamento /studiare

t: trattamento che questa linea di dati corrisponde a (numerate da 1: 3)

b: Indica quale treatme nt è la linea di base a cui vengono confrontati gli altri (impostato su 1 per ogni riga).

bs: trattamento che è il braccio di controllo dello studio

bt: Trattamento che è il braccio di ricerca di questo studio

WinBUGS codice (inclusi i dati):

#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

In definitiva, non sono riuscito a far funzionare correttamente il modello in WinBUGS. Tuttavia, sono riuscito a portare il modello su STAN e ottenere una corrispondenza molto ravvicinata. Vedi sotto per il codice STAN:

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]); 
    } 
} 

Si noti che a causa delle differenze tra STAN/BUGS, confusione sulla precisione/varianza può rendere confuso da immettere per i dati. Per R e prec2, ho caricato:

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

Alcuni suggerimenti per te, che spero ti saranno d'aiuto:

1. Penso che la tua domanda potrebbe essere meglio risposta in https://stats.stackexchange.com/; Sebbene tu voglia spostare la domanda, dovresti riscriverla invece di pubblicare un dump del codice.

2. WinBUGS ha diversi anni, ed è 8 anni dal suo ultimo aggiornamento! Penso che dovresti dimenticarlo, perché ci sono molte alternative migliori.

3. Puoi provare virtualmente lo stesso codice in Jags o Stan, in cui entrambi sono utilizzabili in R via rJags e RStan. Stan è particolarmente importante, perché usa HCMC che converge in molte situazioni che WinBUGS non possiede.

4. Penso che si dovrebbe leggere il libro facile: Doing Bayesian Data Analysis 2e da John K. Kruschke per essere in grado di capire e fare l'analisi dei dati bayesiana da soli.