Практикум по Вероятности и Статистика с R - Записки

Точкови оценки

Ако имаме наблюдения \(x_1, ..., x_n\) над случайна величина \(X \sim \mathcal{N}(\mu, \sigma^2)\), то можем да мислим на всяко наблюдение \(x_i\) като на наблюдение от случайна величина \(X_i\) разпределена като \(X\) и \(X_i\) са независими. Тоест, имаме независими и еднакво разпределени (н.е.р.) случайни величини \(X_1, ..., X_n \overset{distr}{\sim} X \sim \mathcal{N}(\mu, \sigma^2)\), като сме получили по едно наблюдение от всяка. Векторът \((X_1,...,X_n)\) наричаме случайна извадка.

Тогава \(\bar{x} := \frac{1}{n}\sum_{i=1}^n{x_i}\) е наблюдение от случайната величина \(\bar{X} := \frac{1}{n}\sum_{i=1}^n{X_i} \sim \mathcal{N}(\mu, \frac{\sigma^2}{n})\), понеже \[ E[\bar{X}] = E[\frac{1}{n}\sum_{i=1}^n{X_i}]=\frac{1}{n}\sum_{i=1}^n{E[X_i]}=\frac{1}{n}n\mu = \mu\\ Var[\bar{X}]=Var[\frac{1}{n}\sum_{i=1}^n{X_i}] = \frac{1}{n^2}n\sigma^2=\frac{\sigma^2}{n} \]

Забелязваме, че \(Var[\bar{X}] \xrightarrow[n\to \infty]{} 0\). Тоест, колкото повече наблюдения имаме, толкова повече \(\bar{x}\) има вероятност да е близо до истинския параметър \(\mu\), понеже ще е наблюдение от нормално разпределена случайна величина с дисперсия \(\frac{\sigma^2}{n}\) и очакване \(\mu\).

Можем да наблюдаваме този факт със следната симулация

# library("tidyverse")
# library("ggplot2")

# брой извадки които взимаме
n <- 1000

# 1000 пъти взимаме средно на 10, 100 и 1000 наблюдения от N(25, 40^2)
c(10, 100, 1000) %>%
  map(function(size) {
    # Взимаме средното на size на брой наблюдения от N(25, 40^2)
    map(1:n, ~mean(rnorm(size, mean = 25, sd = 40))) %>%
      unlist() %>%
      data.frame(x = ., size = as.character(size))
  }) %>%
  bind_rows() %>%
  ggplot() +
  geom_density(mapping = aes(x = x, fill = size)) +
  facet_wrap(~size) +
  ggtitle("Плътност на извадково средно при n=10, 100, 1000") +
  xlab("извадково средно") +
  ylab("емпирична плътност на извадково средно") +
    # theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  scale_fill_pomological() + 
  scale_color_pomological()

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Още една визуализация

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

## Picking joint bandwidth of 0.721

## Warning: This manual palette can handle a maximum of 9 values. You have supplied
## 10.

\(\bar{X}\) наричаме точкова оценка за параметъра \(\mu\). Това, че \(Е(\bar{X}) = \mu\) наричаме неизместеност на оценката. Когато една точкова оценка \(f_{\theta}(\overrightarrow{X})\) приближава по вероятност своя параметър \(\theta\), тоест \[ \lim_{n\to \infty}Pr(|f_{\theta}(\overrightarrow{X}) - \theta| > \epsilon) = 0 \text{ за } \forall \epsilon > 0 \] казваме, че точковата оценка \(f_{\theta}(\overrightarrow{X})\) е асимптотично консистента.

Стандартното отклонение на оценката, обикновено наричаме стандартна грешка.

Интервални Оценки

Ако имаме наблюдения над случайна величина \(X \sim \mathcal{N}(\mu^?, \sigma^2)\) за която \(\mu^?\) е неизвестен параметър, но \(\sigma^2\) е известна (в практиката рядко ще знаем който и да е параметър), то знаем, че

\[ Z := \frac{\bar{X} - \mu^?}{\sigma/\sqrt{n}} \sim \mathcal{N}(0, 1) \]

и следователно, можем да използваме определен квантил q от \(\mathcal{N}(0, 1)\), така че

\[ P(-q < \frac{\bar{X} - \mu^?}{\sigma/\sqrt{n}} < q) = \gamma \]

Така, с вероятност \(\gamma\) (ниво на достоверност) можем да сме сигурни, че

\[ \bar{X} - q \frac{\sigma}{\sqrt{n}} < \mu^? < \bar{X} + q \frac{\sigma}{\sqrt{n}} \]

В R за 95% интервал на достоверност, това ще стане по следния начин

# брой елементи
n <- 1000

# сигма ни е известно, но мю не ни е и се опитваме да го оценим
known_sd <- 10
x <- rnorm(n, 14.55, known_sd)

emp_mean <- mean(x)
q <- qnorm(0.975)

c(emp_mean - q * known_sd / sqrt(n), emp_mean + q * known_sd / sqrt(n))

## [1] 13.60031 14.83990

В повечето случаи обаче, \(\sigma\) не ни е известно. Тогава, ако заместим с извадковото стандартно отклонение S, случайната величина \(\frac{\bar{X} - \mu}{S/\sqrt{n}}\) има \(\mathcal{T}\) разпределение с \(n - 1\) степени на свобода. В този случай, трябва просто да използваме квантилите от \(\mathcal{T}\) разпределение.

n <- 1000
x <- rnorm(n, 25, 7.654)

emp_mean <- mean(x)
emp_sd <- sd(x)

q <- qt(0.975, df = n - 1)

c(emp_mean - q * emp_sd / sqrt(n), emp_mean + q * emp_sd / sqrt(n))

## [1] 24.59308 25.54813

Или с t.test функцията

t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 103.03, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  24.59308 25.54813
## sample estimates:
## mean of x 
##  25.07061

Получената оценка за \(\mu\) се нарича интервална оценка. Ето как изглеждат интервалите на достоверност при Нормално и \(\mathcal{T}\) разпределение за случайната величина Z дефиниране по-горе.

# library("ggridges")

n <- 100000

list("normal" = rnorm(n), "t" = rt(n, 100000 - 1)) %>%
  imap(~data.frame(x = .x, distribution = .y)) %>%
  bind_rows() %>%
  ggplot(aes(x = x, y = distribution, fill = factor(stat(quantile)))) +
  stat_density_ridges(
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    quantiles = c(0.025, 0.975)
  ) +
  scale_fill_manual(
    name = "Probability", values = c("#FF0000A0", "#62b5b2", "#FF0000A0"),
    labels = c("(0, 0.025]", "(0.025, 0.975]", "(0.975, 1]")
  ) +
  xlab("Z") +
  theme_minimal()

## Picking joint bandwidth of 0.09

Нека симулираме различни извадки от нормално разпределена случайна величина с \(\mu=3\) и \(\sigma^2 = 2\) и да видим в колко от случаите ще попадаме в интервала на достоверност.

simulate_t <- function(n, mu, sigma) {
  x <- rnorm(n, mu, sigma)
  b1 <- mean(x) - qt(0.975, df = n - 1) * (sd(x) / sqrt(n))
  b2 <- mean(x) + qt(0.975, df = n - 1) * (sd(x) / sqrt(n))
  between(mu, b1, b2)
}

# да преброим колко елементи са в интервала на достоверност
simulations <- replicate(10000, simulate_t(40, 4, 4))
mean(simulations)

## [1] 0.9538

Виждме, че наистина в около 95% от случаите попадаме в интервала на достоверност.

Тестване на хипотези върху една извадка

Хипотеза наричаме какво да е твърдение за функцията на разпределение \(F_X\) на случайната величина \(Х\) върху която имаме извадка \((X_1, ..., X_n)\).

Нека наблюдаваме данни идващи от нормално разпределение, примерно от \(\mathcal{N}(176, 6^2)\) за които не знаем истинската стойност на \(\mu\) (за момент забравяме, че сме ги генерирали с \(\mu=176\))

heights <- rnorm(100, 176, 6)

Ако примем, че истинското математическо очакване \(\mu\) е по-голямо или равно на 190, то случайната величина

\[ Z := \frac{\bar{X} - 190}{\sigma / \sqrt{n}} \] ще бъде нормално разпределена с \(\mu=0\) и \(\sigma^2=1\)

# z-statistic
z <- (mean(heights) - 190) / (6 / sqrt(100))
z

## [1] -24.37082

Можем да видим вероятността да имаме стойност z или по-малка използвайки функцията pnorm

pvalue <- pnorm(z)
# Вероятността да сме наблюдавали z или по-крайна стойност на z
pvalue

## [1] 1.744098e-131

Виждаме, че ако хипотезата, че \(\mu \ge 190\) e вярна, то е изключително малко вероятно да наблюдаваме \(z\)-статистиката, която имаме (или по-крайна). Тогава, можем да отхвърлим хипотезата, че \(\mu \ge 190\). Но колко трябва да е малко нашето pvalue, така че да отхвърлим хипотезета \(H_0\) че \(\mu \ge 190\) и да приемем алтернативата \(H_1\) че \(\mu < 190\). Ако изберем някакво \(\alpha\), примерно \(\alpha = 0.05\) (ниво на значимост) можем да отхвърляме хипотезата \(H_0\) само когато \(pvalue < \alpha := 0.05\) и вероятността \(P(\text{Отхвърляне}|H_0=True):= \alpha\) да отхвърлим хипотезата \(H_0\), когато тя е вярна е 0.05 (Грешка от първи род). Областта \(\mathcal{C}\) в която избираме да отхвърлим нулевата хипотеза, когато вектора на наблюденията \((x_1, ..., x_n)\) попадне в нея, наричаме критична област. В нашия случай, критичната област е

\[ \mathcal{C}:= \{(x_1, ..., x_n): \frac{\bar{x}-\mu_{H}}{\sigma/\sqrt{n}} < q^{\mathcal{N}(0,1)}_{\alpha}\} \]

Ако не успеем да отхвърлим хипотезата \(H_0\), а тя се окаже грешна, правиме грешка от втори род. Нейната вероятност \(P(\text{Приемане}|H_0=False)\) бележим с \(\beta\). Сила на теста \(\pi:= 1-\beta\) наричаме вероятността да отхвърлим \(H_0\) при условие, че \(H_0\) наистина е грешна и е вярна \(H_1\).

Вярно	Приемаме \(H_0\)	Отхвърляме \(H_0\)
\(H_0\)	Вярно	Грешка от 1ви род с вероятност \(\alpha\)
\(H_1\)	Грешка от 2ри род с вероятност \(\beta\)	Вярно със сила \(\pi := 1-\beta\)

---

Когато \(\sigma\) също е неизвестна, то

\[ Т := \frac{\bar{X} - 190}{S / \sqrt{n}} \sim \mathcal{T}(n-1) \] и за да видим вероятността да наблюдаваме такава стойност или по-крайна (при услувие, че нулевата хипотеза е вярна), трябва да използваме функцията pt и да заместим \(\sigma\) с \(S\) във формулата за \(Z\)

t <- (mean(heights) - 190) / (sd(heights) / sqrt(100))
pvalue <- pt(t, length(heights) - 1)

# вероятността да сме наблюдавали z или по-крайна стойност от z, при условие, че
# нулевата хипотеза е вярна.
pvalue

## [1] 3.131185e-44

Отново, можем да използваме функцията t.test, която да ни даде цялата информация.

t.test(heights, mu = 190, alternative = "less")

## 
##  One Sample t-test
## 
## data:  heights
## t = -24.575, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is less than 190
## 95 percent confidence interval:
##      -Inf 176.3655
## sample estimates:
## mean of x 
##  175.3775

---

За разлика от подхода при интервални оценки, при тестване на хипотези, центрираме тестовата статистика \(Z\) около хипотетичен параметър и отхвърляме нулевата хипотеза \(H_0\) когато наблюдаваната (или по-крайна) тест статистика \(Z\) се случва с много малък шанс.

---

• Когато тестваме \(H_0:= \mu \ge \mu_{хипотеза}\) срещу \(H_1 := \mu < \mu_{хипотеза}\) правим едностранен тест с лява опашка. В този случай \(pvalue := P(Z \le t |H_0)\), където \(t\) е стойността която сме наблюдавали за тестовата статистика \(T\).

## Picking joint bandwidth of 0.0899

• Когато тестваме \(H_0:= \mu \le \mu_{хипотеза}\) срещу \(H_1 := \mu > \mu_{хипотеза}\) правим едностранен тест с дясна опашка. В този случай \(pvalue := P(T \ge z |H_0)\), където \(t\) е стойността която сме наблюдавали за тестовата статистика \(T\).

# Тестване на хипотезата, че E(X) < 130
t.test(heights, mu = 130, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  heights
## t = 76.263, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is greater than 130
## 95 percent confidence interval:
##  174.3896      Inf
## sample estimates:
## mean of x 
##  175.3775

## Picking joint bandwidth of 0.0901

---

• Когато тестваме \(H_0:= \mu = \mu_{хипотеза}\) срещу \(H_1 := \mu \neq \mu_{хипотеза}\) правим двустранен тест с две опашки. В този случай ще трябва да разделим нивото на значимост \(\alpha\) на две за двете опашки така, че да правим грешка от първи род с \(\alpha\) вероятност. При тест с две опашки, дефинираме \(pvalue := 2\min(P(Z \ge t|H_0), P(Z \le t|H_0))\). Aко разпределението е симетрично около нулата, \(pvalue := P(abs(T) \ge abs(t) |H_0)\), където отново, \(t\) е стойността която сме наблюдавали за тестовата статистика \(T\).

# Тестване на хипотезата, че E(X) = 165
t.test(heights, mu = 130, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  heights
## t = 76.263, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 130
## 95 percent confidence interval:
##  174.1969 176.5581
## sample estimates:
## mean of x 
##  175.3775

Ако не искаме да правим грешка от първи род, можем да намалим областта на \(\alpha\), но така ще отхвърляме прекалено лесно и ще качим вероятността за грешка от втори род.

Нека симулираме \(T\)-статистиката при вярна нулева хипотеза.

simulate_t <- function(n, th_mean, th_std, hyp_mean) {
  x <- rnorm(n, th_mean, th_std)
  (mean(x) - hyp_mean) / (sd(x) / sqrt(n))
}

# да преброим колко елементи са в критичната област
simulations <- replicate(10000, simulate_t(40, 4, 4, 4))
q <- qt(p = 0.05, df = 39)

mean(simulations <= q)

## [1] 0.0488

Виждаме, че дори при вярна нулева хипотеза, в 1 на 20 случаи ще попадаме в критичната област. Тогава, ако попаднем в тази област при вярна нулева хипотеза, ще отхвърлим нулевата хипотеза като грешна и ще направим грешка от първи род.

Нека симулираме данни от две различни разпределение ; \(X_1 \sim \cal{N(\mu_1, \sigma^2)}\) и \(X_2 \sim \cal{N(\mu_2, \sigma^2)}\) с цел да видим как се променят грешките от втори род \(\beta\) при увеличаване на броя на елементите \(n\) с които изграждаме тестовата статистика.

simulate_t <- function (n, m, std = 10) {
  observations <- rnorm(n, m, std)
  mean(observations) / (sd(observations) / sqrt(n))
}

simulate_ht <- function(n, n_sim, mean1, mean2, std = 10) {
  set.seed(1)

  bind_rows(
    data.frame(
      x = replicate(n_sim, simulate_t(n, mean1, std)),
      distribution = paste0("средно=", mean1),
      n = as.factor(n)
    ),
    data.frame(
      x = replicate(n_sim, simulate_t(n, mean2, std)),
      distribution = paste0("средно=", mean2),
      n = as.factor(n)
    )
  )
}

c(5, 10, 20, 30, 70, 140) %>%
  imap(~simulate_ht(.x, 10000, 5, 8, 6)) %>%
  bind_rows() %>%
  ggplot(mapping = aes(x = x, fill = distribution)) +
  geom_density(alpha = 0.55) +
  xlim(c(-5, 20)) +
  facet_wrap(~n) +
  xlab("Т-статистика") +
  ylab("Вероятностни плътности") +
  ggtitle("Промяна на силата при увеличаване на броя на елементи") +
  theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  scale_fill_pomological() + 
  scale_color_pomological() +
  theme(legend.position = "bottom")

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

## Warning: Removed 48 rows containing non-finite values (stat_density).

Виждаме, че при растящ брой елементи \(n\), областта в която бихме могли да направим грешка от втори род намалява, понеже дисперсията на T-статистиката под нулевата хипотеза намалява и разпределението се центрира около хипотетичния си параметър. По този начин, \(\beta\) намалява, докато силата на теста \(\pi:= 1-\beta\) се увеличава.

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Забелязваме, че силата на теста \(\pi\) расте при увеличаване на броя на наблюденията, но при по-голямо стандартно отклонение, скоростта на растеж на \(\pi\) намалява.

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

При качване на нивото на значимост \(\alpha\), качваме силата на теста (\(\pi\)) и намаляваме вероятността за грешка от втори род \(\beta\), но по този начин увеличаваме вероятността за грешка от първи род \(\alpha\).

---

Тестване на хипотези върху две извадки

Ако имаме еднакъв брой наблюдения от две случайни величини \(X\) и \(Y\), които са нормално разпределени с еднакви дисперсии - съответно, \(X \sim \mathcal{N}(\mu_1, \sigma^2)\) и \(Y \sim \mathcal{N}(\mu_2, \sigma^2)\), можем да зададем нулева хипотеза \(H_0 := \mu_1 = \mu_2\). Ако нулевата хипотеза е вярна, то \(\bar{X} - \bar{Y} \sim \mathcal{N}(0, \frac{2\sigma^2}{n})\). Тогава

\[ Т := \frac{\bar{X}-\bar{Y}}{\sqrt{\frac{S_1^2}{n} + \frac{S_2^2}{n}}} \sim \mathcal{T}(2n - 2) \]

Тогава, можем да използваме досегашния подход при тестване на хипотези за да проверим дали наистина \(\mu_1=\mu_2\).

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Aко имаме \(pvalue < \alpha\), където \(\alpha\) е нивото на значимост, а t e наблюдаваната стойност за T-статистиката, можем да отхвърлим нулевата хипотеза, че \(\mu_1=\mu_2\).

Пример:

# брой елементи в двете случайни извадки
n <- 40

# Симулации на случайни извадки от с.в. X ~ N(3, 9) и Y ~ N(4, 9)
x <- rnorm(n, 3, 3)
y <- rnorm(n, 4, 3)

# Извадкови вариации
v1 <- var(x)
v2 <- var(y)

# Степени на свобода за T-разпределението
degrees_of_freedom <- 2 * n - 2

# Наблюдаваната стойност от T-статистиката
t <- (mean(x) - mean(y)) / sqrt(v1 / 40 + v2 / 40)

# Изчисление за pvalue при двустранен тест
pval <- 2 * min(
  pt(q = t, df = degrees_of_freedom, lower.tail = TRUE),
  pt(q = t, df = degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pval, "df" = degrees_of_freedom)

##          t     pvalue         df 
## -1.9004828  0.0610646 78.0000000

Или ако просто използваме t.test

t.test(x, y, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = -1.9005, df = 78, p-value = 0.06106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.67067950  0.06201836
## sample estimates:
## mean of x mean of y 
##  2.835611  4.139941

---

Когато работим с предположението, че дисперсиите са еднакви, но имаме различен брой наблюдения от двете случайни извадки, за да оценим по-точно общата дисперсия посредством извадъчните дисперсии, използваме формулата за комбинирана извадъчна дисперсия и съответно - фoрмулата за комбинирано стандартно отклонение.

\[ S_p^2 = \frac{(n - 1)S_1^2 + (m - 1)S_2^2}{n + m - 2} \]

където \(n\) и \(m\) са броят на елементите в двете извадки, а \(S_i\) са извадъчните дисперсии на двете извадки. В този случай T-статистиката е

\[ Т := \frac{\bar{X}-\bar{Y}}{S_p\sqrt{\frac{1}{n} + \frac{1}{m}}} \sim \mathcal{T}(n + m - 2) \]

Пример:

# Формула за смесено извадъчно стандартно отклонение
pooled_sd <- function (x, y) {
  n1 <- length(x)
  n2 <- length(y)
  sqrt(
    ((n1 - 1) * var(x) + (n2 - 1) * var(y))
    /
    (n1 + n2 - 2)
  )
}

# броя на елементите в двете случайни извадки които ще симулираме
n <- 30
m <- 40

# Симулации с различен брой над с.в.
x <- rnorm(n, 3, 2)
y <- rnorm(m, 4, 2)

# Степени на свобода
degrees_of_freedom <- n + m - 2

# Наблюдавана стойност t на Т-статистиката
t <- (mean(x) - mean(y)) / (pooled_sd(x, y) * sqrt((1 / n) + (1 / m)))

# Изчисление за pvalue при двустранен тест
pvalue <- 2 * min(
  pt(t, degrees_of_freedom, lower.tail = TRUE),
  pt(t, degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pvalue, "df" = degrees_of_freedom)

##           t      pvalue          df 
## -1.86249033  0.06685454 68.00000000

Или с t.test

t.test(x, y, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = -1.8625, df = 68, p-value = 0.06685
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.80818636  0.06232572
## sample estimates:
## mean of x mean of y 
##  2.937953  3.810884

---

Когато работим с предположението, че дисперсиите на двете случайни величини \(X\) и \(Y\) не са еднакви, тоест \(\sigma_X^2 := D(X) \neq D(Y) =: \sigma_Y^2\), за степените на свобода използваме формулата за комбинирани степени на свобода.

\[ df_{c} := \frac{(\frac{S_1^2}{n} + \frac{S_2^2}{m})^2}{\frac{(S_1^2/n)^2}{n-1} + \frac{(S_2^2/m)^2}{m-1}} \]

В този случай Т-статистиката ще е

\[ Т := \frac{\bar{X}-\bar{Y}}{\sqrt{\frac{S_1^2}{n} + \frac{S_2^2}{m}}} \sim \mathcal{T}(df_{c}) \]

Пример:

# брой наблюдения над симулираните случайни извадки
n <- 30
m <- 40

# Симулации с различен брой наблюдения над нормално разпределени с.в.
x <- rnorm(n, 3, 2)
y <- rnorm(m, 4, 3)

# Извадъчни дисперсии
v1 <- var(x)
v2 <- var(y)

# Смесени степени на свобода
degrees_of_freedom <- (
  (v1 / n + v2 / m)^2
  /
  ((v1 / n)^2 / (n - 1) + (v2 / m)^2 / (m - 1)))

t <- (mean(x) - mean(y)) / sqrt(v1 / n + v2 / m)

# Изчисление за pvalue при двустранен тест
pvalue <- 2 * min(
  pt(t, degrees_of_freedom, lower.tail = TRUE),
  pt(t, degrees_of_freedom, lower.tail = FALSE)
)

# резултат
c("t" = t, "pvalue" = pvalue, "df" = degrees_of_freedom)

##         t    pvalue        df 
## -1.151663  0.253680 64.971366

Или за по-кратко с t.test

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

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -1.1517, df = 64.971, p-value = 0.2537
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.9534532  0.5245233
## sample estimates:
## mean of x mean of y 
##  3.077375  3.791840

Сдвоени и несдвоени данни (Paired vs Unpaired)

Сдвоени данни (наблюдения) от две извадки наричаме, такива при при които имаме зависимост между извадките и няблюденията идват от едни и същи индивиди, но в различен етап от време. Тоест \(x_i\) и \(y_i\) са наблюдения върху един и обект. Пример за това би бил ефекта върху нивото на холестерола преди и след приемането на лекарство от едни и същи хора. При несдвоени данни, имаме наблюдения от независими извадки \(X\) и \(Y\). Пример за това би бил, кръвното на група която не приема определено лекарство, и кръвното на група която е приела лекарството. При t.test в този случай, при нулева хипотеза \(H_0 := \mu_0=0\), тест статистиката \(Т\) би била \[ T := \frac{(\bar{X} - \bar{Y})-(\mu_1-\mu_2)}{st.dev(\bar{X}-\bar{Y})} = \frac{\bar{X}_D-\mu_0}{S_D/\sqrt{n}} \sim \mathcal{T}(n - 1) \] където \(\bar{X}_D\) и \(S_D\) са извадъчното средно и стандартно отклонение върху разликите между всички наблюдения.

Пример:

# Брой елементи от които ще симулираме извадка
n <- 30

# Симулираме наблюдения от с.в.
x <- rnorm(n, 80, 4)
y <- rnorm(n, 84, 4)

# Тест статистика при сдвоен тест
t <- mean(x - y) / (sd(x - y) / sqrt(30))

# Изчисление на pvalue
pvalue <- 2 * min(pt(t, n - 1), pt(t, n - 1, lower.tail = FALSE))

# резултат
c("t" = t, "pvalue" = pvalue, "df" = n - 1)

##           t      pvalue          df 
## -2.38593558  0.02378466 29.00000000

Отново, можем да ползваме t.test

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

## 
##  Paired t-test
## 
## data:  x and y
## t = -2.3859, df = 29, p-value = 0.02378
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.5515669 -0.3499634
## sample estimates:
## mean of the differences 
##               -2.450765

Кога ползваме t.test? Тест на Уилкоксън

Т-тестът работи с предположението, че имаме извадки от нормално разпределени случайни величини (популации) и в този случай, трябва да проверим дали нашите данни отговарят на това условие. Поради Централната Гранична Теорема, стига да имаме достатъчно наблюдения (като правило се счита над 30), можем да ползваме Т-тест и без данните ни да следват нормално разпределение, защото

\[ \frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \to_d \mathcal{N}(0, 1) \]

Когато тези условия не са спазени, можем да ползваме Непараметричен тест на Уилкоксън.

Тест на Уилкоксън

Тестът на Уилкоксън е непараметричен тест, който ни позволява да да сравним две извадки, без да правим преположение за вероятностното разпределение което стои зад тях. Предположенията върху които той стъпва са следните

• Наблюденията са независими едно от друго (Случайна Извадка)
• Наблюденията са числени или поне ординални (можем да ги сравним).

При теста на Уилкоксън, работим с нулевата хипотеза \(H_0 := P(x_i > y_i) = \frac{1}{2}\), срещу алтернативата \(H_1 := P(x_i > y_i) \neq \frac{1}{2}\) и следната \(U\)-статистика.

Нека имаме \(X_1, ...,X_n\) е независима случайна извадка от \(X\) и \(Y_1, ...,Y_n\) е независима случайна извадка от \(Y\) и

\[ U := \sum_{i=1}^n\sum_{j=1}^mS(X_i, Y_j) \]

където

\[ S(X, Y) := \begin{cases} 1,& \text{if } Y < X\\ 1/2, & \text{if } Y = X \\ 0, & \text{if } Y > X \\ \end{cases} \] Разпределението на \(U\) е добре известно и го има във всеки справочник. В R функцията на разпределението му е pwilcox. Виждаме, че \(U\) брои броя на инверсиите между Y и Х. Тоест, ако имаме изместеност (Y и X не са с един център), наблюдаваната стойност за \(U\)-статистиката ще бъде в някаква крайна област на разпределението.

Пример:

s <- function(x, y) (y < x) + .5 * (y == x)

x <- c(1, 5, 7,  2, 5,  9, 3, 2, 1, 12)
y <- c(8, 3, 13, 11, 6, 12, 15, 10, 6, 8)

U <- sum(unlist(map(x, ~sum(s(.x, y)))))

c("U" = U, "pvalue" = 2 * pwilcox(U, n = length(x), m = length(y)))

##           U      pvalue 
## 18.00000000  0.01468964

Отново, има функция която изпълнява целия тест

# Дава леко по-различни резултати поради вътрешни оптимизации
wilcox.test(x, y)

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

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

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Виждаме, че разпределението прилича на Нормално разпределение при увеличение на броя на наблюдения \(n\). При достатъчен брой наблюдения, разпределението на тест статистиката може да се апроксимира с нормално.

Тестът на Уилкоксън може да бъде изпълнен при всички от горните случаи

• само върху една извадка - wilcox.test(x, mu = 10, alternative = “two.sided”)
• върху две извадки - wilcox.test(x, y, alternative = “two.sided”)
• сдвоени данни - wilcox.test(x, y, alternative = “two.sided”, paired = TRUE)

---

\(\chi^2\)-тест на Пиърсън

Интуиция за две групи

Ако имаме случайна извадка \((Y_1, ..., Y_n)\) от биномно разпределена сл.в. \(Y \sim Bin(n, p)\), можем да обозначим броя наблюдавани успехи (елементи от едната група) с \(О_S\), а броя наблюдавани неуспехи с \(О_F\) (елементи от другата група), и съответно техните очаквани стойности с \(E_S\) и \(E_F\).

Тогава, от Централната гранична теорема

\[ X := \frac{Y-np}{\sqrt{np(1-p)}} \to_d \mathcal{N}(0, 1) \]

Следователно \(X^2 \to_d \chi^2(1)\) и ако имаме достатъчно наблюдения \(n\), \(X^2\) ще има приблизително \(\chi^2\) разпределение с една степен на свобода.

Oт това, че \(\frac{1}{p(1-p)}=\frac{1}{p} + \frac{1}{n-p}\), и това, че \((Y-np)^2=[(n-Y)-n(1-p)]^2\), можем да представим \(X^2\) по следния начин \[ X^2 := \frac{(Y-np)^2}{np(1-p)} = \frac{(Y-np)^2}{np}+\frac{(Y-np)^2}{n(1-p)}=\\ = \frac{(Y-np)^2}{np}+\frac{[(n-Y)-n(1-p)]^2}{n(1-p)} =\\ = \frac{(O_S-E_S)^2}{E_S}+\frac{(O_F-E_F)^2}{E_F} \] По този начин, при \(n\) наблюдавани стойности, можем да тестваме нулева хипотеза, че очакваният брой успехи (елементи от едната група) \(E_S = np_S\) е определено число, а очаквания брой неуспехи ( елементи от другата група) \(E_F =np_F\)е друго число. Също, можем да тестваме и хипотези за самите вероятности на двете групи \(p_S\) и \(p_F\). Ако нулевата хипотеза не е вярна, наблюдаваната \(X^2\) статистика ще е прекалено надясно от \(0\).

Както при биномно разпределена случайна величина, така и при такава с много групи (мултиномно разпределена), съществува асимптотична апроксимация с нормално разпределение. В такъв случай, можем да направим сходни разсъждения за \(k\) на брой различни групи и \[ X^2 := \sum_{i=1}^k\frac{(O_i-E_i)^2}{E_i} \]

ще има, отново асимптотично приблизително \(\chi^2\) разпределение, но този път с \(k-1\) степени на свобода, където \(О_i\) и \(E_i\) са наблюдавания и очаквания брой елементи под \(H_0\) от \(i\)-тата група.

	Група 1	Група 2	…	Група k
Наблюдавани	\(O_1 := n\hat{p_1}\)	\(O_2 := n\hat{p_2}\)	…	\(O_k := n\hat{p_k}\)
Очаквани под \(H_0\)	\(E_1 := np_1\)	\(E_2 := np_2\)	…	\(E_k := np_k\)

Нека симулираме разпределенията на \(X^2\) при вярна нулева хипотеза и различен брой групи, с цел да видим, че наистина имат приблизително \(\chi^2\) разпределения.

simulate_chsq <- function(n = 60, th_probs, hyp_probs) {
  k = length(hyp_probs)
  # Правим извадка от к на брой групи с теоретични вероятности за всяка - th_probs
  x <- sample(1:k, prob = th_probs, size = n, replace = TRUE)

  # изграждаме брой елементи във всяка група
  observed <- unlist(map(1:k, ~sum(x == .x)))

  # очакваните бройки елементи ако нулевата хипотеза е вярна
  expected <- n * hyp_probs

  # връщаме X^2 статистиката
  sum((observed - expected)^2 / expected)
}

# общ брой елементи от всички групи
n <- 120
# Брой симулации на X^2-статистиката
s <- 2000

# за степени на свобода от 4 до 10 генерираме X^2 статистики
5:11 %>%
  imap(
    ~data.frame(
      # Задаваме истинските вероятности за нулева хипотеза, така че да видим разпределенията
      x = replicate(s, simulate_chsq(n, th_probs = rep(1 / .x, .x), hyp_probs = rep(1 / .x, .x))),
      df = as.factor(.x)
    )
  ) %>%
  # сливаме всички data.frames в един data.frame с цел да визуализираме по-лесно
  bind_rows() %>%
  ggplot(mapping = aes(x = x, fill = df)) +
  geom_density(alpha = 0.25) +
  ggtitle("X^2 статистиката с Хи-квадрат разпределение") +
  xlab("X^2") +
  ylab("Емпирична Плътност") +
  theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  scale_fill_pomological() + 
  scale_color_pomological()

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Нека направим тест с четири групи симулирани с теоретични вероятности 1/4 за всяка и \(n=60\). Тогава, ако тестваме нулевата хипотеза \(H_0 := \{p_1=1/2, p_2 = p_3=p_4=1/6\}\), нашата \(X^2\) статистика би трябвало да е прекалено надясно от нулата.

xsquared <- simulate_chsq(n = 60, th_probs = rep(1 / 4, 4), hyp_probs = c(1/2, 1/6, 1/6, 1/6))

pvalue <- pchisq(xsquared, df = 3, lower.tail = FALSE)

c("X^2" = xsquared, "pvalue" = pvalue, "df" = 3)

##          X^2       pvalue           df 
## 15.800000000  0.001246226  3.000000000

	Група 1	Група 2	Група 3	Група 4
Истински Очаквания	15	15	15	15
\(O_i := n\hat{p_i}\) := Наблюдавани	14	13	17	16
\(E_i := np_i\) := Очаквани под \(H_0\)	30	10	10	10
\(\frac{(O_i-E_i)^2}{E_i}\)	\(\frac{(14-30)^2}{30}=8.53\)	\(\frac{(13-10)^2}{10}=0.90\)	\(\frac{(17-10)^2}{10}=4.90\)	\(\frac{(16-10)^2}{10}=3.60\)

Виждаме, че вероятността да наблюдаваме стойност по-голяма или равна от \(17.93\) на \(X^2\)-статистиката е много малка (\(pvalue := P(X^2\ge17.93) < 0.005\)). В такъв случай, можем да отхвълим нулевата хипотеза, че очакваните вероятности са \(\{p_1=1/2, p_2 = p_3=p_4=1/6\}\).

Критичните области и техните pvalue-та биха изглеждали така

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

## Picking joint bandwidth of 1.71

В R за улеснение, използваме функцията chisq.test(x, p) където x ще е наблюдаван брой елементи в различните групи, а p е вероятностите под нулевата хипотеза, за различните групи.

Приложения на \(\chi^2\)-тест

Проверка за произволно дискретно разпределение

Ако имаме дискретни данни, бихме могли да тестваме хипотезата, че наблюдаваната извадка идва от дискретна случайна величина с избрано от нас разпределение.

Пример

Нека имаме случайна извадка \(X_1, ..., X_{100} \sim^d X \sim Poi(2)\). Искаме да тестваме хипотезата \(H_0 := \lambda = 2.5\) срещу алтернативата \(H_1 := \lambda \neq 2.5\)

# брой наблюдения които имаме
n <- 100

set.seed(2)
# Стойностите които примерно сме наблюдавали
observed_props <- rpois(n, lambda = 2) %>% 
  table() %>%
  prop.table() %>%
  unname() %>%
  as.vector()

# Хипотетични вероятности под нулева хипотеза
hypothet_props <- dpois(0:6, lambda = 2.5)

# Калибрация с цел по-лесно тестване. chisq.test изисква вероятностите да се сумират до 1
hypothet_props <- hypothet_props + (1 - sum(hypothet_props)) / 7

# Визуализация на пропорциите
bind_rows(
  data.frame(x = 0:6, y = hypothet_props, type = "Хипотетични"),
  data.frame(x = 0:6, y = observed_props, type = "Наблюдавани")
)  %>%
  ggplot(mapping = aes(x = x, y = y, fill = type)) +
  geom_bar(position = "dodge", stat = "identity") +
  ggtitle("Наблюдавани и хипотетични вероятности") +
  xlab("X") +
  ylab("Плътност") +
  theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  scale_fill_pomological() + 
  scale_color_pomological()

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Така, ако нулевата хипотеза е вярна,

\[ X^2 := \sum_{i=0}^6\frac{(n\hat{p_i}-np_i)^2}{np_i} \sim^d_{aprox} \chi^2_{6} \]

observed <- n * observed_props
expected <- n * hypothet_props

xsquared <- sum((observed - expected)^2 / expected)
pvalue <- pchisq(xsquared, lower.tail = FALSE, df = 6)

c("X^2" = xsquared, "df" = 6, "pvalue" = round(pvalue, 6))

##       X^2        df    pvalue 
## 16.964442  6.000000  0.009415

Така виждаме, че е много малко вероятно да наблюдаваме съответната тестова статистика, ако истинското разпоределение беше \(Poi(2.5)\). В такъв случай, бихме могли да отхвърлим нулевата хипотеза при избрано от нас ниво на съгласие \(\alpha\).

В R за по-лесно това би станало така

chisq.test(observed, p = hypothet_props)

## Warning in chisq.test(observed, p = hypothet_props): Chi-squared approximation
## may be incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 16.964, df = 6, p-value = 0.009415

---

Проверка за произволно непрекъснато разпределение

Ако имаме непрекъснати данни и искаме да тестваме хипотезата, че те идват от случайна величина с дадено непрекъснато разпределение, бихме могли да дискретизираме данните в под интервали и да проверим дали честотата с която се падат е същата като при същата дискретизация на плътността на тестваното от нас разпределение.

Пример

Нека имаме \(2000\) на брой наблюдения \(X_1, ..., X_{2000}\) от случайна величина \(X \sim \Gamma(12, 4)\) и нека тестваме нулева хипотеза с истинските параметри \(\alpha = 12, \beta = 4\). В този случай, не би трябвало да можем да отхвърлим нулевата хипотеза.

# Брой наблюдения
n <- 2000

# Параметри под нулевата хипотеза за разпределението
hyp_shape <- 12; hyp_rate <- 4;

# seed с цел възпроизвеждане на резултатите
set.seed(16)

# Симулирани данни от Гама разпределена случайна величина
x <- rgamma(n, 12, 4)

# Дискретизация
divisions <- seq(0.5, 6.5, by = 0.5)

# Хипотетични честоти ако нулевата хипотеза е вярна
hyp_props <- pgamma(divisions, hyp_shape, hyp_rate) - pgamma(divisions - 0.5, hyp_shape, hyp_rate)

# таблица с брой наблюдения в интервалите
observed_table <- x %>%
  cut(c(0, divisions)) %>%
  table()

# Наблюдавани честоти
observed <- as.vector(unname(observed_table))

# Визуализация на честотите (подобно е на хистограма)
bind_rows(
  data.frame(x = names(observed_table), y = observed, type = "Наблюдавани"),
  data.frame(x = names(observed_table), y = n * hyp_props, type = "Хипотетични"),
) %>%
  ggplot(mapping = aes(x = x, y = y, fill = type)) +
  geom_bar(position = "stack", stat = "identity") +
  coord_flip() +
  xlab("Интервали на Дискретизация") +
  ylab("Честота на срещане в интервала") +
  ggtitle("Честота на Дискретизирани данни") +
  theme_pomological_nobg(base_family = "Pacifico", base_size = 16) +
  scale_fill_pomological() + 
  scale_color_pomological()

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

# добавяме остатъка за да се сумират до 1
hyp_props <- hyp_props + (1 - sum(hyp_props)) / length(hyp_props)

# изпълнение на теста
chisq.test(observed, p = hyp_props)

## Warning in chisq.test(observed, p = hyp_props): Chi-squared approximation may be
## incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 7.6662, df = 12, p-value = 0.8106

Не успяваме да отхвърлим нулевата хипотеза, че данните са наблюдения от случайна величина \(X \sim \Gamma (12, 4)\).

---

Проверка за независимост на две случайни величини

Ако имаме две случайни извадки \(X_1, ..., X_n \sim^d X\) и \(Y_1, ..., Y_n \sim^d Y\), бихме могли отново да изградим таблици с пропорции \(\hat{p}_{ij} :=\hat{P}(X = i, Y = j)\) от наблюдаваните стойности, и съответно маргиналните пропорции (емпирични вероятности)
\(\hat{p}_{i, \forall}:= \hat{P}(X=i)=\sum_{j = 1}^c\hat{p}_{ij}\) и \(\hat{p}_{\forall, j}:= \hat{P}(Y=j)=\sum_{i = 1}^r\hat{p}_{ij}\), където \(r\) и \(c\) са бройките елементи за \(X\) и \(Y\).

В такъв случай, ако случайните величини \(X\) и \(Y\) са независими, то \(P(X = x_i, Y = y_j) = P(X = x_i)\cdot P(Y=y_j) := p_{i, \forall}p_{\forall,j}\) за всяко \(i, j\) и би трябвало да наблюдаваме емпирични вероятности които приблизително притежават това свойство.

Тогава

\[ X^2 := \sum_{i=1}^r\sum_{j=1}^c \frac{(O_{i,j}-E_{i,j})^2}{E_{i,j}} \underset{appr}{\overset{distr}{\sim}} \chi^2_{(r-1)(c-1)} \]

Пример

Нека симулираме две извадки с големина \(n=600\) от независими случайни величини \(X \sim Bin(8, 0.5), Y \sim Poi(2.3)\) и изпълним \(\chi^2\)-тест за независимост, като евентуално след това, бихме могли да добавим зависимост между част от наблюденията с цел да видим как се променят резултатите от \(\chi^2\) теста.

# генериране на две независими случайни величини
x <- rbinom(600, size = 8, prob = 1 / 2)
y <- rpois(600, 2.3)

# Вкарване на зависимост
y[400:600] <- y[400:600] + x[400:600]

# таблица брояща наблюдавани комбинации от вида (X, Y)
observed <- table(x, y)

# брой на всички елементи, N = n * m
N <- sum(observed)

# Маргинални емпирични вероятности
x_marginal <- prop.table(table(x))
y_marginal <- prop.table(table(y))

# Очаквани бройки при предположение за независимост
expected <- N * (x_marginal %*% t(y_marginal))

# X^2 статистика
xsquare <- sum((observed - expected)^2 / expected)

# Степени на свобода df = (r - 1) * (s - 1)
degrees_of_freedom <- prod(dim(observed) - 1)

# pvalue за съответната наблюдавана Х^2 статистика
pvalue <- pchisq(xsquare, df = degrees_of_freedom, lower.tail = FALSE)

# връщане на резултат
c("X^2" = xsquare, "df" = degrees_of_freedom, "pvalue" = pvalue)

##          X^2           df       pvalue 
## 1.461809e+02 1.040000e+02 4.054662e-03

Или за кратко с функцията chisq.test

chisq.test(observed)

## Warning in chisq.test(observed): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  observed
## X-squared = 146.18, df = 104, p-value = 0.004055

Функцията chisq.test автоматично разбира, че става въпрос за тест за независимост, когато ѝ се подаде таблица.

---

Проста Линейна Регресия

Често се налага да обясним една непрекъсната променлива чрез друга, примерно

• заплата на спрямо години в образование
• скорост на кола на базата на конски сили
• резултат от изпит спрямо часовете вкарани в учене
• продажби на продукт спрямо приходите вложени в различни реклами

Нека разгледаме данни от последния пример

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Виждаме, че има добра линейна връзка между печалбите и приходите вложени в ТВ реклами. Данните като цяло се движат около права линия около която има лек случаен шум. Ако моделираме печалбата като случайна величини \(Y\) и \(X\), бихме могли да кажем, че

\[ Y = \beta_0 + \beta_1X + \epsilon, \] където \(\epsilon\) е нашия случаен елемент, който за леснота можем да считаме като независими наблюдения oт нормално разпределена случайна величина, тоест

\[ \epsilon_1, ..., \epsilon_n \sim^d \epsilon \sim \mathcal{N}(0, \sigma^2) \]

За леснота, можем да считаме, че \(Х\) са наблюдавани константи и техните вероятностни свойства не ни интересуват.

По този начин,

\[ E[Y] = E(\beta_0 + \beta_1x + \epsilon) = \beta_0 + \beta_1x\\ Var[Y] = Var(\beta_0 + \beta_1x + \epsilon) = \sigma^2\\ Y \sim \mathcal{N}(\beta_0+\beta_1x, \sigma^2) \]

и математическото очакването за \(Y\) ще е права линия около която наблюденията варират.

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Един начин по който бихме могли да оценим математическото очакване на \(Y\) е метода на най-малките квадрати.

Понеже знаем, че върху математическото очакване може да се гледа като на център на тежестта, естествен начин за неговата оценка е да изберем минималното разстояние между вектора на наблюденията \(\overrightarrow{y} = (y_1, ..., y_n)'\) и вектора-линия \(\beta_0 + \beta_1\overrightarrow{x}=(\beta_0 + \beta_1x_1, ..., \beta_0 + \beta_1x_n)'\), тоест за \(\beta_0\) и \(\beta_1\) да изберем такива коефициенти, че

\[ (\hat{\beta_0}, \hat{\beta_1}) := \underset{\beta_0, \beta_1}{argmin}\text{ }\rho(\overrightarrow{y}, \beta_0 + \beta_1 \overrightarrow{x})^2=\\ =\underset{\beta_0, \beta_1}{argmin}\text{ }\left\lVert \overrightarrow{y}-\beta_0 - \beta_1 \overrightarrow{x} \right\rVert^2_2=\\ =\underset{\beta_0, \beta_1}{argmin}\text{ }\sum_{i=1}^n(y_i-\beta_0-\beta_1 x_i)^2 \]

Тоест, избираме коефициенти \(\beta_0\) и \(\beta_1\), такива че нашата линия \(f(x) := \beta_0 + \beta_1x\) максимално да приближава наблюденията \((y_1, ..., y_n)\) които имаме. Функцията която минимизираме по-горе, наричаме сума на квадратите на грешките и най-често обозначаваме с RSS,

която като функция на \(\beta_0\) и \(\beta_1\) би изглеждала така

Един начин по който можем да намерим минимума на горната функция е като намерим първите частни производни и ги положим на 0, т.е.

\[ \frac{\partial{(RSS)}}{\beta_0}=0\\ \frac{\partial{(RSS)}}{\beta_1}=0 \]

и да намерим инфлексната точка на функцията. Решавайки за \(\beta_0\) и \(\beta_1\) получаваме

\[ \hat{\beta_0} = \bar{y}-\beta_1\bar{x}\\ \]

и съответно

\[ \hat{\beta_1} = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}\\ \] В R това би станало така

x <- advertisement$tv
y <- advertisement$sales

# Оценка за коефициента β1
beta1_estimate <- sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2)

# Оценка за коефициента β0
beta0_estimate <- mean(y) - beta1_estimate * mean(x)

# връщане на резултат
c("β0 Оценка" = beta0_estimate, "β1 Оценка" = beta1_estimate)

##   β0 Оценка   β1 Оценка 
## 111.6637544   0.8243548

Или просто с функцията lm

lm(sales ~ tv, data = advertisement)

## 
## Call:
## lm(formula = sales ~ tv, data = advertisement)
## 
## Coefficients:
## (Intercept)           tv  
##    111.6638       0.8244

Забелязваме, че оценката на \(\beta_1\) може да се представи и така

\[ \hat{\beta_1} = \langle\frac{(\overrightarrow{x}-\bar{x})}{||\overrightarrow{x}-\bar{x}||_2},\frac{(\overrightarrow{y}-\bar{y})}{||\overrightarrow{x}-\bar{x}||_2} \rangle \]

Тоест, оценката за \(\beta_1\) ни дава колко наблюдаваните вектори \(\overrightarrow{x}\) и \(\overrightarrow{y}\) (центрирани и нормирани спрямо \(\overrightarrow{x}\)) са в една посока, или колко \(\overrightarrow{y}\) се променя спрямо \(\overrightarrow{x}\) (скаларно произведение).

По този начин, получената оценка за \(EY\) e

\[ \hat{Y} :=: \hat{f}(x) := \hat{\beta_0}+\hat{\beta_1}x \]

---

Свойства на \(\hat{\beta}_0\) и \(\hat{\beta}_1\)

Нека разгледаме свойствата на оценките \(\hat{\beta_0}\) и \(\hat{\beta_1}\).

\[ Е[\hat{\beta}_1] = E[\frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2}] =\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})E[Y_i]}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})(\beta_0 + \beta_1x_i)}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})\beta_0}{\sum_{i=1}^n(x_i-\bar{x})^2}+\frac{\beta_1\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x}+\bar{x})}{\sum_{i=1}^n(x_i-\bar{x})^2}=\\ 0+\frac{\beta_1\sum_{i=1}^n(x_i-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}+0=\\ =\beta_1 \]

и така, оценката \(\hat{\beta_1}\) за \(\beta_1\) е неизместена оценка.

\[ Е[\hat{\beta_0}]=E[\bar{y}-\hat{\beta_1}\bar{x}]=E[\beta_0+\beta_1\bar{x}+\bar{\epsilon}-\hat{\beta_1}\bar{x}]=\beta_0 \]

и следователно \(\hat{\beta}_0\) също е неизместена оценка за \(\beta_0\).

\[ Var[\hat{\beta}_1] = Var[\frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2}]=\\ =\frac{\sum_{i=1}^n(x_i-\bar{x})^2\sigma^2}{[\sum_{i=1}^n(x_i-\bar{x})^2]^2}=\frac{\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

За дисперсията на \(\hat{\beta}_0\) получаваме

\[ Var[\hat{\beta}_0]=Var[\bar{y}-\hat{\beta_1}\bar{x}]=\\ =\frac{\sigma^2}{n}+\bar{x}^2Var[\hat{\beta_1}]+Cov[\bar{y}, \hat{\beta_1}]=\\ =\frac{\sigma^2}{n}+\frac{\bar{x}^2\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

понеже \(Cov[\bar{y}, \hat{\beta_1}] = 0\) (Проверете сами!).

За ковариацията между \(\hat{\beta}_0\) и \(\hat{\beta}_1\) получаваме

\[ Cov(\hat{\beta}_0, \hat{\beta}_1) = E[(\bar{y}-\hat{\beta}_1\bar{x})\hat{\beta}_1] - \beta_0\beta_1= E[\bar{y}]E[\hat{\beta}_1]-\bar{x}E[\hat{\beta}_1^2]=\\ =(\beta_0+\beta_1\bar{x})\beta_1-\bar{x}Var[\hat{\beta}_1]-\bar{x}\beta_1^2-\beta_0\beta_1=\\ =-\bar{x}Var[\hat{\beta}_1]=-\frac{\bar{x}\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \]

---

Оценка за \(\sigma^2\)

Нека разгледаме минимизираната стойност на \(RSS\) като случайна величина

\[ \hat{RSS} := \rho(\overrightarrow{Y},\overrightarrow{\hat{Y}})^2 = ||\overrightarrow{Y}- \overrightarrow{\hat{Y}}||_2^2 = \sum_{i=1}^n{(Y_i-\hat{Y}_i)^2} = \sum_{i=1}^n{(Y_i-\hat{\beta}_0 -\hat{\beta}_1x_i)^2} \]

Тогава

\[ \frac{\hat{RSS}}{\sigma^2} \sim \chi^2_{n-2} \] и ако вземем математическото очакване на \(\hat{RSS}\)

\[ E[\frac{\hat{RSS}}{\sigma^2}] = n-2 \implies E[\frac{\hat{RSS}}{n-2}] = \sigma^2 \]

Така получаваме неизместена оценка \(\hat{\sigma}^2 := \frac{\hat{RSS}}{n-2}\)

Получената оценка можем да използваме в гореполучените дисперсии на \(\hat{\beta_0}\) и \(\hat{\beta_1}\) за да изградим оценки за техните дисперсии.

# брой елементи които наблюдаваме
n <- nrow(advertisement)

# Оценка за очакването на Y(x)
y_fitted <- beta0_estimate + beta1_estimate * x

# Оценка за дисперсията
var_estimate <- sum((y_fitted - y)^2) / (n - 2)

# Оценка за стандартното отклонение
sd_estimate <- sqrt(var_estimate)

# Оценка за дисперсията на оценката на β0
beta0_var_est <- var_estimate / n + (mean(x)^2 * var_estimate) / sum((x - mean(x))^2)

# Оценка за дисперсията на оценката на β1
beta1_var_est <- var_estimate / (sum((x - mean(x))^2))

# Оценка за стандартното отклонение на оценката на β0
beta0_sd_est <- sqrt(beta0_var_est)

# Оценка за стандартното отклонение на оценката на β1
beta1_sd_est <- sqrt(beta1_var_est)

# резултат
data.frame(
  estimation = c(
    "Оценка за стандартното отклонение",
    "Оценка за стандартното отклонение на оценката за β0",
    "Оценка за стандартното отклонение на оценката за β1"
  ),
  value = c(sd_estimate, beta0_sd_est, beta1_sd_est)
)

##                                            estimation       value
## 1                   Оценка за стандартното отклонение 37.80484374
## 2 Оценка за стандартното отклонение на оценката за β0  7.63869509
## 3 Оценка за стандартното отклонение на оценката за β1  0.04595558

Или написано по-кратко чрез lm и summary

lm(sales ~ tv, advertisement) %>%
  summary() %>%
  {data.frame(
    estimation = c(
      "Оценка за стандартното отклонение",
      "Оценка за стандартното отклонение на оценката за β0",
      "Оценка за стандартното отклонение на оценката за β1"
    ),
    value = c(.$sigma, .$coefficients[1, 2], .$coefficients[2, 2])
  )}

##                                            estimation       value
## 1                   Оценка за стандартното отклонение 37.80484374
## 2 Оценка за стандартното отклонение на оценката за β0  7.63869509
## 3 Оценка за стандартното отклонение на оценката за β1  0.04595558

---

Интервали на достоверност в Проста Линейна Регресия

Тъй като \(\hat{\beta}_0 \sim \mathcal{N}(\beta_0, \frac{\sigma^2}{n}+\frac{\bar{x}^2\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2})\) и следователно

\[ Z_0 := \frac{\hat{\beta}_0 - \beta_0}{\sqrt{Var(\hat{\beta}_0)}}:= \frac{\hat{\beta}_0 - \beta_0}{\sqrt{\frac{\sigma^2}{n}+\frac{\bar{x}^2\sigma^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}} \sim \mathcal{N}(0, 1) \]

и

\[ T_0 := \frac{\hat{\beta_0}-\beta_0}{\sqrt{\hat{Var}(\hat{\beta}_0)}} = \frac{\hat{\beta}_0 - \beta_0}{\sqrt{\frac{\hat{\sigma}^2}{n}+\frac{\bar{x}^2\hat{\sigma}^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}} \sim \mathcal{T}(n - 2) \]

то можем да изградим \(\gamma\)-доверителни интервали, както и да тестваме различни хипотези на базата на горните статистики.

Аналогични разсъждения са валидни и за \(\beta_1\), тъй като

\[ Z_1 := \frac{\hat{\beta}_1-\beta_1}{\sqrt{Var(\hat{\beta_1})}}= \frac{(\hat{\beta}_1-\beta_1)\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2}}{\sigma} \sim \mathcal{N}(0, 1) \] и

\[ T_1 := \frac{\hat{\beta}_1-\beta_1}{\sqrt{\hat{Var}(\hat{\beta}_1)}} = \frac{(\hat{\beta}_1-\beta_1)\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2}}{\hat{\sigma}} \sim \mathcal{T}(n - 2) \]

• 95% доверителен интервал за \(\beta_0\)

# Квантил от Т-разпределение с n-2 степени на свобода
q <- pt(0.975, df = n - 2)

# 95% доверителен интервал за параметъра β0
c(beta0_estimate - q * beta0_sd_est, beta0_estimate + q * beta0_sd_est)

## [1] 105.2929 118.0346

• 95% доверителен интервал за \(\beta_1\)

c(beta1_estimate - q * beta1_sd_est, beta1_estimate + q * beta1_sd_est)

## [1] 0.7860270 0.8626827

• Тестване на хипотеза че \(\beta_0=0\)

# Хипотетичен параметър под който т статистиката има Т-разпределение
hypothetical_beta0 <- 0

# наблюдавана т-статистика
t0 <- (beta0_estimate - hypothetical_beta0) / beta0_sd_est

# pvalue за т-статистиката
pval0 <- 2 * min(pt(t0, df = n - 2, lower.tail = TRUE), pt(t0, df = n - 2, lower.tail = FALSE))

# резултат
c("t" = t0, "pvalue" = pval0)

##            t       pvalue 
## 1.461817e+01 2.305746e-26

• Тестване на хипотеза че \(\beta_1=0\)

# Параметър под нулева хипотеза
hypothetical_beta1 <- 0

# тест статистика
t1 <- (beta1_estimate - hypothetical_beta1) / beta1_sd_est

# pvalue за тест статистиката
pval1 <- 2 * min(pt(t1, df = n - 2, lower.tail = TRUE), pt(t1, df = n - 2, lower.tail = FALSE))

# резултат
c("t" = t1, "pvalue" = pval1)

##            t       pvalue 
## 1.793808e+01 1.009070e-32

Отново, резултатът го има чрез функцията summary

lm(sales ~ tv, advertisement) %>%
  summary() %>%
  .$coefficients

##                Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 111.6637544 7.63869509 14.61817 2.305746e-26
## tv            0.8243548 0.04595558 17.93808 1.009070e-32

Всъщност, въпросът дали \(\beta_1=0\) може да се преведе в това дали \(X\) има ефект върху \(Y\) и дали участва в обясняването на \(Y\).

---

Прогнозиране

Основен въпрос, след като имаме изграден модел, е това да направим прогноза върху ненаблюдавана досега точка \(x_0\). В нашия случай, колко очакваме да спечелим, ако инвестираме 144лв в реклама?

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

## `geom_smooth()` using formula 'y ~ x'

Използвайки досега напаснатия модел, оценката за очакването на y върху определена точка \(x_0\) е

\[ \hat{y}_{x_0} := \hat{E}[Y|X = x_0] := \hat{\beta}_0 + \hat{\beta}_1 x_0 \]

Като точкова оценка, то притежава свойството неизместеност, защото

\[ Е[\hat{y}_{x_0} - y_{x_0}] = E[\hat{\beta}_0 + \hat{\beta}_1 x_0 - \beta_0 - \beta_1x_0] =\\ = E[\hat{\beta}_0] + E[\hat{\beta}_1 x_0] - \beta_0 - \beta_1x_0 =\\ = \beta_0 + \beta_1x_0 - \beta_0 - \beta_1x_0 = 0 \]

Дисперсията на тази оценка е сумата от дисперсиите на оценките на коефициентите

\[ Var[\hat{y}_{x_0}] = Var[\hat{\beta}_0 + \hat{\beta}_1 x_0] = \\ = Var[\hat{\beta}_0] + Var[\hat{\beta}_1 x_0] + 2Cov[\hat{\beta}_0, \hat{\beta}_1 x_0] =\\ =\frac{\sigma^2}{n}+\frac{\bar{x}^2\sigma^2}{S_{xx}} + \frac{x_0^2\sigma^2}{S_{xx}} - 2\frac{x_0\bar{x}\sigma^2}{S_{xx}} = \\ = \sigma^2(\frac{1}{n} + \frac{(\bar{x}-x_0)^2}{S_{xx}}) \]

където, за леснота сме положили \(S_{xx} := \sum_{i=1}^n{(x_i-\bar{x})^2}\).

Ако правим оценка за дисперсията на \(\hat{y}_{x_0}\), необходимо е в горната формула просто да заместим дисперсията \(\sigma^2\) с оценката за дисперсията \(\hat{\sigma}^2\).

\[ \hat{Var}[\hat{y}_{x_0}] = \hat{\sigma}^2(\frac{1}{n} + \frac{(\bar{x}-x_0)^2}{S_{xx}}) \]

Oсвен точкова оценка за очакваната стойност в \(x_0\), можем да изградим и интервална такава посредстом следната \(T_{x_0}\)-статистика

\[ Т_{x_0} := \frac{\hat{y}_{x_0} - E[y|x_0]}{\sqrt{\hat{Var}[\hat{y}_{x_0}]}} \sim \mathcal{T}(n - 2) \]

# Точка която да прогнозираме
x0 <- 300

# Оценка за очакването в тази точка
fit_at_x0 <- beta0_estimate + beta1_estimate * x0

Sxx <- sum((x - mean(x))^2)

# Дисперсия и стандартно отклонение на у0
y0_sigma_sq_estimate <- var_estimate * ((1 / n) + ((x0 - mean(x))^2 / Sxx))
y0_sigma_estimate <- sqrt(y0_sigma_sq_estimate)

# смятане на квантил от Т-разпределение с цел да сметнем доверителен интервал
q <- qt(0.975, df = n - 2)

# оценка и 95% доверителен интервал
c(
  "оценка" = fit_at_x0, 
  "долен край" = fit_at_x0 - q * y0_sigma_estimate,
  "горен край" = fit_at_x0 + q * y0_sigma_estimate
)

##     оценка долен край горен край 
##   358.9702   342.9216   375.0188

Което може да стане накратко така

predict.lm(
  lm(sales ~ tv, advertisement),
  data.frame(tv = c(300)),
  interval = "confidence",
  level = 0.95
)

##        fit      lwr      upr
## 1 358.9702 342.9216 375.0188

Точност на модела

RSE

След като сме отхвърлили хипотезата, че \(\beta_1 = 0\), тоест това че \(x\) няма ефект върху \(y\), бихме могли да попитаме каква е разликата между напаснатия модел \(\hat{Y}\) и вектора \(Y\), тоест

\[ RSE := \sqrt{\frac{\hat{RSS}}{n-2}} \]

sd_estimate

## [1] 37.80484

summary(lm(sales ~ tv, data = advertisement))

## 
## Call:
## lm(formula = sales ~ tv, data = advertisement)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -89.730 -22.101   0.514  26.104  83.086 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 111.66375    7.63870   14.62   <2e-16 ***
## tv            0.82435    0.04596   17.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37.8 on 98 degrees of freedom
## Multiple R-squared:  0.7665, Adjusted R-squared:  0.7642 
## F-statistic: 321.8 on 1 and 98 DF,  p-value: < 2.2e-16

Това дали наблюдаваното \(RSE\) е допустима грешка зависи от мащаба на y. Ако разделим \(RSE\) на \(\bar{y}\)

sd_estimate / mean(y)

## [1] 0.1638493

Тоест, правим около 16% грешка. RSE e оценка за липса на напасване на модела спрямо данните.

Коефициент на детерминация \(R^2\)

\[ R^2 := \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS} := 1 - \frac{\rho(Y, \hat{Y})^2}{\rho(Y, \bar{Y})^2}=\\= 1 - \frac{||Y - \hat{Y}||^2_2}{||Y - \bar{Y}||^2_2} = 1-\frac{\sum_{i=1}^n(y_i-\hat{\beta}_0-\hat{\beta}_1 x_i)^2}{\sum_{i=1}^n(y_i-\bar{y})^2} \]

1 - sum((y - y_fitted)^2) / sum((y - mean(y))^2)

## [1] 0.7665414

Виждаме, че колкото по-добре е напаснат модела спрямо \(Y\), толкова повече \(R^2\) ще приближава 1, защото и толкова по малко ще е разстоянието между векторите \(Y\) и \(\hat{Y}\).

Регресия с няколко променливи

Простата линейна регресия е полезен подход при интерпретиране на връзка между променливи и прогнозиране на ненаблюдавани точки, но повечето случаи обаче, ни се налага да обясним Y като комбинация на повече от една независима променлива.

Примерно, да разширим горния случай с обясняване на печалбите не само с реклами в телевизия, но и с реклами в радио, вестници, тн. В този случай нашия модел за Y e

\[ y_i = \beta_0 + \beta_1 x_{i1} + ... + \beta_p x_{ip} + \epsilon_i \]

където i варира от 1 до n (брой наблюдения), \(\epsilon_i\) са независими наблюдения от нормално разпределена случайна грешка \(\epsilon \sim \mathcal N(0, \sigma^2)\), a \(\beta_j\) e връзката между \(Y\) и независимата променлива \(Х_j\). Интерпретираме \(\beta_j\) като средния ефект върху \(Y\) при една единица промяна в \(X_j\), държейки другите променливи фиксирани.

Целта на задачата се измества от това да оценим \(\beta_0, \beta_1\) и \(\sigma^2\), към това да оценим \(\beta_0, ..., \beta_p\) и \(\sigma^2\).

Отново, функцията на минимизиране е

\[ RSS = \sum_{i=1}^n(y_i-\hat{y_i})^2 = \sum_{i=1}^n(y_i-\beta_0 - \beta_1x_{i1} - ... - \beta_px_{ip})^2 \]

При \(p = 2\) целта е да напаснем вместо линия, цяла равнина, която да е максимално близо до наблюдаваните точки \(Y\).

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

При прилагане на Линейна регресия едни от въпросите които ни интересуват са следните

• Има ли поне една от променливите \(X_1, ..., X_p\) която е полезна за напасването на модела?
• Всички независими променливи (предиктори) \(X_1,...,X_p\) ли са полезни в това да обяснят \(Y\)?
• Колко добре е напаснат модела?

1. Има ли връзка между отговора (зависимата променлива) \(Y\) и независимите променливи (предиктори) \(X_1, ..., X_p\)?

В простата линейна регресия, за да проверим дали \(X\) има ефект върху \(Y\), беше необходимо просто да проверим хипотезата дали \(\beta_1 = 0\). Когато имаме повече променливи, въпросът се променя в това да проверим хипотезата дали \(H_0 := \beta_1 = \beta_2 = ... = \beta_p = 0\) срещу алтернативата \(H_1 := \exists j : \beta_j \neq 0\).

Нека първо дефинираме няколко помощни метрики

Където \[ TSS := \rho(Y, \bar{Y})^2 = ||Y - \bar{Y}||_2^2 = \sum_{i=1}^n(y_i-\bar{y})^2\\ RSS := \rho(Y, \hat{Y})^2 = ||Y - \hat{Y}||_2^2 = \sum_{i=1}^n(y_i-\hat{y})^2\\ SSR := \rho(\hat{Y}, \bar{Y})^2 = ||\hat{Y} - \bar{Y}||_2^2 = \sum_{i=1}^n(\hat{y_i}-\bar{y})^2 \]

• TSS (Total Sum of Squares) описва цялостната вариация на данните \(Y = (y_1, ..., y_n)\)
• RSS (Residual sum of Squares) описва вариацията определена от напаснатия модел
• SSR (Sum of Squares Regression) описва вариацията която не е определена от модела

Може да се наблюдава следното правило

\[ TSS = RSS + SSR \] Нека дефинираме следната F-статистика

\[ F = \frac{(TSS - RSS) / p}{RSS / (n - p - 1)} = \frac{(SSR) / p}{RSS / (n - p - 1)} = \frac{[\sum_{i=1}^n(y_i-\bar{y})^2 - \sum_{i=1}^n(y_i-\hat{y_i})^2] / p}{\sum_{i=1}^n(y_i-\hat{y_i})^2 / (n - p - 1)} \]

Така, ако нулевата хипотеза, че \(\beta_1 = \beta_2 = ... = \beta_p = 0\), би трябвало да имаме \(F\)-статистика близка до единица, защото вариацията обяснена от модела не би трябвало да се различава особено от тази която не е обяснена от модела (модела не носи стойност и \(\beta_1 = \beta_2 = ... = \beta_p = 0\)).

Също така, може да се покаже, че ако предположенията на линейната регресия са верни,

\[ Е[RSS / (n - p - 1)] = \sigma^2 \\ Е[SSR / p] = \sigma^2 \] което отново показва, че ако нулевата хипотеза е вярна, \(F\)-статистиката би била близка до 1.

Освен това, \((TSS - RSS)/\sigma^2 \sim \chi^2_{p}\) и \(RSS/\sigma^2 \sim \chi^2_{n-p-1}\), откъдето виждаме, че \(F\)-статистиката е деление на две нормирани Хи-квадрат разпраделени случайни величини. Полученото разпределение за \(F\) се нарича \(\mathcal{F}\)-разпределение с параметри \(df_1\) и \(df_2\) които са степените на свобода на хи-квадрат разпределените случайни величини в знаменател и числител.

По този начин, можем да измерим вероятността да наблюдаваме определена \(F\)-статистика (или по-крайна) с функцията pf(q, df1, df2). Ако съответната вероятност е прекалено малка (примерно по-малка от \(\alpha := 0.05\)), можем да отхвърлим нулевата хипотеза \(H_0 := \beta_1 = \beta_2 = ... = \beta_p = 0\) и да примем алтернативата, че съществува променлива \(X_j\) която подобрява модела.

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

## Warning: Removed 1134 rows containing non-finite values (stat_density).

# изграждане на модел с много променливи
multiple_model <- lm(sales ~ tv + radio + newspaper + facebook, data = advertisement2)

# напаснати данни
y_fitted <- fitted(multiple_model)

# истински данни
y <- advertisement2$sales

# брой елементи
n <- 100

# брой предиктори
p <- 4

# Total sum of squares (TSS)
tss <- sum((y - mean(y))^2)
# Residual sum of squares (RSS)
rss <- sum((y_fitted - y)^2)

# изчисление на F-статистика
f_statistic <- ((tss - rss) / p) / (rss / (n - p - 1))

# изчисление на pvalue за F-статистиката
pvalue <- pf(f_statistic, p, n - p - 1, lower.tail = FALSE)

# резултат
c("f_statistic" = f_statistic, "pvalue" = pvalue, "df1" = p, "df2" = n - p - 1)

##  f_statistic       pvalue          df1          df2 
## 3.023998e+02 4.076768e-53 4.000000e+00 9.500000e+01

Или просто със summary

multiple_model %>%
  # обобщение на модела
  summary() %>%
  # взимане на f-статистиката от обобщението
  .$fstatistic

##    value    numdf    dendf 
## 302.3998   4.0000  95.0000

2. Намиране на важни променливи

Един начин за намиране на важни и обясняващи \(Y\) променливи \(Х_1, ..., Х_п\) e първо да пробваме с празен модел (само с единичен вектор умножен по \(\beta_0\))

\[ Y = \beta_0 + \epsilon \] и след това да напаснем p на брой линейни модела за всяка променлива \(X_p\)

\[ Y = \beta_0 + \beta_1 X_j + \epsilon \] И да изберем \(X_j\) такова, че \(R^2\) да увеличи максимално (или RSЕ да се намали максимално). Повтаряме процеса, като вече добавяме недобавяна досега променлива към модела \(Y = \beta_0 + \beta_1 X_j + \epsilon\), тоест \[ Y = \beta_0 + \beta_1 X_j + \beta_1 X_k +\epsilon \] Продължаваме докато се изпълни някакво зададено от нас терминално условие.

F-тестове и Анализ на Вариациите

Ако тестваме линейна регресия без независима променлива \(X\), тоест \(Y = \beta_0 + \epsilon\), забелязваме че \[ \hat{\beta}_0 = \bar{y} \] а това да тестваме дали \(\beta_0 = 0\), се транслира до \[ T := \frac{\hat{\beta}_0 - 0}{\hat{sd}(\hat{\beta_0})} = \frac{\bar{y} - 0}{\hat{\sigma}/\sqrt{n}} \sim \mathcal{T}(n-1) \]

Така виждаме, че всъщност тестването на хипотезата дали \(\mu = 0\) в една нормално разпределена популация, е всъщност частен случай на линейна регресия без променливи (само с единичен вектор).

# общ брой елементи
n <- 100

# генериране на извадка със средно 0.4
y <- rnorm(n, mean = 0.4)

# изпълняване на т-тест с нулева хипотеза Но := µ = 0
t.test(y, mu = 0)

## 
##  One Sample t-test
## 
## data:  y
## t = 4.7652, df = 99, p-value = 6.476e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2786152 0.6761941
## sample estimates:
## mean of x 
## 0.4774047

Или същото с линеен модел

# линеен модел без променливи х (оценява само средното на у)
lm(y ~ 1) %>%
  # извикване на обобщение
  summary() %>%
  # взимане на т-статистиката и нейното pvalue за хипотезата, че ß0 = 0
  {.$coefficients[, c(3, 4)]}

##      t value     Pr(>|t|) 
## 4.765215e+00 6.475705e-06

Ако искаме да тестваме за произволно средно, трябва просто да извадим това средно от данните

# линеен модел без променливи х (оценява само средното на у с изместване единица)
lm(y - 1 ~ 1) %>%
  # обобщение
  summary() %>%
  # взимане на т-статистиката и нейното pvalue за хипотезата, че ß0 = 0 (което е дали средното е 1 като вземем предвид изместването)
  {.$coefficients[, c(3, 4)]}

##       t value      Pr(>|t|) 
## -5.216286e+00  1.003722e-06

Или с t.test

t.test(y, mu = 1)

## 
##  One Sample t-test
## 
## data:  y
## t = -5.2163, df = 99, p-value = 1.004e-06
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
##  0.2786152 0.6761941
## sample estimates:
## mean of x 
## 0.4774047

Но какво правим, когато имаме две нормални извадки, всяка с размер n и искаме да направим тест дали идват от популации с еднакви средни?

Ако разгледаме следния модел

\[ Y = \beta_0 + \beta_1X + \epsilon \]

където Y e конкатенация на двете извадки, а \(X = \unicode{x1D7D9}_{Y \in \text{извадка 2}}\), тоест \(x_i = 0\) за \(i = 1,..., n\) и \(x_i = 1\) за \(i = n + 1, ...., 2n\)

В горното уравнение на линия, отново \(\beta_1\) представлява наклона който имаме. В този смисъл \[ \beta_1 = E[\frac{\Delta y}{\Delta x}] = E[\frac{\Delta y}{1}] = E[\Delta y] \]

{ height: 100px; width: 100px; }

Така, ако теглиме от първата извадка, оценката за \(Y\) става \(\hat{Y} = \hat{\beta}_0 = \bar{Y}_{\text{извадка 1}}\), а ако теглим от втората извадка \(\hat{Y} = \hat{\beta}_0 + \hat{\beta}_11 = \hat{\beta}_0 + \hat{\beta}_1\) и така \(\hat{\beta}_1\) e оценка за разликата в средните. Ако тестваме дали \(\beta_1 = 0\), реално тестваме дали има разлика в средните на популациите от които произлизат двете извадки.

Оценката за \(\beta_1\) в случая ще бъде

\[ \hat{\beta_1} = \frac{\sum_{i=1}^{2n}(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^{2n}(x_i-\bar{x})^2} = \\ = \frac{\sum_{i=1}^{n}(0 - 0.5)(y_i-\bar{y}) + \sum_{i=n+1}^{2n}(1 - 0.5)(y_i-\bar{y})}{\sum_{i=1}^{2n}(0.5)^2} = \\ = \frac{0.5(\sum_{i=n+1}^{2n}y_i - \sum_{i=1}^{n}y_i)}{2n(0.5)^2} = \\ \frac{\sum_{i=n+1}^{2n}y_i - \sum_{i=1}^{n}y_i}{n} = \bar{y}_{\text{извадка 1}} - \bar{y}_{\text{извадка 2}} \]

Виждаме че всъщност, оценката \(\hat{\beta}_1\) за \(\beta_1\) е разликата в средните между двете извадки.

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Така, ако тестваме дали \(\beta_1 = 0\), получаваме следната T-статистика

\[ Т := \frac{\bar{y}_{\text{извадка 1}} - \bar{y}_{\text{извадка 2}}}{\hat{sd}(\hat{\beta}_1)} \]

# генериране на извадка 1
y1 <- rnorm(50, mean = 10)

# генериране на извадка 2
y2 <- rnorm(50, mean = 15)

# конкатениране на цялостния вектор
y <- c(y1, y2)

# индикаторна х променлива
x <- rep(c(0, 1), each = 50)
# може да се дефинира и така
x <- 0 + (1:100 > 50)

# изпълнение на модела обясняващ y срямо индикаторна променлива
lm(y ~ x) %>%
  # обобщение
  summary() %>% 
  # Взимане на Т-статистиката и нейното pvalue за коефициента ß1
  {.$coefficients[2, c(3, 4)]}

##      t value     Pr(>|t|) 
## 2.316026e+01 1.567575e-41

Което е същото като да направим Т-тест за това дали двете извадки са с еднакво средно.

t.test(y2, y1)

## 
##  Welch Two Sample t-test
## 
## data:  y2 and y1
## t = 23.16, df = 97.837, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  4.398213 5.222580
## sample estimates:
## mean of x mean of y 
##  14.82267  10.01228

ANOVA (Анализ на Вариациите)

В горния пример, бихме могли да тестваме по още един начин, дали променливата Х (тоест дали двете извадки се различават значително) има ефект върху модела. Ако изпълниме F-тест, той на практика би проверил точно това - дали \(\beta_1 = 0\), тоест дали втората извадка статистически значително променя цялостното средно (има статистически значима разлика между двете средни).

\[ H_0 := \{ \beta_1 = 0 \iff \bar{Y} = \beta_0 \iff \text{Извадка 2 e същата като Извадка 1 спрямо математическото очакване} \} \\ H_1 := \{ \beta_1 \neq 0 \iff \text{Извадка 2 се различава значително от Извадка 1 спрямо средното} \} \]

# изпълнение на модела обясняващ y срямо индикаторна променлива
lm(y ~ x) %>%
  # обобщение от което да рагледаме F-статистиката и нейното pvalue
  summary()

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8880 -0.7973 -0.1254  0.6250  2.3124 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.0123     0.1469   68.17   <2e-16 ***
## x             4.8104     0.2077   23.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.039 on 98 degrees of freedom
## Multiple R-squared:  0.8455, Adjusted R-squared:  0.8439 
## F-statistic: 536.4 on 1 and 98 DF,  p-value: < 2.2e-16

Ако разширим модела си с повече от две извадки, изпълнявайки \(F\)-тест, всъщност ще питаме дали има извадка чиято популация драстично се различава от тази на другите спрямо средните.

\[ H_0 := \{ \beta_1 = \beta_2 = ... = \beta_p = 0 \iff Y = \beta_0 \iff \text{Всички извадки идват от една и съща популация} \} \\ H_1 := \{ \exists j : \beta_j \neq 0 \iff \text{Има извадка чието разпределение драстично се различава от останалите} \} \]

# общ брой елементи
n <- 250

# дисперсия (трябва да е обща спрямо модела)
sigma <- 6

# петте различни извадки
y1 <- rnorm(n, mean = 13, sd = sigma)
y2 <- rnorm(n, mean = 16, sd = sigma)
y3 <- rnorm(n, mean = 9, sd = sigma)
y4 <- rnorm(n, mean = 6, sd = sigma)
y5 <- rnorm(n, mean = 13, sd = sigma)

# конкатениране на y променливите в eдно y
y <- c(y1, y2, y3, y4, y5)

# индикаторни променливи които ни казват в коя извадка сме
x1 <- rep(c(0, 1, 0, 0, 0), each = n) # извадка 2
x2 <- rep(c(0, 0, 1, 0, 0), each = n) # извадка 3
x3 <- rep(c(0, 0, 0, 1, 0), each = n) # извадка 4
x4 <- rep(c(0, 0, 0, 0, 1), each = n)  # извадка 5

x <- rep(1:5, each = n)

groups_data <- data.frame(
  y = y,
  x1, x2, x3, x4,
  sample = as.factor(rep(paste0("Извадка ", 1:5), each = n))
)

## Warning: Font 'Pacifico' isn't in the extrafont font list (but it may still
## work). If recently installed, you can try running `extrafont::font_import()`.

Ако го разглеждаме от гледна точка на линейна регресия с много променливи (индикатори за съответните групи), оценките за \(\beta_1, .., \beta_p\) ще са разликите между средните на групи 2,3,…, р и първата група (визуално това са наклони в отделни линии). В долната графика, пресечната линия е средното на първата група (оценката за _0).

За да изпълним ANOVA тест който проверява дали някоя от групите значително се различава от останалите (нейния коефициент значително обяснява модела), пускаме линейна регресия обясняваща y спрямо съответните индикаторни променливи \(x_1, ..., x_4\).

lm(y ~ x1 + x2 + x3 + x4, data = groups_data) %>%
  summary()

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = groups_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.0861  -3.9651   0.2251   3.9667  18.0560 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.2239     0.3676  35.975  < 2e-16 ***
## x1            2.6733     0.5198   5.142 3.15e-07 ***
## x2           -4.0030     0.5198  -7.700 2.75e-14 ***
## x3           -7.3097     0.5198 -14.061  < 2e-16 ***
## x4           -0.5039     0.5198  -0.969    0.333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.812 on 1245 degrees of freedom
## Multiple R-squared:  0.2633, Adjusted R-squared:  0.261 
## F-statistic: 111.3 on 4 and 1245 DF,  p-value: < 2.2e-16

# Проверка за коефициентите и нашите разсъждения
c(
  "beta0" = mean(y1),
  "beta1" = mean(y2 - y1),
  "beta2" = mean(y3 - y1),
  "beta3" = mean(y4 - y1),
  "beta4" = mean(y5 - y1)
)

##      beta0      beta1      beta2      beta3      beta4 
## 13.2238764  2.6732943 -4.0029761 -7.3096923 -0.5038667

Предположенията на ANOVA теста си остават същите като при линейна регресия и Т-тест - данни с нормално разпределение и еднаква дисперсия между групите.

Бихме могли да изпълним теста и така

# Функцията stack създава помощна факторна променлива и конкатенира y променливите в една
anova_df <- stack(data.frame(y1, y2, y3, y4, y5))

anova_df %>%
  # Линеен модел
  lm(values ~ ind, data = .) %>%
  # подаваме линейния модел на anova функцията
  anova()

## Analysis of Variance Table
## 
## Response: values
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## ind          4  15034  3758.5  111.27 < 2.2e-16 ***
## Residuals 1245  42056    33.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Можем и с функцията oneway.test

oneway.test(values ~ ind, anova_df, var.equal = TRUE)

## 
##  One-way analysis of means
## 
## data:  values and ind
## F = 111.27, num df = 4, denom df = 1245, p-value < 2.2e-16

Така получения тест бихме могли да използваме за да проверим дали определен фактор (група) променя значително цялата съвкупност.

Примери

• Проверка дали различните методи на учене играят роля в оценката
• Проверка дали ученето в различни университети играе роля в заплащането
• Проверка дали различния тип гориво играе роля върху това колко изразходва колата на 100км
• Проверка дали различни раси получават еднакво заплащане

МANOVA (Многопроменлив Анализ на Вариациите)

Често се налага обаче да проверим дали факторите (групите) имат ефект върху не една променлива, a две променливи. Пример за това би бил въпроса дали различни видове лекарства имат значителен ефект както върху хемоглобина, така и върху хематокрита.

MANOVA Пример

Изпълнението на MANOVA теста става с функцията manova, като като аргумент подаваме матрица от двата отклика manova(cbind(y1, y2) ~ x, data = some_data)

two_response_groups_data %>%
  # изпълняваме функцията manova със Y = cbind(y1, y2)
  manova(cbind(y1, y2) ~ x, data = .) %>%
  # обобщение на модела
  summary()

##             Df Pillai approx F num Df den Df    Pr(>F)    
## x            4 1.1047   384.06      8   2490 < 2.2e-16 ***
## Residuals 1245                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Виждаме, че F-статистиката е далеч над 1 и pvalue-то е изключително малко, откъдето отсъждаме, че определено различните групи имат ефект върху двата отклика.

За да проверим, колко точно групите афектират двата отклика \(y_1\) и \(y_2\), можем да използваме функцията summary.aov за да изпълним поединична ANOVA

two_response_groups_data %>%
  # изпълняваме функцията manova със Y = cbind(y1, y2)
  manova(cbind(y1, y2) ~ x, data = .) %>%
  # обобщение на модела експлицитно за двата отклика
  summary.aov()

##  Response y1 :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## x              4  15034  3758.5  111.27 < 2.2e-16 ***
## Residuals   1245  42056    33.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response y2 :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## x              4 22993.3  5748.3  5946.2 < 2.2e-16 ***
## Residuals   1245  1203.6     1.0                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA с корелирани извадки (разширение на т-тест)

При ANOVA с корелирани извадки, отново тестваме за това дали извадките са еднакви или се различават по тяхното средно, но този път разглеждаме случая когато бихме могли да имаме зависимост между отделните групи.

Пример за това би бил тестване за времето на действие на различни лекарства върху няколко различни пациенти.

patients_df <- data.frame(
  patient = as.factor(rep(1:5, each = 4)),
  drug = as.factor(rep(1:4, times = 5)),
  response_time = c(
    41, 39, 27, 45,
    25, 29, 21, 33,
    35, 31, 29, 41,
    49, 45, 31, 55,
    37, 39, 25, 41
  )
)

patients_df

##    patient drug response_time
## 1        1    1            41
## 2        1    2            39
## 3        1    3            27
## 4        1    4            45
## 5        2    1            25
## 6        2    2            29
## 7        2    3            21
## 8        2    4            33
## 9        3    1            35
## 10       3    2            31
## 11       3    3            29
## 12       3    4            41
## 13       4    1            49
## 14       4    2            45
## 15       4    3            31
## 16       4    4            55
## 17       5    1            37
## 18       5    2            39
## 19       5    3            25
## 20       5    4            41

# Error(patients) задава променливата върху която гледаме двустранно
# преди и слсед лекарството
# У drug 
aov(response_time ~ drug + Error(patient), data = patients_df) %>%
  # Обобщение на модела
  summary()

## 
## Error: patient
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  4  680.8   170.2               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## drug       3  698.2   232.7   24.76 1.99e-05 ***
## Residuals 12  112.8     9.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1