Adattamento regressione polinomiale in R

. (È possibile segnalare problema circa il contenuto in questa pagina qui)Vuoi condividere i tuoi contenuti su R-blogger? clicca qui se hai un blog, o qui se non lo fai.

Una relazione lineare tra due variabili x e y è una delle ipotesi più comuni, efficaci e facili da fare quando si cerca di capire la loro relazione. A volte, tuttavia, la vera relazione sottostante è più complessa di quella, e questo è quando la regressione polinomiale arriva per aiutare.

Vediamo un esempio dall’economia: Supponiamo che si desidera acquistare una certa quantità q di un determinato prodotto. Se il prezzo unitario è p, allora si dovrebbe pagare un importo totale y. Questo è un tipico esempio di relazione lineare. Il prezzo totale e la quantità sono direttamente proporzionali. Per tracciarlo scriveremmo qualcosa del genere:

p <- 0.5q <- seq(0,100,1)y <- p*qplot(q,y,type='l',col='red',main='Linear relationship')

La trama sarà simile a questa:
linear-relationship

Ora, questa è una buona approssimazione della vera relazione tra y e q, tuttavia al momento dell’acquisto e della vendita potremmo prendere in considerazione alcune altre informazioni rilevanti, come: Acquistando quantità significative è probabile che possiamo chiedere e ottenere uno sconto, o acquistando sempre più di un certo bene potremmo spingere il prezzo verso l’alto.
Questo può portare a uno scenario come questo in cui il costo totale non è più una funzione lineare della quantità:

y <- 450 + p*(q-10)^3plot(q,y,type='l',col='navy',main='Nonlinear relationship',lwd=3)

Rplot02
Con la regressione polinomiale possiamo adattare modelli di ordine n > 1 ai dati e provare a modellare relazioni non lineari.

Come adattare una regressione polinomiale

Innanzitutto, ricorda sempre use to set.seed(n) quando generi numeri pseudo casuali. In questo modo, il generatore di numeri casuali genera sempre gli stessi numeri.

set.seed(20)

Predittore (q). Usa seq per generare sequenze equidistanti veloce

q <- seq(from=0, to=20, by=0.1)

Valore da prevedere (y):

y <- 500 + 0.4 * (q-10)^3

Un po ‘ di rumore viene generato e aggiunto al segnale reale (y):

noise <- rnorm(length(q), mean=10, sd=80)noisy.y <- y + noise

Trama del segnale rumoroso:

plot(q,noisy.y,col='deepskyblue4',xlab='q',main='Observed data')lines(q,y,col='firebrick1',lwd=3)

Questa è la trama dei nostri dati osservati simulati. I datapoints simulati sono i punti blu mentre la linea rossa è il segnale (segnale è un termine tecnico che viene spesso utilizzato per indicare la tendenza generale che siamo interessati a rilevare).
Rplot03
Il nostro modello dovrebbe essere qualcosa del genere: y = a*q + b*q2 + c*q3 + costo

Adattiamolo usando R. Quando si adattano i polinomi è possibile utilizzare

model <- lm(noisy.y ~ poly(q,3))

O

model <- lm(noisy.y ~ x + I(X^2) + I(X^3))

Tuttavia, si noti che q, I(q^2) e I(q^3) saranno correlate e le variabili correlate possono causare problemi. L’uso di poly() ti consente di evitare questo producendo polinomi ortogonali, quindi userò la prima opzione.

summary(model)Call:lm(formula = noisy.y ~ poly(q, 3))Residuals: Min 1Q Median 3Q Max -212.326 -51.186 4.276 61.485 165.960 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 513.615 5.602 91.69 <2e-16 ***poly(q, 3)1 2075.899 79.422 26.14 <2e-16 ***poly(q, 3)2 -108.004 79.422 -1.36 0.175 poly(q, 3)3 864.025 79.422 10.88 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 79.42 on 197 degrees of freedomMultiple R-squared: 0.8031,Adjusted R-squared: 0.8001 F-statistic: 267.8 on 3 and 197 DF, p-value: 0 

Utilizzando la funzione confint() possiamo ottenere gli intervalli di confidenza dei parametri del nostro modello.
Intervalli di confidenza per i parametri del modello:

confint(model, level=0.95) 2.5 % 97.5 %(Intercept) 502.5676 524.66261poly(q, 3)1 1919.2739 2232.52494poly(q, 3)2 -264.6292 48.62188poly(q, 3)3 707.3999 1020.65097

Trama di fitted vs residui. Nessun modello chiaro dovrebbe mostrare nella trama residua se il modello è una buona misura

plot(fitted(model),residuals(model))

Rplot04
Nel complesso il modello sembra una buona misura come indica la R al quadrato di 0,8. I coefficienti dei termini del primo e del terzo ordine sono statisticamente significativi come ci aspettavamo. Ora possiamo usare la funzione predict() per ottenere i valori montati e gli intervalli di confidenza per tracciare tutto contro i nostri dati.

Valori previsti e intervalli di confidenza:

predicted.intervals <- predict(model,data.frame(x=q),interval='confidence', level=0.99)

Aggiungi linee alla trama esistente:

lines(q,predicted.intervals,col='green',lwd=3)lines(q,predicted.intervals,col='black',lwd=1)lines(q,predicted.intervals,col='black',lwd=1)

Aggiungi una legenda:

legend("bottomright",c("Observ.","Signal","Predicted"), col=c("deepskyblue4","red","green"), lwd=3)

Ecco la trama:
Rplot05
Possiamo vedere che il nostro modello ha fatto un lavoro decente nel montare i dati e quindi possiamo essere soddisfatti.

Una parola di cautela: i polinomi sono strumenti potenti ma potrebbero ritorcersi contro: in questo caso sapevamo che il segnale originale era generato usando un polinomio di terzo grado, tuttavia quando analizziamo i dati reali, di solito ne sappiamo poco e quindi dobbiamo essere cauti perché l’uso di polinomi di ordine elevato (n > 4) può portare a un eccesso di adattamento. L’over-fitting si verifica quando il modello rileva il rumore anziché il segnale: anche se il modello sta migliorando sempre di più nel montare i dati esistenti, questo può essere negativo quando si sta tentando di prevedere nuovi dati e portare a risultati fuorvianti.

Un gist con il codice completo per questo esempio può essere trovato qui.

Grazie per aver letto questo post, lascia un commento qui sotto se hai qualche domanda.

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.