2_stats.Rmd

---
editor_options: 
  chunk_output_type: console
---

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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

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

```{r, message=FALSE}
library(tidyverse)
library(lme4)
library(lmerTest)
# это стиль для ggplot (можно игнорировать эту команду)
theme_set(theme_bw()+
            theme(text = element_text(size = 16))) 
```

В качестве примера я использую данные из работы [Huttenlocher, Vasilyeva, Cymerman, Levine 2002], в которой авторы проанализировали 46 пар матерей и детей (возрастом от 47 до 59 месяцев, средний возраст ... а пасчитайте сами!). Они записали и затранскрибировали 2 часа каждого дня ребенка. Потом они посчитали количество именных групп на предложение у матерей и ребенка. Данные можно скачать командой:
`np_acquisition <- read.csv("https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/Huttenlocher.csv")`

```{r}
np_acquisition <- read.csv("https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/Huttenlocher.csv")
np_acquisition
```


## 1 Линейная регрессия

### 1.1 Простой случай

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

```{r}
np_acquisition %>% 
  ggplot(aes(mother, child))+
  geom_point()+
  labs(caption = "данные [Huttenlocher, Vasilyeva, Cymerman, Levine 2002]")
```

Простейшая статистическая модель будет выглядеть следующим образом:

$$y_i = \beta_0 + e_i$$

* $y_i$ --- множество ответов
* $\beta_0$ --- некоторая константа
* $e_i$ --- остатки/ошибки регресиионной модели

Так как любая регрессионная модель стремиться улучшить свои предсказания, лучшим коэфециентом $k$ будет тот, при котором $e_i$ будут минимальны. В нашей простой модели --- это понятное дело будет среднее:

```{r}
np_acquisition %>% 
  ggplot(aes(mother, child))+
  geom_hline(aes(yintercept = (mean(mother))), color = "blue")+
  geom_point()+
  labs(caption = "данные [Huttenlocher, Vasilyeva, Cymerman, Levine 2002]")
```

$$child_j = mean(child) + e_j$$

Если мы хотим включить переменную `mother` в нашу регрессию, то мы это можем сделать используя формулу прямой:

$$y_j = \beta_0 + \beta_1\times x_j + e_j$$

```{r}
np_acquisition %>% 
  ggplot(aes(mother, child))+
  geom_smooth(method = "lm", se = FALSE)+
  geom_point()+
  labs(caption = "данные [Huttenlocher, Vasilyeva, Cymerman, Levine 2002]")
```

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

$$child_j = \beta_0 + \beta_1\times mother_j + e_j$$

В R это делается при помощи функции `lm`:

```{r}
lm(child~mother, data = np_acquisition)
```

Теперь мы можем полностью записать формулу:

$$child_j = 0.7038 + 0.5452 \times mother_j + e_j$$

Результаты регрессии можно записать в переменную, а потом посмотреть даже статистическую значимость каждого из полученных коэффициентов:

```{r}
fit <- lm(child~mother, data = np_acquisition)
summary(fit)
```

На что смотреть:

* на статистическую значимость коэффициентов
* Adjusted R-квадрат
* на p-value всей модели

Теперь мы можем предсказывать! Правдо не обязательно все руками вбивать в формулу, например, для значения `mother = 1.69`:

```{r}
predict(fit, newdata = data.frame(mother = 1.69))
```


> 🤔 Скачайте [датасет с параметрами рассказов А. П. Чехова](https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/chekhov.csv): количество уникальных слов, длина рассказа. Постройте регрессионную модель предсказывающую количество уникальных слов на основе длины рассказа. Какой получился интерсепт (с точностью до 4 знаков после запятой)?

<form name="FormOne" onsubmit="return validateFormOne()" method="post">
<input type="text" name="answerOne">
<input type="submit" value="check">
</form><br>

> 🤔 Какой получился коэффициент при переменной длина рассказа (с точностью до 4 знаков после запятой)?

<form name="FormTwo" onsubmit="return validateFormTwo()" method="post">
<input type="text" name="answerTwo">
<input type="submit" value="check">
</form><br>

> 🤔 Что полученная модель предсказывает для рассказа размером 855 слов  (с точностью до 4 знаков после запятой)?

<form name="FormThree" onsubmit="return validateFormThree()" method="post">
<input type="text" name="answerThree">
<input type="submit" value="check">
</form><br>

### 1.2 Категориальные переменные

А что если одна из переменных у нас категориальная? В таком случае вводятся *dummy-переменные*. Рассмотрим наш вчерашний пример [из [Hau 2007]](https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/Hau.2007.csv): 

```{r}
homo <- read.csv("https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/Hau.2007.csv")
```

Попробуем предсказать длительность s на основании ориентации спикера. Для этого в регрссию вводят так-называемые *dummy-переменные*. Dummy-переменные принимают лишь два значения либо 1, либо 0. В нашем случае 1 --- гомосексуал, 0 --- гетеросексуал.

```{r}
homo %>% 
  ggplot(aes(orientation, s.duration.ms))+
  geom_point()
```

В таком случае наша формула примет следующий вид:

$$y_j = \beta_0 + \beta_1\times x + e_j$$

$$s\_duration_j = \beta_0 + \beta_1\times orientation_j + e_j$$

При описании гетеросексуалов переменная $orientation = 0$, тогда модель принимает вид: 
$$s\_duration_j = \beta_0 + \beta_1\times orientation_j + e_j = \beta_0 + \beta_1\times 0 + e_j = \beta_0 + e_j$$


При описании гомосексуалов переменная $orientation = 1$, тогда модель принимает вид:
$$s\_duration_j = \beta_0 + \beta_1\times orientation_j + e_j = \beta_0 + \beta_1\times 1 + e_j = \beta_0 + \beta_1 + e_j$$

```{r}
fit <- lm(s.duration.ms~orientation, data = homo)
summary(fit)
```

Таким образом эта модель возвращает всего два числа: гетеросексуалы --- $58.466$, гомосексуалы --- $58.466+5.517=63.983$.

**ОЧЕНЬ ВАЖНО:** dummy-переменных всегда (n-1). Т. е. если значений категориальной переменной 7, то dummy-переменных будет 6 и т. д.

**ОЧЕНЬ ВАЖНО:** сколько бы не было значений категориальных, регрессия выбирает одну категорию (reference level) и проводит сравнение ее со всеми. Т. е. категории не в интерсепте не сравниваются, однако можно поменять reference level при помощи типа переменных `factor`.

```{r}
homo$orientation <- factor(homo$orientation, levels = c("homo", "hetero"))
fit <- lm(s.duration.ms~orientation, data = homo)
summary(fit)
```

```{r, include=FALSE}
homo$orientation <- factor(homo$orientation, levels = c("hetero", "homo"))
```


> 🤔 Скачайте датасет, который содержит [данные по количеству согласных в языках мира](https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/phoible_ejectives.csv). Постройте регрессию предсказывающую количество согласных на основании наличия абруптивных. Приведите интерсепт модели (с точностью до 4 знаков после запятой):

```{r, include=FALSE}
ejectives <- read_csv("https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/phoible_ejectives.csv")
fit <- lm(total~have_ejectives, data = ejectives)
summary(fit)
```

<form name="FormFour" onsubmit="return validateFormFour()" method="post">
<input type="text" name="answerFour">
<input type="submit" value="check">
</form><br>

> 🤔 Приведите коэффициент при переменной `have_ejectives` (с точностью до 4 знаков после запятой):

<form name="FormFive" onsubmit="return validateFormFive()" method="post">
<input type="text" name="answerFive">
<input type="submit" value="check">
</form><br>


### 1.3 Множественная регрессия

В целом, вообще-то можно в предсказании использовать не одну переменную, а сразу много. Это сложнее визуализировать, но все остальное выглядит так же:

$$y_j = \beta_0 + \beta_1 \times x_{1j} + \dots + \beta_k \times x_{kj} + e_j$$

Кроме того на практике, исследователи сравнивают разные модели, выбирая модели с статистически значимыми предикторами или используя информационные критерии (самый распространенный --- критерий Акаике, который реализован в функции `AIC`)
  
> 🤔 В датасете про гомосексуалов попробуйте предсказать восприятие носителей как гомосексуалов на основании переменных длительность s,  средняя частота основного тона и размер диапозона частоты основного тона. Првиедите $R^2$ получившейся модели (с точностью до трех знаков после запятой):

```{r, include=FALSE}
summary(lm(perceived.as.homo.percent~s.duration.ms+average.f0.Hz+f0.range.Hz, data = homo))
# 0.238
```

<form name="FormSix" onsubmit="return validateFormSix()" method="post">
<input type="text" name="answerSix">
<input type="submit" value="check">
</form><br>

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

### 2.1 Введение
Мы хотим чего-то такого:
$$\underbrace{y}_{[-\infty, +\infty]}=\underbrace{\mbox{β}_0+\mbox{β}_1\cdot x_1+\mbox{β}_2\cdot x_2 + \dots +\mbox{β}_k\cdot x_k +\mbox{ε}_i}_{[-\infty, +\infty]}$$
Вероятность — (в классической статистике) отношение количества успехов к общему числу событий:
$$p = \frac{\mbox{# успехов}}{\mbox{# неудач} + \mbox{# успехов}}, \mbox{область значений: }[0, 1]$$
Шансы — отношение количества успехов к количеству неудач:
$$odds = \frac{p}{1-p} = \frac{p\mbox{(успеха)}}{p\mbox{(неудачи)}}, \mbox{область значений: }[0, +\infty]$$
Натуральный логарифм шансов:
$$\log(odds), \mbox{область значений: }[-\infty, +\infty]$$

Но, что нам говорит логарифм шансов? Как нам его интерпретировать?

```{r}
tibble(n = 10,
       success = 1:9,
       failure = n - success,
       prob.1 = success/(success+failure),
       odds = success/failure,
       log_odds = log(odds),
       prob.2 = exp(log_odds)/(1+exp(log_odds)))
```

Как связаны вероятность и логарифм шансов:
$$\log(odds) = \log\left(\frac{p}{1-p}\right)$$
$$p = \frac{\exp(\log(odds))}{1+\exp(\log(odds))}$$

Как связаны вероятность и логарифм шансов:

```{r}
data_frame(p = seq(0, 1, 0.001),
           log_odds = log(p/(1-p))) %>% 
  ggplot(aes(log_odds, p))+
  geom_line()+
  labs(x = latex2exp::TeX("$log\\left(\\frac{p}{1-p}\\right)$"))+
  theme_bw()
```

### 2.2 Почему не линейную регрессию?
```{r}
lm_0 <- lm(as.integer(have_ejectives)~1, data = ejectives)
lm_1 <- lm(as.integer(have_ejectives)~total, data = ejectives)
lm_0
lm_1
```
Первая модель:
$$ejectives = 1.316 \times consonants$$
Вторая модель:
$$ejectives = 0.4611 + 0.0353 \times consonants$$

```{r}
ejectives %>% 
  ggplot(aes(total, as.integer(have_ejectives)))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()+
  labs(y = "ejectives (yes = 2, no = 1)")
```

### 2.3 Логит: модель без предиктора
Будьте осторожны, `glm` не работает с тибблом.
```{r}
logit_0 <- glm(have_ejectives~1, family = "binomial", data = ejectives)
summary(logit_0)
logit_0$coefficients
table(ejectives$have_ejectives)
log(111/872) # β0
111/(111+872) # p
exp(log(111/872))/(1+exp(log(111/872))) # p
```

### 2.4 Логит: модель c одним числовым предиктором
```{r}
logit_1 <- glm(have_ejectives~total, family = "binomial", data = ejectives)
summary(logit_1)
logit_1$coefficients

ejectives %>% 
  mutate(have_ejectives = as.integer(have_ejectives)) %>% 
  ggplot(aes(total, have_ejectives)) +
  geom_point()+
  theme_bw()+
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"),
              se = FALSE)
```

Какова вероятность, что в языке с 29 согласными есть абруптивные?
```{r}
logit_1$coefficients
```

$$\log\left({\frac{p}{1-p}}\right)_i=\beta_0+\beta_1\times consinants_i + \epsilon_i$$
$$\log\left({\frac{p}{1-p}}\right)=-12.1123347 + 0.4576095 \times 29 = 1.158341$$
$$p = \frac{e^{1.158341}}{1+e^{1.158341}} = 0.7610311$$

```{r}
# log(odds)
predict(logit_1, newdata = data.frame(total = 29))
# p
predict(logit_1, newdata = data.frame(total = 29), type = "response")
```

### 2.5 Логит: модель c одним категориальным предиктором
```{r}
logit_2 <- glm(have_ejectives~area, family = "binomial", data = ejectives)
summary(logit_2)
logit_2$coefficients
table(ejectives$have_ejectives, ejectives$area)
log(1/6) # Eurasia
log(3/1) # North America
```

### 2.6 Логит: множественная регрессия
```{r}
logit_3 <- glm(have_ejectives~total+area, family = "binomial", data = ejectives)
summary(logit_3)
```

### 2.7 Логит: сравнение моделей
```{r}
AIC(logit_0)
AIC(logit_1)
AIC(logit_2)
AIC(logit_3)
```

Для того, чтобы интерпретировать коэффициенты нужно проделать трансформацю:
```{r}
(exp(logit_1$coefficients)-1)*100
```
Перед нами процентное изменние шансов при увеличении независимой переменной на 1.

Было предложено много аналогов R$^2$, например, McFadden's R squared:
```{r}
pscl::pR2(logit_1)
```

> 🤔 В датасете про гомосексуалов попробуйте предсказать ориентацию носителей  на основании переменных длительность s,  средняя частота основного тона и размер диапозона частоты основного тона. Приведите коэффициент при предикторе `average.f0.Hz` получившейся модели (с точностью до 3 знака после запятой).

```{r, include=FALSE}
homo_fit <- glm(orientation~s.duration.ms+average.f0.Hz+ f0.range.Hz, family = "binomial", data = homo)
summary(homo_fit)
# 0.025
```

<form name="FormSeven" onsubmit="return validateFormSeven()" method="post">
<input type="text" name="answerSeven">
<input type="submit" value="check">
</form><br>

> 🤔 Посчитайте вероятность быть гомосексуалом с условием `s.duration.ms = 62, average.f0.Hz =120, f0.range.Hz = 60. Ответ приведите с точностью до 3 знаков после запятой.

```{r, include=FALSE}
predict(homo_fit, newdata = data.frame(s.duration.ms = 62,average.f0.Hz =120, f0.range.Hz = 60), type = "response")
```


<form name="FormEight" onsubmit="return validateFormEight()" method="post">
<input type="text" name="answerEight">
<input type="submit" value="check">
</form><br>

## 3 Регрессия со случайными эффектами

In this dataset we have number of consonants and vowels in 402 languages collected from UPSID database (http://www.lapsyd.ddl.ish-lyon.cnrs.fr/lapsyd/). There is an variable of the area based on Glottolog (http://glottolog.org/). In this part we will try to make models that predict number of vowels by number of consonants.

```{r, warning= FALSE}
upsid <- read_csv("https://raw.githubusercontent.com/agricolamz/2019_data_analysis_for_linguists/master/data/upsid.csv")
upsid

lingtypology::map.feature(upsid$language, 
                          features = upsid$area,
                          label = upsid$language,
                          label.hide = TRUE)
```


```{r}
upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()
```

Обведем наблюдения по каждому спикеру:
```{r}
upsid %>% 
  ggplot(aes(consonants, vowels, color = area))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()+
  stat_ellipse()
```

Построим простую регрессию и добавим ее на график:
```{r}
fit1 <- lm(vowels~consonants, data = upsid)
summary(fit1)

upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()+
  geom_line(data = fortify(fit1), aes(x = consonants, y = .fitted), color = "blue")

fit2 <- lmer(vowels ~ consonants + (1|area), data = upsid)
summary(fit2)

upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw() +
  geom_line(data = fortify(fit2), aes(x = consonants, y = .fitted, color = area))
```

Если мы предполагаем скоррелированность свободных эффектов:
```{r}
fit3 <- lmer(vowels ~ consonants + (1+consonants|area), data = upsid)
summary(fit3)

upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()+
  geom_line(data = fortify(fit3), aes(x = consonants, y = .fitted, color = area))

fit4 <- lmer(vowels ~ consonants + (0+consonants|area), data = upsid)
summary(fit4)

upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()+
  geom_line(data = fortify(fit4), aes(x = consonants, y = .fitted, color = area))
```


Если мы не предполагаем скоррелированность свободных эффектов:
```{r}
fit5 <- lmer(vowels ~ consonants + (1|area) + (0+consonants|area), data = upsid)
summary(fit5)

upsid %>% 
  ggplot(aes(consonants, vowels))+
  geom_point()+
  labs(x = "number of consonants",
       y = "number of vowels",
       caption = "data from LAPSyD")+
  theme_bw()+
  geom_line(data = fortify(fit5), aes(x = consonants, y = .fitted, color = area))

fit0 <- lmer(vowels ~ 1 + (1|area) + (0+consonants|area), data = upsid)
summary(fit0)

anova(fit5, fit4, fit3, fit2, fit1, fit0)
```

### 3.1 Task
Pavel Duryagin ran an experiment on perception of vowel reduction in Russian language. The dataset shva includes the following variables:

* time1 - reaction time 1<\br>
* duration - duration of the vowel in the stimuly (in milliseconds, ms)<\br>
* time2 - reaction time 2<\br>
* f1, f2, f3 - the 1st, 2nd and 3rd formant of the vowel measured in Hz (for a short introduction into formants, see here)<\br>
* vowel - vowel classified according the 3-fold classification (A - a under stress, a - a/o as in the first syllable before the stressed one, y (stands for shva) - a/o as in the second etc. syllable before the stressed one or after the stressed syllable, cf. g[y]g[a]t[A]l[y] gogotala ‘guffawed’). In this part, we will ask you to analyse correlation between f1, f2, and duration.

```{r}
shva <- read_tsv("https://raw.githubusercontent.com/agricolamz/2019_PhonDan/master/data/duryagin_ReductionRussian.txt")
shva %>% 
  ggplot(aes(x=f2, y=f1, color=vowel)) +
  geom_point(show.legend = FALSE) +
  scale_x_reverse()+
  scale_y_reverse()+
  labs(x = "f2",
       y = "f1",
       title = "f2 and f1 of the reduced and stressed vowels",
       caption = "Data from Duryagin 2018")
```

Постройте регрессию со смешанными эффектами, предсказывая  f2 на основе f1 и гласный -- в случайных эффектах. В ответе приведите коэффициент при f1 (с точностью до 4 знаков после запятой:

```{r, include=FALSE}
fit2 <- lmer(f2~f1+(1|vowel), data = shva)
```

<form name="FormNine" onsubmit="return validateFormNine()" method="post">
<input type="text" name="answerNine">
<input type="submit" value="check">
</form><br>


<script>
function validateFormOne() {
    var x = document.forms["FormOne"]["answerOne"].value;
    if (x != "77.9269") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormTwo() {
    var x = document.forms["FormTwo"]["answerTwo"].value;
    if (x != "0.4938") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormThree() {
    var x = document.forms["FormThree"]["answerThree"].value;
    if (x != "500.1245") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormFour() {
    var x = document.forms["FormFour"]["answerFour"].value;
    if (x != "23.4771") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormFive() {
    var x = document.forms["FormFive"]["answerFive"].value;
    if (x != "10.0455") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormSix() {
    var x = document.forms["FormSix"]["answerSix"].value;
    if (x != "0.238") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormSeven() {
    var x = document.forms["FormSeven"]["answerSeven"].value;
    if (x != "0.025") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

function validateFormEight() {
    var x = document.forms["FormEight"]["answerEight"].value;
    if (x != "0.463") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}
function validateFormNine() {
    var x = document.forms["FormNine"]["answerNine"].value;
    if (x != "0.0689") {
        alert("У меня другой ответ...");
        return false;
    } else {
        alert("Да, все правильно");
        return false;
    }
}

</script>