R Language
Modelli lineari (regressione)

Sintassi

• lm (formula, dati, sottoinsieme, pesi, na.action, method = "qr", modello = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasti = NULL, offset, .. .)

Parametri

Regressione lineare sul set di dati mtcars

Il frame dati mtcars integrato contiene informazioni su 32 auto, tra cui peso, efficienza del carburante (in miglia al gallone), velocità, ecc. (Per saperne di più sul set di dati, utilizzare help(mtcars) ).

Se siamo interessati alla relazione tra efficienza del carburante ( mpg ) e peso ( wt ) possiamo iniziare a tracciare quelle variabili con:

plot(mpg ~ wt, data = mtcars, col=2)

I grafici mostrano una relazione (lineare) !. Quindi se vogliamo eseguire la regressione lineare per determinare i coefficienti di un modello lineare, useremo la funzione lm :

fit <- lm(mpg ~ wt, data = mtcars)

Il ~ qui significa "spiegato da", quindi la formula mpg ~ wt significa che stiamo predicendo mpg come spiegato da wt. Il modo più utile per visualizzare l'output è con:

summary(fit)

Che dà l'output:

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Questo fornisce informazioni su:

• la pendenza stimata di ciascun coefficiente ( wt e l'intercetta y), che suggerisce la predizione migliore di mpg è 37.2851 + (-5.3445) * wt
• Il valore p di ciascun coefficiente, che suggerisce che l'intercetta e il peso non sono probabilmente dovuti al caso
• Stime complessive di adattamento come R ^ 2 e R ^ 2 aggiustato, che mostrano quanto della variazione di mpg è spiegata dal modello

Potremmo aggiungere una riga al nostro primo grafico per mostrare il mpg previsto:

abline(fit,col=3,lwd=2)

È anche possibile aggiungere l'equazione a quella trama. Innanzitutto, ottieni i coefficienti con coef . Quindi usando paste0 collassiamo i coefficienti con le variabili appropriate e +/- , per costruire l'equazione. Infine, lo aggiungiamo alla trama usando mtext :

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

Il risultato è:

Tracciare la regressione (base)

Continuando sull'esempio di mtcars , ecco un modo semplice per produrre una trama della regressione lineare potenzialmente adatta alla pubblicazione.

Innanzitutto, adattare il modello lineare e

fit <- lm(mpg ~ wt, data = mtcars)

Quindi traccia le due variabili di interesse e aggiungi la linea di regressione all'interno del dominio di definizione:

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

Quasi lì! L'ultimo passaggio consiste nell'aggiungere al grafico, all'equazione di regressione, alla rsquare e al coefficiente di correlazione. Questo viene fatto usando la funzione vector :

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

Si noti che è possibile aggiungere qualsiasi altro parametro come RMSE adattando la funzione vettoriale. Immagina di voler una leggenda con 10 elementi. La definizione del vettore sarebbe la seguente:

rp = vector('expression',10)

e dovrai definire r[1] .... in r[10]

Ecco l'output:

ponderazione

A volte vogliamo che il modello dia più peso ad alcuni punti o esempi di dati rispetto ad altri. Ciò è possibile specificando il peso per i dati di input durante l'apprendimento del modello. Esistono generalmente due tipi di scenari in cui è possibile utilizzare i pesi non uniformi sugli esempi:

• Pesi analitici: riflettono i diversi livelli di precisione delle diverse osservazioni. Ad esempio, se si analizzano i dati in cui ogni osservazione è il risultato medio di un'area geografica, il peso analitico è proporzionale all'inverso della varianza stimata. Utile quando si affrontano le medie dei dati fornendo un peso proporzionale dato il numero di osservazioni. fonte

• Sampling Weights (Inverse Probability Weights - IPW): una tecnica statistica per il calcolo di statistiche standardizzate a una popolazione diversa da quella in cui sono stati raccolti i dati. I disegni di studio con una popolazione di campionamento disparata e la popolazione di inferenza sui target (popolazione target) sono comuni nell'applicazione. Utile quando si tratta di dati con valori mancanti. fonte

La funzione lm() esegue la ponderazione analitica. Per i pesi di campionamento, il pacchetto di survey viene utilizzato per costruire un oggetto di design del sondaggio ed eseguire svyglm() . Per impostazione predefinita, il pacchetto di survey utilizza pesi di campionamento. (NOTA: lm() e svyglm() con gaussian() famiglia gaussian() generano tutte le stesse stime di punto, poiché entrambi risolvono i coefficienti riducendo al minimo i minimi quadrati pesati e differiscono nel modo in cui vengono calcolati gli errori standard).

Dati di test

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

Pesi analitici

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

Produzione

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

Sampling Weights (IPW)

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

Produzione

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

Controllo della non linearità con regressione polinomiale

A volte, quando si lavora con la regressione lineare, è necessario verificare la non linearità nei dati. Un modo per farlo è quello di adattare un modello polinomiale e verificare se si adatta meglio ai dati rispetto a un modello lineare. Ci sono altri motivi, come quelli teorici, che indicano di adattarsi a un modello di ordine quadratico o superiore perché si ritiene che la relazione delle variabili sia intrinsecamente polinomiale in natura.

mtcars un modello quadratico per il set di dati mtcars . Per un modello lineare vedere Regressione lineare sul set di dati mtcars .

Per prima cosa facciamo un diagramma a dispersione delle variabili mpg (Miles / gallon), disp (Displacement (cu.in.)) e wt (Weight (1000 lbs)). La relazione tra mpg e disp appare non lineare.

plot(mtcars[,c("mpg","disp","wt")])

Un adattamento lineare mostrerà che disp non è significativo.

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

Quindi, per ottenere il risultato di un modello quadratico, abbiamo aggiunto I(disp^2) . Il nuovo modello sembra migliore quando si guarda a R^2 e tutte le variabili sono significative.

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

Poiché abbiamo tre variabili, il modello adattato è una superficie rappresentata da:

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

Un altro modo per specificare la regressione polinomiale è l'utilizzo di poly con parametro raw=TRUE , altrimenti verranno considerati i polinomi ortogonali (consultare l' help(ploy) per maggiori informazioni). Otteniamo lo stesso risultato usando:

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

Infine, cosa succede se abbiamo bisogno di mostrare una trama della superficie stimata? Bene ci sono molte opzioni per fare grafici 3D in R Qui usiamo Fit3d dal pacchetto p3d .

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

Valutazione della qualità

Dopo aver creato un modello di regressione, è importante controllare il risultato e decidere se il modello è appropriato e funziona bene con i dati a portata di mano. Questo può essere fatto esaminando la trama dei residui e altri grafici diagnostici.

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

Questi grafici verificano due ipotesi formulate durante la costruzione del modello:

1. Che il valore atteso della variabile prevista (in questo caso mpg ) è dato da una combinazione lineare dei predittori (in questo caso wt ). Ci aspettiamo che questa stima sia imparziale. Quindi i residui dovrebbero essere centrati attorno alla media per tutti i valori dei predittori. In questo caso vediamo che i residui tendono ad essere positivi alle estremità e negativi nel mezzo, suggerendo una relazione non lineare tra le variabili.
2. Che l'effettiva variabile prevista è normalmente distribuita attorno alla sua stima. Pertanto, i residui dovrebbero essere distribuiti normalmente. Per i dati normalmente distribuiti, i punti in un normale diagramma QQ dovrebbero trovarsi sopra o vicino alla diagonale. C'è una certa quantità di inclinazione alle estremità qui.

Utilizzando la funzione 'Prevedi'

Una volta che un modello è stato costruito, predict è la funzione principale da testare con i nuovi dati. Il nostro esempio userà il set di dati incorporato mtcars per regredire miglia per gallone contro lo spostamento:

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Se avessi una nuova fonte di dati con dislocamento, potrei vedere le miglia stimate per gallone.

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

La parte più importante del processo è creare un nuovo frame di dati con gli stessi nomi di colonna dei dati originali. In questo caso, i dati originali avevano una colonna etichettata disp , ero sicuro di chiamare i nuovi dati con lo stesso nome.

Attenzione

Diamo un'occhiata ad alcune trappole comuni:

1. non utilizzare un data.frame nel nuovo oggetto:

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. non usare gli stessi nomi nel nuovo data frame:

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

Precisione

Per verificare l'accuratezza della previsione sono necessari i valori y effettivi dei nuovi dati. In questo esempio, newdf avrà bisogno di una colonna per "mpg" e "disp".

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148