R Language
Modèles linéaires (régression)
Recherche…
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
Paramètre | Sens |
---|---|
formule | une formule en notation Wilkinson-Rogers ; response ~ ... où ... contient des termes correspondant à des variables dans l'environnement ou dans le bloc de données spécifié par l'argument de data |
Les données | bloc de données contenant les variables de réponse et de prédicteur |
sous-ensemble | un vecteur spécifiant un sous-ensemble d'observations à utiliser: peut être exprimé comme un énoncé logique en termes de variables dans les data |
poids | poids analytiques (voir la section Poids ci-dessus) |
na.action | comment gérer les valeurs manquantes ( NA ): voir ?na.action |
méthode | comment effectuer le montage. Seuls les choix sont "qr" ou "model.frame" (ce dernier renvoie le cadre du modèle sans ajustement du modèle, identique à la spécification du model=TRUE ) |
modèle | s'il faut stocker le cadre du modèle dans l'objet ajusté |
X | s'il faut stocker la matrice de modèle dans l'objet ajusté |
y | s'il faut stocker la réponse du modèle dans l'objet ajusté |
qr | s'il faut stocker la décomposition QR dans l'objet ajusté |
singular.ok | s'il faut autoriser des ajustements singuliers , des modèles avec des prédicteurs colinéaires (un sous-ensemble des coefficients sera automatiquement défini sur NA dans ce cas |
contrastes | une liste de contrastes à utiliser pour des facteurs particuliers du modèle; voir l'argument contrasts.arg de ?model.matrix.default . Les contrastes peuvent également être définis avec les options() (voir l'argument contrasts ) ou en attribuant les attributs de contrast d'un facteur (voir ?contrasts ) |
décalage | utilisé pour spécifier un composant connu a priori dans le modèle. Peut également être spécifié dans la formule. Voir ?model.offset |
... | arguments supplémentaires à transmettre aux fonctions d'ajustement de niveau inférieur ( lm.fit() ou lm.wfit() ) |
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 est37.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)
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:
- 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. - 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:
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
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