Ismétlés és bevezetés

Egymintás két-oldali t-próba

Legyen az \(X_1, X_2, \ldots, X_n\) minta háttéreloszlása: \(X \sim \mathcal{N}(\underbrace{\mu}_{\text{ismeretlen}}, \underbrace{\sigma}_{\text{ismeretlen}})\). Tehát az u-próbától annyiban különbözik, hogy most a szórás is ismeretlen.

Nullhipotézis és alternatív hipotézis: \[ H_0: \mu = m_0 \quad\text{vs.}\quad H_1: \mu \neq m_0 \]

A próbastatisztika

Az egymintás \(t\)-próba próbastatisztikája: \[t = \frac{\overline{X}_n - m_0}{s^*_n / \sqrt{n}} \, \overset{H_0}{\sim} \,t\text{-eloszlás } n-1 \text{ szabadságfokkal} ,\]

ahol

  • \(\bar{X}_n\) az \(n\)-elemű minta átlaga,
  • \(s^*_n\) a minta szórásának torzítatlan becslése (ha \(n\) elég nagy akkor használhatjuk \(s_n\)-t is),
  • \(m_0\) pedig előre megadott érték, amelyhez az átlagot viszonyítjuk.

A \(H_0\) hipotézis fennállása esetén az \(n\)-elemű mintán számolt \(t\) próbastatisztika \(n-1\) szabadságfokú Student t-eloszlású.

Lab 7 – t-próba folytatás

Végezzük el a második napi koleszterines adatokon a t-próbát. Ehhez először töltsük be az adatokat:

cholestr <- read.delim("./cholestr.txt")

A t-próba elvégzése R-ben nagyon egyszerű: a t.test függvényt kell hozzá meghívni. Tegyük fel (\(H_0\)), hogy az elméleti várhatóértéke a második napi koleszterin szintnek 260.

t.test(cholestr$Day2, mu=260)
## 
##  One Sample t-test
## 
## data:  cholestr$Day2
## t = -0.67337, df = 27, p-value = 0.5064
## alternative hypothesis: true mean is not equal to 260
## 95 percent confidence interval:
##  235.4284 272.4288
## sample estimates:
## mean of x 
##  253.9286

Ennek az eredményében több mindent is látunk, például a konfidencia intervallumot is, ami [235.4284, 272.4288], tehát nyugodtan elfogadhatjuk a nullhipotézist, azaz hogy habár a mintaátlag 253.9286, ez lehet a véletlen következménye is és igaz lehet, hogy valójában 260 az átlagos második napi koleszterinszint (a t.test függvény alapértelmezett konfidencia szintje 0.95).

A hipotézis vizsgálatok esetén a nullhipotézis elfogadásának/elutasítása a \(p\)-értéken és a konfidencia szinten múlik. A \(p\)-érték itt azt mondja meg, hogy ha igaz a \(H_0\), azaz a második napi koleszterin szintnek a váratóértéke \(\mu = 260\), akkor mi annak a valószínűsége, hogy a \(t\) tesztstatisztika (\(t=\frac{\bar{X}-\mu }{\widehat{\sigma}/\sqrt{n} }\)) a megfigyelt -0.673 értéket vagy ettől extrémebb értéket vesz fel. Láthatjuk, hogy ebben az esetben a \(p\)-érték = 0.5064. Ha az előírt elsőfajú hiba \(\varepsilon = \alpha = 0.05\), akkor az azt jelenti, hogy akkor utasítanánk el a \(H_0\)-t, ha a \(p\) érték legfeljebb 0.05 lenne.

P-érték

Eredete: “probability-value” (\(p\)-value)

A \(p\)-érték az az elsőfajú hibaszint (szignifikanciaszint) amikor épp határon lenne az elfogadás és az elutasítás.

Ha van egy \(S_n\) statisztikánk (mint val.változó), aminek ismerjük az eloszlását \(H_0\) mellett, és egy adott mintára nézve ennek a statisztikának a kiszámolt értéke \(s\), akkor a \(p\)-érték:

Jelentése

A \(p\)-érték azt mondja meg, hogy ha feltesszük, hogy a \(H_0\) igaz, akkor mi a valószínűsége annak, hogy a tesztstatisztika értéke \(s\) vagy \(s\)-nél is extrémebb értéket vesz fel.

Döntés: \(H_0\) vagy \(H_1\)?

Döntés a \(p\)-érték alapján adott \(\alpha\) szignifikancia-szint (0.05, 0.01, 0.001, stb.) mellett:

Miért van erre szükség, mikor vizsgálhatjuk a szerkesztett elfogadási intervallumot is? A két irány teljesen ekvivalens eredményt ad. A p-érték előnye a következő:

Kérjük le a múltkori egymintás kétoldali t-próba p-értékét, és számoljuk ki manuálisan is.

Ehhez használjuk a már jól ismert koleszterinszinteket tartalmazó adatot.

cholestr <- read.delim("./cholestr.txt") 
eredmenyek <- t.test(cholestr$Day2, mu=260)
eredmenyek
## 
##  One Sample t-test
## 
## data:  cholestr$Day2
## t = -0.67337, df = 27, p-value = 0.5064
## alternative hypothesis: true mean is not equal to 260
## 95 percent confidence interval:
##  235.4284 272.4288
## sample estimates:
## mean of x 
##  253.9286

A számolt \(p\)-érték külön is lekérhető mivel a próba eredményeit egy listaként kaptuk vissza:

eredmenyek$p.value
## [1] 0.5064333

Két oldali tesztet végeztünk, és a Student t-eloszlás szimmetrikus, így visszakapjuk a \(p\)-értéket a következő formulával: \[ p = 2\left(1-F_{t_{27}}(|s|)\right), \] ahol \(F_{t_{27}}\) a 27 szabadságfokú t-eloszlás eloszlásfüggvénye, és \(s\) a számolt T statisztika értéke.

s <- eredmenyek$statistic
dof <- eredmenyek$parameter 

sajatp <- 2*(1-pt(abs(s), dof))

names(sajatp) <- 'p-érték'
sajatp
##   p-érték 
## 0.5064333
sajatp == eredmenyek$p.value
## p-érték 
##    TRUE

t-eloszlás sűrűségfüggvénye

Ismételjük meg az 5. labor 4. feladatát, azaz hasonlítsuk össze a (standard) normális eloszlás és a Student t-eloszlás sűrűségfüggvényét különböző szabadságfokok esetén.

Ahhoz hogy ábrázolni tudjuk a sűrűségfüggvényeket, először szükségünk van egy vektorra, ami az \(x\)-tengely pontjait tartalmazza, ahol kiértékeljük majd a sűrűségfüggvényeket. Ezt a vektort a seq függvénnyel tudjuk létrehozni, aminek meg kell adni a kezdőpontot, végpontot és a lépésközöket. Például:

x = seq(-4, 4, 0.5)

x
##  [1] -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0

Ha az x pontokban szeretnénk ábrázolni egy függvényt, akkor ki kell értékelnünk a függvényt az x vektor elemein. Ezt megtehetjük úgy is, hogy egy for ciklussal végigiterálunk az x elemein, vagy használhatunk erre apply függvényeket is, amiknek megadunk egy vektort/mátrixot/dataframe-et és egy függvényt, és megkapjuk a függvény értékét elemenként/soronként/oszloponként. Itt megnézhetitek hogyan működnek az apply(), lapply(), sapply(), tapply() függvények.

(A függvények beépített súgóját/dokumentációját elérhetitek az RStudio Help fülén, és ott a függvény nevére keresve; vagy a függvény nevére a kattintva a kódon belül az F1 billentyű nyomvatartása mellett.)

De szerencsére a sűrűségfüggvényeknek vektort is meg lehet adni, ami kiszámolja nekünk a vektor minden pontjára a függvény értékét.

R plot függvényes megoldás

A plot függvényről ezen az oldalon találtok egy jó leírást, példákkal.

x <- seq(-4.5, 4.5, 0.01)

plot(x, dnorm(x), type="l", col="blue", ylab="density", main="Standard normális eloszlás sűr. fv.-e")

Ha az előbbi ábrához hozzá szeretnénk adni a \(t\)-eloszlás sűrűségfüggvényét is, akkor azt a lines paranccsal tehetjük meg, illetve a legend-del lehet feliratot adni az ábrához:

x <- seq(-4.5, 4.5, 0.01)

plot(x, dnorm(x), type="l", col="blue", ylab="density", main="Standard normális vs t-eloszlás")
lines(x, dt(x, df=3), col="orange")
legend("topleft", c("Standard normal", "Student's t (df=3)"), fill=c("blue", "orange"))

Megoldás ggplot2 használatával

library(ggplot2)

x <- seq(-4.5, 4.5, 0.01)

df <- data.frame(x_tengely = x, 
                 t_suruseg = dt(x, df=5),
                 norm_suruseg = dnorm(x))

ggplot(data=df, aes(x=x_tengely)) + 
  geom_line(aes(y=norm_suruseg, color="Standard normális eloszlás")) + 
  geom_line(aes(y=t_suruseg, , color="Student t-eloszlás (df = 5)")) + 
  ylab("Sűrűségfüggvény")+
  xlab("X-tengely")+
  theme(legend.position = "bottom")

Írjunk egy függvényt, aminek egy k paramétere van és az outputja az olyan mint az előbbi ábra csak t-eloszlásnak a szabadsági foka legyen k.

t_sur <- function(k) {
  x <- seq(-4.5, 4.5, 0.01)
  df <- data.frame(x_tengely = x, 
                   t_suruseg = dt(x, df=k),
                   norm_suruseg = dnorm(x))
  
  ggplot(data=df, aes(x=x_tengely)) + 
    geom_line(aes(y=norm_suruseg, color="Standard normális eloszlás")) + 
    geom_line(aes(y=t_suruseg, , color=paste(c("t-eloszlás (df = ", k, ")"), collapse=""))) + 
    ylab("Sűrűségfüggvény")+
    xlab("X-tengely")+
    theme(legend.position = "bottom")
  
}
t_sur(2)

t_sur(30)

DataCamp kurzusok és segédanyagok

ggplot2

t-eloszlás és t-próba

A t-eloszlásról, t-próbáról, p-value-ról van egy jó fejezet a DataCamp Inference for Numerical Data in R című kurzusában (Introducing the t-distribution), illetve ugyanebben a kurzusban az Inference for difference in two parameters fejezet témája a hipotézis vizsgálat t-próbával.

T-próba variánsai

Egymintás egyoldali t-próba

Korábban már volt róla szó, hogy nem csak egyenlőséget tesztelhetünk. Egyoldali esetben a hipotéziseink így módosulnak: \[ H_0: \mu \le \mu_0 \quad\text{vs.}\quad H_1: \mu > \mu_0 \] vagy: \[ H_0: \mu \ge \mu_0 \quad\text{vs.}\quad H_1: \mu < \mu_0 \]

Fontos, hogy ez nem fog változtatni a kiszámolt statisztika értékén, csak az elfogadási intervallum fog kibővülni az egyik irányba (az első esetben balra \(-\infty\) felé, míg a második esetben jobbra \(\infty\) felé). Hiszen azzal, hogy az egyik irányú kilengés \(\mu_0\)-hoz képest átkerül a \(H_0\)-ba, ezek a kilengések már nem számítanak nekünk “extrémnek”, csak a másik irányba esők. Ez a \(p\)-érték fenti egyoldali definícióiból is látszik.

Most végezzük el a t-próba egyoldali variánsát arra \(H_0\)-ra, hogy a várható érték \(\le 260\). Ezt az “alternative” paraméter állítja, ide azt szeretnénk, hogy a 260-nál nagyobb értékek essenek.

egyoldali <- t.test(cholestr$Day2, mu=260, alternative = "greater")
egyoldali
## 
##  One Sample t-test
## 
## data:  cholestr$Day2
## t = -0.67337, df = 27, p-value = 0.7468
## alternative hypothesis: true mean is greater than 260
## 95 percent confidence interval:
##  238.571     Inf
## sample estimates:
## mean of x 
##  253.9286
egyoldali <- t.test(cholestr$Day2, mu=260, alternative = "less")
egyoldali
## 
##  One Sample t-test
## 
## data:  cholestr$Day2
## t = -0.67337, df = 27, p-value = 0.2532
## alternative hypothesis: true mean is less than 260
## 95 percent confidence interval:
##      -Inf 269.2862
## sample estimates:
## mean of x 
##  253.9286

Amit még észre kell venni a fenti két példában, hogy a konfidenciaintervallum hogyan változott. A konfidencia-intervallum az alternatív hipotézis miatt alakult így. Tradicionálisan az alternatív hipotézisbe (\(H_1\)) kerül(nek) az(ok) eset(ek), amikre “számítunk”. Ezt követi az R is. Ennek ellenére a \(H_0\) alapján is építhető a konfidencia-intervallum, csak akkor a másik irányba kell nyitni a végtelen felé.

A \(p\)-értékeket visszakaphatjuk a kétoldali esetben számoltból, kihasználva a Student t-eloszlás szimmetriáját.

Kétmintás kétoldali t-próba

Független kétmintás t-próba

Először azt a változatát nézzük a t-próbának, ahol két egymástól független, akár különböző nagyságú, mintánk van, és ezek várható értékeit akarjuk összehasonlítani.

Azonos szórású eloszlások

A következő teszt/próba akkor használható, ha feltehetjük, hogy a két minta háttéreloszlásának a szórása azonos.

Tekintsük a következő háttéreloszlásokat és az azokból vett \(n_1\) és \(n_2\) elemű mintákat:

  • \(X \sim \mathcal{N}(\mu_X, \sigma)\), minta: \(X_1, \ldots, X_{n_1}\)
  • \(Y \sim \mathcal{N}(\mu_Y, \sigma)\), minta: \(Y_1, \ldots, Y_{n_2}\)

Fontos, hogy a minták egymástól is függetlenek.

Hipotézisek

\[ H_0: \mu_X = \mu_Y \quad \text{vs.} \quad H_1: \mu_X \neq \mu_Y \]

Itt a \(t\) próbastatisztika a két minta átlag különbségén és úgynevezett pooled (közös, összevont) variancián alapszik, lásd ezen az oldalon.

Próbastatisztika \[t = \frac{\overline{X}_{n_1} - \overline{Y}_{n_2}}{ \sqrt{\frac{(n_1-1)s_X^2 + (n_2-1)s_Y^2}{n_1 + n_2 -2}} } \sqrt{\frac{n_1n_2}{n_1+n_2}}\]

Ez a próbastatisztika is t-eloszlású, \(df = n_1+n_2-2\).

Különböző szórású eloszlások (Welch-próba)

Ebben az esetben nem tesszük fel, hogy az \(X_1,\dots,X_{n_1}\) és \(Y_1,\dots, Y_{n_2}\) minták szórása megegyezik. Ezt a próbát Welch-próbának is nevezik.

Ekkor a következő háttéreloszlásaink és független mintáink vannak:

  • \(X \sim \mathcal{N}(\mu_X, \sigma_X)\), \(X_1, \ldots, X_{n_1}\)
  • \(Y \sim \mathcal{N}(\mu_Y, \sigma_Y)\), \(Y_1, \ldots, Y_{n_2}\)

A Welch-próba próbastatisztikája szintén a minta átlagok különbségén alapszik, de itt már nem használhatunk összevont varianciát:

Próbastatisztika \[t = \frac{\overline{X}_{n_1} - \overline{Y}_{n_2}}{\sqrt{\frac{s_X^2}{n_1} + \frac{s_Y^2}{n_2}}}\]

Ez a próbastatisztika is t-eloszlású, a szabadsági foka komplikáltabb, lásd itt, ahol lenn az f jelöli a szabadságfokot.

R példa

Nézzük meg hogy a betegek negyedik napi koleszterinszintje várhatóértékben különbözik-e az egészségesek koleszterinszintjétől. Ehhez töltsük le erről a honlapról az első file-t (tab-delimited text, general format) és töltsük be az R-be.

cholestg <- read.delim("./cholestg.txt") 

names(cholestg)
## [1] "patient" "group"   "day"     "cholest"

Ennek a data frame-nek a következő oszlopai/változói vannak:

  • patient: a páciensek azonosítója,
  • group: az 1-es csoport jelöli a betegeket, a 2-es csoport jelöli az egészségeseket,
  • day: szívinfarktus után hány nappal mérték a koleszterinszintet,
  • cholest: koleszterinszint.
View(cholestg)

Kérjük le a betegek negyedik napi és az egészségesek koleszterinszintjét:

egeszsegesek = cholestg[cholestg$group == 2, "cholest"]

betegek = cholestg[cholestg$group == 1 & cholestg$day == 4, "cholest"]

Itt ugye a két minta természetesen függetlennek tekinthető egymástól. A szórásaikról nincs információnk, így biztonságosabb azt mondani, hogy különböznek, és automatikusan a Welch-tesztet alkalmazni.

(Egy másik út lenne a szórások esetleges azonosságának előzetes statisztikai vizsgálata, amiről a következő heti anyagban lesz szó. Ha meggyőznénk magunkat, hogy a háttérváltozók szórásai azonosak, használhatnánk a sima független kétmintás kétoldali t-próbát. Ezt a módszert több kritika is érte az irodalomban, de ez messzire vezetne most minket, így most nem térünk ki rá.)

A Welch-próba elvégzéséhez csak meg kell adni két változót a t.test függvénynek. Azért elég ezt megtennünk, mert a t.test függvény, ami a mindenféle t-teszt változatokért felel, alapértelmezetten a var.equal = FALSE paraméter szerepel.

welch_proba = t.test(egeszsegesek, betegek)
welch_proba
## 
##  Welch Two Sample t-test
## 
## data:  egeszsegesek and betegek
## t = -3.8411, df = 37.991, p-value = 0.0004512
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -57.27838 -17.74067
## sample estimates:
## mean of x mean of y 
##  193.1333  230.6429
welch_proba$p.value
## [1] 0.0004511934

A \(p\)-érték alapján azonnal látjuk, hogy 0.05, 0.01, de még 0.001 szignifikanciaszint mellett is elutasítjuk a nullhipotézist, és az alternatíva mellett döntünk.

Kétmintás egyoldali t-próba

A kétmintás t-próba mindkét fenti változatára létezik egyoldali változat is. A próbastatisztikák értéke ismét nem fog változni, csak az elfogadási/kritikus intervallum, és ezáltal a \(p\)-érték. Az iménti betegek vs. egészséges példára, először legyen: \[ H_0: \mu_X \le \mu_Y \quad \text{vs.} \quad H_1: \mu_X > \mu_Y \] Vagyis azt tesztelnénk, hogy az egészségeseknél magasabb a koleszterinszint. Megjegyezzük, hogy nem ez a természetes feltételezés a példában, de ettől függetlenül tesztelhető.

Nehéz lehet eldönteni, hogy mit is kellene módosítani a fv.-híváson, hogy a megfelelő tesztet végezzük el. Ilyenkor hívjuk segítségül a fv. dokumentációját. Ebben a következő el: alternative = "greater" is the alternative that x has a larger mean than y. Vagyis az alternative = "greater" hívás tesztelné, amit akarunk.

welch_proba_g = t.test(egeszsegesek, betegek, alternative = "greater")
welch_proba_g
## 
##  Welch Two Sample t-test
## 
## data:  egeszsegesek and betegek
## t = -3.8411, df = 37.991, p-value = 0.9998
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -53.97338       Inf
## sample estimates:
## mean of x mean of y 
##  193.1333  230.6429

A fenti összegzés az alábbi alakban tartalmazza a hipotézist: \[ H_0: \mu_X-\mu_Y\le 0 \quad \text{vs.} \quad H_1: \mu_X - \mu_Y>0 \]

Ez természetesen ekvivalens azzal, amit tesztelni akartunk. A magas \(p\)-érték miatt nem tudjuk elutasítani a nullhipotézist, amellett döntünk.

Most nézzük a másik irányú tesztet. Vagyis azt tesztelnénk, hogy a betegeknél magasabb a koleszterinszint. Ez a természetes feltételezés a példában. \[ H_0: \mu_X \ge \mu_Y \quad \text{vs.} \quad H_1: \mu_X < \mu_Y \]

welch_proba_l = t.test(egeszsegesek, betegek, alternative = "less")
welch_proba_l
## 
##  Welch Two Sample t-test
## 
## data:  egeszsegesek and betegek
## t = -3.8411, df = 37.991, p-value = 0.0002256
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf -21.04567
## sample estimates:
## mean of x mean of y 
##  193.1333  230.6429

A fenti összegzés az alábbi alakban tartalmazza a hipotézist: \[ H_0: \mu_X-\mu_Y\ge 0 \quad \text{vs.} \quad H_1: \mu_X - \mu_Y<0 \]

Az alacsony \(p\)-érték miatt elutasíthatjuk a nullhipotézist, és az alternatíva mellett döntünk, amit tesztelni akartunk. Itt is láthatjuk, hogy a (megfelelő, azaz kisebb) egyoldali \(p\)-érték éppen fele a két-oldali tesztnél látott \(p\)-értéknek.

Megjegyzés a végére

Fontos, hogy a fenti kétoldali tesztek független minták esetén értelmesek. A következő heti anyagban megnézzük a fenti kétmintás tesztek összetartozó mintás változatát, ami valójában egymintás próbaként is értelmezhető. Ezen kívül kitérünk a szórással kapcsolatosan fenn megjegyzett problémára is, és nézünk tesztet a szórások egyezőségére.