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
Guarda quello che dici è tutto giusto, anche ben spiegato. Una piccola precisazione: quando dici che devi usare il MLE quella è una maniera di stimare, puoi farlo anche in altre maniere. Comunque quella di massima verosimiglianza quando possibili è una delle migliori.
Come ti ho detto: quando possibile.
Infatti richiede la conoscenza della pdf. Da quello che scrivi l'unica cosa che abbiamo è $f(x)$, deve quindi venire da la. Non è che è proprio quella? Hai provato
ad integrarla per vedere se in un certo dominio la f ha integrale 1 (o conosci il range che x può assumere?)
Come ti ho detto: quando possibile.
Infatti richiede la conoscenza della pdf. Da quello che scrivi l'unica cosa che abbiamo è $f(x)$, deve quindi venire da la. Non è che è proprio quella? Hai provato
ad integrarla per vedere se in un certo dominio la f ha integrale 1 (o conosci il range che x può assumere?)
Grazie della risposta DajeForte! Per farla breve, faccio il medico 
, quel modello rappresenta una curva dose risposta di un tumore ad una data dose di radiazioni utilizzate per distruggerlo: $x$ rappresenta la dose, $f(x)$ è l'effetto, che può variare fra 0 ed 1 (ossia fra l'assenza di risposta e la percentuale della risposta completa 100%). In letteratura ci sono pubblicazioni che hanno già dimostrato come il metodo della MLE sia il più affidabile nel ricostruire i fenomeni osservati descritti mediante quel modello. Sto lavorando ad un sistema che implementi la possibilità di utilizzare quel genere di modello nella pratica clinica corrente (ad oggi in realtà la modellistica radiobiologica è stata sempre abbastanza al di fuori dai circuiti di utilizzo "quotidiano") e quindi, passo dopo passo, mi sto cimentando con la risoluzione del "problema", anche se non ho la preparazione universitaria necessaria per fare tutto da me
Quanto alla tua ipotesi, ossia che la funzione $f(x)$ sia la mia $PDF$ è un'idea che avevo considerato. Il dominio della $x$ è $0
Grazie ancora della risposta.
, quel modello rappresenta una curva dose risposta di un tumore ad una data dose di radiazioni utilizzate per distruggerlo: $x$ rappresenta la dose, $f(x)$ è l'effetto, che può variare fra 0 ed 1 (ossia fra l'assenza di risposta e la percentuale della risposta completa 100%). In letteratura ci sono pubblicazioni che hanno già dimostrato come il metodo della MLE sia il più affidabile nel ricostruire i fenomeni osservati descritti mediante quel modello. Sto lavorando ad un sistema che implementi la possibilità di utilizzare quel genere di modello nella pratica clinica corrente (ad oggi in realtà la modellistica radiobiologica è stata sempre abbastanza al di fuori dai circuiti di utilizzo "quotidiano") e quindi, passo dopo passo, mi sto cimentando con la risoluzione del "problema", anche se non ho la preparazione universitaria necessaria per fare tutto da me Quanto alla tua ipotesi, ossia che la funzione $f(x)$ sia la mia $PDF$ è un'idea che avevo considerato. Il dominio della $x$ è $0
Grazie ancora della risposta.
La situazione non è semplice. Ti dico quello che mi viene in mente.
Considera che non ci sia aleatorietà. Se lanci con x=50 avrai una reazione di f(50). Innanzitutto anche qua non è determinato il valore ma dipende dai parametri.
In quest'ottica quello che dovresti fare è interpolazione dei dati, ovvero se hai delle osservazioni $(x_1,f(x_1))$,...,$(x_n,f(x_n) ) $ devi trovare i parametri in maniera tale che la curva (x,f(x)) fitti, al meglio secondo un certo criterio, i dati.
Se c'è aleatorietà (diciamo che non c'è sulla quantità di radiazioni, ovvero se decido di sparare 50 sparo effettivamente 50), sparando 50 il risultato sarà un qualsiasi valore tra 0 e 1 e questo avrà una distribuzione. Ora come fare MLE ancora non mi risulta chiaro, ovvero come detrminare la PDE o magari cercare di fare MLE direttamente su f (ma non saprei ancora bene come). Te hai scritto che:"n letteratura ci sono pubblicazioni che hanno già dimostrato come il metodo della MLE sia il più affidabile nel ricostruire i fenomeni osservati descritti mediante quel modello.". Da questi non se ne ricava nulla?
Magari sentiamo anche il parere di qualcun'altro se si vuole aggiungere alla discussione perchè il problema è interessante.
Considera che non ci sia aleatorietà. Se lanci con x=50 avrai una reazione di f(50). Innanzitutto anche qua non è determinato il valore ma dipende dai parametri.
In quest'ottica quello che dovresti fare è interpolazione dei dati, ovvero se hai delle osservazioni $(x_1,f(x_1))$,...,$(x_n,f(x_n) ) $ devi trovare i parametri in maniera tale che la curva (x,f(x)) fitti, al meglio secondo un certo criterio, i dati.
Se c'è aleatorietà (diciamo che non c'è sulla quantità di radiazioni, ovvero se decido di sparare 50 sparo effettivamente 50), sparando 50 il risultato sarà un qualsiasi valore tra 0 e 1 e questo avrà una distribuzione. Ora come fare MLE ancora non mi risulta chiaro, ovvero come detrminare la PDE o magari cercare di fare MLE direttamente su f (ma non saprei ancora bene come). Te hai scritto che:"n letteratura ci sono pubblicazioni che hanno già dimostrato come il metodo della MLE sia il più affidabile nel ricostruire i fenomeni osservati descritti mediante quel modello.". Da questi non se ne ricava nulla?
Magari sentiamo anche il parere di qualcun'altro se si vuole aggiungere alla discussione perchè il problema è interessante.
E' vero, è un problema molto interessante.
Ho cercato di informarmi leggendo questo documento: http://www-naweb.iaea.org/nahu/dmrp/pdf ... pter14.pdf
(in particolare i paragrafi da 14.6 a 14.9), in cui compare la curva dose-risposta a forma di sigmoide.
Anche la funzione $f(x)=1/(1+(\theta/x)^\gamma)$ assume una forma a sigmoide.
Ne ho disegnato il grafico per diversi valori di $theta$ (a sinistra, con $gamma=4$) e di $gamma$ (a destra, con $theta=50$):
Se ho capito, l'ordinata rappresenta la probabilità di controllo del tumore (FFBF=Freedom from Biochemical Failure), in qualche modo definita.
Provo a spiegare quello che ho compreso, magari serve per avere conferma.
In pratica somministriamo una dose di raggi (Gray) e in risposta a tale dose c'è una certa probabilità di successo.
Diciamo che possono verificarsi solo due casi: successo o insuccesso, da stabilire in base a qualche criterio (una soglia?).
In pratica credo che la stima di massima verosimiglianza sia per un modello tipo quello di regressione logistica.
Premetto che ho scarse conoscenze di analisi multivariata perciò mi aspetto (benvengano!) critiche su quanto segue
Diciamo di avere $m$ gruppi di pazienti ($1,2,...,i,...,m$), ciascuno di numerosità $n_i$.
Ad un certo gruppo somministriamo una dose $x_i$ ed osserviamo successo in $y_i$ pazienti ed insuccesso nei restanti $n_i-y_i$.
Assumiamo che la variabile [tex]$Y_i\sim B(n_i,p_i)$[/tex] sia distribuita come una binomiale (pazienti indipendenti).
Ora occorre stabilire quale legge assumiamo per $p_i$, e qui entra in gioco la funzione $f(x)$ assegnata.
Nel modello di regressione logistica (con una sole variabile esplicativa) si assume $"logit"(p)=ln(p/(1-p))=a+bx$, da cui $p=1/(1+e^(-(a+bx)))$
Ne risulta una funzione di verosimiglianza $L(a,b)=prod_{i=1}^{m}P(Y_i=y_i)=\prod_{i=1}^{m}((n_i),(y_i))p_i^(y_i)(1-p_i)^(n_i-y_i)$
dove $y_i$ e $n_i$ sono i dati osservati e $p_i$ dipende da $x_i$ e dai parametri $a$ e $b$.
Ho letto che la massimizzazione della logL si può fare con un algoritmo di Newton-Raphson modificato (detto Fischer score), ma ciò si può demandare senz'altro a pacchetti statistici (pensavo ad R, open source, a meno che nickwing abbia confidenza con altri programmi).
Ho dato uno sguardo a questi articoli (i primi che ho trovato, giusto per farmi un'idea):
http://www.redjournal.org/article/S0360 ... 8/fulltext
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2770596/
e mi sembra di poter confermare l'uso del modello di regressione logistica.
Problema: il modello dato da nickwing $p=1/(1+(\theta/x)^\gamma)$ risulta un logit ?
Direi di no, mi risulta (salvo errori): $"logit"(p)=ln(p/(1-p))=-gamma*ln(theta)+gamma*ln(x)=alpha+beta*ln(x)$ (non viene una funzione lineare nella x)
Credo che con la stessa metodologia si possa utilizzare comunque anche quel modello, anzi... si potrebbero anche confrontare i due
Penso sia utile dare anche uno sguardo alle pubblicazioni citate da nickwing, per avere una panoramica più chiara.
@Daje: passo la palla, che ne pensi ?
Ho cercato di informarmi leggendo questo documento: http://www-naweb.iaea.org/nahu/dmrp/pdf ... pter14.pdf
(in particolare i paragrafi da 14.6 a 14.9), in cui compare la curva dose-risposta a forma di sigmoide.
Anche la funzione $f(x)=1/(1+(\theta/x)^\gamma)$ assume una forma a sigmoide.
Ne ho disegnato il grafico per diversi valori di $theta$ (a sinistra, con $gamma=4$) e di $gamma$ (a destra, con $theta=50$):
Se ho capito, l'ordinata rappresenta la probabilità di controllo del tumore (FFBF=Freedom from Biochemical Failure), in qualche modo definita.
Provo a spiegare quello che ho compreso, magari serve per avere conferma.
In pratica somministriamo una dose di raggi (Gray) e in risposta a tale dose c'è una certa probabilità di successo.
Diciamo che possono verificarsi solo due casi: successo o insuccesso, da stabilire in base a qualche criterio (una soglia?).
In pratica credo che la stima di massima verosimiglianza sia per un modello tipo quello di regressione logistica.
Premetto che ho scarse conoscenze di analisi multivariata perciò mi aspetto (benvengano!) critiche su quanto segue
Diciamo di avere $m$ gruppi di pazienti ($1,2,...,i,...,m$), ciascuno di numerosità $n_i$.
Ad un certo gruppo somministriamo una dose $x_i$ ed osserviamo successo in $y_i$ pazienti ed insuccesso nei restanti $n_i-y_i$.
Assumiamo che la variabile [tex]$Y_i\sim B(n_i,p_i)$[/tex] sia distribuita come una binomiale (pazienti indipendenti).
Ora occorre stabilire quale legge assumiamo per $p_i$, e qui entra in gioco la funzione $f(x)$ assegnata.
Nel modello di regressione logistica (con una sole variabile esplicativa) si assume $"logit"(p)=ln(p/(1-p))=a+bx$, da cui $p=1/(1+e^(-(a+bx)))$
Ne risulta una funzione di verosimiglianza $L(a,b)=prod_{i=1}^{m}P(Y_i=y_i)=\prod_{i=1}^{m}((n_i),(y_i))p_i^(y_i)(1-p_i)^(n_i-y_i)$
dove $y_i$ e $n_i$ sono i dati osservati e $p_i$ dipende da $x_i$ e dai parametri $a$ e $b$.
Ho letto che la massimizzazione della logL si può fare con un algoritmo di Newton-Raphson modificato (detto Fischer score), ma ciò si può demandare senz'altro a pacchetti statistici (pensavo ad R, open source, a meno che nickwing abbia confidenza con altri programmi).
Ho dato uno sguardo a questi articoli (i primi che ho trovato, giusto per farmi un'idea):
http://www.redjournal.org/article/S0360 ... 8/fulltext
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2770596/
e mi sembra di poter confermare l'uso del modello di regressione logistica.
Problema: il modello dato da nickwing $p=1/(1+(\theta/x)^\gamma)$ risulta un logit ?
Direi di no, mi risulta (salvo errori): $"logit"(p)=ln(p/(1-p))=-gamma*ln(theta)+gamma*ln(x)=alpha+beta*ln(x)$ (non viene una funzione lineare nella x)
Credo che con la stessa metodologia si possa utilizzare comunque anche quel modello, anzi... si potrebbero anche confrontare i due
Penso sia utile dare anche uno sguardo alle pubblicazioni citate da nickwing, per avere una panoramica più chiara.
@Daje: passo la palla, che ne pensi ?
Hai colto nel segno Cenzo . Ora sono al mare Domani o lunedì posto qualcuno di quei lavori.
Domani o lunedì posto qualcuno di quei lavori.
Tornando al problema: in radiobiologia sono sostanzialmente due i modelli DOSE/RISPOSTA che sono stati formalizzati (ed ampiamente utilizzati) per descrivere il controllo di malattia (o l'insorgenza di effetti collaterali) a seguito della somministrazione di una data dose. Se $x$ è la dose il primo modello è il seguente:
(1) $ f(x)=1/(1+(theta/x)^rho) $ detto modello generalizzato di Niemierko (valido sia per tumori che per organi sani)
il secondo è il seguente:
(2) $ f(x)=1/sqrt(2pi) int_(-oo )^(t) e^{(-x^2/2)}dx $ detto di Lyman Kutcher Burman (valido per gli organi sani)
con (3) $ t=(D-TD50(v))/(m*TD50(v)) $ laddove $TD50$ è la dose che relativamente al volume $v$ dell'organo irradiato produce nel 50% dei pazienti l'effetto considerato, e questi ultimi sono osservazioni sperimentali. In sostanza questo termine serve a definre il limite alto dell'integrale che è funzione della dose nella formula (2), infine la $m$ è una costante che serve ad indicare come varia la risposta in funzione del volume irradiato. Più tardi vi allego una pubblicazione in cui viene utilizzata l'equazione (2) per trovare questi parametri ed un'altra in cui si parla genericamente dei metodi per "fittare le costanti" alle curve dose / risposta, laddove però le funzioni dose / risposta considerate sono funzioni sigmoidi diverse da quelle da me indicate in (1) e (2).
(1) $ f(x)=1/(1+(theta/x)^rho) $ detto modello generalizzato di Niemierko (valido sia per tumori che per organi sani)
il secondo è il seguente:
(2) $ f(x)=1/sqrt(2pi) int_(-oo )^(t) e^{(-x^2/2)}dx $ detto di Lyman Kutcher Burman (valido per gli organi sani)
con (3) $ t=(D-TD50(v))/(m*TD50(v)) $ laddove $TD50$ è la dose che relativamente al volume $v$ dell'organo irradiato produce nel 50% dei pazienti l'effetto considerato, e questi ultimi sono osservazioni sperimentali. In sostanza questo termine serve a definre il limite alto dell'integrale che è funzione della dose nella formula (2), infine la $m$ è una costante che serve ad indicare come varia la risposta in funzione del volume irradiato. Più tardi vi allego una pubblicazione in cui viene utilizzata l'equazione (2) per trovare questi parametri ed un'altra in cui si parla genericamente dei metodi per "fittare le costanti" alle curve dose / risposta, laddove però le funzioni dose / risposta considerate sono funzioni sigmoidi diverse da quelle da me indicate in (1) e (2).
"nickwing":
Più tardi vi allego una pubblicazione in cui viene utilizzata l'equazione (2) per trovare questi parametri ed un'altra in cui si parla genericamente dei metodi per "fittare le costanti" alle curve dose / risposta, laddove però le funzioni dose / risposta considerate sono funzioni sigmoidi diverse da quelle da me indicate in (1) e (2).
Ciao nick,
ad un primo rapidissimo sguardo direi che le sigmoidi da te indicate sono entrambe considerate nel secondo articolo (Bentzen):
La (1) corrisponde al modello logistico con variabile esplicativa il logaritmo della dose (paragrafo 3.4 Suit 1965).
Infatti se vedi la tabella 1, riga "Logistic", colonna "Logarithm of dose as covariate", compare proprio la tua prima $f(x)$.
La (2) dovrebbe corrispondere ad un modello probit (paragrafo 3.3 in cui si cita Lyman 1985). In particolare mi sembra che la tua (2) corrisponda, in tabella 1, alla riga "Probit", colonna "Dose as covariate".
Ci hai dato un bel po' di roba da leggere
.. Tra i due modelli, quale sei maggiormente orientato ad utilizzare ?
Al momento sono orientato ad utilizzare il modello (1) per un semplice motivo: ho realizzato un software che utilizza questo modello per rappresentare le curve dose risposta. Il calcolo delle funzioni (uso Delphi per la programmazione) è decisamente più semplice che fare una integrazione numerica quale sarebbe richiesta per risolvere l'equazione (2). Lo step successivo (per la pubblicazione) sarebbe quello di avere la possibilità di metterci dentro dei dati clinici su serie reali di pazienti per ottenere le costanti da applicare nel modello. Insomma, una sorta di "tool" da mettere in mano ai medici che notoriamente (io appartenendo a tale categoria) sbadigliano non poco quando c'è da vedere formule e numeri. Una sorta di "coltellino svizzero", non so se mi sono spiegato
In effetti rileggendo il lavoro hai ragione, i due modelli sono già citati nella pubblicazione di Bentzen, cambia solo la forma. Poi sai com'è, chi come me non ha molta dimistichezza con i numeri non riconosce immediatamente certe "similitudini"...
Ti allego una schermata del software tanto per capire. La costante $theta$ nel programma è indicata come $TD50$, $gamma$ invece compare all'esponente della formula originaria come $4*gamma$. Il modello nella sua definizione originale è il seguente: $NTCP=1/(1+((TCD50)/(EUD))^(4*gamma50))$, laddove $EUD=x$ e $TCD50=theta$.
In effetti rileggendo il lavoro hai ragione, i due modelli sono già citati nella pubblicazione di Bentzen, cambia solo la forma. Poi sai com'è, chi come me non ha molta dimistichezza con i numeri non riconosce immediatamente certe "similitudini"...
Ti allego una schermata del software tanto per capire. La costante $theta$ nel programma è indicata come $TD50$, $gamma$ invece compare all'esponente della formula originaria come $4*gamma$. Il modello nella sua definizione originale è il seguente: $NTCP=1/(1+((TCD50)/(EUD))^(4*gamma50))$, laddove $EUD=x$ e $TCD50=theta$.
"nickwing":
Il modello nella sua definizione originale è il seguente: $NTCP=1/(1+((TCD50)/(EUD))^(4*gamma50))$, laddove $EUD=x$ e $TCD50=theta$.
Bene, allora restiamo su questo modello.
Come già visto nel mio primo intervento quel modello equivale ad un logit in funzione del logaritmo della dose $x$:
$"logit"(p)=log(p/(1-p))=\alpha+\beta*log(x)$
e si hanno le seguenti uguaglianze: $\{(\gamma_{50}=\beta/4),("TCD"_{50}=e^(-\alpha/\beta)):}$
La curva fittata passa per il punto $G=("TCD"_{50},1/2)$.
In tale punto G ha una pendenza: $(dp)/(dx)=\gamma_{50}/("TDC"_{50})$
La stima dei parametri $\alpha$ e $\beta$, sviluppando l'espressione della verosimiglianza riportata nel mio primo intevento, si ottiene massimizzando la seguente espressione:
$\alpha*\sum_{i=1}^{m}y_i+\beta*\sum_{i=1}^{m}log(x_i)*y_i-\sum_{i=1}^{m}n_i*log(1+e^{\alpha+\beta*log(x_i)}) \ \(1) $
Ricordo che abbiamo $m$ gruppi di osservazioni; al gruppo $i$, di numerosità $n_i$, è stata somministrata una dose $x_i$, osservando $y_i$ successi e $n_i-y_i$ insuccessi.
La massimizzazione della (1) andrebbe condotta ponendo uguale a zero le due derivate parziali rispetto ai due paramertri $\alpha$ e $\beta$.
Ne esce un sistema non lineare di due equazioni in due incognite, risolvibile col metodo noto come "Fischer scoring".
Però, secondo me, non ti conviene imbarcarti nell'implementazione di quel metodo.
Ti conviene, invece, forse, sfruttare dei pacchetti che hanno delle routine dedicate.
Ad esempio puoi sfruttare il pacchetto di statistica open source R (click).
Segue un esempio, all'osso, di quello che potresti fare (*):
# Leggo e preparo i dati
data (menarche, package = "MASS")
attach(menarche)
out <- cbind(y = Menarche, noty = Total-Menarche)
# Modello
mo.menar <- glm(out ~ log(Age), family = binomial)
# Parametri modello
tcd50 <- exp(-coef(mo.menar)[1]/coef(mo.menar)[2])
gamma50 <- coef(mo.menar)[2]/4 # FINE
# Qui faccio vedere alcuni risultati numerici
> tcd50
(Intercept)
12.96293
> gamma50
log(Age)
5.333216
> summary(mo.menar)
Call:
glm(formula = out ~ log(Age), family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5649 -0.6548 -0.2702 0.7318 1.7758
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -54.6568 1.9757 -27.66 <2e-16 ***
log(Age) 21.3329 0.7684 27.76 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3693.884 on 24 degrees of freedom
Residual deviance: 22.077 on 23 degrees of freedom
AIC: 110.13
Number of Fisher Scoring iterations: 4
# Plotto grafici
plot(Age, Menarche/Total, pch=19)
x <- data.frame(Age = seq(9, 18, len=100))
lines(x$Age, predict(mo.menar, x, type = 'response'), col='red')
plot(Age, log((Menarche+0.5)/(Total-Menarche+0.5)), pch = 19, xlab = 'Age', ylab = 'Menarca (logit)')
lines(x$Age, predict(mo.menar, x), col='red')
Ed ecco i grafici risultanti:
La mia idea è che potresti, eventualmente, cercare di fare un porting tra Delphi ed R.
Ti passo questi due link che potrebbero, forse, servire alla causa:
http://www.stats.uwo.ca/faculty/murdoch ... ascal.html
http://r.789695.n4.nabble.com/R-Calling ... 13930.html
(*) I dati sul Menarca e il codice utilizzato nell'esempio di R sono presi da qui;
(trovi anche il codice da usare per eseguire altri test statistici sulla regressione)
PS: Nel grafico che hai allegato per una dose di 45 Gy leggo una NTCP di quasi il 90%.
Non è un po' alto ? Forse non ho compreso il ruolo del frazionamento della dose in 25 frazioni da 1.8 Gy?
Ciao Cenzo! Ti ringrazio della risposta, nel frattempo in verità ero andato avanti e sono arrivato ad una conclusione, la seguente:
considerato il modello di curva dose/rispota un modello logit ho posto questa uguaglianza:
(1) $ 1/(1+e^(-x))=1/(1+((TCD50)/(EUD))^(4*gamma))$
da qui ottengo che : (2) $ e^(-x)=((TCD50)/(EUD))^(4*gamma) $
nella mia funzione $f(x)$ io vado ad individuare i parametri $TCD50$ e $gamma$, per farlo ho provato a considerare la funzione log-likelihood della distribuzione logit classica che sarebbe:
$ ln L=sum_(m = 1)^(N)y_m ln (1/(1+e^(-x_m)))+ (1-y_m) ln (1 - (1/(1+e^(-x_m)))) $ e considerando la (2) ho ottenuto la seguente equazione log-likelihood relativa al mio modello:
$ ln L=sum_(m = 1)^(N)y_m ln (1/(1+((TCD50)/(EUD_m))^(4*gamma)))+ (1-y_m) ln (1 - (1/(1+((TCD50)/(EUD_m))^(4*gamma)))) $
La soluzione al mio problema è trovare il massimo di questa equazione variando rispettivamente i due parametri $TCD50$ e $gamma$. Facile a dirsi... ma, sono andato avanti anche qui. Intanto ho creato una serie di dati in excel e ho usato il solver per risolvere il problema. E qui sembra funzionare tutto come vedi in questa figura dove la curva "fitta" i dati sperimentali dopo che il solver ha trovato i due parametri:
Per fare tutto in Delphi ho scaricato una libreria che si chiama DMath, realizzata da un francese (che sia benedetto!), all'interno ho trovato preconfezionate (da solo mai sarei riuscito ovviamente) delle routine per ricavare il minimo numerico in funzioni a più parametri. Quella che ho utilizzato impiega il metodo BFGS (Broyden-Fletcher-Goldfarb-Shanno), altrimenti detto Quasi-Newton. Memore dell'articolo (che ho anche uploadato) in cui si diceva che in MatLab avevano utilizzato questo metodo per cercare il massimo della funzione Likelihood ho semplicemente scritto un programmino che utilizzando gli stessi dati del foglio Excel ricava il minimo della funzione inversa Log-Likelihood. Questi sono i parametri ricavati dal risolutore di excel:
e questo e lo screenshot del programma fatto ad hoc per ricavare lo stesso calcolo della MLE a partire dagli stessi dati:
Diciamo che coincidono fino alla 5a cifra decimale per la $TCD50$ e fino alla sesta per il $gamma$. E per gli usi che ne devo fare va fin troppo bene così
Che ne pensi?
considerato il modello di curva dose/rispota un modello logit ho posto questa uguaglianza:
(1) $ 1/(1+e^(-x))=1/(1+((TCD50)/(EUD))^(4*gamma))$
da qui ottengo che : (2) $ e^(-x)=((TCD50)/(EUD))^(4*gamma) $
nella mia funzione $f(x)$ io vado ad individuare i parametri $TCD50$ e $gamma$, per farlo ho provato a considerare la funzione log-likelihood della distribuzione logit classica che sarebbe:
$ ln L=sum_(m = 1)^(N)y_m ln (1/(1+e^(-x_m)))+ (1-y_m) ln (1 - (1/(1+e^(-x_m)))) $ e considerando la (2) ho ottenuto la seguente equazione log-likelihood relativa al mio modello:
$ ln L=sum_(m = 1)^(N)y_m ln (1/(1+((TCD50)/(EUD_m))^(4*gamma)))+ (1-y_m) ln (1 - (1/(1+((TCD50)/(EUD_m))^(4*gamma)))) $
La soluzione al mio problema è trovare il massimo di questa equazione variando rispettivamente i due parametri $TCD50$ e $gamma$. Facile a dirsi... ma, sono andato avanti anche qui. Intanto ho creato una serie di dati in excel e ho usato il solver per risolvere il problema. E qui sembra funzionare tutto come vedi in questa figura dove la curva "fitta" i dati sperimentali dopo che il solver ha trovato i due parametri:
Per fare tutto in Delphi ho scaricato una libreria che si chiama DMath, realizzata da un francese (che sia benedetto!), all'interno ho trovato preconfezionate (da solo mai sarei riuscito ovviamente) delle routine per ricavare il minimo numerico in funzioni a più parametri. Quella che ho utilizzato impiega il metodo BFGS (Broyden-Fletcher-Goldfarb-Shanno), altrimenti detto Quasi-Newton. Memore dell'articolo (che ho anche uploadato) in cui si diceva che in MatLab avevano utilizzato questo metodo per cercare il massimo della funzione Likelihood ho semplicemente scritto un programmino che utilizzando gli stessi dati del foglio Excel ricava il minimo della funzione inversa Log-Likelihood. Questi sono i parametri ricavati dal risolutore di excel:
e questo e lo screenshot del programma fatto ad hoc per ricavare lo stesso calcolo della MLE a partire dagli stessi dati:
Diciamo che coincidono fino alla 5a cifra decimale per la $TCD50$ e fino alla sesta per il $gamma$. E per gli usi che ne devo fare va fin troppo bene così
Che ne pensi?
Ciao nick,
sono un po' perplesso sui passaggi con cui hai ottenuto la log-likelihood del tuo modello:
però l'espressione finale mi sembra corretta, a parte questo dubbio:
la quantità $(1-ym)$, in cui $ym$ presumo sia il numero di successi nel gruppo $m$,
non dovrebbe essere invece $(n_m-ym)$, dove $n_m$ è la numerosità del gruppo $m$ ?
Voglio dire, quella parentesi dovrebbe contenere il numero di insuccessi (totale-successi) del gruppo.
Io avevo seguito un'altra strada che mi permettesse di sfruttare la funzione glm del programma R.
Sono convinto che alla fine i risultati sarebbero identici.
Se vuoi avere un'ulteriore conferma posso fittare il tuo modello con R, passandomi i dati che hai usato.
L'unica cosa che aggiungo a favore di R è che avresti anche la possibilità di avere i p-value dei coefficienti, il test di nullità del coefficiente di regressione, la tavola di analisi della devianza, la bontà di adattamento e un intervallo di confidenza dei parametri.
Però, se ritieni non servano, allora direi che hai praticamente finito.
I miei complimenti per tutto il lavoro che hai fatto
sono un po' perplesso sui passaggi con cui hai ottenuto la log-likelihood del tuo modello:
"nickwing":
ho ottenuto la seguente equazione log-likelihood relativa al mio modello:
$ ln L=sum_(m = 1)^(N)ym ln (1/(1+((TCD50)/(EUDm))^(4*gamma)))+ (1-ym) ln (1 - (1/(1+((TCD50)/(EUDm))^(4*gamma)))) $
però l'espressione finale mi sembra corretta, a parte questo dubbio:
la quantità $(1-ym)$, in cui $ym$ presumo sia il numero di successi nel gruppo $m$,
non dovrebbe essere invece $(n_m-ym)$, dove $n_m$ è la numerosità del gruppo $m$ ?
Voglio dire, quella parentesi dovrebbe contenere il numero di insuccessi (totale-successi) del gruppo.
Io avevo seguito un'altra strada che mi permettesse di sfruttare la funzione glm del programma R.
Sono convinto che alla fine i risultati sarebbero identici.
Se vuoi avere un'ulteriore conferma posso fittare il tuo modello con R, passandomi i dati che hai usato.
L'unica cosa che aggiungo a favore di R è che avresti anche la possibilità di avere i p-value dei coefficienti, il test di nullità del coefficiente di regressione, la tavola di analisi della devianza, la bontà di adattamento e un intervallo di confidenza dei parametri.
Però, se ritieni non servano, allora direi che hai praticamente finito.
I miei complimenti per tutto il lavoro che hai fatto
Ho considerato il modello logit nella forma in cui viene espresso così:
$P=e^(x)/(1+e^(x))$ che poi ho ridotto a questo $P=1/(1+e^(-x))$, $P$ esprime la probabilità compresa fra 0 e 1 che si verifichi l'evento. Nel modello che ho inserito ovviamente $m$ andrebbe scritto a pedice (l'editor di formule non mi consente di farlo quindi $m$ non è un fattore ma un indice ovviamente). La $y$ esprime la probabilità osservata. Mi spiego facendo un esempio pratico. Ho un set di 4 pazienti, il primo presenta l'effeto con dose = 30, il secondo con dose = 40 il terzo con dose = 50, il quarto invece non presenta l'effetto (e mettiamo che abbia preso dose pari a 25). La $P$ al primo paziente è 0.25, al secondo paziente (essendo la dose maggiore) è pari al 0.5, al terzo paziente è pari a 0.75. Il quarto paziente incide riducendomi il valore di probabilità osservata che si presenti l'evento. Quest'ultimo è l'unico passaggio che mi lascia perplesso (non so cioè se metodologicamente è corretto). I valori che inserisco nel modello al posto di $y$ sono quindi questi: 0.25, 0.50 e 0.75. In effetti in questa riduzione non considero i pazienti che non presentano affatto l'effetto, se non nel darmi il valore complessivo della probabilità, ma è un dato possibile, potendo spesso queste curve non raggiungere nel range terapeutico di dose il valore del 100% (e meno male altrimenti significherebbe che esiste comunque un valore di dose al quale tutti i pazienti presentano l'effetto collaterale). Matematicamente poi il modello logit "non accetta" il 100%, nel senso che:
$ lim_(x -> oo ) e^(x)/(1+e^(x))=1 $ anche se poi nella pratica salendo con la dose oltre certi livelli terapeutici ad un certo punto tutti presenteranno l'effetto inevitabilmente.
La serie di esempio che ho utilizzato è la seguente, ovviamente la mia $x$ è la EUD (la dose) mentre la probabilità osservata è indicata sotto la colonna "p".
EUD
20
34
40
42
45
50.2
55
68
72
p
0.02
0.4
0.7
0.78
0.85
0.92
0.95
0.98
0.985
[size=150]Grazie mille davvero per l'aiuto!!![/size]
$P=e^(x)/(1+e^(x))$ che poi ho ridotto a questo $P=1/(1+e^(-x))$, $P$ esprime la probabilità compresa fra 0 e 1 che si verifichi l'evento. Nel modello che ho inserito ovviamente $m$ andrebbe scritto a pedice (l'editor di formule non mi consente di farlo quindi $m$ non è un fattore ma un indice ovviamente). La $y$ esprime la probabilità osservata. Mi spiego facendo un esempio pratico. Ho un set di 4 pazienti, il primo presenta l'effeto con dose = 30, il secondo con dose = 40 il terzo con dose = 50, il quarto invece non presenta l'effetto (e mettiamo che abbia preso dose pari a 25). La $P$ al primo paziente è 0.25, al secondo paziente (essendo la dose maggiore) è pari al 0.5, al terzo paziente è pari a 0.75. Il quarto paziente incide riducendomi il valore di probabilità osservata che si presenti l'evento. Quest'ultimo è l'unico passaggio che mi lascia perplesso (non so cioè se metodologicamente è corretto). I valori che inserisco nel modello al posto di $y$ sono quindi questi: 0.25, 0.50 e 0.75. In effetti in questa riduzione non considero i pazienti che non presentano affatto l'effetto, se non nel darmi il valore complessivo della probabilità, ma è un dato possibile, potendo spesso queste curve non raggiungere nel range terapeutico di dose il valore del 100% (e meno male altrimenti significherebbe che esiste comunque un valore di dose al quale tutti i pazienti presentano l'effetto collaterale). Matematicamente poi il modello logit "non accetta" il 100%, nel senso che:
$ lim_(x -> oo ) e^(x)/(1+e^(x))=1 $ anche se poi nella pratica salendo con la dose oltre certi livelli terapeutici ad un certo punto tutti presenteranno l'effetto inevitabilmente.
La serie di esempio che ho utilizzato è la seguente, ovviamente la mia $x$ è la EUD (la dose) mentre la probabilità osservata è indicata sotto la colonna "p".
EUD
20
34
40
42
45
50.2
55
68
72
p
0.02
0.4
0.7
0.78
0.85
0.92
0.95
0.98
0.985
[size=150]Grazie mille davvero per l'aiuto!!![/size]
In effetti una qualche misura della validità dei coefficienti trovati mi sarebbe molto d'aiuto nel dare consistenza al lavoro. Diciamo che per ora l'unica cosa che riesco a verificare (o meglio, che potrei riuscire a verificare) è la funzione "profile likelihood" che identifica i limiti del 95% dell'intervallo di confidenza. Per le altre misure devo "studiare" ancora...
Per quanto riguarda il pedice nelle formule puoi usare il carattere underscore: \$y_m\$ diventa $y_m$.
Il modello logit classico dovrebbe essere questo: $p=(e^(a+bx))/(1+e^(a+bx))=1/(1+e^(-(a+bx)))$ (i parametri sono due: a e b).
Per quanto riguarda il resto, vediamo se ho capito cosa hai fatto.
In pratica hai costruito una ECDF, cioè la funzione di distibuzione cumulativa empirica. Hai ragionato così:
Per $x<=25$ nessun paziente dei 4 presenta effetto, quindi $p(25)=0$
Per $x<=30$ un paziente dei 4 presenta effetto, quindi $p(30)=1/4=0.25$
Per $x<=40$ due pazienti dei 4 presentano effetto, quindi $p(40)=2/4=0.50$
Per $x<=50$ tre pazienti dei 4 presentano effetto, quindi $p(50)=3/4=0.75$
Resta in effetti il dubbio sul paziente sottoposto alla dose 25...
Forse si poteva anche ragionare escludendo tale paziente:
Per $x<=25$ nessun paziente presenta effetto, quindi $p(25)=0$
Per $x<=30$ un paziente dei 3 presenta effetto, quindi $p(30)=1/3=0.33$
Per $x<=40$ due pazienti dei 3 presentano effetto, quindi $p(40)=2/3=0.66$
Per $x<=50$ tre pazienti su3 presentano effetto, quindi $p(50)=1$
Quindi cerchi di fittare la ECDF con la curva data dal modello logistico ?
Così faresti solo un fitting, ma non applichi il modello logistico.
Forse però non ti ho capito bene, ed è solo un modo di rappresentare i dati sul grafico (vei dopo).
Le ipotesi del modello sono queste:
- hai $m$ gruppi di pazienti (per esempio 6 gruppi): ciascun gruppo è costituito anche da diverse unità, diciamo $n_i$ unità.
Ad esempio il primo gruppo è costituito da $n_1=8$ pazienti. Ogni paziente di questo gruppo è soggetto alla stessa dose $x_1=25$ Gy.
Di questi 8 nessuno presenta effetto: $y_1=0$.
Il secondo gruppo ($i=2$) è costituito da $n_2=5$ pazienti. Ogni paziente di questo gruppo è soggetto alla stessa dose $x_2=30$ Gy.
Di questi 5 uno solo presenta l'effetto: $y_2=1$.
Eccetera.. fino all'ultimo gruppo $i=m=6$.
Il modello logistico ipotizza che la distribuzione dei pazienti che presentano l'effetto sia binomiale.
Quindi la probabilità di osservare $y_i$ effetti in un gruppo di $n_i$ pazienti è $P(y_i)=((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i)$
La verosimiglianza è allora:
$L=\prod_{i=1}^{m}P(y_i)=\prod_{i=1}^{m}((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i)$
Il logaritmo della verosimiglianza:
$log(L)=log(\prod_{i=1}^{m}((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i))=\sum_{i=1}^m log((n_i),(y_i))+sum_{i=1}^m y_i log(p_i)+(n_i-y_i) log(1-p_i)$
Dovendo massimizzare la $log(L)$, quella prima somma la possiamo trascurare in quanto da un contributo costante che non dipende dai parametri del modello, che compaiono invece nelle probabilità $p_i$. La funzione da massimizzare è quindi:
$sum_{i=1}^m y_i log(p_i)+(n_i-p_i) log(1-p_i)$
Ora basta inserire al posto della $p_i$ l'espressione che si usa nel modello logit, nel tuo caso $p_i=1/(1+(\theta/x_i)^(4\gamma))$
(questo tuo modello equivale a dire $"logit"(p_i)=\alpha+\beta*log(x_i)$, vedi anche l'articolo di Brentzen )
Faccio la sostituzione e ottengo questa funzione da massimizzare:
$sum_{i=1}^m y_i log(1/(1+(\theta/x_i)^(4\gamma)))+(n_i-y_i) log(1-1/(1+(\theta/x_i)^(4\gamma)))$
Come puoi notare è simile a quella che avevi postato tu stesso qualche post addietro.
Quello che cambia è che figurano la numerosità $n_i$ di ciascun gruppo e il numero di "successi" $y_i$ contati in ciascun gruppo.
Questo, credo, dovrebbe essere il modo corretto di procedere.
Non mi sembra, però tu abbia i dati in questa forma.
Cioè, tu hai sottoposto ad ogni dose un solo paziente, giusto ?
E su questo paziente singolo puoi dire se hai registrato l'effetto o meno. Insomma, se ho capito, hai dei dati in forma binaria.
E' possibile adattare il modello logistico: il numero di successi questa volta si assume distribuito come una Bernoulliana invece che una Binomiale.
Ogni "gruppo" ha numerosità $n_i=1$ e il numero di successi $y_i$ è una variabile binaria. E' un solo paziente, quindi due soli casi:
$y_i=0$ non vedo effetto nel paziente (sottoposto a dose $x_i$)
$y_i=1$ vedo effetto nel paziente (sottoposto a dose $x_i$)
Il numero di gruppi $m$ coincide quindi col numero complessivo di pazienti.
La log(L) che si ottiene è proprio quella che avevi postato prima:
$sum_{i=1}^m y_i log(1/(1+(\theta/x_i)^(4\gamma)))+(1-y_i) log(1-1/(1+(\theta/x_i)^(4\gamma)))$
I dati di input che servono (anche per provare con R) sono gli $x_i$ per ogni paziente e il corrispondente valore di $y_i$, cioè 1 o 0 a seconda se a quella dose si è verificato l'effetto (successo) oppure no.
Il grafico che esce è di questo tipo.
Perdona lo sproloquio
"nickwing":
Ho considerato il modello logit nella forma in cui viene espresso così: $P=e^(x)/(1+e^(x))$
Il modello logit classico dovrebbe essere questo: $p=(e^(a+bx))/(1+e^(a+bx))=1/(1+e^(-(a+bx)))$ (i parametri sono due: a e b).
Per quanto riguarda il resto, vediamo se ho capito cosa hai fatto.
In pratica hai costruito una ECDF, cioè la funzione di distibuzione cumulativa empirica. Hai ragionato così:
Per $x<=25$ nessun paziente dei 4 presenta effetto, quindi $p(25)=0$
Per $x<=30$ un paziente dei 4 presenta effetto, quindi $p(30)=1/4=0.25$
Per $x<=40$ due pazienti dei 4 presentano effetto, quindi $p(40)=2/4=0.50$
Per $x<=50$ tre pazienti dei 4 presentano effetto, quindi $p(50)=3/4=0.75$
Resta in effetti il dubbio sul paziente sottoposto alla dose 25...
Forse si poteva anche ragionare escludendo tale paziente:
Per $x<=25$ nessun paziente presenta effetto, quindi $p(25)=0$
Per $x<=30$ un paziente dei 3 presenta effetto, quindi $p(30)=1/3=0.33$
Per $x<=40$ due pazienti dei 3 presentano effetto, quindi $p(40)=2/3=0.66$
Per $x<=50$ tre pazienti su3 presentano effetto, quindi $p(50)=1$
Quindi cerchi di fittare la ECDF con la curva data dal modello logistico ?
Così faresti solo un fitting, ma non applichi il modello logistico.
Forse però non ti ho capito bene, ed è solo un modo di rappresentare i dati sul grafico (vei dopo).
Le ipotesi del modello sono queste:
- hai $m$ gruppi di pazienti (per esempio 6 gruppi): ciascun gruppo è costituito anche da diverse unità, diciamo $n_i$ unità.
Ad esempio il primo gruppo è costituito da $n_1=8$ pazienti. Ogni paziente di questo gruppo è soggetto alla stessa dose $x_1=25$ Gy.
Di questi 8 nessuno presenta effetto: $y_1=0$.
Il secondo gruppo ($i=2$) è costituito da $n_2=5$ pazienti. Ogni paziente di questo gruppo è soggetto alla stessa dose $x_2=30$ Gy.
Di questi 5 uno solo presenta l'effetto: $y_2=1$.
Eccetera.. fino all'ultimo gruppo $i=m=6$.
Il modello logistico ipotizza che la distribuzione dei pazienti che presentano l'effetto sia binomiale.
Quindi la probabilità di osservare $y_i$ effetti in un gruppo di $n_i$ pazienti è $P(y_i)=((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i)$
La verosimiglianza è allora:
$L=\prod_{i=1}^{m}P(y_i)=\prod_{i=1}^{m}((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i)$
Il logaritmo della verosimiglianza:
$log(L)=log(\prod_{i=1}^{m}((n_i),(y_i))*p^(y_i)*(1-p_i)^(n_i-y_i))=\sum_{i=1}^m log((n_i),(y_i))+sum_{i=1}^m y_i log(p_i)+(n_i-y_i) log(1-p_i)$
Dovendo massimizzare la $log(L)$, quella prima somma la possiamo trascurare in quanto da un contributo costante che non dipende dai parametri del modello, che compaiono invece nelle probabilità $p_i$. La funzione da massimizzare è quindi:
$sum_{i=1}^m y_i log(p_i)+(n_i-p_i) log(1-p_i)$
Ora basta inserire al posto della $p_i$ l'espressione che si usa nel modello logit, nel tuo caso $p_i=1/(1+(\theta/x_i)^(4\gamma))$
(questo tuo modello equivale a dire $"logit"(p_i)=\alpha+\beta*log(x_i)$, vedi anche l'articolo di Brentzen )
Faccio la sostituzione e ottengo questa funzione da massimizzare:
$sum_{i=1}^m y_i log(1/(1+(\theta/x_i)^(4\gamma)))+(n_i-y_i) log(1-1/(1+(\theta/x_i)^(4\gamma)))$
Come puoi notare è simile a quella che avevi postato tu stesso qualche post addietro.
Quello che cambia è che figurano la numerosità $n_i$ di ciascun gruppo e il numero di "successi" $y_i$ contati in ciascun gruppo.
Questo, credo, dovrebbe essere il modo corretto di procedere.
Non mi sembra, però tu abbia i dati in questa forma.
Cioè, tu hai sottoposto ad ogni dose un solo paziente, giusto ?
E su questo paziente singolo puoi dire se hai registrato l'effetto o meno. Insomma, se ho capito, hai dei dati in forma binaria.
E' possibile adattare il modello logistico: il numero di successi questa volta si assume distribuito come una Bernoulliana invece che una Binomiale.
Ogni "gruppo" ha numerosità $n_i=1$ e il numero di successi $y_i$ è una variabile binaria. E' un solo paziente, quindi due soli casi:
$y_i=0$ non vedo effetto nel paziente (sottoposto a dose $x_i$)
$y_i=1$ vedo effetto nel paziente (sottoposto a dose $x_i$)
Il numero di gruppi $m$ coincide quindi col numero complessivo di pazienti.
La log(L) che si ottiene è proprio quella che avevi postato prima:
$sum_{i=1}^m y_i log(1/(1+(\theta/x_i)^(4\gamma)))+(1-y_i) log(1-1/(1+(\theta/x_i)^(4\gamma)))$
I dati di input che servono (anche per provare con R) sono gli $x_i$ per ogni paziente e il corrispondente valore di $y_i$, cioè 1 o 0 a seconda se a quella dose si è verificato l'effetto (successo) oppure no.
Il grafico che esce è di questo tipo.
Perdona lo sproloquio
"cenzo":
Cioè, tu hai sottoposto ad ogni dose un solo paziente, giusto ?
ESATTO!
"cenzo":
E su questo paziente singolo puoi dire se hai registrato l'effetto o meno. Insomma, se ho capito, hai dei dati in forma binaria.
E' possibile adattare il modello logistico: il numero di successi questa volta si assume distribuito come una Bernoulliana invece che una Binomiale.
Ogni "gruppo" ha numerosità $n_i=1$ e il numero di successi $y_i$ è una variabile binaria. E' un solo paziente, quindi due soli casi:
$y_i=0$ non vedo effetto nel paziente (sottoposto a dose $x_i$)
$y_i=1$ vedo effetto nel paziente (sottoposto a dose $x_i$)
Il numero di gruppi $m$ coincide quindi col numero complessivo di pazienti.
La log(L) che si ottiene è proprio quella che avevi postato prima:
$sum_{i=1}^m y_i log(1/(1+(\theta/x_i)^(4\gamma)))+(1-y_i) log(1-1/(1+(\theta/x_i)^(4\gamma)))$
I dati di input che servono (anche per provare con R) sono gli $x_i$ per ogni paziente e il corrispondente valore di $y_i$, cioè 1 o 0 a seconda se a quella dose si è verificato l'effetto (successo) oppure no.
E' proprio come hai descritto in queste righe. Ogni paziente diverso ha valori di dose differenti. Quindi la formula che ho impiegato dovrebbe essere quella giusta
No?"cenzo":
Perdona lo sproloquio
E di che... grazie mille a te della collaborazione!
"nickwing":
Quindi la formula che ho impiegato dovrebbe essere quella giusta
Si, purchè alle $y_i$ ci metti non le probabilità ma 1 o 0 a seconda se hai successo o no.
E di che... grazie mille a te della collaborazione!
Figurati, anche per me è stata l'occasione per approfondire la regressione logistica
"cenzo":
Si, purchè alle $y_i$ ci metti non le probabilità ma 1 o 0 a seconda se hai successo o no.
Io però nei valori $y_i$ ho inserito le probabilità cumulative, come cercavo di spiegare sopra, non 0 o 1...
perchè la NTCP, ossia la $y$ è la probabilità cumulativa di avere un effetto alla data dose...
"nickwing":
Io però nei valori $y_i$ ho inserito le probabilità cumulative, come cercavo di spiegare sopra, non 0 o 1...
Non mi sembra corretto, per lo meno non è coerente al modello logistico con dati binari.
Dovresti mettere 1 o 0.
Esempio dei dati:
ID_Paziente x(Gy) y 1 20 0 2 50 1 3 35 1 4 36 0 ...
La probabilità prevista dal modello la valuti dopo (aver trovato i due parametri) con la solita formula $p=1/(1+(\theta/x)^(4\gamma))$
"cenzo":
[quote="nickwing"]Io però nei valori $y_i$ ho inserito le probabilità cumulative, come cercavo di spiegare sopra, non 0 o 1...
Non mi sembra corretto, per lo meno non è coerente al modello logistico con dati binari.
Dovresti mettere 1 o 0.
Esempio dei dati:
ID_Paziente x(Gy) y 1 20 0 2 50 1 3 35 1 4 36 0 ...
La probabilità prevista dal modello la valuti dopo (aver trovato i due parametri) con la solita formula $p=1/(1+(\theta/x)^(4\gamma))$[/quote]
ed io "brutalmente" volevo vedere i miei dati rappresentati sul medesimo grafico della curva interpolante. Un pò come nell'articolo di Bentzen (grafico 5 a pagina 537, tieni conto che il modello si applica sia per il controllo tumorale -TCP- che per gli effetti collaterali -NTCP- ). Mi sa che devo ricominciare daccapo...
quindi credo che rappresentarli come probabilità cumulative per ogni paziente (con ogni paziente invariabilmente posizionato sopra il precedente) non sia formalmente corretto... 
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo