Introduzione ai modelli non lineari e approfondimenti con R

Questo tutorial è un prezioso contributo della nostra community ed è stato revisionato da DataCamp per chiarezza e accuratezza.

Ti piacerebbe condividere la tua esperienza? Ci farebbe piacere sentirti! Invia i tuoi articoli o idee tramite il nostro Modulo per i contributi della community.

I modelli non lineari non sono più popolari come un tempo. Oggi molti libri tralasciano ogni discussione sui modelli non lineari o li spiegano in modo limitato. Apparentemente, ci sono buoni motivi per questo. Primo, si può adattare una curva non lineare usando la regressione lineare, sia introducendo termini polinomiali di ordine superiore sia applicando una funzione base (in altre parole, applicando una trasformazione log o sqrt alle variabili esplicative). Secondo, si può preferire la semplicità rispetto a modelli complessi, dato che i modelli lineari sono facili da interpretare e spiegare.

Inoltre, per applicare un modello non lineare, è importante visualizzare i dati e avere un'idea della traiettoria, cosa non possibile in un contesto multivariato. Infine, i modelli non lineari sono più adatti a casi meccanicistici in cui i fenomeni sono puramente fisici o deterministici, il che non vale per i big data, dove le attività degli utenti sono registrate dalle macchine.

Tuttavia, nonostante tutti questi buoni motivi, la modellazione non lineare è una scienza affascinante che merita di essere esplorata e confrontata con la sua controparte lineare.

Modello di regressione non lineare

Le tecniche di stima per un modello non lineare sono esplicitamente iterative. Inoltre, c'è una differenza fondamentale nel modo in cui si applica una formula non lineare, che non è la stessa di una formula lineare. Le tre equazioni/formule seguenti sono tutte considerate lineari perché la relazione tra la variabile y e B₀, B₁ e B₂ è una linea retta.

Ciò che rende una formula non lineare è qualcosa di questo tipo:

L'uso di B₀ o B₁ è arbitrario. Tuttavia, le ultime due equazioni non hanno una relazione lineare tra la variabile y e B₀.

Sfortunatamente, questo rende i calcoli più difficili e non è garantita una soluzione numerica di successo al problema di stima. Quindi, i modelli non lineari arrivano con un avvertimento: ottenere risposte potrebbe non essere intuitivo o accurato.

Modello specializzato

Di solito, quando si applica un modello non lineare, si usano tipi specializzati. Questo rende l'adattamento di un modello non lineare più comodo o numericamente più semplice, ma deve essere concordato da un esperto del dominio e allineato ai dati. Per tipi specializzati intendiamo l'uso di un'equazione matematica predefinita. Ecco l'elenco delle equazioni non lineari più usate (alcune sono state omesse per brevità).

Strategia per un modello di regressione non lineare

1. Visualizza i dati per capire se una funzione matematica specializzata può spiegarli al meglio.
2. Osservando il grafico a occhio o usando una self-starter function (spiegata sotto), seleziona i valori iniziali dei parametri.
3. Adatta il modello ed esaminalo.
4. Esegui un'analisi statistica sul modello adattato per capire l'intervallo di confidenza dei parametri, l'errore residuo e la variazione spiegata.
5. Prova ad aggiungere più parametri al modello se aiuta a ridurre l'errore residuo.
6. Adatta sempre un modello lineare e confronta i risultati dei modelli non lineari e lineari.

Spiegare il modello non lineare con un esempio

Andremo ad adattare un modello non lineare specializzato su un dataset famoso, Puromycin, che riporta la velocità di reazione e l'informazione sulla concentrazione per una reazione enzimatica che coinvolge cellule trattate con il farmaco Puromycin. I dati includono due classi di cellule: una trattata con Puromycin e l'altra non trattata.

names(Puromycin)Puromycin'conc''rate''state'Offerto da

Le tre variabili in questo dataframe sono la concentrazione del substrato, la velocità iniziale della reazione e un indicatore di trattato o non trattato.

Tracciando la velocità rispetto alla concentrazione e etichettando i due livelli del trattamento, possiamo osservare un andamento generale. Le curve sono asintotiche e piuttosto distinte per i gruppi di cellule trattate e non trattate.

plot(Puromycin$conc, Puromycin$rate, type="n", xlab = "conc", ylab = "rate")text(Puromycin$conc, Puromycin$rate,      ifelse(Puromycin$state == "treated", "T", "U"))Offerto da

Poiché questi dati mostrano una forte relazione di Michaelis-Menten tra la velocità di reazione e la concentrazione, gli sperimentatori si aspettano di adattare la seguente equazione matematica.

Qui E è l'errore sperimentale. Usando la stessa analogia del modello lineare, con un modello non lineare specializzato cerchiamo di stimare i parametri Vmax e K che minimizzano la somma dei quadrati dei residui (tuttavia, ci sono altri metodi che si possono usare per stimare questi parametri):

Adattare un modello non lineare

Anche se conosciamo la relazione, adattare un modello non lineare non è così immediato perché il modello richiede ipotesi iniziali per V_max e K. Se non si sceglie una buona ipotesi iniziale, il modello può adattare male i dati. Ci sono due opzioni:

• Stimare a occhio il grafico per trovare i valori iniziali.
• Oppure usare una self-starter function (indicata di seguito).

Per ogni tipo di modello, abbiamo una corrispondente self-starter function che può essere usata per un'ipotesi iniziale.

Adattare un modello usando un'ipotesi iniziale

Usando il valore iniziale di Vmax = 160, K = 0,05 stimato a occhio dal grafico, si può usare la funzione R nls() per adattare i dati. Di solito, per la regressione lineare non è necessario specificare i parametri Vm o K, ma per un modello non lineare è diverso. Ogni algoritmo iterativo ha bisogno di un buon punto di partenza. Altrimenti, potrebbe non convergere. Ecco come adattiamo un modello non lineare su Puromycin con valori iniziali.

Purboth_1 <- nls(rate ~ (Vm)*conc/(K+conc), Puromycin, list(Vm=160, K=0.05))summary(Purboth_1)Formula: rate ~ (Vm) * conc/(K + conc)Parameters:Estimate Std. Error t value Pr(>|t|)Vm 190.80624 8.76459 21.770 6.84e-16 ***K 0.06039 0.01077 5.608 1.45e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 18.61 on 21 degrees of freedomNumber of iterations to convergence: 5Achieved convergence tolerance: 5.121e-06Offerto da

Calcolare l'errore residuo

Calcoliamo l'errore residuo come segue:

sse <- Purboth_1$m$deviance()sse7276.54697931423Offerto da

Calcolare la varianza spiegata

Per trovare la varianza spiegata da questo modello (in pratica l'R-quadro nella terminologia del modello lineare), dobbiamo trovare la varianza totale della velocità adattando un modello nullo, stimando la media complessiva ed estraendo la somma totale dei quadrati, sst, dal sommario del modello nullo per poi usarla nel calcolo della percentuale di varianza spiegata.

sse <- Purboth_1$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained85.3488323994681Offerto da

Questo modello spiega l'85,34% delle variazioni nei dati.

Adattare usando il modello con self-start function

Poiché stiamo usando il modello di Michaelis-Menten, si può usare la self-start function SSmicmen() al suo posto. Tuttavia, si ottengono gli stessi risultati.

Purboth_self <- nls(rate~SSmicmen(conc,Vm,K), Puromycin)summary(Purboth_self)Formula: rate ~ SSmicmen(conc, Vm, K)Parameters:Estimate Std. Error t value Pr(>|t|)Vm 190.80667 8.76462 21.770 6.84e-16 ***K 0.06039 0.01077 5.608 1.45e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 18.61 on 21 degrees of freedomNumber of iterations to convergence: 4Achieved convergence tolerance: 6.745e-06Offerto da

Calcolare l'errore residuo

sse <- Purboth_self$m$deviance()sse7276.54697947673Offerto da

Calcolare la varianza spiegata

Si può verificare la varianza spiegata dal modello.

sse <- Purboth_self$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained85.348832399141Offerto da

In generale, l'equazione adattata, cioè la relazione tra velocità e concentrazione, appare così:

Visualizzare l'adattamento

Si può visualizzare come il modello si comporta tracciando una linea adattata sui dati, che appare non lineare.

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state == "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state != "treated",]), col=2)Offerto da

Selezione di un modello migliore

Si può ottimizzare ulteriormente il modello aggiungendo più parametri. Osservando il grafico tra velocità e concentrazione si può notare che classi diverse (trattate vs non trattate) seguono traiettorie differenti, e la differenza di velocità (asse y) tra le due traiettorie può essere approssimata a 40.

plot(Puromycin$conc, Puromycin$rate, type="n", xlab = "conc", ylab = "rate")text(Puromycin$conc, Puromycin$rate,      ifelse(Puromycin$state == "treated", "+", "*"))Offerto da

In sostanza, stiamo cercando di adattare un modello condizionato sulla variabile state, che porta alla seguente equazione non lineare.

Il modello funziona meglio con un parametro aggiuntivo delV. Assegniamo a delV un valore iniziale di 40 e adattiamo il modello usando nls con tre parametri.

Purboth <- nls(rate ~ (Vm + delV*(state=="treated"))*conc/(K+conc), Puromycin,                list(Vm=160, delV=40, K=0.05))summary(Purboth)Formula: rate ~ (Vm + delV * (state == "treated")) * conc/(K + conc)Parameters:      Estimate Std. Error t value Pr(>|t|)    Vm   166.60397    5.80742  28.688  < 2e-16 ***delV  42.02591    6.27214   6.700 1.61e-06 ***K      0.05797    0.00591   9.809 4.37e-09 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 10.59 on 20 degrees of freedomNumber of iterations to convergence: 5 Achieved convergence tolerance: 8.93e-06Offerto da

Il sommario mostra che tutti e tre i parametri sono statisticamente significativi.

Calcolare l'errore residuo

Per questo modello, l'errore residuo si riduce in modo sostanziale:

sse <- Purboth$m$deviance()sse2240.89143883885Offerto da

Calcolare la varianza spiegata

Questo modello spiega una quota maggiore di varianza.

sse <- Purboth$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained95.4880142822744Offerto da

Verifica AIC

Si può anche eseguire il controllo AIC su entrambi i modelli, anche se oggi la pratica comune è usare la cross-validation; tuttavia, questa in generale richiede più dati. In generale, più basso è l'AIC, migliore è il modello.

AIC(Purboth)AIC(Purboth_1)178.591273238391203.680274951742Offerto da

Visualizzare l'adattamento

Si può visualizzare come il modello ottimizzato si comporta tracciando una linea adattata sui dati. Come previsto, il modello ottimizzato fornisce un adattamento migliore per le diverse classi.

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state == "treated",]$conc, predict(Purboth, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc, predict(Purboth, Puromycin[Puromycin$state != "treated",]), col=3)legend(0, 200, legend = c("treated", "untreated"), fill = c(2,3))Offerto da

Intervallo di confidenza dei parametri

Non ci si può basare solo sui p-value; l'intervallo di confidenza al 95% fornisce i valori dei coefficienti per ciascun parametro significativo al 2,5% e 97,5%. Nessuno di essi include 0.

confint(Purboth)Offerto da

Usando l'intervallo di confidenza, si può affermare statisticamente che c'è il 95% di probabilità che il vero valore del parametro, ad esempio Vmax, cada nell'intervallo riportato sopra, che per Vmax è (154,62; 179,25).

Profilazione

Si può profilare ulteriormente il modello per comprendere l'incertezza nella stima dei parametri esaminando direttamente la funzione obiettivo. Nel nostro modello abbiamo tre parametri, Vmax, delV, K; possiamo profilarli:

• Per ottenere l'intervallo di confidenza.
• E usare la funzione di profilo t, simile alla statistica t usata dal modello lineare. I grafici della funzione t di profilo forniscono gli intervalli di verosimiglianza per i singoli parametri e mostrano quanto non lineare sia la stima.

Tracciare l'intervallo di confidenza

Questo traccia l'intervallo di confidenza dei parametri per diversi valori di (che corrispondono a diversi intervalli di confidenza: 99%, 95%, 90%, 80% e 50%). È simile all'intervallo di confidenza (95%) fornito dalla funzione confint(object, parm, level = 0.95, ...)function.

Purbothpf = profile(Purboth)par( mfrow= c(2,2) )plot(Purbothpf)Offerto da

Tracciare _T vs parametri

Tracciare i tre parametri Vmax, delV, K rispetto a _T dà un'idea se l'adattamento è lineare o non lineare. Si osserva che la linea appare lineare per i parametri delV e Vm, e leggermente curva per K. Questo è atteso perché il parametro K mostra una relazione non lineare.

dataVm = data.frame(Purbothpf$Vm$par.vals)datadelV = data.frame(Purbothpf$delV$par.vals)dataK = data.frame(Purbothpf$K$par.vals)par( mfrow= c(2,2) )plot(dataVm$Vm,Purbothpf$Vm$tau, type="l", xlab = "Vm", ylab = "Tau")plot(datadelV$delV,Purbothpf$delV$tau, type="l", xlab = "delV", ylab = "Tau")plot(dataK$K,Purbothpf$K$tau, type="l", xlab = "K", ylab = "Tau")Offerto da

Adattare un modello lineare

Si può usare la stessa funzione della libreria R nls() per adattare un modello lineare come mostrato sotto. È una buona abitudine provare sempre un modello lineare per interpretabilità e confronto.

PurbothRealLinear <- nls(rate ~ i + m*conc, Puromycin, list(i=160, m=0.05) )summary(PurbothRealLinear)Formula: rate ~ i + m * concParameters:  Estimate Std. Error t value Pr(>|t|)    i    93.92       8.00   11.74 1.09e-10 ***m   105.40      16.92    6.23 3.53e-06 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 28.82 on 21 degrees of freedomNumber of iterations to convergence: 2 Achieved convergence tolerance: 2.571e-09Offerto da

Sebbene i parametri risultino significativi, l'adattamento è scarso. L'errore standard residuo è elevato.

Visualizzare l'adattamento

Il modello lineare sembra comportarsi male; lo si può visualizzare tracciando una linea adattata sui dati e confrontandola con la sua controparte non lineare.

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin$conc, predict(PurbothRealLinear), col=4)Offerto da

Confronto analitico tra modello non lineare e lineare

Dal grafico sopra, sappiamo che il modello lineare fa un lavoro scarso, ma come confrontiamo questi modelli? Sarebbe come confrontare mele con arance. I risultati AIC non possono essere usati direttamente.

Bontà di adattamento

Cosa costituisce un buon adattamento? Per la regressione, che il valore predetto e quello reale siano il più possibile vicini. Quindi, se ispezioniamo la qualità di ciò che un modello ha predetto usando una funzione R lm(), li possiamo confrontare. Il modello più accurato avrà i punti vicini alla retta ideale. Una retta ideale rappresenta un modello accurato al 100%: il predetto è uguale al valore reale per ogni punto dati. Pertanto, tutti i punti giaceranno su una linea retta.

Adattamento non lineare

Per ogni punto dati in input, usiamo il modello per predire il valore target e poi applichiamo la funzione lm() come mostrato. Tracciando uno scatter plot tra la velocità reale (asse x) e la velocità predetta (asse y) e disegnando una retta ideale, si può osservare quanto i punti siano vicini alla linea.

nonLinear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(Purboth))fit_nonLinear = lm(nonLinear_df$rate ~ nonLinear_df$pred_rate)AIC(fit_nonLinear)plot(nonLinear_df$rate, nonLinear_df$pred_rate, xlab="actual_rate",      ylab = "pred_rate",      main = "Comparing Actual vs Predicted for Non Linear Model")lines(Puromycin$rate, Puromycin$rate)174.711198508104Offerto da

Adattamento lineare

Il grafico non lineare sopra è chiaramente un adattamento migliore rispetto al grafico lineare (sotto).

linear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(PurbothRealLinear))fit_linear = lm(linear_df$rate~linear_df$pred_rate)AIC(fit_linear)plot(linear_df$rate, linear_df$pred_rate, xlab="actual_rate", ylab = "pred_rate",      main= "Comparing Actual vs Predicted using lm()")lines(Puromycin$rate, Puromycin$rate)223.783416979551Offerto da

Migliorare il modello lineare

Sembriamo ingiusti col modello lineare perché vediamo chiaramente che i dati non seguono una linea retta, quindi possiamo provare ad adattare il modello usando dati trasformati. Scegliamo sqrt() per trasformare conc e adattiamo un modello lineare come mostrato sotto.

PurbothLinearSqrt = nls(rate ~ i + m*sqrt(conc), Puromycin, list(i=160, m=0.05) )summary(PurbothLinearSqrt)Formula: rate ~ i + m * sqrt(conc)Parameters:  Estimate Std. Error t value Pr(>|t|)    i   60.961      8.758   6.961 7.11e-07 ***m  139.134     15.675   8.876 1.50e-08 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 22.31 on 21 degrees of freedomNumber of iterations to convergence: 2 Achieved convergence tolerance: 1.838e-09Offerto da

In effetti, c'è un calo dell'errore standard residuo, ma inferiore a quanto ci saremmo aspettati.

Tracciare l'adattamento

Si vede chiaramente che usare sqrt() su conc migliora il modello lineare. Tuttavia, rispetto al modello non lineare Purboth non è migliore. L'errore standard residuo del modello PurbothLinearSqrt è quasi il doppio di quello di Purboth, quindi non analizzeremo ulteriormente questo modello. Si può provare ad adattare termini polinomiali di ordine superiore (usando la regressione polinomiale), ma lo sforzo richiesto per mettere a punto quel modello è molto maggiore, mentre con un modello non lineare si possono ottenere risultati migliori.

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state != "treated",]$conc, predict(PurbothLinearSqrt, Puromycin[Puromycin$state != "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc, predict(PurbothLinearSqrt, Puromycin[Puromycin$state != "treated",]), col=3)Offerto da

Rendere lineare l'equazione non lineare

Si può applicare un trucco per "srotolare" una funzione non lineare e adattare un modello lineare. Questa tecnica veniva usata in passato quando i calcoli si facevano a mano, non al computer. Tuttavia, tali modelli non riescono a fare un lavoro migliore.

Convertire la forma non lineare in lineare

Come esempio, si può modificare l'equazione di Michaelis-Menten prendendone l'inversa, come mostrato sotto, e usare la forma modificata per applicare un modello lineare.

Possiamo usare la stessa nls() per adattare la nuova equazione.

PurbothPseudoLinear <- nls((1/rate) ~ K /(Vm*conc) + 1/Vm , Puromycin,                            list(Vm=160, K=0.05))summary(PurbothPseudoLinear)Formula: (1/rate) ~ K/(Vm * conc) + 1/VmParameters:    Estimate Std. Error t value Pr(>|t|)    Vm 1.674e+02  1.362e+01  12.287 4.70e-11 ***K  3.900e-02  6.172e-03   6.318 2.89e-06 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.001786 on 21 degrees of freedomNumber of iterations to convergence: 3 Achieved convergence tolerance: 2.997e-09Offerto da

L'errore standard residuo è basso, ma i dati usati non permettono quindi un confronto con i modelli non lineari. La soluzione migliore è visualizzare come entrambi i modelli si comportano sul dataset originale.

Tracciare l'adattamento

Confrontiamo il modello Purboth_1, che usa due parametri Vmax e K, con PurbothPseudoLinear. Il grafico mostra chiaramente che un modello non lineare fa un lavoro migliore.

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate",      main = "Scatterplot conc vs rate")lines(Puromycin[Puromycin$state == "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state != "treated",]), col=2)lines(Puromycin[Puromycin$state == "treated",]$conc,       1/predict(PurbothPseudoLinear, Puromycin[Puromycin$state == "treated",])      , col=4)lines(Puromycin[Puromycin$state != "treated",]$conc,       1/predict(PurbothPseudoLinear, Puromycin[Puromycin$state != "treated",])      , col=4)legend(0, 200, legend = c("nls", "nls-turned-linear"), fill = c(2,4))Offerto da

Bontà di adattamento

Si può confrontare la bontà del modello.

nonLinear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(Purboth_1))fit_nonLinear = lm(nonLinear_df$rate ~ nonLinear_df$pred_rate)AIC(fit_nonLinear)203.017860461858Offerto da

linear_df = data.frame("rate"=Puromycin$rate,                        "pred_rate"= 1/predict(PurbothPseudoLinear))fit_linear = lm(linear_df$rate~linear_df$pred_rate)AIC(fit_linear)206.999538212925Offerto da

Questo mostra che un modello che usa una funzione non lineare funziona molto meglio.

par( mfrow= c(1,2) )plot(nonLinear_df$rate, nonLinear_df$pred_rate, xlab="actual_rate",      ylab = "pred_rate",      main = "Non-Linear")lines(Puromycin$rate, Puromycin$rate)plot(linear_df$rate, linear_df$pred_rate, xlab="actual_rate", ylab = "pred_rate",      main = "Turned Linear")lines(Puromycin$rate, Puromycin$rate)Offerto da

Conclusione

Il modello non lineare funziona meglio della controparte lineare, soprattutto nei casi in cui abbiamo dati meccanicistici. Tuttavia, un modello non lineare non garantisce una soluzione numerica al problema di stima. Si può tracciare e selezionare una funzione matematica che spieghi al meglio la relazione tra i dati. In tali casi, il modello non lineare può superare la sua controparte lineare.

Se hai trovato illuminante questo articolo e vuoi approfondire il mondo dei modelli di regressione, dai un'occhiata al corso di DataCamp Intermediate Regression in R.