Regressione delle componenti principali in r (passo dopo passo)

Dato un insieme di variabili predittive p e una variabile di risposta, la regressione lineare multipla utilizza un metodo noto come minimo quadrato per ridurre al minimo la somma residua dei quadrati (RSS):

RSS = Σ(y _i – ŷ _i ) ²

Oro:

• Σ : simbolo greco che significa somma
• y _i : il valore di risposta effettivo per l’ ^i-esima osservazione
• ŷ _i : il valore di risposta previsto basato sul modello di regressione lineare multipla

Tuttavia, quando le variabili predittive sono altamente correlate, la multicollinearità può diventare un problema. Ciò può rendere inaffidabili le stime dei coefficienti del modello e presentare una varianza elevata.

Un modo per evitare questo problema è utilizzare la regressione delle componenti principali , che trova M combinazioni lineari (chiamate “componenti principali”) dei predittori p originali e quindi utilizza i minimi quadrati per adattare un modello di regressione lineare utilizzando le componenti principali come predittori.

Questo tutorial fornisce un esempio passo passo di come eseguire la regressione dei componenti principali in R.

Passaggio 1: caricare i pacchetti necessari

Il modo più semplice per eseguire la regressione dei componenti principali in R consiste nell’utilizzare le funzioni nel pacchetto pls .

 #install pls package (if not already installed)
install.packages(" pls ")

load pls package
library(pls)

Passaggio 2: modificare il modello PCR

Per questo esempio, utilizzeremo il set di dati R integrato chiamato mtcars che contiene dati su diversi tipi di auto:

 #view first six rows of mtcars dataset
head(mtcars)

                   mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3,460 20.22 1 0 3 1

Per questo esempio, adatteremo un modello di regressione delle componenti principali (PCR) utilizzando hp come variabile di risposta e le seguenti variabili come variabili predittive:

• mpg
• Schermo
• merda
• peso
• qsec

Il codice seguente mostra come adattare il modello PCR a questi dati. Nota i seguenti argomenti:

• scale=TRUE : indica a R che ciascuna delle variabili predittive deve essere ridimensionata per avere una media pari a 0 e una deviazione standard pari a 1. Ciò garantisce che nessuna variabile predittiva abbia troppa influenza nel modello se viene misurata in unità diverse. .
• validation=”CV” : indica a R di utilizzare la convalida incrociata k-fold per valutare le prestazioni del modello. Tieni presente che per impostazione predefinita vengono utilizzate k=10 pieghe. Tieni inoltre presente che puoi specificare “LOOCV” invece di eseguire la convalida incrociata Leave-One-Out .

 #make this example reproducible
set.seed(1)

#fit PCR model
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=mtcars, scale= TRUE , validation=" CV ")

Passaggio 3: scegli il numero di componenti principali

Una volta adattato il modello, dobbiamo determinare quante componenti principali vale la pena mantenere.

Per fare ciò, basta guardare l’errore quadratico medio della radice del test (test RMSE) calcolato mediante convalida k-cross:

 #view summary of model fitting
summary(model)

Data: 
	Y dimension: 32 1
Fit method: svdpc
Number of components considered: 5

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept) 1 comp 2 comps 3 comps 4 comps 5 comps
CV 69.66 44.56 35.64 35.83 36.23 36.67
adjCV 69.66 44.44 35.27 35.43 35.80 36.20

TRAINING: % variance explained
    1 comp 2 comps 3 comps 4 comps 5 comps
X 69.83 89.35 95.88 98.96 100.00
hp 62.38 81.31 81.96 81.98 82.03

Nel risultato ci sono due tabelle interessanti:

1. CONVALIDA: RMSEP

Questa tabella ci dice il test RMSE calcolato mediante convalida incrociata k-fold. Possiamo vedere quanto segue:

• Se utilizziamo solo il termine originale nel modello, l’RMSE del test è 69,66 .
• Se aggiungiamo la prima componente principale, il test RMSE scende a 44,56.
• Se aggiungiamo la seconda componente principale, il test RMSE scende a 35,64.

Possiamo vedere che l’aggiunta di ulteriori componenti principali si traduce effettivamente in un aumento dell’RMSE del test. Sembra quindi che sarebbe ottimale utilizzare solo due componenti principali nel modello finale.

2. FORMAZIONE: % di varianza spiegata

Questa tabella ci dice la percentuale di varianza nella variabile di risposta spiegata dalle componenti principali. Possiamo vedere quanto segue:

• Utilizzando solo la prima componente principale, possiamo spiegare il 69,83% della variazione della variabile di risposta.
• Aggiungendo la seconda componente principale, possiamo spiegare l’ 89,35% della variazione della variabile di risposta.

Si noti che saremo comunque in grado di spiegare una maggiore varianza utilizzando più componenti principali, ma possiamo vedere che l’aggiunta di più di due componenti principali in realtà non aumenta di molto la percentuale della varianza spiegata.

Possiamo anche visualizzare il test RMSE (insieme al test MSE e R-squared) in funzione del numero di componenti principali utilizzando la funzione validationplot() .

 #visualize cross-validation plots
validationplot(model)
validationplot(model, val.type="MSEP")
validationplot(model, val.type="R2") 

In ogni grafico, possiamo vedere che l’adattamento del modello migliora aggiungendo due componenti principali, ma tende a peggiorare quando aggiungiamo più componenti principali.

Pertanto, il modello ottimale include solo le prime due componenti principali.

Passaggio 4: utilizzare il modello finale per fare previsioni

Possiamo utilizzare il modello PCR finale a due componenti principali per fare previsioni su nuove osservazioni.

Il codice seguente mostra come suddividere il set di dati originale in un set di training e un set di test e utilizzare il modello PCR con due componenti principali per fare previsioni sul set di test.

 #define training and testing sets
train <- mtcars[1:25, c("hp", "mpg", "disp", "drat", "wt", "qsec")]
y_test <- mtcars[26: nrow (mtcars), c("hp")]
test <- mtcars[26: nrow (mtcars), c("mpg", "disp", "drat", "wt", "qsec")]
    
#use model to make predictions on a test set
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=train, scale= TRUE , validation=" CV ")
pcr_pred <- predict(model, test, ncomp= 2 )

#calculate RMSE
sqrt ( mean ((pcr_pred - y_test)^2))

[1] 56.86549

Vediamo che l’RMSE del test risulta essere 56.86549 . Questa è la deviazione media tra il valore hp previsto e il valore hp osservato per le osservazioni del set di test.

L’utilizzo completo del codice R in questo esempio può essere trovato qui .