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

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

Интервали поверења за средњу вредност, дисперзија непозната

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

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

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

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

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

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

confidence_interval <- function(x, alfa) {
  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),mean(x) + C * sd(x) / sqrt(n))
}

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

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

## [1] -0.5762474  0.0442271

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

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

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

##            [,1]       [,2]       [,3]       [,4]       [,5]
## [1,] -0.1784162 -0.2659102 -0.3496114 -0.2783862 -0.3566179
## [2,]  0.3241724  0.3655665  0.2506831  0.2109375  0.1952336

Као резултат добијамо матрицу са 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.9485

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

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

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

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.9473

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

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

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

где је \(S_n^2\) (поправљена) узорачка дисперзија.

Имајући у виду ову расподелу, интервал поверења нивоа \(\alpha\) за непознати параметар \(\sigma^2\) се лако изводи и добија се да је једнак \[\left(\frac{(n-1)S_n^2}{c_2},\frac{(n-1)S_n^2}{c_1}\right)\]

где су \(c_1=F^{-1}_{\chi^2_{n-1}}\left(\frac{1-\alpha}{2}\right)\) и \(c_2=F^{-1}_{\chi^2_{n-1}}\left(\frac{1+\alpha}{2}\right)\).

Нека је дат узорак из нормалне расподеле обима 20 за који је израчунато да је узорачка дисперзија 21.12. Одредити \(99%\)-тни двострани интервал поверења за непознати параметар \(\sigma^2\).

n <- 20 # obim uzorka
sn2 <- 21.12 # popravljena uzoracka disperzija
alfa <- 0.99 # nivo poverenja
c1 <- qchisq( (1-alfa)/2 , 19)
c2 <- qchisq( (1+alfa)/2 , 19)
interval_poverenja <- c((n-1)*sn2/c2, (n-1)*sn2/c1)
interval_poverenja

## [1] 10.40064 58.63262

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

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

За ово тестирање можемо да користимо тест статистику: \[\frac{\bar{X}-m_0}{\sigma_0}\sqrt{n} \sim \mathcal{N}(0,1).\]

Претпоставимо да имамо узорак из нормалне \(\mathcal{N}(m,10)\) расподеле, где је \(m\) непознати параметар. Са нивоом значајности \(0.95\) испитати да ли је средња вредност узорка већа од 2.

uzorak <- c(5.6884562, 2.0308148, -0.2508697, -2.5062107, 5.1422191, 1.1849337, 3.6343267,
 4.6528487, 1.5760364, 5.1431170, 1.8864374, 4.4475159, 3.5619849, 4.7377123,
 -1.2994996, 2.4134766, 0.2744983, 5.2431267, -0.1152874, 0.8769407, 2.1222095,
 1.8561701, 3.0470218, 0.4738757, 2.1079643, 2.5249362, -2.2613439, 1.0116844,
 2.1796532, 1.6629011, 0.6026096, 3.0568876, 0.1365397, -5.0585758, 2.3148987,
 4.0924715, -0.9189906, 7.6984261, -1.9539183, 2.7631545, 5.5492733, -1.5838495,
 -2.3246922, 0.7752239, 1.3790237, -1.1731306, 4.8165865, 2.9824243, 6.0295380,
 0.1998773, 9.0725547, -1.2977859, 1.9417684, -3.5287460, 2.5369642, -1.3349065,
 -1.7820453, 4.2251318, 2.9737008, 1.8158107)

# install.packages("BSDA")
library(BSDA)

## Warning: package 'BSDA' was built under R version 4.1.3

## Loading required package: lattice

## 
## Attaching package: 'BSDA'

## The following object is masked from 'package:datasets':
## 
##     Orange

z.test(uzorak, sigma.x = sqrt(10), mu = 2, alternative = "greater")

## 
##  One-sample z-Test
## 
## data:  uzorak
## z = -0.52852, p-value = 0.7014
## alternative hypothesis: true mean is greater than 2
## 95 percent confidence interval:
##  1.112723       NA
## sample estimates:
## mean of x 
##  1.784231

Испитати да ли је дисперзија узорка \(x\) мања од 8.

# instalirati prvo paket EnvStats
library(EnvStats)

## Warning: package 'EnvStats' was built under R version 4.1.3

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

x <- c(5.6884562, 2.0308148, -0.2508697, -2.5062107, 5.1422191, 1.1849337, 3.6343267,
 4.6528487, 1.5760364, 5.1431170, 1.8864374, 4.4475159, 3.5619849, 4.7377123,
 -1.2994996, 2.4134766, 0.2744983, 5.2431267, -0.1152874, 0.8769407, 2.1222095,
 1.8561701, 3.0470218, 0.4738757, 2.1079643, 2.5249362, -2.2613439, 1.0116844,
 2.1796532, 1.6629011, 0.6026096, 3.0568876, 0.1365397, -5.0585758, 2.3148987,
 4.0924715, -0.9189906, 7.6984261, -1.9539183, 2.7631545, 5.5492733, -1.5838495,
 -2.3246922, 0.7752239, 1.3790237, -1.1731306, 4.8165865, 2.9824243, 6.0295380,
 0.1998773, 9.0725547, -1.2977859, 1.9417684, -3.5287460, 2.5369642, -1.3349065,
 -1.7820453, 4.2251318, 2.9737008, 1.8158107)
varTest(x, alternative = "less", sigma.squared = 8)

## $statistic
## Chi-Squared 
##    57.67336 
## 
## $parameters
## df 
## 59 
## 
## $p.value
## [1] 0.4754733
## 
## $estimate
## variance 
## 7.820117 
## 
## $null.value
## variance 
##        8 
## 
## $alternative
## [1] "less"
## 
## $method
## [1] "Chi-Squared Test on Variance"
## 
## $data.name
## [1] "x"
## 
## $conf.int
##      LCL      UCL 
##  0.00000 10.89737 
## attr(,"conf.level")
## [1] 0.95
## 
## attr(,"class")
## [1] "htestEnvStats"

Студентов \(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.3594977

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

## 
##  One Sample t-test
## 
## data:  x
## t = 0.925, df = 49, p-value = 0.3595
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
##  0.7287182 1.7340199
## sample estimates:
## mean of x 
##  1.231369

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

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

confidence_interval(x, 0.95)

## [1] 0.7287182 1.7340199

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

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

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

## 
##  One Sample t-test
## 
## data:  x
## t = 0.925, df = 49, p-value = 0.8203
## alternative hypothesis: true mean is less than 1
## 95 percent confidence interval:
##      -Inf 1.650721
## sample estimates:
## mean of x 
##  1.231369

Ако тестирамо \(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 = 0.925, df = 49, p-value = 0.1797
## alternative hypothesis: true mean is greater than 1
## 95 percent confidence interval:
##  0.8120169       Inf
## sample estimates:
## mean of x 
##  1.231369

У наставку ћемо користити уграђене тестове у 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}-\bar{Y}_{n_2}}{\sqrt{\frac{\bar{S}^2_{n_1}}{n_1}+\frac{\bar{S}^2_{n_2}}{n_2}}},\]

која има Студентову \(t_{n_1+n_2-1}\) расподелу при \(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.3405, df = 37.235, p-value = 0.001911
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.8461314 -0.9423001
## sample estimates:
##   mean of x   mean of y 
## -0.04642389  2.34779186

Овде добијамо малу \(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.3405, df = 37.235, p-value = 0.999
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -3.603204       Inf
## sample estimates:
##   mean of x   mean of y 
## -0.04642389  2.34779186

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

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

## Warning: package 'ISwR' was built under R version 4.1.3

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

Упарени тест

Ако обележја \(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 = 1.8317, df = 49, p-value = 0.07308
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.002245168  0.048474191
## sample estimates:
## mean of the differences 
##              0.02311451

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

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

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

Тест Колмогоров-Смирнова служи за проверу да ли неки узорак \(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(x, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test

## 
##  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(cars$speed, function(x) pnorm(x, 15, 5)): ties should not be
## present for the Kolmogorov-Smirnov test

## 
##  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)\) расподеле.

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

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

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

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.15, p-value = 0.6473
## alternative hypothesis: two-sided

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

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

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.5, p-value = 1.448e-05
## 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)

## Warning: package 'MASS' was built under R version 4.1.3

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:EnvStats':
## 
##     boxcox

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\) вреднсот мала, па одбацујемо нулту хипотезу.