R Language
Modèles linéaires (régression)

Syntaxe

• lm (formule, données, sous-ensemble, poids, na.action, méthode = "qr", modèle = VRAI, x = FAUX, y = FAUX, qr = VRAI, singular.ok = VRAI, contrastes = NULL, décalage, .. .)

Paramètres

Régression linéaire sur le jeu de données mtcars

Le bloc de données intégré mtcars contient des informations sur 32 voitures, notamment leur poids, leur consommation de carburant (en miles par gallon), leur vitesse, etc. (Pour en savoir plus sur le jeu de données, utilisez help(mtcars) ).

Si nous nous intéressons à la relation entre l'efficacité énergétique ( mpg ) et le poids ( wt ), nous pouvons commencer à tracer ces variables avec:

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

Les graphiques montrent une relation (linéaire)! Alors, si nous voulons effectuer une régression linéaire pour déterminer les coefficients d’un modèle linéaire, nous utiliserons la fonction lm :

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

Le ~ signifie ici "expliqué par", donc la formule mpg ~ wt signifie que nous prédisons mpg comme expliqué par wt. Le moyen le plus utile de visualiser le résultat est le suivant:

summary(fit)

Ce qui donne la sortie:

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

Cela fournit des informations sur:

• la pente estimée de chaque coefficient ( wt et l'ordonnée à l'origine), ce qui suggère que la meilleure prévision de mpg est 37.2851 + (-5.3445) * wt
• La valeur p de chaque coefficient, ce qui suggère que l'interception et le poids ne sont probablement pas dus au hasard
• Les estimations globales de l'ajustement telles que R ^ 2 et R ^ 2 ajusté, qui montrent combien de variation du mpg est expliqué par le modèle

Nous pourrions ajouter une ligne à notre premier tracé pour montrer le mpg prévu:

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

Il est également possible d'ajouter l'équation à cette parcelle. Premièrement, obtenez les coefficients avec coef . Ensuite, en utilisant paste0 nous paste0 les coefficients avec les variables appropriées et +/- , pour construire l'équation. Enfin, nous l’ajoutons au graphique en utilisant 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) 

Le résultat est:

Tracé de la régression (base)

En mtcars exemple de mtcars , voici un moyen simple de produire un tracé de votre régression linéaire potentiellement adapté à la publication.

D'abord adapter le modèle linéaire et

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

Tracez ensuite les deux variables d’intérêt et ajoutez la ligne de régression dans le domaine de définition:

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))))))

Presque là! La dernière étape consiste à ajouter à l'intrigue, l'équation de régression, le carré ainsi que le coefficient de corrélation. Ceci est fait en utilisant la fonction 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')

Notez que vous pouvez ajouter tout autre paramètre tel que le RMSE en adaptant la fonction vectorielle. Imaginez que vous vouliez une légende avec 10 éléments. La définition du vecteur serait la suivante:

rp = vector('expression',10)

et vous devrez définir r[1] .... à r[10]

Voici la sortie:

Pondération

Parfois, nous voulons que le modèle donne plus de poids à certains points de données ou exemples que d’autres. Cela est possible en spécifiant le poids des données d'entrée lors de l'apprentissage du modèle. Il existe généralement deux types de scénarios dans lesquels nous pourrions utiliser des poids non uniformes sur les exemples:

• Poids analytiques: reflètent les différents niveaux de précision des différentes observations. Par exemple, si l'analyse des données pour lesquelles chaque observation correspond à la moyenne des résultats d'une zone géographique, le poids analytique est proportionnel à l'inverse de la variance estimée. Utile pour traiter les moyennes des données en fournissant un poids proportionnel compte tenu du nombre d'observations. La source

• Poids d'échantillonnage (poids de probabilité inverse - IPW): technique statistique de calcul de statistiques normalisées à une population différente de celle dans laquelle les données ont été collectées. Les plans d'étude comportant une population d'échantillonnage disparate et une population d'inférence cible (population cible) sont courants dans les applications. Utile lorsque vous traitez des données qui ont des valeurs manquantes. La source

La fonction lm() effectue la pondération analytique. Pour l'échantillonnage des pondérations, le progiciel d' survey est utilisé pour créer un objet de conception d'enquête et exécuter svyglm() . Par défaut, le package d' survey utilise des poids d'échantillonnage. (NOTE: lm() et svyglm() avec la famille gaussian() produiront toutes les mêmes estimations de points, car elles résolvent toutes les deux les coefficients en minimisant les moindres carrés pondérés. La différence entre les erreurs standards est calculée.)

Données de 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))

Poids analytiques

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

Sortie

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

Poids d'échantillonnage (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)

Sortie

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

Vérification de la non-linéarité avec la régression polynomiale

Parfois, lorsque nous travaillons avec une régression linéaire, nous devons vérifier la non-linéarité des données. L'une des manières d'y parvenir consiste à adapter un modèle polynomial et à vérifier s'il correspond mieux aux données qu'un modèle linéaire. Il y a d'autres raisons, telles que théoriques, qui indiquent l'adéquation d'un modèle quadratique ou d'un ordre supérieur, car on pense que la relation entre les variables est intrinsèquement polynomiale.

mtcars un modèle quadratique au mtcars données mtcars . Pour un modèle linéaire, voir Régression linéaire sur le jeu de données mtcars .

Nous commençons par faire un nuage de points des variables mpg (Miles / gallon), disp (Déplacement (cu.in.)) et wt (Weight (1000 lbs)). La relation entre mpg et disp semble non linéaire.

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

Un ajustement linéaire montrera que disp n'est pas significatif.

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

Ensuite, pour obtenir le résultat d'un modèle quadratique, nous avons ajouté I(disp^2) . Le nouveau modèle semble mieux regarder R^2 et toutes les variables sont significatives.

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 

Comme nous avons trois variables, le modèle ajusté est une surface représentée par:

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

Une autre façon de spécifier la régression polynomiale consiste à utiliser poly avec le paramètre raw=TRUE , sinon les polynômes orthogonaux seront pris en compte (voir l' help(ploy) pour plus d'informations). Nous obtenons le même résultat en utilisant:

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

Enfin, que faire si nous avons besoin de montrer un tracé de la surface estimée? Eh bien, il existe de nombreuses options pour créer 3D parcelles 3D dans R Ici, nous utilisons Fit3d du package p3d .

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

Évaluation de la qualité

Après avoir construit un modèle de régression, il est important de vérifier le résultat et de décider si le modèle est approprié et fonctionne bien avec les données disponibles. Cela peut être fait en examinant la parcelle de résidus ainsi que d'autres parcelles de diagnostic.

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

Ces tracés vérifient deux hypothèses qui ont été faites lors de la construction du modèle:

1. Que la valeur attendue de la variable prédite (dans ce cas mpg ) est donnée par une combinaison linéaire des prédicteurs (dans ce cas, wt ). Nous nous attendons à ce que cette estimation soit impartiale. Les résidus doivent donc être centrés autour de la moyenne de toutes les valeurs des prédicteurs. Dans ce cas, nous voyons que les résidus ont tendance à être positifs aux extrémités et négatifs au milieu, suggérant une relation non linéaire entre les variables.
2. La variable prédite réelle est normalement distribuée autour de son estimation. Ainsi, les résidus doivent être normalement distribués. Pour les données normalement distribuées, les points d'un graphique QQ normal doivent se trouver sur la diagonale ou à proximité de celle-ci. Il y a une certaine quantité de biais aux extrémités ici.

Utiliser la fonction "prédire"

Une fois qu'un modèle est construit, predict est la fonction principale à tester avec de nouvelles données. Notre exemple utilisera le jeu de données intégré mtcars pour régresser les miles par gallon par rapport au déplacement:

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

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

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Si j'avais une nouvelle source de données avec déplacement, je pouvais voir les miles estimés par gallon.

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 partie la plus importante du processus consiste à créer un nouveau bloc de données avec les mêmes noms de colonne que les données d'origine. Dans ce cas, les données d'origine avaient une colonne étiquetée disp , j'étais sûr d'appeler les nouvelles données du même nom.

Mise en garde

Regardons quelques pièges courants:

1. ne pas utiliser un data.frame dans le nouvel objet:

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

2. n'utilisant pas les mêmes noms dans le nouveau bloc de données:

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

Précision

Pour vérifier la précision de la prédiction, vous aurez besoin des valeurs y réelles des nouvelles données. Dans cet exemple, newdf aura besoin d'une colonne pour 'mpg' et '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