Embedded Files
Rappel de la deuxième année de Licence: l'intérêt de la régression linéaire est de savoir si une variable dépend linéairement d'une autre. Par exemple, considérez l'étude réalisée dans l'article "Allelic Variation in Gene Expression Is Common in the Human Genome" par H. Shuen Lo, Zhining Wang, Ying Hu, et al., Genome Res. 2003 13: 1855-1862.
Pour la figure 1(A), les auteurs on effectués des expériences a partir de materiel de type DNA (a gauche) ou (droite) créé à partir de RNA de rein. Chaque petit cercle représente une paire de valeur d'intensité entre deux probes (sondes ?) pour chaque SNP. En abscisse, on représente la valeur de cette différence pour l'expérience 1 et en ordonnée, la valeur de cette différence pour l'expérience 2. On voudrait savoir s'il y a un lien entre l'expérience 1 et l'expérience 2. Un moyen pour le savoir est de chercher si les valeurs pour l'expérience 2 ne pourraient pas être exprimées comme une fonction linéaire des valeurs de l'expérience 1 (plus des fluctuations propres à l'exprérience 2).
Imaginons pour simplifier qu'on ait les valeurs suivantes
> exp1 <- c(10.42 , 7.18 , 7.2 , 5.3 , 5.4)
> exp2 <- c(9.34 , 8.50 , 7.62 , 6.93 , 6.60)
Affichons ces données à l'aide d'un scatter plot et remarquons tout de suite que la relation semble affine (linéaire). Pour confirmer notre découverte empirique, nous calculons le coefficient de corrélation:
> plot(exp1,exp2,
main="exp 2 vs exp 1",
sub="Allelic Variation in Gene Expression Is Common in the Human Genome par H. Shuen Lo,
Zhining Wang, Ying Hu, et al., Genome Res. 2003 13: 1855-1862")
Le coefficient de corrélation est obtenu en faisant:
> cor(exp1,exp2)
[1] -0.9880813
La commande pour déterminer la droite de régression par la méthode des moindres carrés est la commande "lm" (comme "linear model" traduction de "modèle linéaire" comme vous pouvez vous en douter). Cette commande a beaucoup d'options. Si vous êtes intéressés, tapez help(lm) pour en apprendre plus. On considère donc dans la suite le modèle:
exp2= (pente) exp1 + (ordonnée a l'origine)
On procède alors comme suit:
> res=lm(exp2 ~ exp1)
> res
Call:
lm(formula = exp2 ~ exp1)
Coefficients:
(Intercept) exp1
4.177 0.510
Quand vous faites l'appel à la fonction "lm" on se voit retourner beaucoup d'informations. Si vous débutez avec la régression linéaire, vous êtes probablement seulement intéressés par deux choses: la pente et l'intercept (ou ordonnée à l'origine). Si vous tapez juste le nom de la variable rendue par "lm" il imprimera ces informations minimales à l'écran.
Pour savoir tout ce qui est contenu dans le résultat de la fonction "lm" sur vos données, il suffit d'utiliser la commande "attributes":
> attributes(res)
$names
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model"
$class
[1] "lm"
Pour connaître les coefficient seulement:
> res$coefficients[1]
(Intercept)
4.177
> res$coefficients[[1]]
[1] 4.177
> res$coefficients[2]
year
0.510
> res$coefficients[[2]]
[1] 0.510
Pour connaître une estimation de ce que donnerait l'expérience 2 avec une valeur de l'expérience 1 égale à 12 il suffit d'utiliser la ligne de commande suivante:
> res$coefficients[[2]]*12+res$coefficients[[1]]
[1] 10.29720
Calculons maintenant les "résidus", c'est à dire les erreurs commises par ce simple modèle linéaire:
> resid <- exp2 - (res$coefficients[[2]]*exp1+res$coefficients[[1]])
> resid
[1] -0.15133768 0.66119668 -0.22900415 0.05007465 -0.33092950
Cette formule est un peu obscure a première vue mais il y a une fonction qui permet de calculer les résidus automatiquement:
> residuals(res)
1 2 3 4 5
-0.15133768 0.66119668 -0.22900415 0.05007465 -0.33092950
Pour afficher la droite de régression en même temps que les données, il suffit d'utiliser la fonction "abline":
> plot(exp1,exp2)
> abline(res$coefficients)
Finalement, on peut tester la validité du modèle au moyen de la fonction summary:
> summary(res)
Call:
lm(formula = exp2 ~ exp1)
Residuals:
1 2 3 4 5
-0.15134 0.66120 -0.22900 0.05007 -0.33093
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.1767 0.8080 5.169 0.0140 *
exp1 0.5100 0.1101 4.632 0.0190 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4563 on 3 degrees of freedom
Multiple R-squared: 0.8773, Adjusted R-squared: 0.8364
F-statistic: 21.45 on 1 and 3 DF, p-value: 0.01896
---------------------------------------------------------------------------
Notez qu'il est parfois pratique lorsqu'on a beaucoup de facteurs, de construire un dataframe et de l'utiliser dans la regression comme dans l'exemple suivant :
y <- c(125,158,207,182,196,175,145,144,160,175,151,161,200,173,175,162,155,230,162,153)
x1 <- c(13,39,52,29,50,64,11,22,30,51,27,41,51,37,23,43,38,62,28,30)
x2 <- c(18,18,50,43,37,19,27,23,18,11,15,22,52,36,48,15,19,56,30,25)
x3 <- c(25,59,62,50,65,79,17,31,34,58,29,53,75,44,27,65,62,75,36,41)
x4 <- c(11,30,53,29,56,49,14,17,22,40,31,39,36,27,20,36,37,50,20,33)
data <- data.frame(cbind(y,x1,x2,x3,x4))
### 2 ###
reg <-lm(y~x1 + x2 + x3 + x4,data=data)
summary(reg)
Google Sites
Report abuse