Previsioni sequenziali delle lunghezze del ciclo...

Previsioni sequenziali delle lunghezze del ciclomestruale

Paola Bortot1 Bruno Scarpa2 Guido Masarotto3

1Universita di BolognaE–mail: [email protected]

2Universita di PaviaE–mail: [email protected]

3Universita di PadovaE–mail: [email protected]

17 febbraio 2005

Paola Bortot, Bruno Scarpa, Guido Masarotto Previsioni sequenziali delle lunghezze del ciclo mestruale

Il problemaIl Modello

Motivazioni

La previsione della lunghezza del ciclo mestruale e delle suefasi e un problema di grande importanza nella regolazionenaturale delle nascite e nella gestione dell’infertilita.

Lo sviluppo di un metodo semplice ed efficace per laprevisione dell’ovulazione accrescerebbe l’efficacia delcosiddetto “metodo del ritmo” per la contraccezione

L’accuratezza di alcune procedure (ad esempio testspost-coitali) e il successo di alcune misure terapeutichedipendono da un adeguata identificazione dell’ovulazione

Motivazioni

Obiettivi

L’obiettivo di questo lavoro e di sviluppare un approcciostatistico al problema

Si consideri una particolare donna e sia yt la lunghezza del suot-esimo ciclo mestruale (t = 1, 2, . . .).Il nostro obiettivo consiste nella stima dell’intera distribuzionepredittiva

Ft+1(x) = P{yt+1 ≤ x |y1, . . . , yt}

La variabilita delle lunghezze di una particolare donna non eirrilevante, per questo siamo piu interessati ad intervalli diprevisione che a previsioni puntuali (Marshall, 1965).

Un modello analogo puo essere adattato per la previsione dialtri parametri legati al ciclo, ad esempio la fase ipotermicanella temperatura basale.

Obiettivi

Ft+1(x) = P{yt+1 ≤ x |y1, . . . , yt}

Obiettivi

Ft+1(x) = P{yt+1 ≤ x |y1, . . . , yt}

I dati

I dati che usiamo sono stati raccolti dal Catholic MarriageAdvisory Council of England and Wales, un centro che forniscegratuitamente servizi di consulenza di educazione senzadiscriminazione di razza, nazionalita o affiliazione religiosa.

Le informazioni raccolte (Miolo et al. 1993) si riferiscono a uncampione di 1798 donne, per ciascuna delle quali sono adisposizione le lunghezze totali di almeno 6 cicli consecutivi,per un totale di 36641 cicli.

La piu lunga sequenza di cicli consecutivi comprende 109misure.

I dati

Il modello

Nel cercare un modello statistico per descrivere le lunghezze dei cicliosservati dobbiamo tener conto che sono presenti due fonti di variabilita,all’interno di ciascuna donna e tra le diverse donne:

le lunghezze dei cicli di una donna non sono costanti nel tempo;

donne diverse spesso hanno diverse lunghezze medie e diversa

variabilita attorno a queste medie.

L’idea base per modellare queste due fonti di variabilita puo esseredescritta come segue:

1 Per spiegare la variabilita intra–donna proponiamo un modelloparametrico, cioe

P{yt+1 ≤ x |y1, . . . , yt , θ} = Gt+1(x |y1, . . . , yt , θ)

dove Gt+1 e completamente specificato e θ e un vettore di parametriignoti.

2 Per descrivere la variabilita tra–donne facciamo variare θ tra ledonne seguendo una distribuzione di probabilita p(θ|ζ).

3 Infine, adottiamo un approccio Bayesiano specificando una

distribuzione a priori per ζ.

Il modello

Un modello per una singola donnaUn modello per tutte le donne

Evidenze empiriche

I grafici di yt contro t e le relative autocorrelazioni delle differenzeprime sottolineano alcune caratteristiche delle sequenze osservatedelle lunghezze dei cicli:

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Autocorrelation of differences

0 5 10 15 20

1 Le lunghezze dei cicli osservate sono discretecon solo un numero limitato di valori. (lelunghezze sono misurate in giorni)

2 Un trend lentamente decrescente e ingenerale osservato per sequenze che copronodiversi anni. (La lunghezza media del ciclopuo variare nel lungo periodo ed e legatoall’eta della donna, vedi Vollman, 1977)

3 Una autocorrelazione negativa a un passosulle differenze prime e presente perparecchie donne. (Le donne possono maleinterpretare i segnali della fine del ciclo,anticipandolo, quindi aggiungono al ciclosuccessivo i giorni appartenenti alprecedente)

4 Alcune osservazioni sono molto diverse dallealtre e possono essere considerate comevalori anomali (In alcuni ci puo essere unaspiegazione biologica, come perdite precocinon osservabili...)

5 Dopo aver considerato gli effetti descrittisopra, nel processo rimane del “rumore”

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Evidenze empiriche

0 20 40 60 80 100

series

0 20 40 60 80 100

0 5 10 15 20

Modellare la variabilita intra–donna

Come tener conto delle caratteristiche sopra descritte all’interno diun modello statistico?

Adottiamo una formulazione state–space

Le lunghezze osservate yt possono essere viste come unamisura discreta della vera lunghezza non osservata del ciclozt , che e una variabile continua:

yt = [zt ].

Le vere non osservate lunghezze dei cicli sono generate dal seguentemodello state-space:

zt = mt + ψt + εt ,

dove mt , ψt e εt sono definite come segue.

Assumendo che le traiettorie sono localmente costanti, adottiamoper mt il modello a passeggiata casuale

mt = mt−1 + ηt ,

dove ηt e una sequenza di variabili i.i.d. N (0, σ2η).

Dopo aver considerato un possibile trend, per ammettere unaautocorrelazione negativa a un passo, specifichiamo per ψt ilmodello

ψt = ρψt−1 + νt ,

dove νt e una sequenza di variabili i.i.d. N (0, σ2ν).

Per considerare la possibile contaminazione di una piccola frazionedi outliers, assumiamo che εt e una sequenza di variabili aleatoriei.i.d. che seguono una distribuzione mistura di componenti N (0, σ2

ε)e N (0, ασ2

ε), con pesi (1− π) e π, rispettivamente.

mt = mt−1 + ηt ,

ε)e N (0, ασ2

mt = mt−1 + ηt ,

ε)e N (0, ασ2

mt = mt−1 + ηt ,

ε)e N (0, ασ2

Il modello di variabilita intra-donna

Il modello completo e

yt = [zt ]

zt = mt + ψt + εt

mt = mt−1 + ηt

ψt = ρψt−1 + νt

εt ∼ (1− π)N (0, σ2ε) + πN (0, ασ2

ε) i.i.d.,

ηt ∼ N (0, σ2η) i.i.d., νt ∼ N (0, σ2

ν) i.i.d.

Specificazione delle distribuzioni a priori

Il vettore dei parametri del modello e θ = (σ2ε , π, α, σ

2η, σ

2ν , ρ).

Assumiamo a priori che le componenti di θ sono mutuamenteindipendentiSpecifichiamo a priori piatte per σ−2

ε , σ−2η e ρ.

Per garantire l’identificabilita scegliamo le distribuzioni a prioriper π e α−1 in modo da assegnare una maggior massa su(0, 0.5) e (0, 1), rispettivamente.Useremo quindi le seguenti distribuzioni

π ∼ Beta(1, 10),

α−1 ∼ Gamma(5, 10),

σ−2ε ∼ Gamma(0.13, 0.13),

σ−2η ∼ Gamma(0.13, 0.13),

σ−2ν ∼ Gamma(0.13, 0.13),

ρ ∼ N (0, 100).

2η, σ

2ν , ρ).

ε , σ−2η e ρ.

Per garantire l’identificabilita scegliamo le distribuzioni a prioriper π e α−1 in modo da assegnare una maggior massa su(0, 0.5) e (0, 1), rispettivamente.

Useremo quindi le seguenti distribuzioni

π ∼ Beta(1, 10),

α−1 ∼ Gamma(5, 10),

σ−2ε ∼ Gamma(0.13, 0.13),

σ−2η ∼ Gamma(0.13, 0.13),

σ−2ν ∼ Gamma(0.13, 0.13),

ρ ∼ N (0, 100).

2η, σ

2ν , ρ).

ε , σ−2η e ρ.

Per garantire l’identificabilita scegliamo le distribuzioni a prioriper π e α−1 in modo da assegnare una maggior massa su(0, 0.5) e (0, 1), rispettivamente.Useremo quindi le seguenti distribuzioni

π ∼ Beta(1, 10),

α−1 ∼ Gamma(5, 10),

σ−2ε ∼ Gamma(0.13, 0.13),

σ−2η ∼ Gamma(0.13, 0.13),

σ−2ν ∼ Gamma(0.13, 0.13),

ρ ∼ N (0, 100).

Inferenza

L’inferenza e fatta usando procedure MCMC.

Abbiamo aumentato lo spazio parametrico per includere lelunghezze dei cicli non osservate zt e delle variabili latenti diBernoulli per identificare i valori anomali. Condizionatamentea questo vettore aumentato di parametri ignoti, il modelloproposto si riduce a un modello state–space lineare Gaussiano.

L’algoritmo Gibbs Sample multi–mosse di Shephard (1994) eCarter and Kohn (1996) puo allora essere applicato per fareinferenza.

Dato il relativamente grande numero di varianze ignote persemplificare la stima poniamo il vincolo

σν/σε = 0.5.

Inferenza

σν/σε = 0.5.

Inferenza

σν/σε = 0.5.

Inferenza

σν/σε = 0.5.

Risultati: Donna I

Mostriamo i risultati ottenuti nell’adattare il modello alla sequenza piulunga di lunghezze di cicli consecutivi da una singola donna. (109osservazioni).

Il processo stimato mt e tracciato sul grafico delle lunghezze osservate e ilprocesso stimato ψt :

0 20 40 60 80 100

Indicatori di sintesi delle stime dei parametri del modello

ση σε ρ π αmediana 1.01 1.50 -0.27 0.07 3.2495% intervallo di credibilita (0.72, 1.36) (0.98, 1.92) (-0.75, 0.71) (0.01, 0.24) (1.2, 11.0)

Risultati: Donna I

Mostriamo i risultati ottenuti nell’adattare il modello alla sequenza piulunga di lunghezze di cicli consecutivi da una singola donna. (109osservazioni).

Il processo stimato mt e tracciato sul grafico delle lunghezze osservate e ilprocesso stimato ψt :

0 20 40 60 80 100

Indicatori di sintesi delle stime dei parametri del modello

ση σε ρ π αmediana 1.01 1.50 -0.27 0.07 3.2495% intervallo di credibilita (0.72, 1.36) (0.98, 1.92) (-0.75, 0.71) (0.01, 0.24) (1.2, 11.0)

Risultati: Donna I

Buone prestazioni dell’algoritmoMCMC. Come esempio mostriamol’output (5000 iterazioni) dellacatena MCMC per ση:

0 1000 2000 3000 4000 5000

iteration

Le assunzioni del modello sembranovalide. A fianco mostriamo il qqplotNormale per i residui del modello:

−2 −1 0 1 2

Normal Q−Q Plot

Theoretical Quantiles

Risultati: Donna I

0 1000 2000 3000 4000 5000

iteration

−2 −1 0 1 2

Normal Q−Q Plot

Risultati: Donna I

0 1000 2000 3000 4000 5000

iteration

−2 −1 0 1 2

Normal Q−Q Plot

Risultati: Donna I

0 1000 2000 3000 4000 5000

iteration

−2 −1 0 1 2

Normal Q−Q Plot

Risultati: Donna II

Mostriamo anche i risultati dell’inferenza per una sequenza di lunghezzedi cicli consecutivi da una seconda donna. (102 osservazioni).

Il processo stimato mt disegnato sul grafico delle lunghezze osservate e ilprocesso ψt stimato:

0 20 40 60 80 100

ση σε ρ π αmediana 0.98 1.40 0.14 0.10 13.495% intervallo di credibilita (0.58, 1.42) (0.98, 1.88) (-0.66, 0.90) (0.03, 0.21) (5.5, 30.1)

Le differenze principali tra le due donne si osservano nella stima di ρ e α.

Tuttavia nell’insieme di tutte le donne si osservano differenze anche su σε.

Risultati: Donna II

0 20 40 60 80 100

Risultati: Donna II

0 20 40 60 80 100

La maggior parte di donne hanno sequenze relativamente brevi:non sonopossibili stime individuali

Una soluzione possibile e l’uso di un modello gerarchico che permette untrasferimento di informazioni tra le donne e la stima di parametri dipopolazione per la previsione su donne non incluse nell’insieme di dati.La formulazione del modello gerarchico segue direttamente dal modellointra–donna, permettendo a θ di variare tra le donne secondo unadistribuzione di probabilita p(θ|ζ).Pare sensato considerare variabili tra le donne solo alcuni parametri: σε

(la varianza degli errori) e π (la proporzione di cicli anomali nellasequenza). Consideriamo, quindi, effetti casuali solo σε e π considerandogli altri parametri come fissi tra le donne.Assumiamo

σ−2ε ∼ Gamma(a1, b1) e π ∼ Beta(a2, a2b2)

con specificazione a priori

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

L’inferenza viene effettuata attraverso un’estensione del Gibbs sampling

con aggiornamento degli stati in blocco (Shephard, 1994; Carter and

Kohn, 1996)

La maggior parte di donne hanno sequenze relativamente brevi:non sonopossibili stime individualiUna soluzione possibile e l’uso di un modello gerarchico che permette untrasferimento di informazioni tra le donne e la stima di parametri dipopolazione per la previsione su donne non incluse nell’insieme di dati.

La formulazione del modello gerarchico segue direttamente dal modellointra–donna, permettendo a θ di variare tra le donne secondo unadistribuzione di probabilita p(θ|ζ).Pare sensato considerare variabili tra le donne solo alcuni parametri: σε

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

Kohn, 1996)

La maggior parte di donne hanno sequenze relativamente brevi:non sonopossibili stime individualiUna soluzione possibile e l’uso di un modello gerarchico che permette untrasferimento di informazioni tra le donne e la stima di parametri dipopolazione per la previsione su donne non incluse nell’insieme di dati.La formulazione del modello gerarchico segue direttamente dal modellointra–donna, permettendo a θ di variare tra le donne secondo unadistribuzione di probabilita p(θ|ζ).

Pare sensato considerare variabili tra le donne solo alcuni parametri: σε

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

Kohn, 1996)

La maggior parte di donne hanno sequenze relativamente brevi:non sonopossibili stime individualiUna soluzione possibile e l’uso di un modello gerarchico che permette untrasferimento di informazioni tra le donne e la stima di parametri dipopolazione per la previsione su donne non incluse nell’insieme di dati.La formulazione del modello gerarchico segue direttamente dal modellointra–donna, permettendo a θ di variare tra le donne secondo unadistribuzione di probabilita p(θ|ζ).Pare sensato considerare variabili tra le donne solo alcuni parametri: σε

(la varianza degli errori) e π (la proporzione di cicli anomali nellasequenza). Consideriamo, quindi, effetti casuali solo σε e π considerandogli altri parametri come fissi tra le donne.

Assumiamo

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

Kohn, 1996)

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

Kohn, 1996)

a1 ∼ Gamma(0.13, 0.13), b1 ∼ Gamma(0.13, 0.13)

a2 ∼ Gamma(0.13, 0.13), b2 ∼ Gamma(1, 10)

Kohn, 1996)Paola Bortot, Bruno Scarpa, Guido Masarotto Previsioni sequenziali delle lunghezze del ciclo mestruale

Alcuni primi risultati: 50 donne

0 1000 2000 3000 4000 5000

a0 1000 3000 5000

0 1000 3000 5000

Indexra

0 1000 3000 5000

ση ρ σν α a1 b1 a2 b2media 2.35 –0.01 0.0013 6.96 1.15 0.26 0.84 18.11

standard deviation 0.14 0.1085 4.47·10−5 1.95 0.20 0.08 0.22 1.45

Buone prestazioni dell’MCMC. Risultati promettenti.Elevata complessita computazionale e uso massiccio dimemoria

Alcuni primi risultati: 50 donne

0 1000 2000 3000 4000 5000

a0 1000 3000 5000

0 1000 3000 5000

Indexra

0 1000 3000 5000

ση ρ σν α a1 b1 a2 b2media 2.35 –0.01 0.0013 6.96 1.15 0.26 0.84 18.11

standard deviation 0.14 0.1085 4.47·10−5 1.95 0.20 0.08 0.22 1.45

Buone prestazioni dell’MCMC. Risultati promettenti.Elevata complessita computazionale e uso massiccio dimemoriaPaola Bortot, Bruno Scarpa, Guido Masarotto Previsioni sequenziali delle lunghezze del ciclo mestruale

Lavori futuriStima dei modelli e applicazioniEstensioni del modello

Stima dei modelli e applicazioni

Stima del modello gerarchico per le lunghezze del ciclo usandotutte le osservazioni disponibili

Stima del modello gerarchico per le lunghezze delle fasi pre– epost– ovulatoria

Applicazione alla previsione della probabilita di concepimentosotto un fissato pattern di rapporti.

Estensioni del modello e dell’approccio inferenziale

Inclusione di covariate (eta della donna, temperatura, tipo dimuco, ...).

Stima dei parametri usando procedure Monte Carlosequenziali, per sviluppare algoritmi dinamici per la previsione.

...non senza fadiga si giunge al...

Previsioni sequenziali delle lunghezze del ciclo...

Documents

Transcript of Previsioni sequenziali delle lunghezze del ciclo...