Интервали поверења

Дефиниција: Нека је \(\theta\) непознат параметар у расподели обележја \(X\) и нека је \((X_1,X_2,\dots,X_n)\) прост случајан узорак обима \(n\) за посматрано обележје. Нека су \(U_n\) и \(V_n\) статистике дефинисане на основу узорка такве да је \(P\{U_n\leq V_n\}=1\) за које важи \(P\{U_n\leq \theta \leq V_n\} = \beta\). Интервал \([U_n, V_n]\) се назива \(\beta\%\) двострани интервал поверења за параметар \(\theta\), а \(\beta\) је ниво поверења.

Претпоставимо да узорак \(X_1,X_2\dots X_n\) потиче из нормалне \(N(m,\sigma^2)\) расподеле, тада знамо да вaжи следеће: \[\frac{\bar{X}_n-m}{S_n}\sqrt{n-1} \sim t_{n-1},\]

где је \(\bar{X}_n\) узорачка средина, а \(\bar{S}_n\) (поправљена) узорачка стандардна девијација.

Имајући у виду ову расподелу, интервал поверења нивоа \(\alpha\) за непознати параметар \(m\) се лако изводи и добија се да је једнак \[\left(\bar{X}_n-C\frac{\bar{S}_n}{\sqrt{n-1}},\, \bar{X}_n+C\frac{\bar{S}_n}{\sqrt{n-1}}\right)\]

где је \(C=F^{-1}_{t_{n-1}}\left(\frac{1+\alpha}{2}\right)\).

Овакав интервал се може користити не само за нормалну расподелу, већ и за друге расподеле ако је узорак довољно велики да може да се примени централна гранична теорема. У том случају тражи се интервал поверења за очекивање, што је у случају нормалне баш \(m\).

• Направимо функцију која за дати узорак враћа овакав интервал поверења.

confidence_interval <- function(x, alfa = 0.95) {
  n <- length(x) # obim uzorka
  # iz formule vrednost C - kvantil t raspodele
  C <- qt((1 + alfa)/2, df = n - 1)
  
  # vracamo interval poverenja
  c(mean(x) - C * sd(x) / sqrt(n-1),mean(x) + C * sd(x) / sqrt(n-1))
}

Примена фуннкције:

x <- rnorm(50)
confidence_interval(x, 0.95)

## [1] -0.1314343  0.3573549

Обратимо пажњу на тренутак на интерпретацију нивоа поверења \(\alpha\). То је вероватноћа да добијени интервал поверења обухвати стварну вредност параметра \(m\). То значи да ако је \(\alpha=0.95\), да ако много пута извућемо узорак, у 95\(\%\) случајева ће инетрвал поверења садржати вредност \(m\). Испитјамо то:

Прво генеришемо 10.000 узорака и одговарајућих интервала поверења.

intervals <- replicate(1e4, {
  x <- rnorm(50)
  confidence_interval(x, 0.95)
})
intervals[, 1:5]

##            [,1]       [,2]        [,3]       [,4]        [,5]
## [1,] -0.2526964 -0.3868946 -0.01446017 -0.4292072 -0.43089978
## [2,]  0.3334463  0.2251439  0.53234242  0.2114943  0.06359259

Као резултат добијамо матрицу са 2 врсте и 10000 колона, где свака колона представља један интервал поверења.

Погледајмо колико од тих интервала поверења садржи нулу, која је била стварна вредност параметра \(m\).

# pravimo logicki vektor koji oznacava da li interval sadrzi nulu
contains_zero <- apply(intervals, 2, function(interval) {
  interval[1] < 0 && 0 < interval[2]
})

# gledamo u koliko intervala je sadrzana nula
mean(contains_zero)

## [1] 0.9484

Видимо да је нула садржана у приближно 95\(\%\) интервала.

Задатак.

Из популације чије обележје \(X\) има нормалну \(\mathcal{N}(m,\sigma^2)\) расподелу извучен је узорак обима \(20\) и добијени су резултати: 12.81, 9.35, 11.03, 12.34, 5.86, -2.00, -0.74, 7.03, 7.93, 14.86, 9.91, 6.61, -5.58, 5.32, -0.69, 13.07, 1.52, 6.60, 3.93, -1.81.

(а) Одредити \(90\%\) интервал поверења за непознат параметар \(m\).

(б) Одредити \(90\%\) једнострани (горњи и доњи) и двострани интервал поверења за непознат параметар \(\sigma^2\), као и за непознат параметар \(\sigma\).

Приметимо да је \(m\) непознато, па ћемо користити случајну величину \(\frac{n\bar{S}^2_n}{\sigma^2}\) која има \(\chi^2_{n-1}\) расподелу. Интервал поверења за непознат параметар \(\sigma^2\) је \[I_{\sigma^2} = \left[\frac{n\bar{S}^2_n}{\chi^2_{n-1;\frac{1-\beta}{2} }}, \, \frac{n\bar{S}^2_n}{\chi^2_{n-1;\frac{1+\beta}{2} }}\right]\]

x <- c(12.81, 9.35, 11.03, 12.34, 5.86, -2.00, -0.74, 7.03, 7.93, 14.86, 9.91, 6.61, -5.58,5.32, -0.69, 13.07,1.52, 6.60, 3.93, -1.81)

confidence_interval_m <- function(x, alfa = 0.90) {
  n <- length(x) 
  C <- qt((1 + alfa)/2, df = n - 1)
  
  c(mean(x) - C * sd(x) / sqrt(n-1),mean(x) + C * sd(x) / sqrt(n-1))
}

confidence_interval_m(x)

## [1] 3.554153 8.180847

confidence_interval_s2 <- function(x, alfa = 0.90) {
  n <- length(x) 
  c1 <- qchisq((1 + alfa)/2, df = n - 1)
  c2 <- qchisq((1 - alfa)/2, df = n - 1)
  c(n * var(x) / c1, n * var(x) / c2)
}

confidence_interval_s2(x)

## [1] 22.56391 67.22891

confidence_interval_s <- function(x, alfa = 0.90) {
  sqrt(confidence_interval_s2(x))
}

confidence_interval_s(x)

## [1] 4.750148 8.199324

Овај инетрвал се може применити и на неке друге расподеле кад је узорак велики (за тражење интервала поверења за очекивање).

Посматрајмо на пример експоненцијалну расподелу.

mean(replicate(1e4, {
  lambda = 2
  x <- rexp(1000, 2)
  interval <- confidence_interval(x, 0.95)
  interval[1] < 1/lambda && 1/lambda < interval[2]
}))

## [1] 0.95

Опет је у приближно 95 случајева стварна вредност очекивања \(\frac{1}{\lambda} = \frac{1}{2}\) упала у интервале поверења.

Статистички тестови

Студентов \(t\) тест

One Sample t-test - тест једног узорка

Ако имамо узорак из нормалне \(N(m, \sigma^2)\) расподеле можемо тестирати хипотезу \(H_0 (m=m_0)\) против неке од алтернативних облика \[H_1(m<m_0), \quad H_1(m\neq m_0), \quad H_1(m> m_0).\]

За ово тестирање можемо да користимо тест статистику: \[t=\frac{\bar{X}-m_0}{\bar{S}}\sqrt{n} \sim t_{n-1}.\] Критична област овог тест, са нивоом поверења \(\alpha\) је облика (за одговарајуће алтернативне) \[W=\left\{t<F^{-1}_{t_{n-1}}(\alpha)\right\}, \quad W=\left\{|t|>F^{-1}_{t_{n-1}}\left(1-\frac{\alpha}{2}\right)\right\}, \quad W=\left\{t>F^{-1}_{t_{n-1}}(1-\alpha)\right\}.\]

У статистичким пакетима тетсирање се обично врши налажењем \(p\) вредности тетса, па се пореди та вредност са нивоом значајности. Угрубо, \(p\) вредност теста се може описати као “количина доказа за нулту хипотезу”, уколико је велика \((p>\alpha)\), онда не одбацујемо нулту хипотезу, а ако је \((p<\alpha)\), онда одбацујемо нулту хипотезу у корист алтернативне.

Ако је алтернативна хипотеза двострана, \(p\) вредност можемо да израчунамо на следећи начин: \[p=P\{|t|>|t_0|\},\] где је \(t_0\) реализована вредност тест статистике на основу узорка.

У случају \(t\) теста \(p\) вредност ће бити једнака \[p=F_{t_{n-1}}(t_0), \quad p=2(1-F_{t_{n-1}}(|t_0|)), \quad p=1-F_{t_{n-1}}(t_0)\] за одговарајуће алтернативне хипотезе, редом.

Имплементирајмо овај тест за један узорак.

pval_t_test <- function(x, m0, alternative){
  n <- length(x)
  stat <- (mean(x) - m0)/sd(x)*sqrt(n)
  
  if(alternative == "less") {
    pval <- pt(stat, df = n - 1)
  } else if(alternative == "two.sided") {
    pval <- 2 * (1 - pt(abs(stat), df = n - 1))
  } else if(alternative == "greater") {
    pval <- 1 - pt(stat, df = n - 1)
  } else {
    stop("Unknown alternative")
  }
  
  return(pval)
}

Уграђена функција у R-у која спроводи \(t\)-тест је t.test(uzorak, mu) - ово су обавезни параметри.

Могуће је додати и алтернативне параметре

• conf.level за ниво поверења,

• alternative за облик алтернативне хипотезе. Можемо поставити на "greater" или "less", ако хоћемо такву алтернативну хипотезу.

Овај тест враћа:

• \(t\) вредност тест статистике,

• \(df\) број степени слободе Студентове расподеле тест статистике,

• \(p\) вредност теста на основу које закључујемо да ли прихватамо нулту хипотезу или не (ако је већа од нивоа значајности прихватамо \(H_0\))

• интервал поверења за \(m\), по default-у је 95%-ни, а ако хоћемо неки други интервал поверења то назначимо параметром conf.level

Упоредимо наше резултате и резултате уграђене функцје.

x <- rnorm(50, 1, 2)

pval_t_test(x, 1, "two.sided")

## [1] 0.05079599

t.test(x, mu=1, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  x
## t = 2.0024, df = 49, p-value = 0.0508
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
##  0.9980766 2.0699389
## sample estimates:
## mean of x 
##  1.534008

Примећујемо да уграђена функција даје више информација од саме \(p\) вредност, а \(p\) вредност се поклапа у нашој и уграђеној функцији.

Уграђена функција нам даје и интервал поверења, па можемо и то проверити:

confidence_interval(x, 0.95)

## [1] 0.9926355 2.0753800

Наравно, вредности се поклапају.

Погледајмо могуће алтернативе кроз уграђену функцију.

# H1: m < m_0
t.test(x, mu=1, alternative = "less")

## 
##  One Sample t-test
## 
## data:  x
## t = 2.0024, df = 49, p-value = 0.9746
## alternative hypothesis: true mean is less than 1
## 95 percent confidence interval:
##      -Inf 1.981125
## sample estimates:
## mean of x 
##  1.534008

Ако тестирамо \(H_1(m<m_0)\), добијамо врло високу \(p\) вредност, па ћемо прихватити нулту хипотезу \(H_0(m=1)\).

Приметимо да је у случају једностраног теста и интервал поверења једностран - лева граница му је \(-\infty\).

Ако тестирамо \(H_1(m>m_0)\) резултат је сличан.

t.test(x, mu=1, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  x
## t = 2.0024, df = 49, p-value = 0.0254
## alternative hypothesis: true mean is greater than 1
## 95 percent confidence interval:
##  1.08689     Inf
## sample estimates:
## mean of x 
##  1.534008

У наставку ћемо користити уграђене тестове у R-у и нећемо имплементирати своје.

Two Sample t-test - тест два узорка

Уколико имамо два независна узорка \(X_1, \dots X_{n_1}\) и \(Y_1,\dots, Y_{n_2}\), из расподела \(N(m_1,\sigma_1^2)\) и \(N(m_2,\sigma_2^2)\) редом, за тестирање хипотезе \(H_0(m_1=m_2)\) користимо тест статистику \[t=\frac{(\bar{X}_{n_1} - m_1)-(\bar{Y}_{n_2}-m_2)}{\sqrt{n_1\bar{S}^2_{n_1} + n_2\bar{S}^2_{n_2}}}\sqrt{\frac{n_1n_2}{n_1+n_2}(n_1+n_2-2)},\]

која има Студентову \(t_{n_1+n_2-2}\) расподелу при \(H_0\).

Алтернативне хипотезе могу бити облика \(H_1:\;(m_1<m_2)\), \(H_1:\;(m_1>m_2)\) или \(H_1:\;(m_1\neq m_2)\).

Овај тест се у R-у извршава додајући још један узорак у позив функције t.test. У овом случају параметри функције су:

• x - први вектор
• y - други вектор
• formula - ако нећемо да прослеђујемо одвојено векторе \(x\) и \(y\) него делимо базу података на два дела преко модел формуле: kolona1~kolona2 (где је kolona2 фактор са два нивоа)
• paired - логички параметар, TRUEако хоћемо упарени t-тест
• var.equal - логички параметар, TRUE је ако можемо да претпоставимо да су дисперзије два узорка једнаке (подразумева се да нису), што се проверава тестом var.test(uzorak1, uzorak2)

Овај тест враћа вредност тест статистике, број степени слободе, \(p\) вредност теста, интервал поверења за разлику средњих вредности, као и оцењене средње вредности за оба узорка.

x <- rnorm(50)
y <- rnorm(35, mean = 2, sd = 4)

t.test(x, y)

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -3.1198, df = 37.769, p-value = 0.003459
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.1914343 -0.6792784
## sample estimates:
##   mean of x   mean of y 
## -0.08989631  1.84546003

Овде добијамо малу \(p\) вредност, што указује на то да треба одбацити нулту хипотезу у корист алтернативне, која је по default-у \(H_1(m_1 \neq m_2)\), што и пише у излазу функције.

Ако бисмо тестирали са алтернативном хипотезом \(H_1(m_1>m_2)\)

t.test(x, y, alternative = "greater")

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -3.1198, df = 37.769, p-value = 0.9983
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -2.981392       Inf
## sample estimates:
##   mean of x   mean of y 
## -0.08989631  1.84546003

добијена је јако велика \(p\) вредност, па не бисмо могли да одбацимо нулту хипотезу. То не значи да је нулта хипотеза тачна, већ да на основу добијеног узорка не можемо одбацити нулту хипотезу у корист алтернативне.

# Varijanta kada uzorke vadimo kao podskupove neke baze podataka
library(ISwR)
data(energy)
attach(energy)
energy

##    expend stature
## 1    9.21   obese
## 2    7.53    lean
## 3    7.48    lean
## 4    8.08    lean
## 5    8.09    lean
## 6   10.15    lean
## 7    8.40    lean
## 8   10.88    lean
## 9    6.13    lean
## 10   7.90    lean
## 11  11.51   obese
## 12  12.79   obese
## 13   7.05    lean
## 14  11.85   obese
## 15   9.97   obese
## 16   7.48    lean
## 17   8.79   obese
## 18   9.69   obese
## 19   9.68   obese
## 20   7.58    lean
## 21   9.19   obese
## 22   8.11    lean

Видимо да се колона expend може поделити на два подскупа у зависности од тога да ли је у колони stature вредност OBESE или LEAN (фактор са два нивоа). Зато можемо користити \(t\) тест са формулом expend~stature. Прво проверавамо да ли су дисперзије једнаке.

var.test(expend~stature)

## 
##  F test to compare two variances
## 
## data:  expend by stature
## F = 0.78445, num df = 12, denom df = 8, p-value = 0.6797
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1867876 2.7547991
## sample estimates:
## ratio of variances 
##           0.784446

Видимо да је \(p\) вредност већа од 0.05, па прихватамо хипотезу о једнакости дисперзија.

t.test(expend~stature, var.equal=T)

## 
##  Two Sample t-test
## 
## data:  expend by stature
## t = -3.9456, df = 20, p-value = 0.000799
## alternative hypothesis: true difference in means between group lean and group obese is not equal to 0
## 95 percent confidence interval:
##  -3.411451 -1.051796
## sample estimates:
##  mean in group lean mean in group obese 
##            8.066154           10.297778

Видимо да је \(p\) вредност мања од нивоа значајности 0.05, па одбацујемо хипотезу о једнакости средњих вредности у ове две групе.

Наравно, исто се добија и овако

x<-expend[stature=="lean"]
y<-expend[stature=="obese"]
t.test(x,y,var.equal=T)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = -3.9456, df = 20, p-value = 0.000799
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.411451 -1.051796
## sample estimates:
## mean of x mean of y 
##  8.066154 10.297778

Задатак.

Мерен је горњи притисак на узорку од 12 мушкараца и добијено је 130, 148, 122, 140, 132, 142, 124, 150, 170, 136, 146, 140, а на узорку од 13 жена добијене су следеће вредности 140, 150, 130, 132, 150, 138, 123, 124, 160, 138, 170, 144, 108. Сматра се да горњи притисак и код мушкараца и код жена има нормалну расподелу. Са прагом значајности 0.1 тестирати хипотезу да су средње вредности притисака мушкараца и жена једнаке против алтернативе да се разликују.

m <- c(130, 148, 122, 140, 132, 142, 124, 150, 170, 136, 146, 140)
z <- c(140, 150, 130, 132, 150, 138, 123, 124, 160, 138, 170, 144, 108)

var.test(m, z, ratio = 1, alternative = "two.sided")

## 
##  F test to compare two variances
## 
## data:  m and z
## F = 0.62684, num df = 11, denom df = 12, p-value = 0.4475
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1887217 2.1498034
## sample estimates:
## ratio of variances 
##          0.6268356

t.test(m, z, var.equal=T)

## 
##  Two Sample t-test
## 
## data:  m and z
## t = 0.16762, df = 23, p-value = 0.8683
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -11.34121  13.34121
## sample estimates:
## mean of x mean of y 
##       140       139

Упарени тест

Ако обележја \(X\) и \(Y\) нису независна, већ имамо узорак парова \((X_1,Y_1), \dots, (X_n,Y_n)\) тестирање хипотезе \(H_0(m_1=m_2)\) се врши упареним \(t\) тестом који је у R-у имплементиран такође у функцији t.test, где се само дода аргумент paired = TRUE.

x <- rnorm(50, mean = 2)
y <- x + rnorm(50, sd = 0.1) # Y nije nezavisno od X nego je X + mali sum

t.test(x, y, paired = TRUE)

## 
##  Paired t-test
## 
## data:  x and y
## t = -0.58645, df = 49, p-value = 0.5603
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.03912975  0.02145071
## sample estimates:
## mean difference 
##     -0.00883952

Као резултат имамо велику \(p\) вредност и не одбацујемо нулту хипотезу.

Непараметарски тестови

Вилкоксонов тест заснован на ранговима и знаковима

Вилкоксонов тест има сличну сврху као \(t\) тест, али претпоставља да је расподела обележја \(X\) симетрична и тестира \(H_0 (m=m_0)\), без претпоставке нормалности. Користи се тест статистика \[T=\sum_{i=1}^n r_i I\{ X_i-m_0\ge 0\},\] где су \(r_i\) рангови елемената \(|X_i-m_0|\) у узорку \(|X_1-m_0|, \dots, |X_n-m_0|\).

У R-у овај тест је имплементиран у функцији wilcox.test. Ова функција има исти интерфејс као функција t.test (аргументи и слично).

Проверимо тест неким узорком из биномне расподеле \(B(10,0.3)\). Видимо да је очекивање 3.

x <- rbinom(50, 10, 0.3)

wilcox.test(x, mu = 3)

## Warning in wilcox.test.default(x, mu = 3): cannot compute exact p-value with
## ties

## Warning in wilcox.test.default(x, mu = 3): cannot compute exact p-value with
## zeroes

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  x
## V = 398, p-value = 0.8744
## alternative hypothesis: true location is not equal to 3

Применимо овај тест и на неки узорак из \(t\) расподеле за алтернативну хипотезу \(H_1(m>0)\).

x <- rt(20, df = 2)
wilcox.test(x, mu = 0, alternative = "greater")

## 
##  Wilcoxon signed rank exact test
## 
## data:  x
## V = 63, p-value = 0.943
## alternative hypothesis: true location is greater than 0

Дакле, уколико имамо претпоставку о нормалној расподели, користимо t.test, а уколико немамо wilcox.test може бити задовољавајући.

Као и t.test и wilcox.test се може применити на два узорка (упарени и неупарени тест).

x <- rexp(50)
y <- rexp(30)
wilcox.test(x, y)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x and y
## W = 696, p-value = 0.5949
## alternative hypothesis: true location shift is not equal to 0

Колмогоров-Смирновљев тест сагласности са расподелом

Тест Колмогоров-Смирнова служи за проверу да ли неки узорак \(X_1,\dots ,X_n\)одговара расподели са функцијом расподеле \(F_0\). Заснован је на тест статистици \[T=\sup_{x} |F_n(x)-F_0(x)|,\]

где је \(F_n\) емпиријска функција расподеле узорка.

У R-у је имплементиран кроз функцију ks.test, а као аргументе прима узорак, као и функцију расподеле (функције обично почињу са p*, pnorm, pexp,…)

Тестирамо да ли колона speed из базе података cars има стандардну нормалну расподелу.

x <- cars$speed
ks.test(x, "pnorm")

## Warning in ks.test.default(x, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.99997, p-value < 2.2e-16
## alternative hypothesis: two-sided

Видимо да је \(p\) вредност практично нула, па одбацујемо хипотезу која каже да узорак има нормалну \(N(0,1)\) расподелу.

Можемо да тестирамо да ли има нормалну \(N(15,25)\) расподелу.

ks.test(cars$speed, function(x) pnorm(x, 15, 5))

## Warning in ks.test.default(cars$speed, function(x) pnorm(x, 15, 5)): ties should
## not be present for the Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  cars$speed
## D = 0.10575, p-value = 0.631
## alternative hypothesis: two-sided

У овом случају \(p\) вредност је 0.6, што је веће од \(\alpha\), па не одбацујемо хипотезу која каже да је узорак из нормалне \(N(15,25)\) расподеле.

set.seed(27)
x <- rnorm(50)
grid <- seq(-4, 4, length.out = 100)

plot(grid, pnorm(grid, mean = 0, sd = 1), type = "l",
     ylim = c(0, 1), ylab = "", lwd = 2, col = "red", xlab = "", main = "ECDFs")
lines(ecdf(x), verticals = TRUE, do.points = FALSE, col.01line = NULL)

# Does 'x' come from a standard normal distribution?
ks.test(x, "pnorm", mean = 0, sd = 1) 

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.15121, p-value = 0.1833
## alternative hypothesis: two-sided

set.seed(27)
y <- rnorm(50)

grid <- seq(-5, 5, length.out = 100)

plot(grid, punif(grid, min = -4, max = 4), type = "l",
     ylim = c(0, 1), ylab = "", lwd = 2, col = "red", xlab = "", main = "ECDFs")
lines(ecdf(y), verticals = TRUE, do.points = FALSE, col.01line = NULL)

# Does 'y' come from a uniform distribution in the interval (-4, 4)?
ks.test(y, "punif", min = -4, max = 4) 

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.23968, p-value = 0.005178
## alternative hypothesis: two-sided

И овај тест се може применити на тестирање о сагласности расподеле два узорка, тј. може да тестира да ли два узорка имају исту расподелу.

На пример, ако имамо два узорка из исте нормалне расподеле, очекујемо велику \(p\) вредност.

x <- rnorm(50)
y <- rnorm(40)

plot(ecdf(x), verticals = TRUE, do.points = FALSE, col.01line = NULL, xlab = "", main = "ECDFs")
lines(ecdf(y), verticals = TRUE, do.points = FALSE, col.01line = NULL, col = 4)

ks.test(x, y)

## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.235, p-value = 0.1477
## alternative hypothesis: two-sided

А ако имамо узорке из различитих расподела очекујемо малу \(p\) вредност.

x <- rnorm(50)
y <- rexp(40)
ks.test(x, y)

## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.52, p-value = 5.31e-06
## alternative hypothesis: two-sided

\(\chi^2\) тест незавиности

Користи се за тестирање независности два обележја \(X\) и \(Y\), а заснован је на тест статистици \[T= \sum_{i,j} \frac{(M_{ij}-n\widehat{p}_{ij})^2}{n\widehat{p}_{ij}}.\]

Користимо податке survey из пакета MASS.

У овом скупу постоје променљиве Smoke и Exer које говоре о томе да ли је студент пушач или непушач и у којој мери, као и о учесталости бављења физичком активношћу. Табелу контингенције добијамо позивом функције table.

library(MASS)
table(survey$Smoke, survey$Exer) 

##        
##         Freq None Some
##   Heavy    7    1    3
##   Never   87   18   84
##   Occas   12    3    4
##   Regul    9    1    7

Испитајмо да ли постоји зависност између чињенице да је студент пушач и нивоа физичке активности. За тестирање нулте хипотезе да су ова два обележја независна, можемо да користимо функцију chisq.test.

chisq.test(survey$Smoke, survey$Exer)

## Warning in chisq.test(survey$Smoke, survey$Exer): Chi-squared approximation may
## be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  survey$Smoke and survey$Exer
## X-squared = 5.4885, df = 6, p-value = 0.4828

Видимо да је \(p\) вредност 0.49 што нам указује да не можемо да одбацимо нулту хипотезу о независности.

Функцији chisq.test можемо да проследимо и табелу (или матрицу) са подацима коју R схвата као табелу контингенције

library(ISwR)
data(juul)
attach(juul)

## The following object is masked from package:MASS:
## 
##     menarche

head(juul)

##    age menarche sex igf1 tanner testvol
## 1   NA       NA  NA   90     NA      NA
## 2   NA       NA  NA   88     NA      NA
## 3   NA       NA  NA  164     NA      NA
## 4   NA       NA  NA  166     NA      NA
## 5   NA       NA  NA  131     NA      NA
## 6 0.17       NA   1  101      1      NA

# Ispitujemo nezavisnost kolona tanner i sex:
chisq.test(tanner,sex) # ovako je kada prosledjujemo vektore odvojeno

## 
##  Pearson's Chi-squared test
## 
## data:  tanner and sex
## X-squared = 28.867, df = 4, p-value = 8.318e-06

# Drugi nacin je da sami napravimo tabelu i prosledimo je
tabela<-table(tanner,sex)
tabela

##       sex
## tanner   1   2
##      1 291 224
##      2  55  48
##      3  34  38
##      4  41  40
##      5 124 204

chisq.test(tabela)

## 
##  Pearson's Chi-squared test
## 
## data:  tabela
## X-squared = 28.867, df = 4, p-value = 8.318e-06

Приметимо да је \(p\) вреднсот мала, па одбацујемо нулту хипотезу.

Задатак за вежбу

1. Учитати базу \(vit2005\) (пакет \(PASWR\)) и испитати њен садржај. Применити \(\chi^2\)-квадрат тест независности за податке из табеле. Да ли се хипотеза о независности прихвата за ниво значајности \(\alpha=0.01\)?