Probability Density Function (PDF) di una funzione data
Salve a tutti! Il mio primo post qui.
Inizio con un problema per me molto complesso che cerco di esplicitare in breve. Ho una curva "dose/risposta" che risponde ad una legge data dalla funzione:
$ f(x):1/(1+(theta/x)^gamma) $
Ho una serie di valori sperimentali funzione di $ x $ che dovrebbero "fittare" tale funzione. In base a questi valori sperimentali devo calcolare il valore di $ theta $ e $gamma $ che sono costanti della funzione.
Dopo essermi informato leggendo di tutto e di più ho trovato che la maniera per risolvere questo problema risiede nella determinazione della cosiddetta "stima di massima verosimiglianza" (Maximum Likelihood Extimation - MLE) che consente di ricavare i valori delle costanti più probabili a determinare la mia serie di dati.
Il primo passo del problema dovrebbe essere, se ho capito bene, ricavare la "funzione di verosimiglianza" (Likelihood function: $ L $) che a sua volta dipende da una una "probability density function" (PDF) per tale distribuzione di valori.
Quindi:
1) Come posso ricavare la PDF a partire dalla mia $ f(x) $ ?
Una volta determinata $ PDF = p(x; theta, gamma) $ posso ricavare la $ L $ in questo modo:
$ L=prod_(i = 1)^(n) p(x;theta,gamma) $
Mi piacerebbe sapere se fino ad ora ho scritto castronerie o sono andato bene. Dopodichè mi servirebbe ottenere la $ PDF $ per iniziare i calcoli...
Grazie
Inizio con un problema per me molto complesso che cerco di esplicitare in breve. Ho una curva "dose/risposta" che risponde ad una legge data dalla funzione:
$ f(x):1/(1+(theta/x)^gamma) $
Ho una serie di valori sperimentali funzione di $ x $ che dovrebbero "fittare" tale funzione. In base a questi valori sperimentali devo calcolare il valore di $ theta $ e $gamma $ che sono costanti della funzione.
Dopo essermi informato leggendo di tutto e di più ho trovato che la maniera per risolvere questo problema risiede nella determinazione della cosiddetta "stima di massima verosimiglianza" (Maximum Likelihood Extimation - MLE) che consente di ricavare i valori delle costanti più probabili a determinare la mia serie di dati.
Il primo passo del problema dovrebbe essere, se ho capito bene, ricavare la "funzione di verosimiglianza" (Likelihood function: $ L $) che a sua volta dipende da una una "probability density function" (PDF) per tale distribuzione di valori.
Quindi:
1) Come posso ricavare la PDF a partire dalla mia $ f(x) $ ?
Una volta determinata $ PDF = p(x; theta, gamma) $ posso ricavare la $ L $ in questo modo:
$ L=prod_(i = 1)^(n) p(x;theta,gamma) $
Mi piacerebbe sapere se fino ad ora ho scritto castronerie o sono andato bene. Dopodichè mi servirebbe ottenere la $ PDF $ per iniziare i calcoli...
Grazie
Risposte
Non penso devi ricominciare da capo, magari occorre solo aggiustare un po' il tiro 
Questo mi sembra una cosa diversa. Sono d'accordo che si vuole anche rappresentare i dati sul grafico.
La cosa importante, per quanto ne ho capito (ma potrei anche sbagliarmi!) è che il fitting, la massimizzazione della log-verosimiglianza, sia fatta con una delle due formule viste in precedenza:
1) Se hai dati binari, con $y_i$ uguale a 0 o 1. E questa potrebbe essere una possibile visualizzazione:
(ogni cerchietto è un paziente, ho inventato dei dati: cerchietto con $y_i=1$, il paziente presenta l'effetto se trattato con la relativa dose $x_i$)
2) Se invece (ma non mi sembra il tuo caso) hai sottoposto un gruppo di $n_i$ pazienti alla stessa dose $x_i$ e hai osservato $y_i$ successi in questo gruppo e quindi $n_i-y_i$ insuccessi, allora usi quall'altra formula della log-verosimiglianza (entrambe riportate verso la fine di questo post).
Il grafico che ne uscirebbe è più del tipo che vorresti.
C'è però una cosa che forse sto/stiamo trascurando.
Leggendo questo articolo http://www.redjournal.org/article/S0360 ... 8/fulltext
mi sembra di capire che hanno 1530 pazienti e ne fanno 4 gruppi, sebbene ogni paziente è trattato con una diverse dose.
Di ogni gruppo considerano la dose mediana e poi hanno fatto una regressione logistica su questi 4 gruppi, tipo quella al punto 2), e ne esce un grafico simil-Bentzen.
Inoltre la faccenda pare complicarsi, perchè per stimare i risultati all'interno di ciascun gruppo (outcome rates) hanno usato un modello "Multivariate Cox proportional hazards regression".
Altra cosa: non ho ben capito di che tipo di dati disponi. Magari può essere utile.
Oltre la dose EUD, cosa misuri/associ come risposta al paziente ?
"nickwing":
ed io "brutalmente" volevo vedere i miei dati rappresentati sul medesimo grafico della curva interpolante.
Questo mi sembra una cosa diversa. Sono d'accordo che si vuole anche rappresentare i dati sul grafico.
La cosa importante, per quanto ne ho capito (ma potrei anche sbagliarmi!) è che il fitting, la massimizzazione della log-verosimiglianza, sia fatta con una delle due formule viste in precedenza:
1) Se hai dati binari, con $y_i$ uguale a 0 o 1. E questa potrebbe essere una possibile visualizzazione:
(ogni cerchietto è un paziente, ho inventato dei dati: cerchietto con $y_i=1$, il paziente presenta l'effetto se trattato con la relativa dose $x_i$)
2) Se invece (ma non mi sembra il tuo caso) hai sottoposto un gruppo di $n_i$ pazienti alla stessa dose $x_i$ e hai osservato $y_i$ successi in questo gruppo e quindi $n_i-y_i$ insuccessi, allora usi quall'altra formula della log-verosimiglianza (entrambe riportate verso la fine di questo post).
Il grafico che ne uscirebbe è più del tipo che vorresti.
C'è però una cosa che forse sto/stiamo trascurando.
Leggendo questo articolo http://www.redjournal.org/article/S0360 ... 8/fulltext
mi sembra di capire che hanno 1530 pazienti e ne fanno 4 gruppi, sebbene ogni paziente è trattato con una diverse dose.
Di ogni gruppo considerano la dose mediana e poi hanno fatto una regressione logistica su questi 4 gruppi, tipo quella al punto 2), e ne esce un grafico simil-Bentzen.
Inoltre la faccenda pare complicarsi, perchè per stimare i risultati all'interno di ciascun gruppo (outcome rates) hanno usato un modello "Multivariate Cox proportional hazards regression".
Altra cosa: non ho ben capito di che tipo di dati disponi. Magari può essere utile.
Oltre la dose EUD, cosa misuri/associ come risposta al paziente ?
Tra l'altro, ora ci faccio caso, vedendo l'articolo di Zhu che hai allegato (in cui usa il modello probit), pag. 451(colonna sinistra in alto), nella log-verosimiglianza compare la variabile misurata $R_i$: "Let Ri = 1 if the patient i experienced AE and Ri = 0 otherwise".
Forse (non capisco) nel grafico che presenta in basso a pag. 452, le osservazioni, rappresentate con dei quadratini, le ottiene aggregando i pazienti in gruppi e valutando in ogni gruppo la % di successi, in modo da non avere un grafico con i cerchietti a y=0 o y=1 ?
Forse (non capisco) nel grafico che presenta in basso a pag. 452, le osservazioni, rappresentate con dei quadratini, le ottiene aggregando i pazienti in gruppi e valutando in ogni gruppo la % di successi, in modo da non avere un grafico con i cerchietti a y=0 o y=1 ?
Allora, con ordine:
il dato clinico che ho è la comparsa di un effetto collaterale valutata in visita. L'esito può essere solo 0 o 1 (ossia assente o presente). Ogni paziente viene sottoposto a livelli di dose che possono quindi, ad un certo punto, suscitare tale effetto. La variabile EUD è una variabile continua, certo potrei fare dei raggruppamenti ma questo imporrebbe ulteriori approssimazioni. Durante il trattamento però gli effetti compaiono in maniera variabile, diversa da paziente a paziente. Se per un paziente è stato preparato un piano di cura in cui, a fine trattamento, la sua EUD sull'organo che sto analizzando è, mettiamo, 60 Gy, e questi manifesta l'effetto a due terzi del trattamento, la dose a cui quel paziente verà valutato positivo (quindi "1") sarà di 40 Gy. Se un paziente non presenta mai tale effetto alla EUD totale del piano valuterò il paziente come "0". Quindi nella tabella che ti inserisco i pazienti valutati "1" hanno presentato l'effetto alla EUD indicata per la positività, i pazienti valutati come "0" a fine trattamento sono arrivati senza manifestare l'effetto. Ma se considero ogni paziente come gruppo a sè, questi indica che tutti quelli che hanno avuto l'effetto per dosi uguali o inferiori rappresentano una parte del tutto, e quindi esprimono la probabilità cumulativa di avere l'effetto per quella data dose. Ogni paziente che presenta l'effetto ad una data dose quindi "cumula" i pazienti precedenti, per questo pensavo di posizionare i miei dati direttamente sulla curva come pallini in cui all'ordinata corrisponde il numero di pazienti che hanno presentato l'effetto a dose inferiore o uguale a quella del paziente osservato. In numeri ecco i dati su 14 pazienti:
Nell'articolo di Zhu hanno definito le "osservazioni", cioè i quadratini nel grafico della curva dose/riposta proprio come "raggruppamenti" di pazienti. In realtà non dicono molto altro ma ritengo che se 100% è il numero di pazienti questi autori abbiano semplicemente raprpesentato gruppi con $EUD<= $ dei livelli opportunamente scelti in base ai piani di cura a loro disposizione per il loro pazienti. L'idea di lavorare in questo modo, anche se non proprio ortodosso, mi è derivata dal fatto che le curve sono proprio curve di probabilità cumulativa, cioè ogni valore della $x$ e cioè ogni livello di $EUD$, include in se tutti i livelli precedenti perchè il paziente per raggiungere tale EUD deve essere passato attraverso tutti i livelli precedenti
Se mi copi in tabella i dati che hai usato per la tua curva provo a vedere che succede nel mio foglio elettronico per il calcolo delle costanti
il dato clinico che ho è la comparsa di un effetto collaterale valutata in visita. L'esito può essere solo 0 o 1 (ossia assente o presente). Ogni paziente viene sottoposto a livelli di dose che possono quindi, ad un certo punto, suscitare tale effetto. La variabile EUD è una variabile continua, certo potrei fare dei raggruppamenti ma questo imporrebbe ulteriori approssimazioni. Durante il trattamento però gli effetti compaiono in maniera variabile, diversa da paziente a paziente. Se per un paziente è stato preparato un piano di cura in cui, a fine trattamento, la sua EUD sull'organo che sto analizzando è, mettiamo, 60 Gy, e questi manifesta l'effetto a due terzi del trattamento, la dose a cui quel paziente verà valutato positivo (quindi "1") sarà di 40 Gy. Se un paziente non presenta mai tale effetto alla EUD totale del piano valuterò il paziente come "0". Quindi nella tabella che ti inserisco i pazienti valutati "1" hanno presentato l'effetto alla EUD indicata per la positività, i pazienti valutati come "0" a fine trattamento sono arrivati senza manifestare l'effetto. Ma se considero ogni paziente come gruppo a sè, questi indica che tutti quelli che hanno avuto l'effetto per dosi uguali o inferiori rappresentano una parte del tutto, e quindi esprimono la probabilità cumulativa di avere l'effetto per quella data dose. Ogni paziente che presenta l'effetto ad una data dose quindi "cumula" i pazienti precedenti, per questo pensavo di posizionare i miei dati direttamente sulla curva come pallini in cui all'ordinata corrisponde il numero di pazienti che hanno presentato l'effetto a dose inferiore o uguale a quella del paziente osservato. In numeri ecco i dati su 14 pazienti:
EUD [Gy] Effetto [assente = 0, presente = 1] 10 0 12 1 15 1 16 0 18 1 19 0 20 1 22 0 26 1 35 1 40 1 55 1 67 1 72 1 NTCP[10 Gy] = 0% NTCP[12 Gy] = 1/14 = 7.1% NTCP[15 Gy] = 2/14 = 14.3% NTCP[18 Gy] = 3/14 = 21.4% NTCP[20 Gy] = 4/14 = 28.6% NTCP[26 Gy] = 5/14 = 35.7% NTCP[35 Gy] = 6/14 = 42.9% NTCP[40 Gy] = 7/14 = 50% NTCP[55 Gy] = 8/14 = 57.1% NTCP[67 Gy] = 9/14 = 64.3% NTCP[72 Gy] = 10/14 = 71.4%
Nell'articolo di Zhu hanno definito le "osservazioni", cioè i quadratini nel grafico della curva dose/riposta proprio come "raggruppamenti" di pazienti. In realtà non dicono molto altro ma ritengo che se 100% è il numero di pazienti questi autori abbiano semplicemente raprpesentato gruppi con $EUD<= $ dei livelli opportunamente scelti in base ai piani di cura a loro disposizione per il loro pazienti. L'idea di lavorare in questo modo, anche se non proprio ortodosso, mi è derivata dal fatto che le curve sono proprio curve di probabilità cumulativa, cioè ogni valore della $x$ e cioè ogni livello di $EUD$, include in se tutti i livelli precedenti perchè il paziente per raggiungere tale EUD deve essere passato attraverso tutti i livelli precedenti
Se mi copi in tabella i dati che hai usato per la tua curva provo a vedere che succede nel mio foglio elettronico per il calcolo delle costanti
Quanto al discorso "Redjournal" l'analisi multivariata credo sia stata resa necessaria per stratificare i loro pazienti in base all'effetto finale che è il controllo di malattia (ti ripeto che queste curve si applicano indifferentemente al controllo tumorale o agli effetti collaterali), però la dose somministrata al tumore in genere è molto omogenea, e per questo i pazienti possono raggiungere solo livelli predefiniti (con limitata variabilità) decisi a priori per controllare la malattia. L'analisi multivariata è stata fatta per verificare quanto la dose impattasse rispetto al controllo di malattia indipendentemente da altri fattori prognostici (PSA, stadio o quant'altro).
Ecco i dati (inventati), controlla dovrebbero essere 200 "osservazioni".
Scusa nick, stasera leggo con calma il tuo messaggio, ora sono "in a hurry"
Ciao
x y 20 0 20 0 20 0 20 0 21 0 21 0 21 0 21 0 21 0 21 0 21 1 22 0 22 0 22 0 22 0 23 0 23 0 24 0 24 0 24 0 24 0 24 0 25 0 25 0 25 0 26 0 26 0 26 0 26 0 26 0 26 0 26 0 26 0 26 0 26 0 26 0 27 0 27 0 27 0 27 0 27 0 27 0 27 1 27 1 27 1 28 0 28 0 28 0 28 1 29 0 29 0 29 0 29 0 29 0 29 0 29 1 29 1 30 0 30 0 30 0 30 0 30 0 31 0 31 0 31 0 31 0 31 0 31 0 32 0 32 0 33 0 33 0 33 0 33 0 33 0 33 0 34 0 34 0 34 0 34 0 34 1 34 1 35 0 35 0 35 0 35 0 35 1 36 0 36 0 36 0 36 0 37 0 37 0 37 0 37 0 37 0 37 0 37 0 37 0 37 0 37 0 38 0 38 0 38 0 39 0 39 0 39 1 39 1 39 1 40 0 40 1 41 0 42 1 42 1 42 1 42 1 43 1 43 1 44 0 44 1 44 1 45 1 45 1 45 1 45 1 45 1 45 1 46 1 46 1 46 1 47 1 47 1 47 1 47 1 48 1 48 1 48 1 48 1 48 1 48 1 49 1 49 1 49 1 50 0 50 1 50 1 50 1 51 1 51 1 52 1 52 1 52 1 52 1 52 1 52 1 53 1 53 1 53 1 53 1 54 1 55 1 55 1 55 1 55 1 55 1 57 1 57 1 58 1 58 1 59 1 59 1 59 1 59 1 59 1 59 1 59 1 59 1 60 1 60 1 61 1 61 1 61 1 61 1 61 1 63 1 63 1 63 1 63 1 64 1 64 1 64 1 64 1 64 1 65 1 66 1 66 1 68 1 68 1 70 1 72 1
Scusa nick, stasera leggo con calma il tuo messaggio, ora sono "in a hurry"
Ciao
Allora ho "fatto i conti" ed ho trovato i tuoi stessi valori. Utilizzando però la funzione nel modo in cui dici tu, ossia associando il valore 1 o 0 alla $y_i$, facendo come dicevo io i risultati semplicemente erano diversi, e quindi i dati visibili come plottati sul grafico non coincidevano per niente. Il problema è che ora non riesco a capire ancora come renderli visibili. Quel che vorrei è che l'utente (o "utonto") medico che utilizza il programma semplicemente veda su schermo i dati così come li ha raccolti a confronto con la curva ottenuta. Come facevano i cinesi nell'articolo di Zhu.
Ho capito!
Nel grafico devo fare un plot per gruppi di pazienti, come è stato fatto qui:
divido il mio set di pazienti in gruppi per intervalli di dose uguali, io ho diviso il tuo set in intervalli per 13 Gy (giusto perchè la differenza fra dose massima e minima era 52 e quindi potevo fare 4 gruppi con ampiezza uguale), dopodichè devo fare un plot della probabilità misurata all'interno del subset di pazienti sulla dose media nel subset. Il risultato lo vedi in questa figura
Che ne dici?
Nel grafico devo fare un plot per gruppi di pazienti, come è stato fatto qui:
divido il mio set di pazienti in gruppi per intervalli di dose uguali, io ho diviso il tuo set in intervalli per 13 Gy (giusto perchè la differenza fra dose massima e minima era 52 e quindi potevo fare 4 gruppi con ampiezza uguale), dopodichè devo fare un plot della probabilità misurata all'interno del subset di pazienti sulla dose media nel subset. Il risultato lo vedi in questa figura
Che ne dici?
P.S. Ecco perché nell' articolo di Zhu raggruppano i pazienti
"nickwing":
Ho capito! [...] Che ne dici?
OK, così sono abbastanza d'accordo.
Fermo restando che la massimizzazione della log-verosimiglianza va fatta con $y_i=0$ o $y_i=1$ e dati scorporati, cioè non raggruppati.
Il raggruppamento lo fai solo per rappresentare i dati sul grafico.
Dicevo abbastanza (d'accordo) perchè mi sembra più coerente con i dati disponibili (che sono binari: 0 o 1) la rappresentazione che avevo proposto qui sopra.
E' vero che graficamente sembra "brutta": la curva appare "lontana" dai dati, ma l'artificio del raggruppamento è solo per un fatto estetico di mostrare che i punti sono vicini alla curva...
Per quanto riguarda l'interpretazione della curva come una probabilità cumulata ho qualche perplessità.
Pensavo: la curva a sigmoide non è l'unica curva dose/risposta come leggevo qui.
Immagino una curva dose/risposta a forma di campana, tipo la curva normale.
Potrebbero cioè esistere fenomeni per cui la risposta massima si ha per certi valori (centrali) della dose e si ha una scarsa risposta per valori inferiori ma anche superiori. In tal caso non ci vedo adatto il discorso della probabilità cumulata.
Questa è la perplessità.. ci dovrei riflettere un po'.
A parte queste cose minori, comunque direi che il resto funziona, giusto ?
PS: con pochi pazienti a disposizione il raggruppamento per il grafico potrebbe essere comunque esteticamente non molto gradevole.
Invece che la media del gruppo puoi anche considerare la mediana, come nell'articolo del Redjournal.
"cenzo":
[quote="nickwing"]Ho capito! [...] Che ne dici?
OK, così sono abbastanza d'accordo.
Fermo restando che la massimizzazione della log-verosimiglianza va fatta con $y_i=0$ o $y_i=1$ e dati scorporati, cioè non raggruppati.
Il raggruppamento lo fai solo per rappresentare i dati sul grafico.
Dicevo abbastanza (d'accordo) perchè mi sembra più coerente con i dati disponibili (che sono binari: 0 o 1) la rappresentazione che avevo proposto qui sopra.
E' vero che graficamente sembra "brutta": la curva appare "lontana" dai dati, ma l'artificio del raggruppamento è solo per un fatto estetico di mostrare che i punti sono vicini alla curva...
[/quote]Sono interamente d'accordo con te. Devo dire l'artificio del raggruppamento è qualcosa che deriva dall'abitudine a leggere pubblicazioni che rappresentano i dati in quella maniera.
"cenzo":
Per quanto riguarda l'interpretazione della curva come una probabilità cumulata ho qualche perplessità.
Pensavo: la curva a sigmoide non è l'unica curva dose/risposta come leggevo qui.
Immagino una curva dose/risposta a forma di campana, tipo la curva normale.
Potrebbero cioè esistere fenomeni per cui la risposta massima si ha per certi valori (centrali) della dose e si ha una scarsa risposta per valori inferiori ma anche superiori. In tal caso non ci vedo adatto il discorso della probabilità cumulata.
Questa è la perplessità.. ci dovrei riflettere un po'.
In effetti "capendo" il discorso nel suo insieme ho realizzato che avevo male interpretato il ruolo dei punti "sperimentali" nei grafici visti nelle pubblicazioni. E dunque anche il discorso della probabilita cumulata va a farsi benedire. Era una mia interpretazione (scorretta).
"cenzo":
Allora con questo hai finito ?
Per ora si. Ma hai un messaggio privato
Ora si aprirebbe un nuovo problema: come valutare l'affidabilità delle costanti calcolate con la MLE?
Ossia quale test posso fare per valutare la differenza fra casi osservati e modello proposto? Esiste un qualche modo per calcolare "una sorta di p-value" che stimi la probabilità che modello e dati siano differenti (o vicecversa che siano uguali)? Chiedo di questo calcolo di probabilità (p-value) perchè conosco l'utenza a cui si indirizza questo progetto (i medici), e per essi il 90% delle conoscenze statistiche si limitano a verificare "se $p<= 0.05$"...
In letteratura ho trovato esempi nei quali si parla di valutazione mediante il $chi^2$ o l'F-test. Mi servirebbe una dritta per implementarlo numericamente. Ho trovato poi anche un altro modo per valutare la "bontà" del fit, ossia l'Akaike information criterion (AIC). Quest'ultimo mi sembra più semplice da implementare ma mi servirebbe solo sapere quale peso dare ai risultati del calcolo che si basa sul valore di massima verosimiglianza e sul numero di parametri stimati secondo questa relazione:
$AIC=2k-2ln(L)$, dove $k$ è il numero di costanti stimate e $L$ è l'ormai nota Maximum LIkelihood Estimation.
Grazie...
Ossia quale test posso fare per valutare la differenza fra casi osservati e modello proposto? Esiste un qualche modo per calcolare "una sorta di p-value" che stimi la probabilità che modello e dati siano differenti (o vicecversa che siano uguali)? Chiedo di questo calcolo di probabilità (p-value) perchè conosco l'utenza a cui si indirizza questo progetto (i medici), e per essi il 90% delle conoscenze statistiche si limitano a verificare "se $p<= 0.05$"...
In letteratura ho trovato esempi nei quali si parla di valutazione mediante il $chi^2$ o l'F-test. Mi servirebbe una dritta per implementarlo numericamente. Ho trovato poi anche un altro modo per valutare la "bontà" del fit, ossia l'Akaike information criterion (AIC). Quest'ultimo mi sembra più semplice da implementare ma mi servirebbe solo sapere quale peso dare ai risultati del calcolo che si basa sul valore di massima verosimiglianza e sul numero di parametri stimati secondo questa relazione:
$AIC=2k-2ln(L)$, dove $k$ è il numero di costanti stimate e $L$ è l'ormai nota Maximum LIkelihood Estimation.
Grazie...
Ciao nick,
con riferimento ai 200 dati inventati più su, dovrebbe risultare una log-verosimiglianza massima pari a $-61.003$, quindi $-2logL=122.01$.
Risulta quindi un $"AIC"=122.01+2*2=126.01$ (nel nostro modello abbiamo $k=2$ parametri da stimare).
L'AIC è utile per confrontare diversi modelli tra loro: risulta preferibile quello con un AIC più basso.
Anche il test F è usato nell'articolo per confrontare i vari modelli con quello "base" di Lyman (il probit).
Ma nel tuo caso hai un unico modello, non vuoi fare un confronto con altri modelli.
Per valutare la bontà di adattamento potresti usare il test $\chi^2$.
La devianza residua non spiegata dal tuo modello (logit binario) è uguale a $-2logL=122.01$.
Potresti confrontare questo valore con una $\chi^2$ con 198 gradi di libertà (200-2, dove 2 sono i parametri del modello).
Risulta \( P(\chi^2_{198} > 122.01)=0.9999951 \)
Valori inferiori alla solita soglia $0.05$ indicano uno scarso adattamento. Nel nostro caso si avrebbe un buon adattamento.
C'è un però...
Nel caso di modello binario (è il nostro caso ma non quello dell'articolo Cancers), non è corretto usare questo test.
Siccome la dimensione dei gruppi e’ 1 (non abbiamo gruppi..), la devianza non è distribuita come una chi-quadrato, quindi non e’ possibile valutare la bonta’ di adattamento con quell'indice. Non vale più la convergenza alla chi-quadrato.
Ho letto che per ovvviare a ciò esistono varie alternative, tra cui il test di Hosmer–Lemeshow, che prevede comunque un raggruppamento dei dati in 10 classi..
Un atro modo potrebbe essere di utilizzare dei "Pseudo R-Squared", sulla falsariga di quello che si usa nel modello dei minimi quadrati.
Ad esempio, con gli stessi dati, il modello nullo ($\beta=0$, solo l'intercetta) ha una devianza di $-2logL=-2*(-138.63)=277.26$ (*)
Si avrebbe quindi (McFadden's): \( \displaystyle R^2=1-\frac{122.01}{277.26}=0.560 \) (più alto è meglio).
(*) Per il modello nullo puoi calcolare:
\( \displaystyle \log L=n \cdot [\bar y \cdot \log (\bar y) + (1-\bar y) \cdot \log (1- \bar y)] = 200 \cdot [0.5 \cdot \log (0.5) + (1-0.5) \cdot \log (1- 0.5)] = -138.63 \)
Dove \( \displaystyle \bar y = \frac {\sum_{i=1}^n y_i}{n} = \frac {100}{200} = 0.5 \)
Anche questo metodo però ha le sue critiche e sono state proposte varie correzioni (come puoi vedere nel link..)
con riferimento ai 200 dati inventati più su, dovrebbe risultare una log-verosimiglianza massima pari a $-61.003$, quindi $-2logL=122.01$.
Risulta quindi un $"AIC"=122.01+2*2=126.01$ (nel nostro modello abbiamo $k=2$ parametri da stimare).
L'AIC è utile per confrontare diversi modelli tra loro: risulta preferibile quello con un AIC più basso.
Anche il test F è usato nell'articolo per confrontare i vari modelli con quello "base" di Lyman (il probit).
Ma nel tuo caso hai un unico modello, non vuoi fare un confronto con altri modelli.
Per valutare la bontà di adattamento potresti usare il test $\chi^2$.
La devianza residua non spiegata dal tuo modello (logit binario) è uguale a $-2logL=122.01$.
Potresti confrontare questo valore con una $\chi^2$ con 198 gradi di libertà (200-2, dove 2 sono i parametri del modello).
Risulta \( P(\chi^2_{198} > 122.01)=0.9999951 \)
Valori inferiori alla solita soglia $0.05$ indicano uno scarso adattamento. Nel nostro caso si avrebbe un buon adattamento.
C'è un però...
Nel caso di modello binario (è il nostro caso ma non quello dell'articolo Cancers), non è corretto usare questo test.
Siccome la dimensione dei gruppi e’ 1 (non abbiamo gruppi..), la devianza non è distribuita come una chi-quadrato, quindi non e’ possibile valutare la bonta’ di adattamento con quell'indice. Non vale più la convergenza alla chi-quadrato.
Ho letto che per ovvviare a ciò esistono varie alternative, tra cui il test di Hosmer–Lemeshow, che prevede comunque un raggruppamento dei dati in 10 classi..
Un atro modo potrebbe essere di utilizzare dei "Pseudo R-Squared", sulla falsariga di quello che si usa nel modello dei minimi quadrati.
Ad esempio, con gli stessi dati, il modello nullo ($\beta=0$, solo l'intercetta) ha una devianza di $-2logL=-2*(-138.63)=277.26$ (*)
Si avrebbe quindi (McFadden's): \( \displaystyle R^2=1-\frac{122.01}{277.26}=0.560 \) (più alto è meglio).
(*) Per il modello nullo puoi calcolare:
\( \displaystyle \log L=n \cdot [\bar y \cdot \log (\bar y) + (1-\bar y) \cdot \log (1- \bar y)] = 200 \cdot [0.5 \cdot \log (0.5) + (1-0.5) \cdot \log (1- 0.5)] = -138.63 \)
Dove \( \displaystyle \bar y = \frac {\sum_{i=1}^n y_i}{n} = \frac {100}{200} = 0.5 \)
Anche questo metodo però ha le sue critiche e sono state proposte varie correzioni (come puoi vedere nel link..)
In effetti come al solito ero andato avanti e quindi avevo omesso di aggiungere che sto provando ad implementare oltre al modello logit (Niemierko) che abbiamo analizzato fino ad adesso anche il modello probit (Lyman) che risponde alla formula:
$ NTCP=1/sqrt(2pi)int_(-oo )^(t) e^((-x^2/2)) dx $
con $t=(EUD-TD_50)/(m*TD_50)$
Per mia fortuna nell'articolo di Zhu c'era la funzione Log(L) del modello di Lyman già bella e pronta, e quindi mi sono limitato a "spalmarla" nel foglio elettronico dopo aver trovato la maniera di risolvere quell'integrale (che ho risolto usando la conversione alla funzione errore -calcolabile in Excel- così come l'ho trovata esplicitata nell'equazione 17 dell'articolo di Cancers).
Analizzando i dati test che abbiamo utilizzato fino ad ora mi ritrovo un AIC per Niemierko pari a 126.01, come dicevi tu, mentre per Lyman ho trovato AIC = 121.13. Il confronto lo vedi nella figura che segue:
Da quel poco che mi pare di aver capito in questo caso Lyman fitta meglio (lo si vede anche "ad occhio" soprattutto per i valori alti del grafico)
Riguardo le altre modalità di verifica del fitting mi sa che mi tocca studiare ancora un pò. Ad ogni modo Cenzo il tuo aiuto è stato insostituibile come al solito, davvero grazie e complimenti!!!
$ NTCP=1/sqrt(2pi)int_(-oo )^(t) e^((-x^2/2)) dx $
con $t=(EUD-TD_50)/(m*TD_50)$
Per mia fortuna nell'articolo di Zhu c'era la funzione Log(L) del modello di Lyman già bella e pronta, e quindi mi sono limitato a "spalmarla" nel foglio elettronico dopo aver trovato la maniera di risolvere quell'integrale (che ho risolto usando la conversione alla funzione errore -calcolabile in Excel- così come l'ho trovata esplicitata nell'equazione 17 dell'articolo di Cancers).
Analizzando i dati test che abbiamo utilizzato fino ad ora mi ritrovo un AIC per Niemierko pari a 126.01, come dicevi tu, mentre per Lyman ho trovato AIC = 121.13. Il confronto lo vedi nella figura che segue:
Da quel poco che mi pare di aver capito in questo caso Lyman fitta meglio (lo si vede anche "ad occhio" soprattutto per i valori alti del grafico)
Riguardo le altre modalità di verifica del fitting mi sa che mi tocca studiare ancora un pò. Ad ogni modo Cenzo il tuo aiuto è stato insostituibile come al solito, davvero grazie e complimenti!!!
"nickwing":
In effetti come al solito ero andato avanti e quindi avevo omesso di aggiungere che sto provando ad implementare oltre al modello logit (Niemierko) che abbiamo analizzato fino ad adesso anche il modello probit (Lyman) che risponde alla formula:
$ NTCP=1/sqrt(2pi)int_(-oo )^(t) e^((-x^2/2)) dx $
con $t=(EUD-TD_50)/(m*TD_50)$
Perfetto
"nickwing":
Per mia fortuna nell'articolo di Zhu c'era la funzione Log(L) del modello di Lyman già bella e pronta, e quindi mi sono limitato a "spalmarla" nel foglio elettronico dopo aver trovato la maniera di risolvere quell'integrale (che ho risolto usando la conversione alla funzione errore -calcolabile in Excel- così come l'ho trovata esplicitata nell'equazione 17 dell'articolo di Cancers).
Va benissimo l'articolo Zhu. Per risovere l'integrale è inutile passare per Cancers e per la "erf". (*)
In Excel puoi usare la funzione =DISTRIB.NORM.ST(t)
Bastano 3 (ma anche 2) colonne. Nella prima ci metti $t$ in funzione dei parametri $m$ e $"TD"_50$ (celle a parte) e della dose $"EUD"_i=x_i$
Nella seconda ci metti $"NTCP"_i=p_i$ calcolata come " =DISTRIB.NORM.ST(t) " (sarebbe la CDF normale standard).
Nella terza ci metti la logL come Zhu, in funzione di $p_i$ e della risposta $R_i=y_i$ binaria (0 o 1).
La somma di quest'ultima colonna è la log-Lik che occorre massimizzare variando i due parametri $m$ e $"TD"_50$.
Infine -è questo forse che ti ha indotto a passare per Cancers...- vale la relazione $\gamma_50=1/(m*sqrt(2*pi))=1/(0.195*sqrt(2*pi))=2.049$
(*) La erf mi sembra solo un inutile ingarbugliamento. Userei solo la CDF normale.
Per produrre il grafico hai la colonna $p_i$ (ordinata) che è corretta dopo la massimizzazione e $x_i$ ascissa.
"nickwing":
Analizzando i dati test che abbiamo utilizzato fino ad ora mi ritrovo un AIC per Niemierko pari a 126.01, come dicevi tu, mentre per Lyman ho trovato AIC = 121.13. [..] Lyman fitta meglio (lo si vede anche "ad occhio" soprattutto per i valori alti del grafico)
Confermo l'AIC del modello Lyman. Noto anche che hai cambiato un po' i gruppi (se ricordo bene, prima erano 4)
Aggiungo che se proviamo il modello logistico classico, in cui compare la variabile dose (e non il suo logaritmo) - vedi Betnzen tabella 1 Logistic Dose as covariate.. esce un AIC=119.07 quindi ancora meglio di Lyman, per lo meno su questi dati inventati.
E questa è un'immagine riassuntiva (anch'io ho aggiornato i gruppi
):
In definitiva, non ci vedo una grossa differenza tra i 3 modelli. Più o meno stanno tutti lì (con questi dati).
PS: ottimo lavoro nick
"cenzo":
(*) La erf mi sembra solo un inutile ingarbugliamento. Userei solo la CDF normale.
Non ho dubbi... Ma io non conoscevo la funzione DISTRIB.NORM.ST di Excel
"cenzo":
PS: ottimo lavoro nick
Non avrei risolto niente senza di te. Grazie mille!!!
"nickwing":
[quote="cenzo"](*) La erf mi sembra solo un inutile ingarbugliamento. Userei solo la CDF normale.
Non ho dubbi... Ma io non conoscevo la funzione DISTRIB.NORM.ST di Excel
[/quote]Epperò.. complimenti per come hai saputo aggirare l'ostacolo!
"nickwing":
Non avrei risolto niente senza di te. Grazie mille!!!
Figurati, alla prossima. Ciao nick
ancora ulteriori evoluzioni...
Sia $ TCP=0.5^(e^{((2*gamma_50)/ln2)*(1-((EUD)/(TCD_50))) }) $ la funzione che descrive la probabilità di controllare un tumore espressa in maniera "poissoniana". Nel senso che questa è la funzione di distribuzione cumulativa poissoniana di controllare il tumore per una data dose EUD somministrata allo stesso. Siamo sempre alle solite... devo ritrovare la funzione di verosimiglianza in maniera da calcolare la MLE di questa funzione relativa ai parametri $TCD50$ e $gamma_50$ (gli stessi che utilizzavo negli altri modelli). La reference che utilizza questo modello è la seguente.
Sia $ TCP=0.5^(e^{((2*gamma_50)/ln2)*(1-((EUD)/(TCD_50))) }) $ la funzione che descrive la probabilità di controllare un tumore espressa in maniera "poissoniana". Nel senso che questa è la funzione di distribuzione cumulativa poissoniana di controllare il tumore per una data dose EUD somministrata allo stesso. Siamo sempre alle solite... devo ritrovare la funzione di verosimiglianza in maniera da calcolare la MLE di questa funzione relativa ai parametri $TCD50$ e $gamma_50$ (gli stessi che utilizzavo negli altri modelli). La reference che utilizza questo modello è la seguente.
Se al paziente soggetto alla dose \(x\) associamo una risposta \(y\) binaria (0/1) si ha sempre una famiglia binomiale \( Y \sim Bin(1,p) \), TCP è la probabilità $p$, la log-Lik è sempre la stessa:
\( \displaystyle \log(L)=\sum_{i=1}^n y_i \log(p_i)+(1-y_i) \log(1-p_i) \)
Quello che cambia è la funzione "link", in queso caso:
Con i soliti dati mi risulta \( \text{TCD}_{50}=37.486 \), \( \gamma_{50}=1.641 \), \( \log(L)=-65.144 \) e \( \text{AIC}=134.29 \).
In pratica questo è il peggiore modello di quelli analizzati finora rispetto ai dati da fittare.
\( \displaystyle \log(L)=\sum_{i=1}^n y_i \log(p_i)+(1-y_i) \log(1-p_i) \)
Quello che cambia è la funzione "link", in queso caso:
"nickwing":
$p=0.5^(e^{((2*gamma_50)/ln2)*(1-((EUD)/(TCD_50))) }) $
Con i soliti dati mi risulta \( \text{TCD}_{50}=37.486 \), \( \gamma_{50}=1.641 \), \( \log(L)=-65.144 \) e \( \text{AIC}=134.29 \).
In pratica questo è il peggiore modello di quelli analizzati finora rispetto ai dati da fittare.
Vero, nel frattempo ci avevo provato anche io ad usare la stessa funzione... con i medesimi risultati
qui i diagrammi...
qui i diagrammi...
Nel grafico però non vedo il modello logistico classico (Dose ad covariate) che ha un AIC=119.07
(la funzione link è riportata su Bentzen).
Inoltre ho provato anche un altro modello - complementary log log (*) - con funzione link (espressa rispetto ai tuoi parametri):
\( \displaystyle p=1- 0.5^{ e^{ \frac{-2 \gamma_{50}}{\log2} \left( 1-\frac{\text{EUD}}{\text{TCD}_{50}} \right) } } \)
Con questo risulta TCD50=40.012, γ50=2.408, log(L)=−55.287 e AIC=114.57 (il più basso finora)
(*) La sua espressione originaria è \( \log(-\log(1-p))=a+bx \)
(la funzione link è riportata su Bentzen).
Inoltre ho provato anche un altro modello - complementary log log (*) - con funzione link (espressa rispetto ai tuoi parametri):
\( \displaystyle p=1- 0.5^{ e^{ \frac{-2 \gamma_{50}}{\log2} \left( 1-\frac{\text{EUD}}{\text{TCD}_{50}} \right) } } \)
Con questo risulta TCD50=40.012, γ50=2.408, log(L)=−55.287 e AIC=114.57 (il più basso finora)
(*) La sua espressione originaria è \( \log(-\log(1-p))=a+bx \)
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo