lunedì 17 agosto 2009

Regressione lineare multivariata

Accolgo volentieri l'invito di Fabio, e mi accingo a cominciare alcuni post sulla statistica multivariata.
Nella regressione lineare semplice, abbiamo immaginato che una certa variabile Y dipendesse dall'andamento di un'altra variabile (X), in maniera lineare con andamento crescente o decrescente. Abbiamo quindi visto come realizzare e disegnare la retta che pone in relazione le due variabili, e come valutare la bontà del modello.
Nella realtà (scientifica, economica, psicometrica, etc.), quasi mai un evento dipende solamente dall'andamento di un certo fattore. Tutti gli eventi (anche i più comuni) sono influenzati da numerosissimi elementi. Risulta pertanto molto più utile formulare un modello che tenga conto di tutte queste influenze (o, come vedremo in seguito, delle influenze maggiori). Ciò si ottiene con lo studio della regressione lineare multipla (o multivariata). In generale si indica con Y la variabile dipendente, e con X seguito da un numero in pedice le variabili indipendenti che si suppone abbiano un effetto. Le X vengono chiamate predittori e la formula generale del modello che cerchiamo è:



è l'intercetta, mentre eccetera, sono i regressori, i quali rappresentano il coefficiente angolare della retta che otterremmo al variare del predittore corrispondente, qualora tutti gli altri predittori fossero costanti.
La rappresentazione grafica dipende dal numero di predittori che si vogliono considerare. Con un solo predittore, si ottiene (come abbiamo visto) una retta; con 2 predittori si ottiene un piano nello spazio a 3 dimensioni; con 3 o più predittori non è possibile la rappresentazione grafica, in quanto occorrerebbe rappresentare uno spazio a più di 3 dimensioni.

Cominciamo a vedere i passaggi da effettuare quando si vuole realizzare un modello di regressione multipla. Utilizzerò qui il dataset disponibile su questo sito. In esso viene raccolto il consumo di gelati pro-capite, considerando il periodo dell'anno, il prezzo del gelato, i guadagni delle famiglie e la temperatura media ambientale. Potremmo ad esempio chiederci: il consumo di gelato aumenta o diminuisce quando la temperatura aumenta? Sul consumo di gelati influisce maggiormente la temperatura o il suo prezzo?
Importiamo i dati. Sul sito è disponibile un file in formato .dat. Possiamo aprirlo con il blocco note per visualizzarlo, e renderci conto di come è strutturato. In R è facile importare questo formato di dati, anche senza la necessità di salvare il file sul nostro computer, utilizzando questo comando:


dati <- read.table("http://www.sci.usq.edu.au/staff/dunn/Datasets/Books/Hand/Hand-R/icecream-R.dat", header = TRUE, row.names = 1)


Adesso in dati è contenuta una tabella, con i valori riportati nel file specificato. Perchè R possa accedere al file .dat, non ci deve essere alcuna limitazione di connessione; se non riesce ad accedere ad internet, conviene controllare il proprio firewall (che spesso blocca la connessione a internet da parte di alcuni programmi). Per comodità trasformiamo la tabella in un dataframe:


datiDF <- data.frame(dati)


Disponiamo quindi di 4 variabili e 30 osservazioni. Per prima cosa cerchiamo di indagare il consumo di gelato in base alla temperatura e al prezzo (ignoriamo le entrate). Le colonne che ci interessano, per ora, sono datiDF$Consumption, datiDF$Price, datiDF$Temp. Cerchiamo un modello di regressione multipla a 3 predittori. La tecnica per stimare i regressori è detta ordinary least square OLS, che ripercorre gli stessi principi della tecnica dei minimi quadrati della regressione semplice.
In R stimiamo il modello in questo modo:


model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
summary(model)

Call:
lm(formula = datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

Residuals:
Min 1Q Median 3Q Max
-0.08226 -0.02051 0.00184 0.02272 0.10076

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.59655 0.25831 2.309 0.0288 *
datiDF$Price -1.40176 0.92509 -1.515 0.1413
datiDF$Temp 0.00303 0.00047 6.448 6.56e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04132 on 27 degrees of freedom
Multiple R-squared: 0.6328, Adjusted R-squared: 0.6056
F-statistic: 23.27 on 2 and 27 DF, p-value: 1.336e-06


Il model contiene la formula che vogliamo indagare. Dapprima si specifica la variabile dipendente (il consumo), quindi si specificano le variabili indipendenti, tra loro unite da un più.
Con la funzione summary() otteniamo una serie di informazioni. Innanzitutto il valore dell'intercetta e dei regressori; in questo caso:



Pertanto il modello di regressione multipla stimato è:



Osserviamo il segno del primo regressore, che descrive il predittore "prezzo": esso è negativo, il che significa che per ogni aumento unitario del prezzo del gelato, il consumo procapite diminuisce di circa una unità e mezzo (1.40); viceversa la temperatura ha un effetto positivo (seppure molto scarso) sul consumo di gelato (all'aumentare della temperatura, aumenta anche il consumo).

Possiamo rappresentare graficamente il modello creato. Installiamo il pacchetto scatterplot3d, e utilizziamo l'unica funzione in esso compreso, specificando le variabili nell'ordine x, y, z (nel nostro caso meglio specificare (X1, X2, Y)):


library(scatterplot3d)
sc <- scatterplot3d(datiDF$Price, datiDF$Temp, datiDF$Consumption, pch=16)




Possiamo anche intersecare il piano che abbiamo trovato, in questo modo:


model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
sc$plane(model, col="red")




Analizziamo adesso il modello. Il coefficiente di determinazione (Multiple R-squared) consente di valutare la bontà del modello, ossia la proporzione di variabilità della Y spiegata dalle variabili esplicative considerate. In questo caso R-squared è pari a 0.63 (il 63% dei dati viene spiegato dalle variabili esplicative). Il valore R-squared adjusted è una correzione del coefficiente di determinazione che tiene conto del numero di predittori utilizzati (secondo alcuni statistici è più preciso quest'ultimo valore per definire la bontà del modello o goodness-of-fit).
Inoltre è opportuno effettuare una verifica dell'intercetta e dei regressori, ossia un t-test con ipotesi nulla H0: beta = 0. Ossia si deve verificare se il valore campionario possa essere esteso all'intera popolazione. Nel nostro caso i p-value associati alle variabili sono:
p-value = 0.0288 (significativo)
p-value = 0.1413 (non significativo - dovremmo ignorarla)
p-value = 6.56e-07 (significativo)

La statistica F ci indica se il modello è da scartare nella sua interessa, oppure se può essere ritenuto valido. Essendo F-statistic maggiore dell'F-tabulato, rifiutiamo l'ipotesi nulla che il modello sia da scartare nella sua interezza.

Possiamo calcolare in questo modo gli intervalli di confidenza dell'intercetta e dei 2 regressori:


confint(model)
2.5 % 97.5 %
(Intercept) 0.066549244 1.126551076
datiDF$Price -3.299896373 0.496375437
datiDF$Temp 0.002066035 0.003994569


Se vogliamo invece calcolare l'intervallo di confidenza di un valore predetto (ossia se vogliamo prevedere il valore che assumerà Y, dati i predittori), utilizziamo la funzione predict():


model <- lm(Consumption ~ Price + Temp, data=datiDF)
predict(model, data.frame(Price=0.280, Temp=30), interval="confidence")

fit lwr upr
1 0.2949663 0.2700105 0.3199220





Proviamo a complicare un pò le cose, e supponiamo di considerare nel nostro modello anche la variabile "entrate". All'aumentare dei guadagni, aumenta anche il consumo di gelato?


model <- lm(Consumption ~ Price + Temp + Income, data=datiDF)
summary(model)

Call:
lm(formula = Consumption ~ Price + Temp + Income, data = datiDF)

Residuals:
Min 1Q Median 3Q Max
-0.065302 -0.011873 0.002737 0.015953 0.078986

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1973151 0.2702162 0.730 0.47179
Price -1.0444140 0.8343573 -1.252 0.22180
Temp 0.0034584 0.0004455 7.762 3.1e-08 ***
Income 0.0033078 0.0011714 2.824 0.00899 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03683 on 26 degrees of freedom
Multiple R-squared: 0.719, Adjusted R-squared: 0.6866
F-statistic: 22.17 on 3 and 26 DF, p-value: 2.451e-07


Il regressore beta3 associato al predittore entrate risulta essere significativo, e di segno positivo. Quindi all'aumentare delle entrate, aumenta anche il consumo di gelati.
Il modello a 3 predittori viene complessivamente accettato (la statistica F risulta significativa), e si assiste ad un miglioramento della goodness-of-fit (è aumentato il valore dell'R quadro aggiustato).

Una domanda che potremmo porci è la seguente: tra le tre variabili considerate, qual è quella che maggiormente influisce sulla variabile consumo? Per far questo possiamo confrontare i valori dei regressori beta, ma solo dopo averli standardizzati (altrimenti non è possibile il confronto). Nel pacchetto QuantPsyc è presenta la funzione lm.beta() che fa al caso nostro:


library(QuantPsyc)
lm.beta(model)

Price Temp Income
-0.1324351 0.8632558 0.3140085


Quelli qui riportati sono i coefficienti beta standardizzati. Osserviamo che la variabile esplicativa che maggiormente influisce sul consumo di gelati è la temperatura.
Se non volete scaricare il pacchetto per utilizzare solo questa funzione, è sufficiente sapere che il calcolo di coefficienti beta standardizzati è pari al prodotto del valore beta per il rapporto tra la deviazione standard della X considerata, e la deviazione standard della Y:


coef(model)[2]*(sd(datiDF$Price))/(sd(datiDF$Consumption))
Price
-0.1324351

coef(model)[3]*(sd(datiDF$Temp))/(sd(datiDF$Consumption))
Temp
0.8632558

coef(model)[4]*(sd(datiDF$Income))/(sd(datiDF$Consumption))
Income
0.3140085


In maniera analoga (ma anche un pò più rischiosa), si potrebbe valutare la dipendenza anche osservando la correlazione tra le variabili. Se calcoliamo il coefficiente di correlazione (con la funzione cor())tra il consumo e le variabili prezzo, entrate, temperatura, osserviamo che il maggior valore è relativo alla variabile temperatura (è maggiore la correlazione tra queste due variabili).

Spesso si finisce col cercare una quantità di possibili variabili esplicative, da rendere poco agevole lo sviluppo di un modello di regressione multipla, che teoricamente dovrebbe contenere il minor numero di predittori possibile. La scelta non si basa solo sulla valutazione "ad occhio" del peso di ciascun predittore (come abbiamo fatto prima), ma viene utilizzato il test F parziale, che in pratica va a testare se c'è un incremento dell'R-quadro nel momento in cui si aggiunge un'altra variabile esplicativa al modello.
In R per eseguire un test F-parziale, si confronta la tabella dell'analisi della varianza dei due modelli. Supponiamo di voler confrontare le bontà dei modelli (Price + Temp) e (Price + Temp + Income):


model1 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
model2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income)
anova(model1, model2)

Analysis of Variance Table

Model 1: datiDF$Consumption ~ datiDF$Price + datiDF$Temp
Model 2: datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27 0.046090
2 26 0.035273 1 0.010817 7.9734 0.008989 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Se rifiutiamo l'ipotesi H0, vuol dire che l'inserimento della variabile esplicativa ha apportato un miglioramento al modello. In questo caso F = 7.97, con p-value < 0.05. Rifiutiamo quindi l'ipotesi nulla: l'aggiunta del predittore entrate migliora il modello.

Possiamo infine calcolare i coefficiente di determinazioni parziali. Fissata una delle variabili esplicative, a quanto ammonta la bontà del modello considerando l'altra variabile?
Riprendiamo a considerare solo 2 predittori (Temp e Price). Costruiamo prima i due modelli lineari semplici, quindi il modello a 2 predittori. Tenendo presente le formule per calcolare il coefficiente di determinazione parziale, abbiamo:


lm.x1 <- lm(datiDF$Consumption ~ datiDF$Price)
lm.x2 <- lm(datiDF$Consumption ~ datiDF$Temp)
lm.x1.x2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

ssr.x1.x2 <- sum(anova(lm.x1.x2)$"Sum Sq"[1:2])
ssr.x1 <- anova(lm.x1)$"Sum Sq"[1]
ssr.x2 <- anova(lm.x2)$"Sum Sq"[1]
ssr.x1.x2 - ssr.x1
sse.x1 <- anova(lm.x1)$"Sum Sq"[2]
sse.x2 <- anova(lm.x2)$"Sum Sq"[2]

(ssr.x1.x2 - ssr.x2) / sse.x2
[1] 0.07837315

(ssr.x1.x2 - ssr.x1) / sse.x1
[1] 0.6062858


Il coefficiente di determinazione parziale corrispondente alla variabile X1, quando X2 è tenuta costante, ci dice che per una data temperatura, il 7.83% della variabilità del consumo di gelati può essere spiegata dalla variabilità del prezzo nei negozi. Il coefficiente di determinazione parziale corrispondente alla variabile X2, quando X1 è tenuta costante, ci dice che per un dato prezzo, il 60.63% della variabilità del consumo di gelati può essere spiegata dalla variabilità della temperatura (sono valori in accordo con il peso definito con il calcolo dei regressori standardizzati).


Queste viste finora solo alcune delle analisi preliminari che dobbiamo effettuare quando si cerca un modello di regressione lineare multipla. Nei post successivi vedremo altri esempi, e andremo ad effettuare altre analisi, come l'analisi dei residui, la scelta dei predittori e il confronto tra più modelli.

17 commenti:

  1. Ottimo articolo. Fabio

    RispondiElimina
  2. Molto molto molto utile....
    Esaustivo!!
    :D

    RispondiElimina
  3. grande articolo, complimentoni

    RispondiElimina
  4. grazie mille..molto utile e chiaro :)

    RispondiElimina
  5. Molto utili gli articoli. Ho un problema. Se i dati di cui voglio fare la regressione sono ripetuti come faccio? Per es.:
    -x = a, a, a, b, b,b, c, c,c, d, d, d, (quindi 4 osservazioni ripetute 3 volte per la variabile indipendente);
    -y = 2, 3, 1, 2, 2,2, 6,7,5, 8, 9, 6 (variabile dipendente)

    Grazie in anticipo

    RispondiElimina
  6. Salve a tutti! Il mio problema è che una volta caricati i dati, ho la colonna relativa ai nomi delle osservazioni che non so come eliminare per procedere con l'analisi.Purtroppo sto studiando da autodidatta...e ho un pò di difficoltà. Vi ringrazio per l'aiuto.

    RispondiElimina
  7. complimenti, un articolo molto ampio!

    RispondiElimina
  8. grazie, forse grazie a te supererò l'esame!

    RispondiElimina
  9. che cos'è F tabulato per favore ?

    ossia avendo il valore F-static come faccio a capire quale sia la bontà del modello? devo usare il p-value finale ?

    RispondiElimina
  10. Molte molte molte grazie, per me è stato importantissimo!

    RispondiElimina
  11. per accettare il modello il t value deve essere tale da essere in modulo >2 ?o sbaglio?

    RispondiElimina
  12. Il link alla risorsa icecream-R.dat ritorna l'errore Forbidden.

    Ho trovato una risorsa alternativa con gli stessi dati: http://www.math.univ-angers.fr/~loustau/icecream-R.dat

    RispondiElimina
  13. io non riesco a fare apparire il piano di regressione.. continua a non uscire!! cosa sbaglio?
    > sc=scatterplot3d(incstrad$n_mort, incstrad$n_veic_contr, incstrad$X._pop_alc, pch=16)
    > model=lm(n_mort~n_veic_contr+X._pop_alc,data=incstrad)
    > summary(model)

    Call:
    lm(formula = n_mort ~ n_veic_contr + X._pop_alc, data = incstrad)

    Residuals:
    Min 1Q Median 3Q Max
    -12.191 -3.711 1.621 3.056 12.779

    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) -494.3771 112.2055 -4.406 0.00313 **
    n_veic_contr 0.4951 0.1241 3.989 0.00526 **
    X._pop_alc 6.1764 1.8095 3.413 0.01124 *
    ---
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Residual standard error: 7.734 on 7 degrees of freedom
    Multiple R-squared: 0.8877, Adjusted R-squared: 0.8556
    F-statistic: 27.66 on 2 and 7 DF, p-value: 0.0004751

    > summary(sc)
    Length Class Mode
    xyz.convert 1 -none- function
    points3d 1 -none- function
    plane3d 1 -none- function
    box3d 1 -none- function
    >
    > sc$plane3d(model, col="red")

    RispondiElimina
  14. Regressione lineare multivariata e multipla sono due cose diverse! La prima individua regressioni effettuate su una variabile risposta Y che non è più un vettore ma un insieme di vettori (ad es., si vuol vedere l'andamento di temperatura e inquinamento in una certa zona), la seconda tratta semplicemente regressioni (univariate o multivariate) con più di una variabile esplicativa. Watch out!

    RispondiElimina
  15. Articolo utilissimo!una domanda: nel mio caso devo valutare se il numero di pescatori stranieri è variato in modo significativo in 3 diversi porti e nel tempo. Considero i pescatori stranieri come variabile dipendente e il tempo come variabile indipendente, come procedo in R?

    RispondiElimina
  16. Nonostante abbia sentito altre volte questo tipo di equivoco, credo che nella letteratura specialistica il termine "multivariato" sia concordemente riferito alle statistiche con più di una variabile dipendente (MANOVA, Correlazione canonica, Regressione multivariata...).

    Il termine "multipla" riferito alla regressione è invece usato per indicare modelli con più di una variabile indipendente.

    Anch'io, come l'utente mbw, per chiarezza non userei i due termini in maniera indiscriminata e mi adeguerei invece alla terminologia più diffusa.

    RispondiElimina
  17. mi dite come si calcola P value?

    RispondiElimina