1. No category

Oversigt
Kursus 02402/02323 Introducerende Statistik
Forelæsning 8: Simpel lineær regression
Per Bruun Brockhoff
DTU Compute, Statistik og Dataanalyse
Bygning 324, Rum 220
Danmarks Tekniske Universitet
2800 Lyngby – Danmark
e-mail: perbb@dtu.dk
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
1 / 43
1
Motiverende eksempel: Højde-vægt
2
Lineær regressionsmodel
3
Mindste kvadraters metode (least squares)
4
Statistik og lineær regression??
5
Hypotesetests og konfidensintervaller for β̂0 og β̂1
6
Konfidensinterval og prædiktionsinterval
Konfidensinterval for linien
Prædiktionsinterval
7
Korrelation
8
Residual Analysis: Model control
Per Bruun Brockhoff (perbb@dtu.dk)
Motiverende eksempel: Højde-vægt
Heights (xi )
Weights (yi )
168
65.5
161
58.3
167
68.1
Introduktion til Statistik, Forelæsning 8
2 / 43
Motiverende eksempel: Højde-vægt
179
85.7
184
80.5
166
63.4
198
102.6
187
91.4
191
86.7
179
78.9
Heights (xi )
Weights (yi )
168
65.5
161
58.3
167
68.1
179
85.7
184
80.5
166
63.4
198
102.6
187
91.4
191
86.7
179
78.9
100
7
100
7
9
80
10
5
70
5
70
10
9
4
Weight
80
8
90
90
8
4
Weight
Foråret 2015
3
3
1
6
60
60
6
2
160
170
180
190
2
160
Height
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
1
170
180
190
Height
Foråret 2015
4 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
5 / 43
Motiverende eksempel: Højde-vægt
168
65.5
161
58.3
167
68.1
179
85.7
184
80.5
166
63.4
198
102.6
187
91.4
191
86.7
179
78.9
Heights (xi )
Weights (yi )
summary(lm(y ~ x))
Call:
lm(formula = y ~ x)
161
58.3
167
68.1
179
85.7
184
80.5
166
63.4
198
102.6
187
91.4
191
86.7
179
78.9
100
7
8
90
3Q
2.234
Max
6.477
9
4
Weight
Residuals:
Min
1Q Median
-5.876 -1.451 -0.608
5
10
70
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -119.958
18.897
-6.35 0.00022 ***
x
1.113
0.106
10.50 5.9e-06 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3
6
60
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
168
65.5
80
Heights (xi )
Weights (yi )
Motiverende eksempel: Højde-vægt
1
2
160
170
180
190
Height
Residual standard error: 3.88 on 8 degrees of freedom
Multiple R-squared: 0.932, Adjusted R-squared: 0.924
F-statistic: 110 on 1 and 8 DF, p-value: 5.87e-06
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
6 / 43
Lineær regressionsmodel
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
7 / 43
Lineær regressionsmodel
Et scatter plot af noget data
Opstil en lineær model
Opstil en lineær model
Vi har n par datapunkter (xi , yi )
800
yi = β0 + β1 xi
600
800
●
●
●
400
●
600
●
●
●
●
●
200
400
y
●
●●
●
●
●
●
●
0
●
●
y
●
●●
200
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
−1
0
1
2
−200
−200
0
●
3
●
●
●
x
−1
0
1
2
data punkter
lineaer model
3
x
men den der mangler noget til at beskrive den tilfældige variation!
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
9 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
10 / 43
Lineær regressionsmodel
Lineær regressionsmodel
Model-illustration
800
Opstil en lineær regressionsmodel
β0 + β1x
600
Opstil den lineære regressionsmodel
Y i = β 0 + β 1 x i + εi
200
y
400
Yi er den afhængige variabel (dependent variable). En stokastisk
variabel.
xi er en forklarende variabel (explanatory variable)
σ
0
εi er afvigelsen (error). En stokastisk variabel.
−200
og vi antager
εi er independent and identically distributed (i.i.d.) og N (0, σ 2 )
−1
0
1
2
3
x
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
11 / 43
Mindste kvadraters metode (least squares)
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
12 / 43
Foråret 2015
15 / 43
Mindste kvadraters metode (least squares)
Illustration af model, data og fit
800
Mindste kvadraters metode
600
σ2
God ide: Minimer variansen
på afvigelsen. Det er på næsten alle
måder det bedste valg i dette setup.
●
400
●
●
●
^εi = ei
●
●
●
●
0
n
X
^ ^
β0 + β1x
●
200
Minimer summen af de kvadrerede afvigelser (Residual Sum of Squares
(RSS ))
●
ε2i
●●
●
σ
●
●
●
−200
i=1
●
●
y
But how!?
●
β0 + β2x
Hvad kan vi gøre for at estimere parametrene β0 og β1 ?
RSS (β0 , β1 ) =
Foråret 2015
●
●
βˆ0 og βˆ1 minimerer RSS
−1
0
1
data punkter
lineaer model
lineaer fit
2
3
x
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
14 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Mindste kvadraters metode (least squares)
Mindste kvadraters metode (least squares)
Least squares estimator
Least squares estimater
Theorem 5.4 (her for estimatorer som i eNoten)
Theorem 5.4 (her for estimater)
The least squares estimators of β0 and β1 are given by
Pn
(Yi − Ȳ )(xi − x̄)
β̂1 = i=1
Sxx
β̂0 =Ȳ − β̂1 x̄
P
where Sxx = ni=1 (xi − x̄)2 .
The least squares estimatates of β0 and β1 are given by
Pn
(yi − ȳ)(xi − x̄)
β̂1 = i=1
Sxx
β̂0 =ȳ − β̂1 x̄
P
where Sxx = ni=1 (xi − x̄)2 .
Tænk ikke længere over det for nu!
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
16 / 43
Mindste kvadraters metode (least squares)
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
17 / 43
Statistik og lineær regression??
R eksempel
Parameter estimaterne er stokastiske variabler
## Simuler en lineær model med normalfordelt afvigelse og estimer parametrene
## FØRST LAV DATA:
## Generer n værdier af input x som uniform fordelt
x <- runif(n=20, min=-2, max=4)
Hvis vi tog en ny stikprøve ville estimaterne β̂0 og β̂1 have samme udfald?
## Simuler lineær regressionsmodel
beta0=50; beta1=200; sigma=90
y <- beta0 + beta1 * x + rnorm(n=length(x), mean=0, sd=sigma)
Nej, de er stokastiske variabler. Tog vi en ny stikprøve så ville vi have en
anden realisation.
## HERFRA ligesom virkeligheden, vi har dataen i x og y:
## Et scatter plot af x og y
plot(x, y)
Hvordan er parameter estimaterne i en lineær regressionsmodel fordelt
(givet normalfordelte afvigelser)?
## Udregn least squares estimaterne, brug Theorem 5.4
(beta1hat <- sum( (y-mean(y))*(x-mean(x)) ) / sum( (x-mean(x))^2 ))
(beta0hat <- mean(y) - beta1hat*mean(x))
Prøv lige at simulere for at se på det...
## Brug lm() til at udregne estimaterne
lm(y ~ x)
## Plot den estimerede linie
abline(lm(y ~ x), col="red")
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
18 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
20 / 43
Statistik og lineær regression??
Statistik og lineær regression??
Estimater af standard afvigelserne på β̂0 og β̂1
Hvordan er parameter estimaterne i en lineær regressionsmodel fordelt
(givet normalfordelte afvigelser)?
De er normalfordelte (for n < 30 brug t-fordeling) og deres varians kan
estimeres:
Theorem 5.7 (første del)
x̄2 σ 2
σ2
+
n
Sxx
σ2
V [β̂1 ] =
Sxx
x̄σ 2
Cov[β̂0 , β̂1 ] = −
Sxx
V [β̂0 ] =
Introduktion til Statistik, Forelæsning 8
Where σ 2 is usually replaced by its estimate (σ̂ 2 ). The central estimator for
σ 2 is
Pn 2
e
RSS(β̂0 , β̂1 )
2
σ̂ =
= i=1 i .
n−2
n−2
When the estimate of σ 2 is used the variances also become estimates and
we’ll refer to them as σ̂β20 and σ̂β21 .
Estimat af standard afvigelserne for β̂0 og β̂1 (ligningerne (5-73))
s
s
1
x̄2
1
σ̂β0 = σ̂
+
; σ̂β1 = σ̂ Pn
2
n Sxx
i=1 (xi − x̄)
Kovariansen Cov[β̂0 , β̂1 ] (covariance) gør vi ikke mere ud af her.
Per Bruun Brockhoff (perbb@dtu.dk)
Theorem 5.7 (anden del)
Foråret 2015
21 / 43
Hypotesetests og konfidensintervaller for β̂0 og β̂1
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
22 / 43
Foråret 2015
25 / 43
Hypotesetests og konfidensintervaller for β̂0 og β̂1
Hypotesetests for parameter estimaterne
Se Eksempel 5.12 for eksempel på hypotesetest.
Vi kan altså udføre hypotesetests for parameter estimater i en lineær
regressionsmodel:
Test om parametrene er signifikant forskellige fra 0
H0,i :
βi = β0,i
H0,i :
βi = 0
H1,i :
βi 6= β1,i
H1,i :
βi 6= 0
Se resultatet i R
Vi bruger de t-fordelte statistikker:
## Hypotesetests om signifikante parametre
## Generer x
x <- runif(n=20, min=-2, max=4)
## Simuler Y
beta0=50; beta1=200; sigma=90
y <- beta0 + beta1 * x + rnorm(n=length(x), mean=0, sd=sigma)
Theorem 5.11
Under the null-hypothesis (β0 = β0,0 and β1 = β0,1 ) the statistics
Tβ0
β̂0 − β0,0
=
;
σ̂β0
Tβ1
## Brug lm() til at udregne estimaterne
fit <- lm(y ~ x)
β̂1 − β0,1
=
,
σ̂β1
## Se summary, deri står hvad vi har brug for
summary(fit)
are t-distributed with n − 2 degrees of freedom, and inference should be based on
this distribution.
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
24 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Hypotesetests og konfidensintervaller for β̂0 og β̂1
Hypotesetests og konfidensintervaller for β̂0 og β̂1
Konfidensintervaller for parametrene
## Lav konfidensintervaller for parametrene
## Antal gentagelser
nRepeat <- 100
## Fangede vi den rigtige parameter
TrueValInCI <- logical(nRepeat)
Method 5.14
(1 − α) confidence intervals for β0 and β1 are given by
## Gentag simuleringen og estimeringen nRepeat gange
for(i in 1:nRepeat){
## Generer x
x <- runif(n=20, min=-2, max=4)
## Simuler y
beta0=50; beta1=200; sigma=90
y <- beta0 + beta1 * x + rnorm(n=length(x), mean=0, sd=sigma)
β̂0 ± t1−α/2 σ̂β0
β̂1 ± t1−α/2 σ̂β1
where t1−α/2 is the (1 − α/2)-quantile of a t-distribution with n − 2
degrees of freedom.
## Brug lm() til at udregne estimaterne
fit <- lm(y ~ x)
husk at σ̂β0 og σ̂β1 findes ved ligningerne (5-74)
## Heldigvis kan R beregne konfidensintervallet (level=1-alpha)
(ci <- confint(fit, "(Intercept)", level=0.95))
i R kan σ̂β0 og σ̂β1 aflæses ved "Std. Error"ved "summary(fit)"
}
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
26 / 43
## Var den rigtige parameterværdi "fanget" af intervallet?
(TrueValInCI[i] <- ci[1] < beta0 & beta0 < ci[2])
Per
(perbb@dtu.dk)
Introduktion
Statistik, Forelæsning 8
##Bruun
HvorBrockhoff
ofte blev
den rigtige
værdi til
"fanget"?
Foråret 2015
27 / 43
sum(TrueValInCI) / nRepeat
Konfidensinterval og prædiktionsinterval
Konfidensinterval for linien
Konfidensinterval og prædiktionsinterval
Method 5.17: Konfidensinterval for β0 + β1 x0
Prædiktionsinterval
Method 5.17: Prædiktionsinterval for β0 + β1 x0 + ε0
Konfidensinterval for β0 + β1 x0 svarer til et konfidensinterval for linien
i punktet x0
Prædiktionsintervallet (prediction interval) for Y0 beregnes med en
værdi x0
Dette gøres før Y0 observeres med
Beregnes med
s
s
(βˆ0 + βˆ1 x0 ) ± tα/2 · σ̂
(βˆ0 + β̂1 x0 ) ± tα/2 · σ̂
x̄)2
1 (x0 −
+
n
Sxx
Introduktion til Statistik, Forelæsning 8
Foråret 2015
1 (x0 − x̄)2
+
n
Sxx
Prædiktionsintervallet vil 100(1 − α)% af gangene indeholde den
observerede y0
Konfidensintervallet vil i 100(1 − α)% af gangene indeholde den
rigtige linie, altså β0 + β1 x0
Per Bruun Brockhoff (perbb@dtu.dk)
1+
Et prædiktionsinterval bliver altså større end et konfidensinterval for
fastholdt α
29 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
30 / 43
Konfidensinterval og prædiktionsinterval
Prædiktionsinterval
Konfidensinterval og prædiktionsinterval
Eksempel med konfidensinterval for linien
Prædiktionsinterval
Eksempel med prædiktionsinterval
●
●
800
800
●
## Eksempel med konfidensinterval for linien
## Eksempel med prædiktionsinterval
●
●
−200
●
●
●
●
●
−2
−1
0
1
2
3
4
600
400
y
●
●
●
0
y
●
−400
## Se lige hvad der kom tilbage
head(PI)
●
●
●
200
400
●
●
●
## Plot data, model og intervaller
plot(x, y, pch=20)
abline(fit)
lines(xval, CI[, "lwr"], lty=2, col="red", lwd=2)
lines(xval, CI[, "upr"], lty=2, col="red", lwd=2)
## Beregn interval for hvert x
PI <- predict(fit, newdata=data.frame(x=xval),
interval="prediction",
level=.95)
●
●
●
●
●
●
## Plot data, model og intervaller
plot(x, y, pch=20)
abline(fit)
lines(xval, PI[, "lwr"], lty=2, col="blue", lwd=2)
lines(xval, PI[, "upr"], lty=2, col="blue", lwd=2)
−200
●
●
0
## Se lige hvad der kom
head(CI)
## Lav en sekvens a x værdier
xval <- seq(from=-2, to=6, length.out=100)
●
200
## Brug predict funktionen
CI <- predict(fit, newdata=data.frame(x=xval),
interval="confidence",
level=.95)
●
−400
600
## Lav en sekvens af x værdier
xval <- seq(from=-2, to=6, length.out=100)
●
Introduktion til Statistik, Forelæsning 8
Foråret 2015
31 / 43
●
●
−2
−1
Max
279.1
3
4
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
32 / 43
Testen er H0,i : βi = 0 vs. H1,i : βi 6= 0
Stjernerne er sat efter p-værdien
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)
51.5
31.1
1.66
0.12
x
216.3
15.2
14.22 3.1e-11 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: XXX on XXX degrees of freedom
εi ∼ N (0, σ 2 ) udskrevet er σ̂ og ν frihedsgrader (brug til
hypotesetesten)
Multiple R-squared: XXX
Forklaret varians r2
Resten bruger vi ikke i det her kursus
Residual standard error: 126 on 18 degrees of freedom
Multiple R-squared: 0.918, Adjusted R-squared: 0.914
F-statistic: 202 on 1 and 18 DF, p-value: 3.14e-11
Introduktion til Statistik, Forelæsning 8
2
Residuals:
Min
1Q
Median
3Q
Max:
Residualernes: Minimum, 1. kvartil, Median, 3. kvartil, Maximum
Coefficients:
Estimate Std. Error t value Pr(>|t|) "stjerner"
Koefficienternes:
Estimat
σ̂βi
tobs
p-værdi
Call:
lm(formula = y ~ x)
Per Bruun Brockhoff (perbb@dtu.dk)
1
summary(lm(y∼x)) wrap up
summary(fit)
3Q
86.6
0
Korrelation
Hvad bliver mere skrevet ud af summary?
Residuals:
Min
1Q Median
-184.7 -96.4 -20.3
●
x
Korrelation
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
●
●
x
Per Bruun Brockhoff (perbb@dtu.dk)
●
●
●
●
Foråret 2015
34 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
35 / 43
Korrelation
Korrelation
Forklaret varians og korrelation
Forklaret varians og korrelation
Forklaret varians af en model er r2 , i summary "Multiple R-squared"
Beregnes med
Korrelationen ρ er et mål for lineær sammenhæng mellem to
stokastiske variable
Estimeret (i.e. empirisk) korrelation
√
ρ̂ = r = r2 sgn(β̂1 )
P
(yi − ŷi )2
2
r = 1 − Pi
2
i (yi − ȳ)
hvor sgn(β̂1 ) er: −1 for β̂1 ≤ 0 og 1 for β̂1 > 0
hvor ŷi = β̂0 + β̂1 xi
Altså:
Andel af den totale varians der er forklaret med modellen
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
Positiv korrelation ved positiv hældning
Negativ korrelation ved negativ hældning
36 / 43
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Korrelation
Foråret 2015
37 / 43
Foråret 2015
39 / 43
Korrelation
Test for signifikant korrelation
## Korrelation
## Generer x
x <- runif(n=20, min=-2, max=4)
## Simuler y
beta0=50; beta1=200; sigma=90
y <- beta0 + beta1 * x + rnorm(n=length(x), mean=0, sd=sigma)
Test for signifikant korrelation (lineær sammenhæng) mellem to
variable
## Scatter plot
plot(x,y)
H0 : ρ = 0
## Brug lm() til at udregne estimaterne
fit <- lm(y ~ x)
H1 : ρ 6= 0
## Den rigtige linie
abline(beta0, beta1)
## Plot fittet
abline(fit, col="red")
er ækvivalent med
## Se summary, deri står hvad vi har brug for
summary(fit)
H0 : β1 = 0
## Korrelation mellem x og y
cor(x,y)
H1 : β1 6= 0
hvor β̂1 er estimatet af hældningen i simpel lineær regressionsmodel
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
38 / 43
## Kvadreret er den "Multiple R-squared" fra summary(fit)
cor(x,y)^2
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Residual Analysis: Model control
Residual Analysis: Model control
Residual Analysis
Residual Analysis in R
fit <- lm(y ~ x)
par(mfrow = c(1, 2))
qqnorm(fit$residuals)
plot(fit$fitted, fit$residuals)
Method 5.26
Check normality assumption with qq-plot.
Check (non)systematic behavior by plotting the residuals ei as a
function of fitted values ŷi
150
150
Normal Q−Q Plot
●
●
●
●
●
●
100
●
100
●
●
●
●
●
50
●
●
●
●
−50
● ●
●
0
fit$residuals
0
●●
●
−50
Sample Quantiles
50
●
●
●●
●
●
●
●
●
●
●
●
−150
−150
●
●
−2
−1
0
Theoretical Quantiles
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
Foråret 2015
41 / 43
Foråret 2015
43 / 43
Outline
Outline
1
Motiverende eksempel: Højde-vægt
2
Lineær regressionsmodel
3
Mindste kvadraters metode (least squares)
4
Statistik og lineær regression??
5
Hypotesetests og konfidensintervaller for β̂0 og β̂1
6
Konfidensinterval og prædiktionsinterval
Konfidensinterval for linien
Prædiktionsinterval
7
Korrelation
8
Residual Analysis: Model control
Per Bruun Brockhoff (perbb@dtu.dk)
Introduktion til Statistik, Forelæsning 8
●
●
−100
−100
●
Per Bruun Brockhoff (perbb@dtu.dk)
1
2
●
−200
0
200
400
600
fit$fitted
Introduktion til Statistik, Forelæsning 8
Foråret 2015
42 / 43