R statisztika

Kategória: Adattudomány R-ben

Adatok eloszlása

Áttekintés

Sokféle eloszlású számokkal lehet 4 féle műveletet végrehajtani. A függvények mindegyike a műveletnek megfelelő betűvel kezdődik:

  • d - sűrűség (density), pl. normális eloszlás esetén ez a haranggörbe
  • p - kumulatív eloszlás (cumulative distribution); ez a sűrűség kumulálásából adódik, 0-ból indul, 1-be ér
  • q - kvantilis (quantile); ez a p inverze
  • r - véletlen szám generálás (random number generation)

Néhány fontosabb eloszlás:

  • norm - normális
  • unif - egyenletes
  • pois - Poisson
  • exp - exponenciális

Pl. a rnorm(5) 5 darab standard normális eloszlású véletlen számot generál.

Normális eloszlás

A normális eloszlás, más néven haranggörbe, a statisztika egyik legfontosabb eloszlása. Olyan adatokat ír le, amelyek a központi érték körül csoportosulnak, és az értékek szimmetrikusan csökkennek a szélek felé. Sok természeti jelenség, mint például az emberi testmagasság vagy a IQ-tesztek eredményei, követik ezt az eloszlást.

Függvények:

  • dnomr(): ennek az intuitív jelentése az, hogy adott értéket mekkora valószínűséggel vehet fel. Matematikailag ez a megfogalmazás nem teljesen precíz, ugyanis pont azt az értéket 0 eséllyel veheti fel; precízebben: a görbe alatti terület összege 1.
  • pnorm(): a kumulatív valószínűséget adja meg, azaz a görbe alatti területet addig a pontig. Azt mutatja meg, hogy mekkora esélye van annak, hogy az érték legfeljebb annyi. A lenti példában a standard normális eloszlás esetén annak az esélye, hogy az érték kisebb mint -1, kb. 15.9%.
  • qnorm(): a pnorm() inverze. Egy adott valószínűséghez tartozó értéket adja meg. Például a qnorm(0.159) eredménye -1, mivel -1 az az érték, amire igaz az az állítás, hogy annak az esélye, hogy az érték nála kisebb, pont 0.159.
  • rnorm(): véletlen szám generálás.
x <- seq(-4, 4, length = 100)
plot(x, dnorm(x), type = "l", col = "blue", lwd = 2, main = "Normális eloszlás", xlab = "Érték (x)", ylab = "Sűrűség")
p_val <- pnorm(-1)
x_shade <- seq(-4, -1, length = 100)
y_shade <- dnorm(x_shade)
polygon(c(-4, x_shade, -1), c(0, y_shade, 0), col = "gray")
text(-3, 0.1, paste("P(X < -1) =", round(p_val, 3)), col = "black", cex = 1)
normal.png

Egyenletes eloszlás

Az egyenletes eloszlás akkor használható, ha egy adott tartományon belül minden lehetséges kimenet valószínűsége azonos. Például egy szabályos dobókocka dobásának eredménye egyenletes eloszlású. Gyakran használják szimulációkban véletlenszámok generálásához.

Az unif folytonos egyenletes eloszlású, és az R-ben nincs is diszkrét egyenletes. Függvények:

  • dunif(): ez a lehetséges értékek szakaszán egy konstans érték, máshol nulla.
  • punif(): ez egy 0-tól 1-ig egyenesen növő szakasz a minimum és a maximum érték között.
  • qunif(): 0 és 1 között értelmezett, ami egyenletesen nő minimumtól maximumig.
  • runif(): véletlen szám generálás.

Mindegyik függvénynek második és harmadik paraméterként megadhatjuk a minimum és maximum értéket, ami alapértelmezésben 0 és 1.

Sűrűségfüggvénye:

x <- seq(0, 10, length = 100)
plot(x, dunif(x, 2, 5), type = "l", col = "blue", lwd = 2, main = "Egyenletes eloszlás", xlab = "Érték (x)", ylab = "Sűrűség")
uniform.png

Az alábbi példa 10 egyenletes eloszlású véletlen számot generál 2 és 5 között:

runif(10, 2, 5)

Egy lehetséges futási eredmény:

 [1] 3.676575 4.860592 4.029691 4.289486 2.725711 2.016264 2.736340 2.547027
 [9] 3.542124 3.403058

Exponenciális eloszlás

Az exponenciális eloszlás gyakran használatos az események közötti idő modellezésére, például a vevők érkezése egy üzletbe vagy a telefonhívások egy call centerben. Az eloszlás kulcsfontosságú jellemzője, hogy az idő múlásával az esemény bekövetkezésének valószínűsége csökken.

Függvények:

  • dexp()
  • pexp()
  • qexp()
  • rexp()

A függvények második paramétere a rate, azaz arány, ami azt mondja meg, hogy egységnyi idő alatt átlagosan hány esemény következik be. Pl. ha a rate értéke 2, akkor azt azt jelenti, hogy az események között átlagosan eltelt idő 0.5.

A sűrűségfüggvénye az alábbi:

x <- seq(0, 10, length = 100)
plot(x, dexp(x), type = "l", col = "blue", lwd = 2, main = "Exponenciális eloszlás", xlab = "Érték (x)", ylab = "Sűrűség")
exponential.png

Poisson

A Poisson eloszlás egy adott idő- vagy térintervallumon belül bekövetkező események számát modellezi. Például a tévesztések száma egy oldalon vagy a balesetek száma egy útszakaszon. Az eloszlást a lambda (λ) paraméter határozza meg, amely az események átlagos előfordulását jelenti.

  • Függvények:
  • dpois()
  • ppois()
  • qpois()
  • rpois()

A második paraméter itt a lambda, ami kötelező, nincs alapértelmezett értéke.

A Poisson eloszlás diszkrét. A sűrűségét emiatt oszlopdiagrammal illusztráljuk. Itt egy-egy oszlop azt jelenti, hogy adott értéket mekkora valószínűséggel vesz fel a λ=5 Poisson eloszlású valószínűségi változó. Tehát ha például az átlagos gépelési hibák száma oldalanként 5, akkor mekkora eséllyel lesz pont 2, pont 5, pont 14 stb.

x <- 0:15
y <- dpois(x, lambda = 5)
barplot(
    y,
    names.arg = x,
    main = "Poisson-eloszlás (lambda = 5)",
    xlab = "Események száma (x)",
    ylab = "Valószínűség",
    col = "lightblue"
)
poisson.png

További eloszlások

A támogatott eloszlásokat a következő paranccsal tudjuk lekérdezni:

?Distributions

Statisztikai próbák

Feltételvizsgáló tesztek

Bizonyos statisztikai próbáknak vannak előfeltételei, például az eloszlással vagy a varianciával kapcsolatosan. Először ezeket a feltételvizsgáló teszteket nézzük meg.

Normalitás teszt (Shapiro-Wilk)

A normalitás teszt azt vizsgálja, hogy egy numerikus változó eloszlása szignifikánsan eltér-e a normális (haranggörbe alakú) eloszlástól. A példában megnézzük, hogy az autók fogyasztása (egészen pontosan azok inverze, azaz az egy gallon üzemanyaggal megtehető mérföldek száma) tekinthető-e normálisnak:

shapiro.test(mtcars$mpg)

Eredmény:

        Shapiro-Wilk normality test

data:  mtcars$mpg
W = 0.94756, p-value = 0.1229

Itt a null-hipotézis az, hogy az adatok normális eloszlásúak. Mivel a p érték elég nagy (nagyobb mint 0.05), nem vethetjük el a null-hipotézist, azaz az adatok tekinthetően normális eloszlásúnak.

Vizuális illusztráció:

hist(mtcars$mpg, freq = FALSE, col = "lightblue", main = "Fogyasztás (mpg) eloszlása", xlab = "mpg")
lines(density(mtcars$mpg), col = "red", lwd = 2)
normalitasteszt.png

Egy másik gyakori illusztrációja a Q-Q plot:

qqnorm(
    mtcars$mpg,
    main = "Q-Q plot: Fogyasztás (mpg)",
    xlab = "Elméleti kvantilisek",
    ylab = "Minta kvantilisek"
)
qqline(mtcars$mpg, col = "red", lwd = 2)
qqplot.png

Variancia homogenitás teszt (Levene)

Ez a teszt azt vizsgálja, hogy két vagy több csoportban a varianciák (szórások) szignifikánsan eltérnek-e egymástól. Az alábbi példában azt vizsgáljuk meg, hogy a fogyasztás varianciája ugyanakkora-e kézi és automata váltó esetén:

var.test(mpg ~ am, data = mtcars)

Eredmény:

        F test to compare two variances

data:  mpg by am
F = 0.38656, num df = 18, denom df = 12, p-value = 0.06691
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.1243721 1.0703429
sample estimates:
ratio of variances 
         0.3865615

Itt is a null-hipotézis az, hogy nincs szignifikáns eltérés. Mivel a p-érték 0.05 feletti, nem tudjuk elvetni a null-hipotézist, azaz a két variancia nem tér el szignifikánsan.

par(mfrow = c(1, 2))
breaks <- seq(min(mtcars$mpg), max(mtcars$mpg), length.out = 8)
hist(mtcars$mpg[mtcars$am == 0],
     main = "Automata váltós autók",
     xlab = "Fogyasztás (mpg)",
     ylab = "Gyakoriság",
     breaks = breaks,
     xlim = c(10, 35),
     col = "lightblue",
     border = "black")
hist(mtcars$mpg[mtcars$am == 1],
     main = "Kézi váltós autók",
     xlab = "Fogyasztás (mpg)",
     ylab = "Gyakoriság",
     breaks = breaks,
     xlim = c(10, 35),
     col = "lightgreen",
     border = "black")
varianciateszt.png

Korrelációs tesztek

Ezek a tesztek két numerikus változó közötti kapcsolatot vizsgálják, azaz azt, hogy van-e kapcsolat a két változó között.

Pearson korrelációs teszt

Azt vizsgálja, hogy van-e lineáris kapcsolat két numerikus változó között. Ez érzékeny a kiugró értékekre.

A példában azt vizsgáljuk, hogy van-e kapcsolat a fogyasztás és a tömeg között:

cor.test(mtcars$mpg, mtcars$wt, method = "pearson")

Eredmény:

        Pearson's product-moment correlation

data:  mtcars$mpg and mtcars$wt
t = -9.559, df = 30, p-value = 1.294e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9338264 -0.7440872
sample estimates:
       cor 
-0.8676594

Ahogy számítani lehetett rá, igen erős kapcsolat áll fenn a két érték között. A kapcsolat itt amiatt negatív, mert nem a fogyasztást, hanem annak inverzét (egységnyi mennyiségű üzemanyag által megtett út) vizsgáljuk.

A korreláció igen szemléletes az alábbi ábrán:

plot(mtcars$wt, mtcars$mpg, pch = 16, col = "blue",
     main = "Fogyasztás és tömeg kapcsolata",
     xlab = "Tömeg (1000 font)", ylab = "Fogyasztás (mpg)")
abline(lm(mpg ~ wt, data = mtcars), col = "red", lwd = 2)
korrelacio.png

Spearman korrelációs teszt

A sorrendi kapcsolatot vizsgálja, azaz azt, hogy ha az elemeket egyik ill. másik érték szerint rendezzük, akkor milyen mértékben kapjuk ugyanazt a sorrendet. Ez nem érzékeny a kiugró értékekre.

cor.test(mtcars$mpg, mtcars$wt, method = "spearman")

Khi-négyzet teszt

A Khi-négyzet teszt két kategorikus változó közötti kapcsolatot vizsgál. Azt teszteli, hogy az egyik változó értékei függetlenek-e a másikétól.

Az alábbi teszt azt vizsgálja, hogy a a hengerek száma (cyl) és a sebességfokozatok száma (gear) független-e.

contingency_table <- table(mtcars$cyl, mtcars$gear)
chisq.test(contingency_table)

Eredmény:

        Pearson's Chi-squared test

data:  contingency_table
X-squared = 18.036, df = 4, p-value = 0.001214

Warning message:
In chisq.test(contingency_table) :
  Chi-squared approximation may be incorrect

A p-érték alapján kijelenthető, hogy a két érték nem független egymástól. Lássunk intuitív módon! Először nézzük meg a kontingencia táblát (contingency_table)!

     3  4  5
  4  1  8  2
  6  2  4  1
  8 12  0  2

Ha függetlenek lennének, akkor az eloszlás sorok és oszlopok szerint is hasonló lenne:

  • Sorok szerint: a sebességfokozatok száma független lenne a hengerek számától. Olyan értelemben független, hogy majdnem minden kombináció létezik, de az eloszlás egészen más. A négyhengeres autók többsége 4 fokozatú, a nyolchengeresek többsége pedig 3 fokozatú.
  • Oszlopok szerint: a háromsebességes autók többség nyolchengeres, a négysebességesek mindegyike négy ill. hathengeres. Öthengeresből akkor még kevés volt, és annak viszonylag egyenletes volt az eloszlása.

Lássuk oszlopdiagram segítségével:

barplot(
    contingency_table,
    beside = TRUE,
    col = c("lightblue", "lightgreen", "lightyellow"),
    main = "Hengerek száma és a sebességfokozatok kapcsolata",
    xlab = "Sebességfokozatok száma",
    ylab = "Autók száma",
    legend.text = rownames(contingency_table),
    args.legend = list(title = "Hengerek száma")
)
khinegyzet.png

Ha függetlenek lennének, akkor az eloszlás nagyjából egyenletes lenne.

A hőtérkép is jól illusztrálja a kapcsolat hiányát:

library(ggplot2)

df <- as.data.frame(contingency_table)
colnames(df) <- c("cyl", "gear", "count")

ggplot(df, aes(x = cyl, y = gear, fill = count)) +
  geom_tile(color = "white") +
  geom_text(aes(label = count), color = "black", size = 6) +
  scale_fill_gradient(low = "lightblue", high = "darkblue", name = "Autók száma") +
  labs(title = "Hengerek és sebességfokozatok kapcsolata",
       x = "Hengerek száma",
       y = "Sebességfokozatok") +
  theme_minimal()
heatmap.png

Összehasonlítás tesztek

Egymintás t-próba

Az egymintás t-próba a minta átlagát hasonlítja össze egy hipotetikus populációs átlaggal. A nullhipotézis azt mondja ki, hogy a minta átlaga megegyezik a populációs átlaggal.

Példaként próbáljuk meg megbecsülni a teljesítmények átlagát úgy, hogy vesszük a legkisebb és a legnagyobb teljesítményű autó átlagát:

avg_hp <- (max(mtcars$hp)+ min(mtcars$hp)) / 2

Eredmény:

193.5

Vajon ez jól közelíti-e a tényleges átlagot?

t.test(mtcars$hp, mu = avg_hp)

Eredmény:

        One Sample t-test

data:  mtcars$hp
t = -3.8623, df = 31, p-value = 0.0005348
alternative hypothesis: true mean is not equal to 193.5
95 percent confidence interval:
 121.9679 171.4071
sample estimates:
mean of x 
 146.6875

A tényleges átlag 147 körüli, azaz közel 50-nel kevesebb, és a teszt ki is mutatta, hogy a két legszélsőségesebb érték átlaga nem jól közelíti a tényleges átlagot ebben az esetben. Ha megnézzük az adatok eloszlását, akkor látni fogjuk, hogy miért:

boxplot(mtcars$hp)
egymintas.png

Egy felül kiugró érték nagyon felhúzza a két szélsőséges érték átlagát.

Egymintás U-próba

Az egymintás U-próba a megadott értéket nem az átlaggal, hanem a mediánnal hasonlítja össze.

wilcox.test(mtcars$hp, mu = avg_hp)

Eredmény:

        Wilcoxon signed rank test with continuity correction

data:  mtcars$hp
V = 89.5, p-value = 0.001134
alternative hypothesis: true location is not equal to 193.5

Warning message:
In wilcox.test.default(mtcars$hp, mu = avg_hp) :
  cannot compute exact p-value with ties

Az eredmény nem tartalmazza, de a tényleges medián 123. Ebben az esetben is szignifikáns az eltérés a tényleges medián és a két legszélső érték átlaga között.

Megjegyzés: a kétmintás U-próba másik neve Wilcox teszt; a függvény neve az egymintás és a páros esetben is wilcox.test().

Kétmintás t-próba

A t-próba azt vizsgálja, hogy két adathalmaz átlaga eltér-e jelentősen egymástól. Vegyük példaként az automata váltós és a kézi váltós autók fogyasztását!

boxplot(
    mpg ~ am,
    data = mtcars,
    col = c("lightgreen", "lightblue"),
    main = "Fogyasztás (mpg) a sebességváltó típusa szerint",
    xlab = "Sebességváltó (0: automata, 1: kézi)",
    ylab = "Mérföld/gallon (mpg)"
)
ketmintas.png

Ebből az ábrából az a sejtésünk, hogy az automata váltós autók fogyasztása magasabb (egységnyi üzemanyaggal kevesebb utat lehet megtenni), mint a kézi váltósé.

t.test(mpg ~ am, data = mtcars)

Az eredmény meggyőző:

        Welch Two Sample t-test

data:  mpg by am
t = -3.7671, df = 18.332, p-value = 0.001374
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -11.280194  -3.209684
sample estimates:
mean in group 0 mean in group 1 
       17.14737        24.39231

Itt érdemes megjegyezni, hogy az alap t-teszt előfeltételei között szerepel az, hogy az összehasonlítandó értékek közel normális eloszlásúak legyenek, és a varianciájuk is hasonló legyen. Ezeknél az adatoknál már megvizsgáltuk, hogy ez (nagyjából) teljesül, viszont az R által alkalmazott Welch kétmintás t-próbának már ez nem előfeltétele.

Kétmintás U-próba

A kétmintás t-próba érzékeny a kiugró értékekre; egy szélsőséges eset is nagyon félre viszi az átlagot és ezáltal a teszt eredményét. A kétmintás U-próba (amit Wilcox-tesztnek is hívnak) a sorrendiséget vizsgálja; jelen esetben tehát azt, hogy ha sorba rendezzük az autókat fogyasztás szerint, akkor a sebességváltó típusa teljesen véletlenszerű lesz-e, vagy lesz benne valamiféle rendszer, pl. előbb jönnek az egyik majd a másik fajták. Itt egyetlen (vagy kis számú) kilógó eset érdemben nem befolyásolja az eredményt.

wilcox.test(mpg ~ am, data = mtcars)

Az eredmény ugyancsak meggyőző:

        Wilcoxon rank sum test with continuity correction

data:  mpg by am
W = 42, p-value = 0.001871
alternative hypothesis: true location shift is not equal to 0

Warning message:
In wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...) :
  cannot compute exact p-value with ties

Páros t-próba

A páros próbák tipikusan olyanok, hogy ugyanazokra a példányokra vannak előtte-utána mérések. Például egy vérnyomáscsökkentő gyógyszer vizsgálata esetén megmérjük a gyógyszer alkalmazása előtt és az alatt. Sajnos az ismert adathalmazokban ilyen jellegű méréseket nem találtam, emiatt itt kivételesen generált adatokkal dolgozunk. Tegyük fel, hogy valaminek az eloszlása normális, a várható értéke 100, a szórása 1. Ezekhez az értékekhez szimuláltan hozzáadunk egy olyan értéket, aminek a várható értéke és a szórása is 1.

Az értékek generálása:

a <- rnorm(10, 100)
diff <- rnorm(10, 1)
b <- a + diff

a:

100.84723 101.51646 100.91806  98.26597  99.44357  99.92761  99.49150  98.15583 100.04470 100.41947

diff:

 0.44451250  2.26588293  0.52425370  1.57982843  0.46865469  1.08064446 -0.95542330  1.35718786  2.32523957 -0.08679268

Itt látható, hogy az értékek alapvetően nőttek, de bizonyos értékek még csökkentek is.

b:

101.29174 103.78235 101.44232  99.84580  99.91223 101.00826  98.53608  99.51301 102.36994 100.33268

A kérdés az az, hogy a b értékek páronként szignifikánsan nagyobbak-e, mint az a értékek, tehát úgy, hogy az adott indexű elemeket hasonlítjuk össze a két adatsorban. A teszt:

t.test(b, a, paired = TRUE)

Az eredmény:

        Paired t-test

data:  b and a
t = 2.7638, df = 9, p-value = 0.02197
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.1634373 1.6373604
sample estimates:
mean difference 
      0.9003988

A p-érték tehát 0.05 alatti, azaz szignifikánsnak tekinthető.

Páros U-próba

A t-teszt ebben az esetben is érzékeny a kiugró értékekre. A páros U-próba a sorrendiséget veszi alapul, ezáltal jól kezeli a szélsőséges értékeket.

wilcox.test(b, a, paired = TRUE)

Az eredmény ebben az esetben is szignifikáns.

        Wilcoxon signed rank exact test

data:  b and a
V = 49, p-value = 0.02734
alternative hypothesis: true location shift is not equal to 0

ANOVA

Az ANOVA teszt több mint két csoport átlagát hasonlítja össze. Ezzel megállapíthatjuk, hogy van-e szignifikáns különbség az átlagok között, de azt nem, hogy melyik csoportpár között van ez a különbség.

Példa: van-e szignifikáns különbség a fogyasztásban a hengerek száma szerint?

anova_model <- aov(mpg ~ as.factor(cyl), data = mtcars)
summary(anova_model)

Eredmény:

               Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(cyl)  2  824.8   412.4    39.7 4.98e-09 ***
Residuals      29  301.3    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Az eredmény szignifikáns. Ezt a következő ábra is jól illusztrálja:

boxplot(
    mpg ~ cyl,
    data = mtcars,
    col = c("lightblue", "lightgreen", "lightyellow"),
    main = "Fogyasztás (mpg) a hengerek száma szerint",
    xlab = "Hengerek száma",
    ylab = "Mérföld/gallon (mpg)"
)
anova.png
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License