Statistick´ e metody v marketingu Ing. Michael Rost, Ph.D. ˇ esk´ Jihoˇ cesk´ a univerzita v C ych Budˇ ejovic´ıch
Regresn´ı anal´ yza doplnˇ en´ı z´ aklad˚ u Vzhledem k poˇ zadavku Vaˇ sich koleg˚ u zaˇ razuji doplˇ nuj´ıc´ı partii o regresn´ı anal´ yze.
c Rost 2007 °
Motivaˇ cn´ı pˇ r´ıklad: V souboru Engel.xls m´ ate k dispozici ´ udaje o roˇ cn´ım disponibiln´ım pˇ r´ıjmu a roˇ cn´ıch v´ ydaj´ıch za j´ıdlo, kter´ e byly ´ daje jsou uvedeny v belgick´ zaznamen´ any u 235 rodin. U ych franc´ıch. Vytvoˇ rte korelaˇ cn´ı pole. Prostˇ rednictv´ım funkce yˆ = βˆ0 + βˆ1x popiˇ ste vztah mezi pˇ r´ıjmem a v´ ydaji na j´ıdlo. Jak lze interpretovat odhadnut´ y regresn´ı koeficient u vysvˇ etluj´ıc´ı promˇ enn´ e? Jak jej naz´ yvaj´ı ekonomov´ e? Je to vhodn´ y model? S
c Rost 2007 °
Pˇ r´ıprava v R
belgie<-read.table("P:/KURZ/Engel.csv",header=TRUE,dec=",",sep=";") belgie[1:5,] 1 2 3 4 5
prijem vydajezajidlo 420.1577 255.8394 541.4117 310.9587 901.1575 485.6800 639.0802 402.9974 750.8756 495.5608
par(mfrow=c(1,2)) plot(prijem,vydajezajidlo,col="blue",pch=20,xlab="Prijem", ylab="Vydaje za jidlo") obal<-chull(prijem,vydajezajidlo) belgie1<-belgie[-obal,] plot(belgie1,col="blue",pch=20,xlab="Prijem",ylab="Vydaje za jidlo")
c Rost 2007 °
1000 3000
Prijem 5000 500 1000 1500 2000 2500 Prijem
c Rost 2007 °
200
400
500
600
800
1000
Vydaje za jidlo
1000
Vydaje za jidlo
1200
1500
1400
1600
2000
Dva ´ ukoly: V pr˚ ubˇ ehu zkoum´ an´ı korelaˇ cn´ı z´ avislosti hled´ ame odpovˇ ed’ na dvˇ e ot´ azky:
• Jak nejl´ epe vystihnout pr˚ ubˇ eh z´ avislosti mezi sledovan´ ymi znaky prostˇ rednictv´ım odpov´ıdaj´ıc´ı matematick´ e funkce? To ˇ reˇ s´ı regresn´ı anal´ yza.
c Rost 2007 °
Dva ´ ukoly: V pr˚ ubˇ ehu zkoum´ an´ı korelaˇ cn´ı z´ avislosti hled´ ame odpovˇ ed’ na dvˇ e ot´ azky:
• Jak nejl´ epe vystihnout pr˚ ubˇ eh z´ avislosti mezi sledovan´ ymi znaky prostˇ rednictv´ım odpov´ıdaj´ıc´ı matematick´ e funkce? To ˇ reˇ s´ı regresn´ı anal´ yza.
• Jak´ y je stupeˇ n (tˇ esnosti, intenzity, s´ıly) z´ avislosti mezi sledovan´ ymi znaky? Odpovˇ ed’ d´ av´ a korelaˇ cn´ı anal´ yza.
c Rost 2007 °
Pˇ redpoklady regresn´ıho modelu
• Stˇ redn´ı hodnota rezidu´ı je nulov´ a. Nebo-li E(²i) = 0 • Rozptyl rezidu´ı je konstantn´ı pro vˇ sechny pozorov´ an´ı, tedy V ar(²i) = σ 2 • Rezidua sleduj´ı norm´ aln´ı rozdˇ elen´ı ²i ∼ N (0, σ 2) • Jednotliv´ e pozorov´ an´ı z´ avisl´ e promˇ enn´ e yi jsou navz´ ajem nez´ avisl´ e. V d˚ usledku toho pak i jednotliv´ e ²i • Jednotliv´ e´ urovnˇ e - hodnoty regresor˚ u jsou pevn´ e, pokud jsou n´ ahodn´ e, pak jsou navz´ ajem nez´ avisl´ e. c Rost 2007 °
Model V pˇ r´ıpadˇ e jednoduch´ e line´ arn´ı regrese vych´ az´ıme z pˇ redpokladu, ˇ ze lze i-t´ e pozorov´ an´ı, i = 1, 2, · · · , n, n ≥ 3 z´ avisle promˇ enn´ e Y, vyj´ adˇ rit prostˇ rednictv´ım nez´ avisle promˇ enn´ e X. Konkr´ etnˇ e jako: yi = β0 + β1xi1 + ²i , a tedy n rovnic pro n pozorov´ an´ı:
c Rost 2007 °
Model V pˇ r´ıpadˇ e jednoduch´ e line´ arn´ı regrese vych´ az´ıme z pˇ redpokladu, ˇ ze lze i-t´ e pozorov´ an´ı, i = 1, 2, · · · , n, n ≥ 3 z´ avisle promˇ enn´ e Y, vyj´ adˇ rit prostˇ rednictv´ım nez´ avisle promˇ enn´ e X. Konkr´ etnˇ e jako: yi = β0 + β1xi1 + ²i , a tedy n rovnic pro n pozorov´ an´ı:
y1 = β0 + β1x11 + ²1 , y2 = β0 + β1x21 + ²2 , ... yn = β0 + β1xn1 + ²n , c Rost 2007 °
Abychom nemuseli vypisovat vˇ sech n rovnic, vyuˇ zijme maticov´ e symboliky:
y
=
y1 y2 ... yn
=
X
1 x11 1 x21 ... ... 1 xn1
"
β=
β0 β1
#
²
=
²1 ²2 ... ²n
Situaci pak m˚ uˇ zeme elegantnˇ e zachytit takto y = Xβ + ². Ot´ azkou je, jak zvolit hodnoty β, tak aby regresn´ı funkce co nejl´ epe vystihovala analyzovan´ a data?
c Rost 2007 °
Metoda nejmenˇ s´ıch ˇ ctverc˚ u Odhady regresn´ıch parametr˚ u β prov´ ad´ıme pomoc´ı metody ˇ. nejmenˇ s´ıch ˇ ctverc˚ u - MNC
Jej´ı podstatou je minimalizace souˇ ctu ˇ ctverc˚ u rezidu´ı.
c Rost 2007 °
Metoda nejmenˇ s´ıch ˇ ctverc˚ u Odhady regresn´ıch parametr˚ u β prov´ ad´ıme pomoc´ı metody ˇ. nejmenˇ s´ıch ˇ ctverc˚ u - MNC
Jej´ı podstatou je minimalizace souˇ ctu ˇ ctverc˚ u rezidu´ı. Zˇ rejmˇ e lze rezidua definovat jako ² = y − Xβ
c Rost 2007 °
Podstata metody
Rezidua: (yi − yˆi)
i = 1, 2, . . . , n .
c Rost 2007 °
Podstata metody
Rezidua: (yi − yˆi)
i = 1, 2, . . . , n .
ˇ tverce rezidu´ı: (yi − yˆi)2 C
i = 1, 2, . . . , n .
c Rost 2007 °
Podstata metody
Rezidua: (yi − yˆi)
i = 1, 2, . . . , n .
ˇ tverce rezidu´ı: (yi − yˆi)2 C
i = 1, 2, . . . , n .
Pn Souˇ cet ˇ ctverc˚ u rezidu´ı: i=1(yi − yˆi)2 .
c Rost 2007 °
Podstata metody
Rezidua: (yi − yˆi)
i = 1, 2, . . . , n .
ˇ tverce rezidu´ı: (yi − yˆi)2 C
i = 1, 2, . . . , n .
Pn Souˇ cet ˇ ctverc˚ u rezidu´ı: i=1(yi − yˆi)2 .
Minimalizace souˇ ctu ˇ ctverc˚ u rezidu´ı:
Pn 2 → min (y − y ˆ ) i i i=1
c Rost 2007 °
Podstata metody
Elegantnˇ e pomoc´ı maticov´ eho z´ apisu S=
n X
²i²i = ²t² → min
i=1
Pokud hodl´ ame minimalizovat funkci S, pak je nutno funkci derivovat a takto derivovanou funkci poloˇ zit rovno nule.
c Rost 2007 °
ˇ Geometrick´ a interpretace MNC Geometricka interpretace metody nejmensich ctvercu
4.0
3.5
y
3.0 ●
●
2.5
● ●
●
●
●
● ●
●
2.0
1.5 −1.0
−0.5
0.0
0.5
1.0
1.5
x
c Rost 2007 °
T´ım je splnˇ en nutn´ y pˇ redpoklad. Odhad jednotliv´ ych sloˇ zek vektoru β tj. regresn´ıch koeficient˚ u z´ısk´ ame takto: S = ²t² = (y − Xβ )t(y − Xβ ) = = yty − ytXβ − (Xβ )ty + (Xβ )tXβ = = yty − 2(Xβ )ty + β tXtXβ Derivaci funkce S poloˇ z´ıme rovnou nule a vyˇ reˇ s´ıme (to je nutn´ a podm´ınka): ∂S = −2Xty + 2XtXβ = 0 ∂β
c Rost 2007 °
Odhad ˆ β Lze tedy ps´ at 2XtXβ = 2Xty
|(XtX)−I
Z´ısk´ ame tak odhad vektoru regresn´ıch koeficient˚ u
c Rost 2007 °
Odhad ˆ β Lze tedy ps´ at 2XtXβ = 2Xty
|(XtX)−I
Z´ısk´ ame tak odhad vektoru regresn´ıch koeficient˚ u
ˆ = (XtX)−I Xty . β
c Rost 2007 °
. . . pokraˇ cov´ an´ı pˇ r´ıkladu attach(belgie) model<-lm(vydajezajidlo~prijem,belgie) summary(model) Residuals: Min 1Q -725.699 -60.239
Median -4.317
3Q 53.411
Max 515.772
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 147.47539 15.95708 9.242 <2e-16 *** prijem 0.48518 0.01437 33.772 <2e-16 *** --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Residual standard error: 114.1 on 233 degrees of freedom Multiple R-Squared: 0.8304, Adjusted R-squared: 0.8296 F-statistic: 1141 on 1 and 233 DF, p-value: < 2.2e-16 plot(prijem,vydajezajidlo,col="blue",pch=20,xlab="Prijem", ylab="Vydaje za jidlo") abline(model,lwd=2,col="red") c Rost 2007 °
Test vˇ sech prediktor˚ u - vysvˇ etluj´ıc´ıch promˇ enn´ ych Jsou vysvˇ etluj´ıc´ı promˇ enn´ e uˇ ziteˇ cn´ e pro predikci z´ avisle promˇ enn´ e? Form´ alnˇ e testujeme hypot´ ezu: H0 : β1 = β2 = . . . = βp−1 = 0 Testovou statistikou je F =
(T SS − RSS)/(p − 1) RSS/(n − p)
F ∼ Fp−1,n−p
Kde RSS a T SS:
ˆ t(y − Xβ) ˆ RSS = (y − Xβ) T SS = (y − y¯)t(y − y¯) Vysok´ e hodnoty F vedou k zam´ıtnut´ı testovan´ e hypot´ ezy. c Rost 2007 °
. . . pokraˇ cov´ an´ı pˇ r´ıkladu
summary(model) Residuals: Min 1Q -725.699 -60.239
Median -4.317
3Q 53.411
Max 515.772
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 147.47539 15.95708 9.242 <2e-16 *** prijem 0.48518 0.01437 33.772 <2e-16 *** --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Residual standard error: 114.1 on 233 degrees of freedom Multiple R-Squared: 0.8304, Adjusted R-squared: 0.8296 F-statistic: 1141 on 1 and 233 DF, p-value: < 2.2e-16
c Rost 2007 °
Wald˚ uv test pro vysvˇ etluj´ıc´ı promˇ ennou Lze j´ım odpovˇ edˇ et na ot´ azku, zda je moˇ zn´ e vyˇ radit pˇ r´ısluˇ snou vysvˇ etluj´ıc´ı promˇ enou z regresn´ıho modelu. Form´ alnˇ e tedy umoˇ zˇ nuje testovat hypot´ ezu H0 : βi = 0.
βˆi ti = s.e.(βˆi)
ti ∼ tn−p
Malou modifikac´ı Waldova testu m˚ uˇ zeme otestovat hypot´ ezu kterou lze vyj´ adˇ rit jako H0 : βi = konst. Testov´ e krit´ erium m´ a pak n´ asleduj´ıc´ı tvar: βˆi − konst. ti = s.e.(βˆi)
ti ∼ tn−p
c Rost 2007 °
Konfidenˇ cn´ı intervaly pro β Je dobr´ e si uvˇ edomit, ˇ ze CI (angl. confidence interval) umoˇ zˇ nuj´ı alternativnˇ e vyj´ adˇ rit nejistotu naˇ sich odhad˚ u! Obecn´ a forma konfidenˇ cn´ıch interval˚ u pro odhady regresn´ıch koeficient˚ u: Odhad ± kritick´ a hodnota × SE odhadu V pˇ r´ıpadˇ e klasick´ eho line´ arn´ıho modelu z´ısk´ ame intervalov´ y odhad pro regresn´ı koeficient βi jako q
βˆi ± t1−α/2,n−pσ ˆ (XtX)−I ii
c Rost 2007 °
Lze sestrojit i simult´ ann´ı konfidenˇ cn´ı oblast pro v´ıce regresn´ıch koefeicent˚ u. Tento pˇ r´ıstup, pokud je umoˇ znˇ en statistick´ ym softwarem, je pochopitelnˇ e preferov´ an. Viz grafick´ e zn´ azornˇ en´ı konfidenˇ cn´ı elipsy. Oblast, resp. 100(1 − α)% konfidenˇ cn´ı region lze vyj´ adˇ rit takto: (βˆ − β)tXtX(βˆ − β) ≤ pˆ σ 2F1−α,p,n−p
c Rost 2007 °
Konfidenˇ cn´ı intervaly pro regresn´ı
0.50 0.45 0.40
β1
0.55
0.60
koeficienty βi
100
120
140
160 β0
180
200
c Rost 2007 °
. . . pokraˇ cov´ an´ı pr´ ace v R
confint(model) 2.5 % 97.5 % (Intercept) 116.0367905 178.913984 prijem 0.4568738 0.513483
c Rost 2007 °
Konfidenˇ cn´ı intervaly pro predikci V podstatˇ e je nutn´ e rozliˇ sit dva v´ yznamovˇ e odliˇ sn´ e pˇ r´ıpady:
♣ Odhad pr˚ umˇ ern´ e hodnoty Y , pˇ resnˇ eji odhad podm´ınˇ en´ e stˇ redn´ı – oˇ cek´ avan´ e – hodnoty veliˇ ciny Y vzhedem ke zvolen´ e kombinaci hodnot vysvˇ etluj´ıc´ı (vysvˇ etluj´ıc´ıch) promˇ enn´ e (promˇ enn´ ych): q
yˆ0 ± t1−α/2,n−pσ ˆ xt0(XtX)−I x0 ♣♣ Odhad konkr´ etn´ı hodnoty Y pˇ ri urˇ cit´ e kombinaci vysvˇ etluj´ıc´ı promˇ enn´ e, ˇ ci urˇ cit´ e kombinaci vysvˇ etluj´ıc´ıch promˇ enn´ ych: q
yˆ0 ± t1−α/2,n−pσ ˆ 1 + xt0(XtX)−I x0
c Rost 2007 °
Konfidenˇ cn´ı intervaly pro predikci v R
attach(belgie) range(prijem) x0<-seq(500,5200,10) prij<-data.frame(prijem=x0) pred.konfid<-predict(model,prij,se=T,interval="confidence") pred.pred<-predict(model,prij,se=T,interval="prediction") pred.konfid$fit[1:5,] fit lwr upr 1 390.0646 370.0255 410.1037 2 394.9164 375.0691 414.7636 3 399.7682 380.1105 419.4258 4 404.6200 385.1496 424.0903 5 409.4717 390.1864 428.7570
c Rost 2007 °
2000 1000
1500
Vydaje za jidlo
1500
500
1000 500
Vydaje za jidlo
2000
2500
2500
3000
Grafy predikˇ cn´ıch interval˚ u
1000
2000
3000 Prijem
4000
5000
1000
2000
3000
4000
5000
Prijem
c Rost 2007 °
Regresn´ı diagnostika Mezi z´ akladn´ı diagnostick´ e prostˇ redky patˇ r´ı pˇ redevˇ s´ım anal´ yza rezidu´ aln´ıch hodnot prostˇ rednictv´ım kvantilov´ ych graf˚ u spolu s diagnostick´ ymi statistikami DF BET AS, DF F IT S, COV RAT IO, Cookovou vzd´ alenost´ı a diagon´ aln´ımi prvky projekˇ cn´ı matice H (angl. leverage), kde H = X(XtX)Xt
c Rost 2007 °
Motivaˇ cn´ı pˇ r´ıklad Dalˇ s´ım datov´ ym souborem vyuˇ zit´ ym pro ilustraci je soubor obsahuj´ıc´ı morfometrick´ e´ udaje∗ nˇ ekolika lini´ı kapra obecn´ eho (tatajsk´ y kapr, amursk´ y sazan, syntetick´ y ˇ supinat´ y kapr, jihoˇ cesk´ y ˇ supinat´ y kapr, litomyˇ slsk´ yˇ supinat´ y kapr). U jednotliv´ ych ryb byla sledov´ ana celkov´ a d´ elka tˇ ela X1, d´ elka tˇ ela X2, d´ elka trupu za ˇ ritn´ı ploutv´ı vˇ cetnˇ e d´ elky hlavy X3, d´ ale d´ elka trupu pˇ red ˇ ritn´ı ploutv´ı vˇ cetnˇ e d´ elky hlavy X4, d´ elka hlavy X5 a v´ yˇ ska tˇ ela X6. Pro lepˇ s´ı n´ azornost jsou tyto charakteristiky zachyceny na obr´ azku .
∗ Vˇ sechny
d´ elky byly mˇ eˇ reny v mm v ose tˇ ela od pˇ redn´ıho okraje rypce. Pouze v´ yˇ ska tˇ ela byla stanovena prostˇ rednictv´ım kolmice spuˇ stˇ en´ e od prvn´ıho paprsku hˇ rbetn´ı ploutve k bˇ riˇ sn´ımu okraji tˇ ela. c Rost 2007 °
Motivaˇ cn´ı pˇ r´ıklad
c Rost 2007 °
Literatura Problematika je diskutov´ ana napˇ r´ıklad v n´ asleduj´ı literatuˇ re: Norman R. Draper, Harry Smith: Applied Regression Analysis,Wiley Series in Probability and Statistics, ISBN 1-58488-425-8 Julian J. Faraway: Linear Models with R, Chapman & Hall/CRC, Boca Raton, 2005, ISBN 1-58488-425-8 John Fox: An R and S-plus Companion to Applied Regression, Sage Publication, Thousand Oaks, 2002, ISBN 0-7619-2280-6
c Rost 2007 °
Dˇ ekuji za pozornost.
c Rost 2007 °
