1. No category

MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. harjoitukset / Tehtävät
Kotitehtävä: 3,4
Tehtävä 2.1.
Jatkoa tietokonetehtävälle 1.2:
(a)
Piirrä aineistosta pisteparvikuvaaja (KULUTUS, SAIRAST) ja siihen estimoitu
regressiosuora. KULUTUS on selittävä muuttuja.
(b)
Määrää estimoidusta mallista sovitteet yˆi ja residuaalit ei ja tallenna ne tiedostoon
”tupakka.txt” muuttujiksi FIT (= sovite) ja RES (= residuaali).
(c)
Piirrä pistediagrammit (SAIRAST, FIT) ja (FIT, RES).
(d)
Tutki havainnon 7 = USA:n poikkeavuutta kohdassa (c) piirretyissä kuviossa.
(e)
Tutki havainnon 7 = USA poikkeavuutta Cookin etäisyyksien avulla. Voisiko USA olla
poikkeava havainto?
(f)
Estimoi malli ilman havaintoa USA ja vertaa tuloksia tehtävän 1.4. tuloksiin.
Tehtävä 2.1. – Mitä opimme?
Tehtävässä tarkastellaan regressioanalyysin tulosten havainnollistamista sekä
poikkeavien havaintojen tunnistamista ja vaikutusta.
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Tehtävä 2.1. – Ratkaisu:
data=read.table(”tupakka.txt”,header=T,sep=”\t”)
tupakka = data[,c(“SAIRAST”,”KULUTUS”)]
malli=lm(SAIRAST~KULUTUS,data=tupakka)
(a)
Pisteparvikuvaaja:
plot(tupakka[,”KULUTUS”],tupakka[,”SAIRAST”])
abline(malli,col=”red”)
USA
Havainnot voidaan merkitä kuvaajaan komennon identify() avulla. (Käyttö ei ihan
triviaalia, tule harjoituksiin jos haluat lisätietoa komennosta.)
(b)
Sovitteet ja residuaalit ovat tallennettuna malliin nimillä ”fit” ja” res”, näihin voidaan
viitata seuraavasti:
malli$fit
malli$res
Tallennetaan sovitteet
yˆi
(fitted values = FIT) ja residuaalit ei (residuals = RES).
tupakka$FIT=malli$fit
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
tupakka$RES=malli$res
Kirjoitetaan vielä data-frame tupakka tiedostoon tupakka.txt
write.table(tupakka,”tupakka2.txt”)
(c)
Pistediagrammi (Selitettävä, Sovite).
Piirretään sovitteet
yˆi
selitettävän muuttujan SAIRAST arvoja vastaan.
plot(tupakka$SAIRAST,tupakka$FIT)
USA
Diagrammi kuvaa mallin hyvyyttä:
( yi , yˆi ) , i = 1, 2, … , n ovat suoraa
–
Malli on sitä parempi, mitä lähempänä pisteet
viivaa.
–
Myös poikkeavat havainnot erottuvat usein selvästi.
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Huomaa, että pisteistä ( yi , yˆi ) , i = 1, 2, … , n määrätty Pearsonin tulomomenttikorrelaatiokertoimen neliö on sama kuin selitysaste:
Cor( y, yˆ )2  R2
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Pistediagrammi (Sovite, Residuaali).
Piirretään residuaalit ei sovitteita
yˆi
vastaan.
plot(malli$fit,malli$res)
USA
Diagrammi kuvaa mallin hyvyyttä:
(d)
( yˆi , ei ) , i = 1, 2, … , n ovat suoraa
–
Malli on sitä parempi, mitä lähempänä pisteet
e = 0.
–
Myös poikkeavat havainnot erottuvat usein selvästi.
Havainnon 7 = USA poikkeavuus näkyy kaikissa yo. kuvissa.
MS-C2128 Ennustaminen ja aikasarja-analyysi
(e)
2. tietokoneharjoitukset
Talletetaan Cookin etäisyydet (DIST) sarakkeeksi ja piirretään niitä kuvaamaan
pistediagrammi, jossa vaaka-akselilla on muuttuja MAA.
cooks.distance(malli)
tupakka$DIST=cooks.distance(malli)
n <- dim(tupakka)[1]
plot(seq(1,n),cooks.distance(malli))
USA
(f)
Estimoidaan malli ilman havaintoa 7 = USA.
malli2=lm(SAIRAST[-7]~KULUTUS[-7],data=tupakka)
summary(malli2)
Call:
lm(formula = SAIRAST[-7] ~ KULUTUS[-7])
Residuals:
Min
1Q
-62.353 -28.923
Coefficients:
Median
-7.861
3Q
35.321
Max
66.919
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.55343
28.26713
0.479
0.644
KULUTUS[-7] 0.35767
0.04547
7.867 4.93e-05 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 44.92 on 8 degrees of freedom
Multiple R-squared: 0.8855, Adjusted R-squared: 0.8712
F-statistic: 61.88 on 1 and 8 DF, p-value: 4.928e-05
Kulmakertoimen estimaatti on kasvanut arvosta 0.23 arvoon 0.36. Tämä todistaa
(jäljelle jääneiden 10 maan joukossa) paljon voimakkaammin tupakoinnin
vaarallisuudesta.
Kysymys: Saako havainnon 7 =USA:n poistaa?
Vastaus:
USA:n poikkeuksellisen pieni keuhkosyöpätapausten suhteellinen lukumäärä
saattoi johtua siitä, että siellä savukkeet olivat erilaisia (tupakan laatu oli
miedompaa ja savukkeissa oli filtterit) kuin muissa tutkimuksen maissa.
Jos näin oli, USA:ta voidaan pitää poikkeavana havaintona, jonka
sulkeminen pois estimoinnista on luvallista. Muista, että havaintoja ei voida
sivuuttaa ilman päteviä perusteluja!
Tehtävä 2.2.
Sementin kovettuessa kehittyy lämpöä, jonka määrä riippuu sementin koostumuksesta.
Tiedostossa HALD.txt on seuraavat tiedot 13 erilaisesta sementtierästä:
HEAT
= lämpömäärä cal/g
CHEM1, CHEM2, CHEM3, CHEM4 = sementin ainesosia (% kuiva-aineesta)
(a)
Estimoi regressiomalli, jossa ovat mukana kaikki selittäjät. Vertaile kertoimien
tilastollista merkitsevyyttä ja tarkastele vastaavien selittäjien varianssin inflaatiotekijöitä.
(b)
Etsi paras selittäjien yhdistelmä käyttäen jotain mallinvalintakriteeriä, esim. Mallowsin
Cp-tunnuslukua.
Tehtävä 2.2. – Mitä opimme?
Tehtävässä näytetään miten erilaiset mallinvalintastrategiat saattavat johtaa erilaisiin
malleihin.
Tehtävä 2.2. – Ratkaisu:
Tavoitteena tehtävässä on selvittää, mitkä selittävistä tekijöistä
CHEM1, CHEM2, CHEM3, CHEM4
vaikuttavat selitettävän muuttujan
HEAT
käyttäytymiseen.
Ladataan aluksi data, ja asennetaan paketti ”car” myöhempää käyttöä varten.
install.packages("car")
library(car)
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
hald=read.table("hald.txt",header=T)
attach(hald)
Komento attach avulla voimme viitata suoraan muuttujan hald sarakkeisiin. Kokeile
kirjoittaa esimerkiksi CHEM1 R-konsoliin.
(a)
Täyden mallin estimointi
Tilanteessa, jossa ei olla selvillä, mitkä selittäjistä vaikuttavat selitettävän muuttujan
käyttäytymiseen, on usein järkevää estimoida ensin ns. täysi malli eli malli, jossa
käytetään selittäjinä kaikkia selittäjäkandidaatteja.
Ennen parhaan selitysmallin etsimistä on syytä tarkastella muuttujien välisiä
korrelaatioita:
cor(hald)
CHEM1
CHEM2
CHEM3
CHEM4
HEAT
SUM
CHEM1 1.00000000 0.2285795 -0.8241338 -0.2454451 0.7307175 0.05010722
CHEM2 0.22857947 1.0000000 -0.1392424 -0.9729550 0.8162526 -0.26044918
CHEM3 -0.82413376 -0.1392424 1.0000000 0.0295370 -0.5346707 -0.11025122
CHEM4 -0.24544511 -0.9729550 0.0295370 1.0000000 -0.8213050 0.32907694
HEAT
0.73071747 0.8162526 -0.5346707 -0.8213050 1.0000000 -0.16458053
SUM
0.05010722 -0.2604492 -0.1102512 0.3290769 -0.1645805 1.00000000
Muuttuja HEAT korreloi melko voimakkaasti kaikkien selittäjäkandidaattien kanssa.
Korrelaatio on positiivinen muuttujien CHEM1 ja CHEM2 kanssa ja negatiivinen
muuttujien CHEM3 ja CHEM4 kanssa. Selittäjäkandidaatit CHEM1 ja CHEM3
korreloivat voimakkaan negatiivisesti keskenään, samoin kandidaatit CHEM2 ja
CHEM4.
Estimoidaan seuraavaksi täysi malli
(1)
HEAT = 0 + 1 CHEM1 + 2 CHEM2 + 3 CHEM3 + 4 CHEM4 + 
taysimalli=lm(HEAT~CHEM1+CHEM2+CHEM3+CHEM4)
summary(taysimalli)
Estimointitulokset mallista (1):
Call:
lm(formula = HEAT ~ CHEM1 + CHEM2 + CHEM3 + CHEM4)
Residuals:
Min
1Q
-3.1750 -1.6709
Median
0.2508
3Q
1.3783
Max
3.9254
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.4054
70.0710
0.891
0.3991
CHEM1
1.5511
0.7448
2.083
0.0708 .
CHEM2
0.5102
0.7238
0.705
0.5009
CHEM3
0.1019
0.7547
0.135
0.8959
CHEM4
-0.1441
0.7091 -0.203
0.8441
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.446 on 8 degrees of freedom
Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736
F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07
Mallin (1) selitysaste on korkea (98.2 %). F-testisuureen arvo nollahypoteesille
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
H0 : 1 = 2 = 3 = 4 = 0
on 111.5 ja sitä vastaava p-arvo on (4:llä desimaalilla) 0.0000, joten malli on
kokonaisuudessaan tilastollisesti erittäin merkitsevä ja ainakin yksi regressiokertoimista
1, 2, 3, 4 poikkeaa nollasta.
Kuitenkaan yksikään mallin (1) selittäjistä ei ole tilastollisesti merkitsevä, jos
merkitsevyyden rajana pidetään (tavanomaista) 5 %:n merkitsevyystasoa. Tämä johtuu
selittäjien multikollineaarisuudesta.
Selittäjien multikollineaarisuutta voidaan mitata VIF-kertoimilla. Jos selittäjän i ja
muiden selittäjien välillä ei ole lineaarista tilastollista riippuvuutta (eli selittäjä i on
ortogonaalinen muiden selittäjiin suhteen),
VIFi = 1
Mitä suurempi on kertoimen VIFi arvo, sitä voimakkaammin selittäjä i riippuu
lineaarisesti muista selittäjistä. Jos
VIFi > 10
jollekin selittäjälle i, multikollineaarisuudesta saattaa olla haittaa.
Vif-kertoimet saadaan paketin car() funktiolla vif()
vif(taysimalli)
CHEM1
CHEM2
38.49621 254.42317
CHEM3
CHEM4
46.86839 282.51286
Mallissa (1) sekä selittäjiä CHEM2 ja CHEM4 vastaavien varianssin inflaatiotekijöiden
arvot ovat suurempia kuin 200, mikä viittaa voimakkaaseen multikollineaarisuuteen.
Tarkastellaan selittäjien multikollineaarisuutta estimoimalla regressiomallit, joissa
selitettävinä muuttujina ovat muuttujat CHEM2 ja CHEM4 ja kummassakin tapauksessa
selittäjinä käytetään kaikkia muita alkuperäisen mallin (1) selittäjiä.
Olkoon mallina
(2)
CHEM2 = 0 + 1 CHEM1 + 3 CHEM3 + 4 CHEM4 + 
Estimointitulokset mallista (2):
mallic2=lm(CHEM2~CHEM1+CHEM3+CHEM4)
summary(mallic2)
Call:
lm(formula = CHEM2 ~ CHEM1 + CHEM3 + CHEM4)
Residuals:
Min
1Q
-2.2494 -0.7280
Median
0.3881
3Q
0.7033
Max
0.9512
Coefficients:
Estimate Std. Error t value
(Intercept) 96.59382
2.16253
44.67
CHEM1
-0.97860
0.10602
-9.23
CHEM3
-1.00350
0.09443 -10.63
CHEM4
-0.97759
0.02111 -46.30
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01
Pr(>|t|)
7.06e-12
6.94e-06
2.15e-06
5.12e-12
***
***
***
***
‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.126 on 9 degrees of freedom
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Multiple R-squared: 0.9961, Adjusted R-squared: 0.9948
F-statistic: 760.3 on 3 and 9 DF, p-value: 3.864e-11
Mallin selitysaste on 99.6 %, joten CHEM2 riippuu hyvin voimakkaasti muista
selittäjistä.
Huomaa, että muuttujan CHEM2 VIF-kerroin mallissa (1) on
VIF2 
1
1  R22
2
jossa R2 on selitysaste mallista (2).
Olkoon mallina
(3)
CHEM4 = 0 + 1 CHEM1 + 2 CHEM2 + 3 CHEM3 + 
Estimointitulokset mallista (3):
Call:
lm(formula = CHEM4 ~ CHEM1 + CHEM2 + CHEM3)
Residuals:
Min
1Q
-2.3264 -0.6836
Median
0.4439
3Q
0.7463
Max
1.0379
Coefficients:
Estimate Std. Error t value
(Intercept) 98.65079
1.94627 50.687
CHEM1
-1.00504
0.10175 -9.878
CHEM2
-1.01865
0.02200 -46.303
CHEM3
-1.02809
0.09187 -11.191
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01
Pr(>|t|)
2.27e-12
3.96e-06
5.12e-12
1.39e-06
***
***
***
***
‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.15 on 9 degrees of freedom
Multiple R-squared: 0.9965, Adjusted R-squared: 0.9953
F-statistic: 844.5 on 3 and 9 DF, p-value: 2.413e-11
Mallin selitysaste on 99.7 %, joten CHEM4 riippuu voimakkaasti muista selittäjistä.
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Huomaa, että muuttujan CHEM4 VIF-kerroin mallissa (1)
VIF4 
1
1  R42
2
jossa R4 on selitysaste mallista (3).
Perussyy multikollineaarisuudelle mallissa (1) on se, että sementti koostuu lähes
kokonaan ainesosista CHEM1, CHEM2, CHEM3, CHEM4: Muuttujien CHEM1,
CHEM2, CHEM3, CHEM4 summa vaihtelee välillä 95 – 99 %. Siten yhden ainesosan
lisäämisen on pakko vähentää joidenkin muiden ainesosien osuutta sementin koostumuksessa. Tämä on selittää sen, miksi muuttujapareilla CHEM1 ja CHEM3 sekä
CHEM2 ja CHEM4 on voimakkaat negatiiviset korrelaatiot.
(b)
Paras selittäjien yhdistelmä
Regressiomallin selittäjien valikointiin voidaan käyttää erilaisia strategioita.
Valittaessa parasta selittäjien yhdistelmää kaikkia mahdollisia mallivaihtoehtoja
verrataan toisiinsa käyttämällä jotakin mallinvalintakriteeriä ja lopulliseksi malliksi
valitaan se, joka on käytetyn kriteerifunktion mielessä optimaalinen.
Tilastotieteellisessä kirjallisuudessa on esitetty lukuisia mallinvalintakriteereitä.
Tunnettuja kriteereitä ovat esim. Akaiken informaatiokriteeri AIC, Schwarzin
bayeslainen informaatiokriteeri SBIC (BIC), Hannanin ja Quinnin kriteeri HQ,
Allenin PRESS ja CAT.
Mallin valinnan kriteerifunktio on muotoa
min
𝑀⊆{1,…,𝑞}
2
𝐶(|𝑀|, 𝜎̂|𝑀|
).
2
Missä M on selittäjäkandidaattien yhdistelmä, ja 𝜎̂|𝑀|
sitä vastaavan mallin
jäännösvarianssin suurimman uskottavuuden estimaattori, ja C näiden suhteen kasvava
funktio. Yleisesti kriteerifunktiolta toivotaan:
(1) Mahdollisimman suurta selitysastetta
(2) Mahdollisimman vähillä selittäjillä.
R:ssä paketin ”leaps” funktio ”regsubsets()” laskee kullekin selittäjämäärälle regression,
joka maksimoi selitysasteen (tai muun halutun suureen). Tällöin saadaan tietyllä tapaa
lista parhaista regressioista. Tehtävän tapauksessa koodi
install.packages(”leaps”)
library(leaps)
osaj=regsubsets(HEAT~CHEM1+CHEM2+CHEM3+CHEM4,data=hald)
summary(osaj)
tulostaa
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
Subset selection object
Call: regsubsets.formula(HEAT ~ CHEM1 + CHEM2 + CHEM3 + CHEM4, data = hald)
4 Variables (and intercept)
Forced in Forced out
CHEM1
FALSE
FALSE
CHEM2
FALSE
FALSE
CHEM3
FALSE
FALSE
CHEM4
FALSE
FALSE
1 subsets of each size up to 4
Selection Algorithm: exhaustive
CHEM1 CHEM2 CHEM3 CHEM4
1 ( 1 ) " "
" "
" "
"*"
2 ( 1 ) "*"
"*"
" "
" "
3 ( 1 ) "*"
"*"
" "
"*"
4 ( 1 ) "*"
"*"
"*"
"*"
Tästä nähdään, että CHEM4 on paras yksittäinen selittäjä muuttujalle HEAT, kun taas
pari (CHEM1,CHEM2) selittää paremmin kuin mikään muu kahden selittäjän
osajoukko. Jos halutaan tutkia mallinvalintakriteereitä, niin niitä voidaan visualisoida
helposti. Esimerkkinä Mallowsin Cp:
plot(osaj, scale=”Cp”)
Kuvaa tulkitaan seuraavasti: Pysty-akselilla on Mallowsin Cp-luku (pienempi on
parempi). Musta laatikko merkitsee, että malli sisältää vaaka-akselin muuttujan. Tällöin
MS-C2128 Ennustaminen ja aikasarja-analyysi
2. tietokoneharjoitukset
tulkitaan, että Cp-luvun mielessä paras malli on malli, joka sisältää selittäjät CHEM1 ja
CHEM2. Vastaava Cp-luku on 2.7.
Kokeile argumentille ”scale” arvoja: ”bic” ja ”adjr2”. Mallille voidaan laskea Akaiken
informaatiokriteeri kirjoittamalla
AIC(malli)
Huomautus:
Jäännösneliösummaa tai selitysastetta ei voida käyttää mallinvalintakriteereinä, koska
sekä jäännösneliösumman minimointi että selitysasteen maksimointi johtavat aina
täyteen (maksimaaliseen) malliin (esimerkin tapauksessa malliin, jossa on selittäjinä
kaikki kandidaatit CHEM1, CHEM2, CHEM3, CHEM4).
Tehtävä 2.3.
Jatkoa tehtävälle 2.2. Käytä mallin valinnassa alaspäin askellusta. Suorita alaspäin
askellus permutaatiotestin avulla. Käytä apunasi luentokalvoja sekä viime viikon
demotehtävää 1. Vertaa tuloksia tehtävän 2.2 (b)-kohtaan.
Tehtävä 2.4.
Lannoiteaineen määrä vaikuttaa vehnän satoon. Määrän vaikutusta tutkittiin vaihtelemalla
lannoiteaineen määrää (11 tasoa) 33:lla koealalla (sama määrä lannoitetta 3:lla koealalla) ja
rekisteröimällä saatu sato. Tiedot kokeesta on annettu tiedostossa VehnanSato.txt.
Muuttujina tiedostossa ovat
Sato
= Sadon määrä (yksikkönä kg/pinta-alayksikkö)
Lannoite = Lannoiteaineen määrä (yksikkönä kg/pinta-alayksikkö)
(a)
Estimoi lineaarinen regressiomalli, jossa selitettävänä on muuttuja Sato ja selittäjänä
muuttuja Lannoite. Tutki mallin hyvyyttä regressiografiikan avulla.
(b)
Estimoi lineaarinen regressiomalli, jossa kohdan (a) malliin on lisätty selittäjäksi
muuttuja
LSqrd = Lannoite×Lannoite
Tutki mallin hyvyyttä regressiografiikan avulla.
(c)
Kumpi malleista on parempi? Miksi?