WEEK3-Quick guide to linear regression

Explore Statistics with R (EDX)

WEEK3-Quick guide to linear regression 笔记

一、引子

Imagine you have measured two variables for each subject in a study.
1. Variables are on interval scale or ratio scale
2. Variables are at least roughly normally distributed

You suspect that the variables may be associated. Maybe a linear model would fit:

y = a + bx + error.

Here we will not go into any depth describing the theory of linear regression. Use the example code below to learn how to add a linear regression line to a scatter plot. If you already know something about linear regression you will find some additional useful lines of code below.

(已经测量了两个变量，它们满足：1. 在区间尺或比例尺上，2. 大致正态分布；它们可能是线性关系哦～)

二、例子

#Create some data

set.seed(278)
x <- rnorm(25, mean=100, sd=10)
y <- 2 * x + 20 + rnorm(25, mean=10, sd=4)

plot(x,y)  #Do you think a linear model would fit?

相关系数R和决定系数R^2

> cor(x,y)  #If you just want the correlation coefficient
[1] 0.9548738
> cor(x,y)^2 #Or the coefficient of determination
[1] 0.911784

> lm.obj <- lm(y~x) # See how models are described in R. y depends on x
> abline(lm.obj)    #We can add the regression line to the scatterplot
> predict(lm.obj) #The predicted y-values for your x-values
       1        2        3        4        5        6        7        8 
226.2604 209.0715 239.2173 210.2221 219.6353 211.7828 202.1819 238.3487 
       9       10       11       12       13       14       15       16 
236.4330 215.9328 243.1476 228.5103 214.0388 241.6755 232.5671 206.4833 
      17       18       19       20       21       22       23       24 
216.5090 212.1111 234.9200 234.7542 239.8300 206.8384 248.6047 215.7893 
      25 
221.9042 
> points(x,predict(lm.obj), col="green") #Add predicted values to the graph
> summary(lm.obj)   #Lets look at the content of lm.obj

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.7166 -2.1476 -0.5456  2.2163  9.9858 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  48.9354    11.4055    4.29 0.000273 ***
x             1.7973     0.1166   15.42 1.28e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.375 on 23 degrees of freedom
Multiple R-squared:  0.9118,	Adjusted R-squared:  0.9079 
F-statistic: 237.7 on 1 and 23 DF,  p-value: 1.284e-13

> str(lm.obj)     #You can retrieve parts of the lm object
List of 12
 $ coefficients : Named num [1:2] 48.9 1.8
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
 $ residuals    : Named num [1:25] -1.659 3.734 0.911 -4.294 1.325 ...
  ..- attr(*, "names")= chr [1:25] "1" "2" "3" "4" ...
 $ effects      : Named num [1:25] -1121.35 67.46 1.89 -4.71 1.36 ...
  ..- attr(*, "names")= chr [1:25] "(Intercept)" "x" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:25] 226 209 239 210 220 ...
  ..- attr(*, "names")= chr [1:25] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:25, 1:2] -5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:25] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.2 1.23
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 23
 $ xlevels      : Named list()
 $ call         : language lm(formula = y ~ x)
 $ terms        :Classes 'terms', 'formula' length 3 y ~ x
  .. ..- attr(*, "variables")= language list(y, x)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "y" "x"
  .. .. .. ..$ : chr "x"
  .. ..- attr(*, "term.labels")= chr "x"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(y, x)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
 $ model        :'data.frame':	25 obs. of  2 variables:
  ..$ y: num [1:25] 225 213 240 206 221 ...
  ..$ x: num [1:25] 98.7 89.1 105.9 89.7 95 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 y ~ x
  .. .. ..- attr(*, "variables")= language list(y, x)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "y" "x"
  .. .. .. .. ..$ : chr "x"
  .. .. ..- attr(*, "term.labels")= chr "x"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(y, x)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
 - attr(*, "class")= chr "lm"
> lm.obj$coefficients
(Intercept)           x 
  48.935353    1.797343 

par(mfrow=c(2,2)) #prepare for a 2x2 layout
plot(lm.obj) #The built in controls for your regression analysis
par(mfrow=c(1,1)) #Restore 1x1 layout