library(readr)
library(vcd)
library(dplyr)
library(tidyverse)
library(BSDA)
library(car)
home = read.csv("C:/future/Statistica/solo/home_data.csv")
gen = read.csv("C:/future/Statistica/solo/genetherapy.csv")
atherosclerosis = read.csv("C:/future/Statistica/solo/atherosclerosis.csv")
birds = read.csv("C:/future/Statistica/solo/birds.csv")
titanic <- read.csv("https://stepic.org/media/attachments/course/524/train.csv")
titanic <- mutate(titanic, 
                  Survived = factor(Survived, labels = c("No", "Yes")), 
                  Pclass = factor(Pclass, labels = c("First", "Second", "Third")), 
                  Sex = factor(Sex, labels = c("Female", "Male")))

STATISTICS

1. Введение

1.1 Выборки

  1. Простая случайная выборка (simple random sample)
  2. Стратифицированная выборка (stratified sample) Перед тем как извлекать данные в выборку, изначально распределеить генеральную совокупность на страты, а только потом из этих стртат случайным образом вытаскивать данные
  3. Групповая выборка (claster sample) Изначально разбиваем на кластеры, а затем из этих групп выбираем (не из всех) случайно данные. Например, при предположении что районы СПБ одинаковы по своим хар-ам, будем выбирать данные из любых нескольких районов

1.2 Типы переменных

  • Количественные
    • Непрерывные
    • Дискретные
  • Номинальные

1.3 Меры центральной тенденции

Мода

Это значение измеряемого признака, которое встречается максимально часто

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)]
}

Mode(home$grade)
## [1] 7
ggplot(data = home) +
  geom_bar(aes(x = factor(grade))) +
  theme_bw()

Медиана

Это значение признака, которое делит упорядоченное множество данных пополам

median(home$price)
## [1] 450000
ggplot(data = home) +
  geom_histogram(aes(x = price/10000), binwidth = 5) +
  scale_x_continuous(breaks = seq(from = 0, to = 1000, by = 45)) +
  geom_vline(xintercept=median(home$price)/10000, color="red",  linetype = "dashed", size = 1) + theme_bw()+
  ggtitle("median")

Среднее значение

Сумма всех значений на число элементов

Если ваше распределение симметрично и не имеет выбросов, можно использовать моду и медиану и среднее значение, они дадут +- то же самое значение. Если есть ассиметрия, среднее значение будет подвержено выбросам

Свойства среднего

  • Каждое значение выборки + C = Среднее + C
  • Каждое значение выборки * C = Среднее + C
  • СУММ(Каждое значение выборки - среднее) = 0

1.4 Меры изменчивости

Размах (Range)

Разность максимального и минимального значения. Выбросы могут искажать размах.

range(home$price)
## [1]   75000 7700000

Дисперсия (Variance)

Посмотреть насколько в среднем наши значения(каждое) отклоняется от среднего значения по выборке - рассчитать дисперсию

## Формула для расчета дисперсии генеральной совокупности
var_pop <- function(x) {
  mean((x - mean(x))^2)
}

set.seed(1234)
sample = sample_n((home %>% select(price)), 3000)


ggplot(data = sample) +
  geom_histogram(aes(x = price/10000), binwidth = 5) + theme_bw() +
  ggtitle("Распределение случайной выборки")

ggplot(data = home) +
  geom_histogram(aes(x = price/10000), binwidth = 5) +
  scale_x_continuous(breaks = seq(from = 0, to = 1000, by = 45)) +
  geom_vline(xintercept=sqrt(var_pop(home$price))/10000, color="red",  linetype = "dashed", size = 1) +
  geom_vline(xintercept=sqrt(var(sample$price))/10000, color="blue", linetype = "dashed", size = 1) +
 theme_bw() +
  ggtitle("Сигма(red) генеральной совокупности и SD(blue")

asf

sd (standart deviation) - сигма для выборки. Используя степени свободы при расчете на выборке, мы получим значение более близкое к значению, которое получилось бы в генеральной совокупности.

Свойства дисперсии:

Квартили распределения

Это три точки, которые делят наши данные на 4 равные части

Box-Plot

test = as.data.frame(quantile(home$price/10000, prob=c(.25)))
colnames(test) = "test1"
ggplot(data = home) +
  geom_boxplot(aes(y = price/10000)) + theme_bw() + scale_y_continuous(breaks = seq(from = 0, to = 800, by = 40)) +
  coord_flip() +
  geom_hline(yintercept=quantile(home$price/10000, prob=c(.25)), color="red",  linetype = "dashed", size = 1) +
  geom_hline(yintercept=quantile(home$price/10000, prob=c(.75)), color="red",  linetype = "dashed", size = 1) +
  geom_hline(yintercept=median(home$price/10000), color="red",  linetype = "dashed", size = 1) +
  geom_hline(yintercept=mean(home$price/10000), color="blue",  linetype = "dashed", size = 1)+
  geom_hline(yintercept=(32.195 - 1.5*32.305), color="green",  linetype = "dashed", size = 1) +
  geom_hline(yintercept=(64.5 + 1.5*32.305), color="green",  linetype = "dashed", size = 1)

### 1.5 Нормальное распределение

Это унимодальное и симметричное распределение.

Z стандартизация - приведение среднего к нулю, а сигму к единице

Оно позволяет нам ответить на вопрос, какой процент наблюдений лежит в абсолютно любом интересующем нас диапазоне.

1.6 Центральная предельная теорема

SE(Standart Error) - стандартное отклонение распределенния средних значений выборки. Показывает, насколько в среднем средние значения выборок отличаются от среднего в генеральной совокупности.

Чем больше размер выборки, тем меньше SE. Чем меньше стандартное отклонение генеральной совокупности, тем меньше SE.

Замечание к ЦПТ: Если число наблюдений выборки >30 и выборка репрезентативна, то в формуле SE вместо стандартного отклонения генеральной совокупности можно использовать стандартное отклонение любой выборки. То есть можно предсказать, какую стандартную ошибку будет иметь распределение выборочного среднего.

Замечание к ЦПТ: Пусть есть признак, распределенный КАК УГОДНО с некоторым средним и некоторым стандартным отклонением. Тогда, если мы будем выбирать из этой совокупности выборки объема n, то их средние тоже будут распределены нормально со средним равным среднему признака в ГС и отклонением стандартным отклонением, se, формулу которого мы уже приводили в прошлых шагах.

Замечание к ЦПТ: Распределение должно обладать конечной дисперсией, да и так бывает, бывают распределения, у которых дисперсия и вовсе не определена!

Замечание к ЦПТ: Распределение должно обладать конечной дисперсией, да и так бывает, бывают распределения, у которых дисперсия и вовсе не определена!

1.7 Доверительные интервалы

Обычно нам неизвестна генеральная совокупность и необходимое среднее ГС приходится искать через доверительные интервалы на выборочных данных. Если бы многократно повторяли наш эксперимент, то средние выборок распределились бы нормальным образом вокруг среднего и со стандартной ошибкой среднего. Если бы мы для каждого выборочного среднего рассчитали, например, интервал 1,96se, то в 95% случаев выборочные средние включили бы среднее по генеральной совокупности.

test = z.test(gen$expr, sigma.x = sqrt(var_pop(gen$expr)), conf.level = 0.95)
test$conf.int
## [1] 94.90672 97.72661
## attr(,"conf.level")
## [1] 0.95

1.8 p-уровень значимости

Сначала допускаем, что верна нулевая гипотеза. Чем меньше p уровень значимости, тем больше оснований отклонить нулевую гипотезу. Односторонний критерий используется, когда отклонение в одну сторону невозможно. p-уровень значимости:

  • Не позволяет сказать, с какой вероятностью верна нулевая гипотеза.

  • Ничего не говорит о силе различий

  • Если p > 0.05 не позволяет отклонить нулевую гипотезу

  • Он ничего не говорит о правильности получаемых результатов

Ошибка 1 рода - отклонили нулевую, хотя она была верна

Ошибка 2 рода не отклонили нулевую, хотя была верна альтернативная

2. Сравнение средних

2.1 t-распределение

Чем меньше выборка, тем больше отклонения от среднего генеральной совокупности.

Если число наблюдений невелико и сигма неизвестна, используется распределение Стьюдента. Основное отличие от нормального распределения - у него более высокие хвосты распределения. Т.е. в отрезке +-2 сигмы будет лежать большее количество наблюдений. Чем больше степеней свободы df(n-1), тем больше данное распределение похоже на нормальное В отличие от нормального, форма t распределение изменяется в зависимости от степеней свободы.

2.2 Сравнение двух средних; t-критерий Стьюдента

Важно, чтобы дисперсии внутри наших групп были примерно одинаковы. Это, так называемое, требование гомогенности дисперсий. Проверить можно, используя критерии Фишера и Ливена. Важно Если n < 30, то очень важным требованием является нормальность распределения выборок. Если n > 30 то все должно быть норм, но лучше проверить.

Таблица для t-значения по df

shiny калькулятор для t-значения еще один

2.3 Проверка распределения на нормальность, QQ-Plot

library(car)
home$price = as.numeric(home$price)
qqPlot(home$sqft_living)

## [1] 12778  7253

Если распределение нормальное, значения должны идти +- по линии. 0 на графике - середина распределения. Спрва от нуля - значения выше ожидаемых. Слева от нуля - значения также выше ожидаемых.

Тесты для проверки нормальности: Колмагорова Shapiro - Wilk

2.4 Однофакторный дисперсионный анализ

Очень часто в экспериментах и исследованиях возникает необходимость сравнить несколько групп между собой. В таком случае мы можем применять однофакторный дисперсионный анализ. Та переменная, которая будет разделять наших испытуемых или наблюдения на группы (номинативная переменная с нескольким градациями) называется независимой переменной. А та количественная переменная, по степени выраженности которой мы сравниваем группы, называется зависимая переменная.

Shiny для F статистики

SST Состоит из SSB(Междугрупповая сумма квадратов) и SSW(Внутригрупповая сумма квадратов)

Если большинство сумм квадратов обеспечивается SSB, то можно говорить о значительных различиях между группами

library("ggpubr")
library(BSDA)
a = aov(expr ~ Therapy, data = gen)
summary(a)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Therapy      3  560.7  186.91   8.037 0.000152 ***
## Residuals   56 1302.3   23.25                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(ggplot2)
ggplot(data = gen) +
  geom_boxplot(aes(x=factor(Therapy), y = expr)) +
  theme_bw()

Основная идея Если бы мы многократно повторяли эксперимент, и нулевая гипотеза об отсутствии различий была верна, то числитель формулы стремился бы к нулю вместе с итоговым зачением F. Тогда большинство F значений были бы маленькими и находились бы в левой стороне распрределения фишера. Значение фишера всегда принимает положительные значения, поэтому вероятность всегда смотрим одностороннюю.

2.5 Множественные сравнения в ANOVA

Используем, когда нужно много групп. Для тестов - используются всякие поправки. Например, поправка Бенферонни предлагает делить alpha на количество парных сравнений между группами.

Поправка Tukey HSD Похожа на t.test, но тут по другому расчитывается standrad error

2.6 Многофакторный ANOVA

summary(aov(expr ~ age + dose + age:dose, data = atherosclerosis))
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## age          1  197.5  197.45   7.450 0.00831 **
## dose         1   16.9   16.91   0.638 0.42755   
## age:dose     1    0.9    0.93   0.035 0.85227   
## Residuals   60 1590.3   26.50                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(var4 ~ hormone + sex + hormone:sex, data = birds))
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## hormone      1    0.8    0.85   0.087 0.76965   
## sex          1    0.1    0.12   0.012 0.91232   
## hormone:sex  1   89.5   89.48   9.136 0.00368 **
## Residuals   60  587.7    9.79                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

При сравнении по двум независимым переменным: SS_TOTAL = SSW + SSB_A + SSB_B + SSB_A * SSB_B

Требования к данным:

  • Нормальность распределения зависмой в каждой из групп

  • Гомогенность дисперсий (примерно одинакова в каждой из групп). Проверять, например, боксплотом. Проверять можно тестом Levene’s test for Homogeneity

При больше 50 наблюдений, ANOVA устойчива к нарушению данных требований

2.6 АБ Тесты

На степе есть ссылки

3. Корреляция и регрессия

3.1 Понятие корреляции

cov > 0 положительная взаимосвязь

Формула - коэфф. корр. Пирсона.

Shared variance - коэффициент детерминации. Очень важное понятие, которое будет в теме регрессиии. Это часть изменчивости одной переменной, которая обусловлена тем, что она взаимосвязана с другой переменной. Показывает, в какой степени дисперсия одной переменной обусловлена “влиянием” другой переменной

Для рассчета p-value для корреляции используем t критерий при df = N-2

3.2 Условия применения коэффициента корреляции

Необходимо, чтобы характер взаимосвязи был линейный и монотонный.

Желательно, чтобы переменные имели нормальное распределение, поскольку сама формула критерия связана со средними значениями.

Если условия не удовлетворяются, можно использовать непараметрические параметры: Spearman correlation

Ошибка коррелляции: она не говорит о причинно-следственной зависимости.

Коэффициент показывает силу и направление взаимосвязи двух количественных переменных.

3.3 Регрессионный анализ

Линия регрессии. На Y - зависимая. X - независимая(предиктор).

Смотрим, как одна переменная помогает предсказать зависимую переменную. Желательно, чтобы наша линия прохоидла через самый центр облака этих точек.

Метод наименьших квадратов - метод нахождения оптимальных параметров линейной регрессии, таких, что сумма квадратов ошибок (остатков) была минимальна.

Можно построить математическую модель, отображающую зависимость переменных

3.4 Гипотеза о значимости взаимосвязи и коэффициент детерминации

Если cor = 0, то линия регрессии проходит горизонтально по Y_bar.

При нулевых гипотезах угол наклона линии равняется нулю.

Здесь и появляется t критерий, готорый говорит о том, что при многократном повторении эксперимента коэффиценты угла наклона распределились бы нормальным образом.

Где SS_total - расстояние от точки до красного пунктира.

Если SS_res = SS_total, то R^2 стремится к нулю и наша модель никак не объясняется изменчивостью зависимой переменной.

3.5 Условия применения линейной регрессии с одним предиктором

Требования:

  • Линейная взаимосвязь X и Y

  • Нормальное распределение остатков

  • Гомоскедастичность - постоянная изменчивость остатков на всех уровнях независимой переменной

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

Пример регрессий со степа

Сначала смотрим распределение остатков, scatter plot и потом решаем о приминении регрессии.

3.6 Применение регрессионного анализа и интерпретация результатов

Когда среднее образование в штате увеличивается на 1%, уровень бедности уменьшается на 0.62%

3.7 Задача предсказания значений зависимой переменной

Иногда регрес. прямую называют линией тренда. Мы можем предсказывать значения, которых у нас на самом деле нет.

3.8 Регрессионный анализ с несколькими независимыми переменными

Можно включать сразу несколько независимых переменных в модель.

Расчитывается исправленый R-squared из за множественности сравнений.

3.9 Выбор наилучшей модели

Когда у нас сильная коррелляция между предикторами, то становится сложно рассчитать необходимые показатели. Когда строим регр. модель, то не факт, что чем больше переменных - тем лучше модель.

4. Анализ номинативных данных

4.1 Расстояние Пирсона

Проверка гипотез между наблюдаемым и теоретическим распределением частот.

Observed - наблюдаемые значения. Expected значения согласно вероятности в нулевой гипотезе.

Чем больше значение ожидаемой частоты, тем большее отклонение от нее можно считать возможными и ожидаемыми. Нужно как то нормировать на общий размер выборки

Распределение данного коэффицента будет подчиняться специальному распределению Chi squared.

4.2 Распределение Pearson Chi-Squared

library(dplyr)

library(ggplot2)

o <-  replicate(10000, sum(sample(x = c(0,1), size = 60, replace = TRUE)))
df <- data.frame(o)

df <- 
  df %>% 
  rowwise() %>% 
  mutate(chis = chisq.test(c(o,60-o))$statistic)

ggplot(df, aes(x=chis)) + 
  geom_histogram(fill="white", color="black",bins=30) + theme_bw()

Если бы многократно повторяли наш эксперимент, то чаще всего отклонения будут незначительны.

Распределение хи-квадрат с K степенями свободы - это распределение суммы квадратов k независимых стандартных нормальных случайных величин.

Чем больше степеней свободы, тем больше кумулятивный эффект от сложения их хи-квадратов. И, поэтому, при очень большом кол-ве степеней свободы хи-квадрат будет стремиться к нормальному распределению.

4.3 Расчет p-уровня значимости

Нулевая гипотеза - что распределение частот распределено не равномерно.

При многократном повторении эксперимента хи-квадраты распределяются по закону хи-квадрата.

### 4.4 Анализ таблиц сопряженности

t <- rbind(c(10, 6), c(5, 15))
chisq.test(t)$expected
##          [,1]      [,2]
## [1,] 6.666667  9.333333
## [2,] 8.333333 11.666667

Здесь используем df = 1, поскольку, зная, сколько всего человек в каждой из групп, можно узнать остальные ячейки, даже если мы знаем только данные одной ячейки.

В R при расчете таблицы 2x2 автоматически используется Поправка Йетса.

students = rbind(c(15,9), c(11,6))
chisq.test(students)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  students
## X-squared = 1.2684e-31, df = 1, p-value = 1

Свойство про больше 5 - если меньше, то меньше сумма нормально распределенных величин и распределение будет совсем другим.

Интерпетация результатов при хи-квадрате - он не рассчитывает никакого показателя про отдельно связанную ячейку. Это чревато тем, что при большой таблице сопряженности результаты не говорят о том, в какой ячейке было отклонение.

t = rbind(c(20,15), c(11,12), c(7,9))
chisq.test(t)
## 
##  Pearson's Chi-squared test
## 
## data:  t
## X-squared = 0.95441, df = 2, p-value = 0.6205
chisq.test(t)$expected
##           [,1]      [,2]
## [1,] 17.972973 17.027027
## [2,] 11.810811 11.189189
## [3,]  8.216216  7.783784
t = rbind(c(18,7), c(6,13))
chisq.test(t)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  t
## X-squared = 5.5772, df = 1, p-value = 0.0182
library(vcd)
a = as.data.frame(rbind(c(18,7), c(6,13)))
colnames(a) <- c("Yes", "No")
row.names(a)<- c("Placebo","Aspirin")
file_path <- "http://www.sthda.com/sthda/RDoc/data/housetasks.txt"
housetasks <- read.delim(file_path, row.names = 1)
dt <- as.table(as.matrix(a))
assoc(dt, gp=shading_max, main= "hz")

4.5 Точный критерий Фишера

Не нужно забывать о требованиях к Хи-Квадрату.

Когда нужен намного меньший объем выборы, для этого существует такой критерий. Работает на экстремально маленьких размерах выборки.

Он решает тест вручную. Считает вероятность получить такие переменные вручную.

С помощью зеленой формулы появляется возможность получить вероятность получить данный набор данных

Рассчитаем вероятность получить такие данные, при условии что нет никакой разницы выздоровления между двумя лекарствами

На скрине показан смысл: алгоритм делает перестановки данных и рассчитывает вероятность получить такие данные случайно. Если p-value меньше необходимого, то отвергаем гипотезу о том, что вероятность вылечиться Лекарством 1 или Лекарством 2 одинаковы.

fisher.test(rbind(c(3,1), c(1,3)))
## 
##  Fisher's Exact Test for Count Data
## 
## data:  rbind(c(3, 1), c(1, 3))
## p-value = 0.4857
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    0.2117329 621.9337505
## sample estimates:
## odds ratio 
##   6.408309

5. Логистическая регрессия

Логистическая регрессия применяется при номинативных переменах и одной независимой.

При ее помощи можно проверять гипотезы о взаимосвязи нескольких номинативных переменных.

Единственное ограничение - зависимая переменная должна иметь две градации, но и это можно будет обойти.

Непараметрические тесты можно использовать и при параметрических данных, но популярность низка.

Основная идея в том, что зависимая п. - номинативная с двумя градациями. Независимые - как номинативные, так и нет с любым кол-во градаций.

В анализе шансов может помочь натуральный логарифм. Значение шансов меньше единицы будут давать отрицательные значение, больше -> положительные.

Выше решали проблему: была математическая несостыковка. В правой части линейная комбинация предикторов, а слева была только две градации номинативной переменной. После применения логита теперь в левой части логарифм шансов.

  • В качестве предикторов могут выступать как количественные, так и номинативные переменные.

  • После преобразования вероятности в логарифм шансов предсказания регрессионной модели могут быть на промежутке от минус до плюс бесконечности.

5.1 Модель без предикторов. Intercept only model

intecept - число, которое максимально удачно спрогнозирует нам модель. Такое число - натуральный логарифм шансов получить положительный исход зависимой переменной для случайно выбранного (пассажира титаника).

Это также можно интерпетировать как проверку нулевой гипотезы о равномерном распределении градаций номинальной зависимой переменной. (Pr(>|z|))

Когда мы говорим про лог. регрессию, все коэффициенты, которые мы получаем - обычно это значения логарифма.

mosaic(~ Sex + Survived, data=titanic) 

simple_fit <- glm(Survived ~ 1, titanic, family = "binomial")
summary(simple_fit)
## 
## Call:
## glm(formula = Survived ~ 1, family = "binomial", data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9841  -0.9841  -0.9841   1.3839   1.3839  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.47329    0.06889   -6.87  6.4e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance: 1186.7  on 890  degrees of freedom
## AIC: 1188.7
## 
## Number of Fisher Scoring iterations: 4

std.Error - если бы многократно повторяли эксперимент и выбирали выборку с такой же номинативной переменной и смотрели бы значимые различия логарифма шансов от нуля. Если верна H0, то мы не должны часто получать отклонения от данных в нулевой гипотезе.

То логарифм шансов распределился бы нормальным образом, поскольку мы бы совершенно случайно получали отклонения в разные стороны.

sd - это и есть стандартная ошибка этого нормального распределения.

ln(1) = 0

Если мы разделим Intercept / sd, то мы ответим на вопрос, сколько стандартных отклонений поместилось в значение от нуля до полученного коэффициента.

И, поэтому, поскольку, это нормальное распределение, можем получить p-value для таких или более выражденных значений intercept.

По сравнению с chi-squared, такой intercept помогает сделать больше интерпретаций вывода.

Если подставить экспоненту - избавимся от логарифма.

При помощи нормального распределения описывается распределение коэффициентов логистической регрессии, если верна нулевая гипотеза

Модель с одним номинативным предиктором

fit1 <- glm(Survived ~ Sex, titanic, family = "binomial")
summary(fit1)
## 
## Call:
## glm(formula = Survived ~ Sex, family = "binomial", data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6462  -0.6471  -0.6471   0.7725   1.8256  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0566     0.1290   8.191 2.58e-16 ***
## SexMale      -2.5137     0.1672 -15.036  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance:  917.8  on 889  degrees of freedom
## AIC: 921.8
## 
## Number of Fisher Scoring iterations: 4
coef(fit1)
## (Intercept)     SexMale 
##    1.056589   -2.513710

Получился только SexMale. Почему? Это так называемое отношение шансов

table(titanic$Survived, titanic$Sex)
##      
##       Female Male
##   No      81  468
##   Yes    233  109

У женщин больше выжило, у мужчин больше погибло.

Шансы выжить для мужчин и женщин:

odds_male <- 93 / 360
odds_female <- 197 / 64

odds_male
## [1] 0.2583333
odds_female
## [1] 3.078125

Посмотрим логарифмы шансов:

log(odds_male)
## [1] -1.353505
log(odds_female)
## [1] 1.124321

Сравним шансы выжить для женщин и для мужчин (поделенные шансы мужчин на шансы женщин):

odds_ratio <- odds_male / odds_female
log(odds_ratio)
## [1] -2.477825

Видно, что их поделенные шансы 0.08392555. Это значит, что у женщин шансы выше.

А логарифм отношения их шансов это и есть SexMale.

Еще раз почему нету SexFemale?

Всегда, когда выполняется регрессионная модель, компьютер за нас переписывает данные в немножко другой формат, где из одного столбика номинативных переменных с градациями делает столбики по каждой градацией.

Можно оставить один столбик, потому что зная значение одного столбика, всегда можно узнать значение другого столбика, поскольку градаций две(или может быть несколько)

Для женщин sex_male будет 0, поэтому при математическом вычислении останется только intercept.

Коэффициент -2.74 означает, что log(odds_m/odds_f) < 0

Лог регрессия позволяет сравнивать шансы.

Intercept обычно выбирается по названию фактора(алфавит и т.д.)

Факторы + intercept - цена перехода, например, человека из 1 класса титаника в третий класс.

Модель с двумя номинативными предикторами

fit2 <- glm(Survived ~ Sex * Pclass, titanic, family = "binomial")
summary(fit2)
## 
## Call:
## glm(formula = Survived ~ Sex * Pclass, family = "binomial", data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6248  -0.5853  -0.5395   0.4056   1.9996  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.4122     0.5868   5.815 6.06e-09 ***
## SexMale               -3.9494     0.6161  -6.411 1.45e-10 ***
## PclassSecond          -0.9555     0.7248  -1.318  0.18737    
## PclassThird           -3.4122     0.6100  -5.594 2.22e-08 ***
## SexMale:PclassSecond  -0.1850     0.7939  -0.233  0.81575    
## SexMale:PclassThird    2.0958     0.6572   3.189  0.00143 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance:  798.1  on 885  degrees of freedom
## AIC: 810.1
## 
## Number of Fisher Scoring iterations: 6

intercept - логарифм шансов выжить для женщин в первом классе

Взаимодействие номинативных признаков

coef(fit2)
##          (Intercept)              SexMale         PclassSecond 
##            3.4122472           -3.9493901           -0.9555114 
##          PclassThird SexMale:PclassSecond  SexMale:PclassThird 
##           -3.4122472           -0.1849918            2.0957553

Взаимодействие между факторами говорит о том, что взаимосвязь между полом и выживаемости проявляется по разному в зависимости от класса билета.

table(titanic$Survived, titanic$Pclass , titanic$Sex)
## , ,  = Female
## 
##      
##       First Second Third
##   No      3      6    72
##   Yes    91     70    72
## 
## , ,  = Male
## 
##      
##       First Second Third
##   No     77     91   300
##   Yes    45     17    47
# (Intercept) 
female_p1_odds <- 82 / 3
log(female_p1_odds)
## [1] 3.308107
# Sexmale  
male_p1_odds <- 40  /  61 
log(male_p1_odds)
## [1] -0.4219944
log(male_p1_odds / female_p1_odds )
## [1] -3.730101
# PclassSecond
female_p2_odds <- 68  /  6 
log(female_p2_odds / female_p1_odds )
## [1] -0.8803587
# PclassThird
female_p3_odds <- 47  /  55 
log(female_p3_odds / female_p1_odds )
## [1] -3.465293
# SexMale:PclassSecond
male_p2_odds <- 15 / 84
log(male_p2_odds / female_p2_odds ) - log(male_p1_odds / female_p1_odds )
## [1] -0.4204135
#Sexmale:factorThird 
male_p3_odds <- 38 / 215
log(male_p3_odds / female_p3_odds ) - log(male_p1_odds / female_p1_odds )
## [1] 2.154235

anova(fit1, fit2, test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: Survived ~ Sex
## Model 2: Survived ~ Sex * Pclass
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       889      917.8                          
## 2       885      798.1  4   119.71 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Survived
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                         890    1186.66              
## Sex         1  268.851       889     917.80 < 2.2e-16 ***
## Pclass      2   90.916       887     826.89 < 2.2e-16 ***
## Sex:Pclass  2   28.791       885     798.10 5.598e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Модель с разными переменными:

fit3 <- glm(Survived ~ Sex + Pclass + Age, titanic, family = "binomial")
anova(fit3, test="Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Survived
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     713     964.52              
## Sex     1  213.816       712     750.70 < 2.2e-16 ***
## Pclass  2   78.269       710     672.43 < 2.2e-16 ***
## Age     1   25.148       709     647.28 5.311e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Когда нужно использовать непараметрические методы и почему?

Чаще всего связано с выбросами

Такие ошибки не будут сильно сказываться на ошибке t-test’a.

Когда мы используем критерий t Стьюдента:

  • Число 30 не обладает магическим свойством.

  • Если распределение отлично от нормального, можно играть со статистикой в испорченный телефон.

Важно смотреть на нормальность распределения!

Если данные отличаются от нормального распределения, дисперсии гомогенны - применяем непараметрический

Критерий Манна-Уитни не сравнивает медианы!

Если одновременно нарушается требование гомогенности дисперсий и нормальности распределения при дисперсионном анализе, следует использовать его непараметрический аналог - Критерий Краскела - Уоллиса

Кластерный анализ и метод главных компонент

Кластерный анализ методом k - средних

Они сами ищут “правильный ответ” в наших данных.

Кластерный анализ бежит по данным и находит группы.

Метод главных компонент смотрит на переменные и ставит вопрос: можно ли сократить и сгруппировать наши данные?

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

library(ggplot2)

ggplot(iris, aes(Sepal.Length, Petal.Length, color = Species)) + 
geom_point()+
ggtitle("Данные о лепестках и чашечках")

Центроид в кластерах - так сказать геометрический центр точек. Ровна середина облака точек.

Центроид - точка, где x,y - средние значения по длине и ширине облака

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