1. Ryhmien jäsenet eivät esiinny muissa kuin omassa ryhmässään (ns. Independence-oletus)
2. Ryhmien arvot on poimittu normaalijakaumasta (tosin mallin residuaalien normaalisuus riittää)
3. Ryhmien varianssit ovat yhtäsuuret
4. Otoskoot ovat samansuuruiset.
Tärkein näistä oletuksista on ensimmäinen, koska jos joidenkin ryhmien jäsenet esiintyvät myös muissa ryhmissä, niin tulokset ovat vääristyneitä, eikä tälle ole muuta korjausta kuin palauttaa oletus toteen. Oletukset 2 ja 3 eivät usein täysin toteudu, mutta esim. varianssin yhtäsuuruus -oletukseen on keksitty ratkaisuja (esim.Welch-correction), joilla tuloksista saadaan luotettavia, vaikka oletus ei pätisi. Mikäli arvoja ei ole poimittu normaalijaukamasta, niin voi käyttää "ei-parametrista" testiä, kuten Kruskall-Wallis -testiä, jossa tätä oletusta ei ole. Toisaalta oleellisinta on, että mallin residuaalit ovat normaalisti jakautuneita, joten ennen Kruskall-Wallis -testin käyttöä on hyvä tutkia residuaalit. Viimeinen ehto on helppo toteuttaa, eli jos jostain ryhmästä on enemmän dataa kuin toisesta, niin käyttää tästä ryhmästä vain toisen ryhmän verran dataa, jolloin otoskoot ovat ryhmille samansuuruiset.
Teen seuraavaksi esimerkin sukupuolten keskipituuksien tutkimisesta. Muodostan datan, jossa kummankin sukupuolen pituudet ovat normaalijakautuneita. Ajatellaan, että data "populaatio" kertoo koko populaation pituudet ja että meillä on mahdollisuus ottaa tuosta populaatiosta 400:n havainnon otos. Tiedämme siis naisten ja miesten oikeat keskipituudet, joten voimme tutkita, pääsemmekö oikeisiin tuloksiin otoksen avulla.
R-koodi:
pituudet <- rnorm(1000,mean=170,sd=20) # Muodostetaan normaalijakautunut populaatio-data
sukupuoli <- factor(rep(c("N","M"),500))
populaatio <- data.frame(pituudet,sukupuoli)
f <- sample(nrow(populaatio),0.2*nrow(populaatio))
otos <- populaatio[f,]
tapply(otos$pituudet,otos$sukupuoli,mean) # Saadaan pituuden keskiarvo ja varianssi sukupuolen mukaan. Keskiarvot: M=170.54 & N= 170.92. Varianssit: M=349.47 & N=432.53
tapply(otos$pituudet,otos$sukupuoli,var)
boxplot(otos$pituudet~otos$sukupuoli) # Boxplot-kuva
library(car) # Seuraavaksi tehdään Levene-testi, jossa tutkitaan varianssien yhtäsuurutta. Nollahypoteesi varianssien yhtäsuuruudesta jää voimaan
leveneTest(otos$pituudet~otos$sukupuoli)
shapiro.test(otos$pituudet[which(otos$sukupuoli=="N")]) # Tehdään Shapiro-testi, joka kertoo, onko "pituudet" jakautunut normaalisti. Testin mukaan "pituudet" on jakautunut normaalisti
# Tehdään histogrammi, josta näemmä naisten ja miesten pituusjakaumat
hist(otos$pituudet[which(otos$sukupuoli=="N")],col="blue",main="Naisten pituudet",xlab="")
hist(otos$pituudet[which(otos$sukupuoli=="M")],col="red",main="Miesten pituudet",xlab="")
R-koodi:
oneway.test(otos$pituudet~otos$sukupuoli,var.equal=T) # P-arvo on 0.85, joten nollahypoteesi keskiarvojen yhtäsuuruudesta jää voimaan. Tämä ei ole yllätys, koska kumpikin sukupuoli sai arvonsa samasta normaalijakaumasta, jossa keskiarvo oli 170 ja keskihajonta 20.
Katsotaan vielä, etteivät mallin residuaalit poikkea hirveästi normaalijakautuneista.
m1 <- lm(otos$pituudet~otos$sukupuoli)
par(mfrow=c(2,1))
res <- resid(m1)
x <- seq(min(res),max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
hist(resid(m1),freq=F)
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)
Kuvista nähdään, että residuaalit ovat normaalijakautuneet. Residuaalien ei tarvitse kuitenkaan olla täysin normaalijakautuneita, koska otoskoon ollessa riittävän suuri Central limit theorem pätee, jolloin luottamusvälit saadaan laskettua luotettavasti. Katso tästä lisää loppukappaleesta.
Yllä olevassa kvantiili-kvantiili -käyrän kuvassa (QQ-plot) residuaalit laitetaan suuruusjärjestykseen ja niitä verrataan normaalijakauman kvantiileihin. Jos havainnot osuvat lähelle suoraa, niin havaintoja voidaan pitää normaalijakautuneina. Katsotaan vielä, miltä kuva näyttää, kun havainnot ovat positiivisesti ja negatiivisesti vinoista jakaumista.
Muodostetaan seuraavaksi uusi data, jossa sukupuolten pituudet tulevat eri normaalijakaumista ja eroavat.
pituudetN <- rnorm(500,mean=165,sd=5)
pituudetM <- rnorm(500,mean=185,sd=5)
sukupuoliN <- factor(rep("N",500))
sukupuoliM <- factor(rep("M",500))
populaatio2 <- data.frame(c(pituudetN,pituudetM),factor(c(as.character(sukupuoliN),as.character(sukupuoliM))))
colnames(populaatio2) <- c("pituudet","sukupuoli")
f <- sample(nrow(populaatio2),0.2*nrow(populaatio2))
otos <- populaatio[f,]
boxplot(otos$pituudet~otos$sukupuoli)
hist(otos$pituudet[which(otos$sukupuoli=="N")],col="blue",main="Naisten ja miesten pituudet",xlab="",xlim=c(150,210))
hist(otos$pituudet[which(otos$sukupuoli=="M")],col="red",xlab="",add=TRUE)
legend("topright",c("Naiset","Miehet"),col=c("blue","red"),lty=c(1,1),lwd=c(5,5))
Saimme siis tietoa, että eri sukupuolten keskipituuksissa on eroa. Emme kuitenkaan aina tiedä, millä ryhmillä nämä erot ovat, koska varianssianalyysin vaihtoehtoinen hypoteesi on, että ainakin kaksi keskiarvoa eroaa toisistaan tilastollisesti merkitsevästi. Esimerkissämme oli vain kaksi ryhmää, joten on helppoa nähdä, että miesten pituudet ovat suurempia. Mikäli ryhmiä on enemmän ja haluamme tietää, mitkä ryhmät eroavat toisistaan ja kuinka paljon, voimme tehdä Tukey:n HSD -testin. Huom. Testiä voi käyttää vain, jos varianssianalyysin nollahypoteesi keskiarvojen yhtäsuuruudesta on ensin hylätty. Lisäksi ryhmien varianssien tulee olla samansuuruisia. Jos näin ei ole, niin voi käyttää Games-Howell -testiä, joka löytyy R:n "userfriendlyscience" -pakkauksesta. Tehdään nyt Tukey:n HSD -testi esimerkkidatalle seuraavasti:
TukeyHSD(aov(otos$pituudet~otos$sukupuoli))
Testi antaa 95%:n luottamusvälin, jolle ryhmien keskiarvojen erotus osuu. Esimerkissämme tämä luottamusväli on naisten ja miesten pituuden erotukselle (N-M) : -21.35, -18.50.
Tarkastetaan vielä, että mallin residuaalit ovat approksimatiivisesti normaalisti jakautuneet. Alla olevasta kuvasta nähdään, että näin myös on.
m2 <- lm(otos$pituudet~otos$sukupuoli)
res <- resid(m2)
x <- seq(-15,max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
par(mfrow=c(2,1),mar=c(3,3,2,2))
hist(res,freq=F,ylim=c(0,0.074))
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)
Muodostan seuraavaksi dataa kolmesta ryhmästä, jotta näemme, kuinka yksisuuntainen varianssianalyysi toimii, kun ryhmiä on yli kaksi.
R-koodi:
pituusR1 <- rnorm(500,mean=165,sd=5)
pituusR2 <- rnorm(500,mean=185,sd=5)
pituusR3 <- rnorm(500,mean=175,sd=5)
ryhmä1 <- factor(rep("R1",500))
ryhmä2 <- factor(rep("R2",500))
ryhmä3 <- factor(rep("R3",500))
populaatio3 <- data.frame(c(pituusR1,pituusR2,pituusR3),factor(c(as.character(ryhmä1),as.character(ryhmä2),
as.character(ryhmä3))))
colnames(populaatio3) <- c("pituudet","ryhmä")
f <- sample(nrow(populaatio3),0.2*nrow(populaatio3))
otos <- populaatio3[f,]
boxplot(otos$pituudet~otos$ryhmä)
library(ggplot2) # Tehdään tiheyskuva
ggplot(otos, aes(pituudet, fill = ryhmä)) + geom_density(alpha = 0.7)
oneway.test(otos$pituudet~otos$ryhmä,var.equal=T) # H0 hylätään
TukeyHSD(aov(otos$pituudet~otos$ryhmä),conf.level=0.95) # Saadaan alla oleva taulukko
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = otos$pituudet ~ otos$ryhmä)
$`otos$ryhmä`
diff lwr upr p adj
R2-R1 20.53805 18.720077 22.356033 0
R3-R1 10.38782 8.537313 12.238317 0
R3-R2 -10.15024 -11.991966 -8.308513 0
R2:n keskipituus on suurempi kuin R1:n
R3:n keskipituus on suurempi kuin R1:n
R2:n keskipituus on suurempi kuin R3:n
Ryhmät ovat keskipituuden mukaisessa järjestyksessä siis R2, R3 ja R1.
Tarkastetaan lopuksi mallin residuaalit ja näemmä, että ne ovat approksimatiivisesti normaalisti jakautuneet.
m3 <- lm(otos$pituudet~otos$ryhmä)
res <- resid(m3)
par(mfrow=c(2,1),mar=c(3,3,2,2))
x <- seq(-15,max(res),length=length(res))
y=dnorm(x,mean=mean(res),sd=sd(res))
hist(res,freq=F)
lines(x,y)
qqnorm(res, main = "Normal Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
qqline(res)
Tarvitseeko residuaalien olla normaalisti jakautuneita?
Seuraavaksi katsotaan, kuinka Central limit theorem:n "lait" pätevät simuloidulle datalle. Simuloin regressiomalleja, joissa residuaalit noudattavat erilaisia jakaumia ja näytän, että poikkeustapauksia lukuunottamatta regressiomallin muuttujan regressiokerroin saa arvonsa normaalijakaumasta, riippumatta virhetermien jakaumasta. On tärkeää, että regressiokertoimet tulevat normaalijakaumasta, koska luottamusvälit ja p-arvot perustuvat vahvasti tähän oletukseen. Kun regressiokertoimet tulevat normaalijakaumasta, voimme esim. sanoa, että n.68% kerroista oikea regressiokertoimen arvo sijaitsee kerroinestimaattimme ja +-1:n keskivirheen sisällä. Jos jakauma on joku toinen, niin luottamusvälimme ja p-arvomme ovat todennäköisesti vääriä.
Regressiomalli on tässä simuloinnissa muotoa y=2+5*x+"virhetermi", eli oikea regressiokerroin on 5, joten katsotaan, kuinka lähelle tuota arvoa päästään, kun simuloidaan 1000 mallia erilaisilla virhetermin jakaumilla.
Ensimmäinen virhetermin jakauma on positiivisesti vino jakauma.
kertoimet = matrix(0,1000,2) # Luodaan matriisi, johon regressiomallin kertoimet tulevat
for(i in 1:1000) # Toistetaan alla oleva 1000 kertaa
{
x = rnorm(100)
y = 2 + 5*x + rbeta(100,1,5)
m = lm(y~x)
kertoimet[i,] = coef(m)
}
par(mfrow=c(2,2))
hist(rbeta(100,1,5),main="Positiivisesti vino jakauma (n=100)")
qqnorm(resid(m),main="Regressiomallin residuaalit QQ Plot")
qqline(resid(m))
qqnorm(kertoimet[,2],main="Kertoimet QQ Plot")
qqline(kertoimet[,2])
hist(kertoimet[,2], main="Kertoimet Histogram")
range(kertoimet[,2]) # Kertoimet vaihtelevat välillä 4.952 - 5.044
Seuraavaksi vastaavat kuvat, kun virhetermi seuraa negatiivisesti vinoa jakaumaa. Kertoimen arvo vaihtelee välillä 4.959-5.043
Yllä näimme, että virhetermien vinous ei vääristä regressiokertoimia. Seuraavaksi katsotaan ongelmallista T-jakaumaa pienillä vapausasteilla. Tällainen virhetermin jakauma hankaloittaa regressiokertoimen estimointia. Alla kuvassa diagnostiikkakuvat T-jakautuneille virhetermeille (df=2.01). Nyt regressiokerroin saa arvoja väliltä: 1.106 - 9.787, eli hännät ovat paksuja.
Seuraavaksi lisätään hieman vapausasteiden määrää (df=3), jolloin kertoimien jakauma alkaa jo muistuttamaan enemmän normaalijakaumaa ja kerroin vaihtelee lähempänä oikeaa: 4.359 - 6.158
Seuraavaksi kuva, kun df=5. Nyt kertoimet vaihtelevat välillä 4.525 - 5.437.
Katsotaan vielä kuvat, kun virhetermit tulevat normaalijakaumasta (mean=0,sd=1). Kertoimien vaihteluväli on 4.656-5.277.
Seuraavaksi lisätään hieman vapausasteiden määrää (df=3), jolloin kertoimien jakauma alkaa jo muistuttamaan enemmän normaalijakaumaa ja kerroin vaihtelee lähempänä oikeaa: 4.359 - 6.158
Seuraavaksi kuva, kun df=5. Nyt kertoimet vaihtelevat välillä 4.525 - 5.437.
Katsotaan vielä kuvat, kun virhetermit tulevat normaalijakaumasta (mean=0,sd=1). Kertoimien vaihteluväli on 4.656-5.277.
Ei kommentteja:
Lähetä kommentti