Приложение. Примеры заданий

Парная регрессия. Метод наименьших квадратов (Курс “Прикладное ПО”)

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

Регрессия – метод, позволяющий обнаружить, как одна или несколько независимых переменных (факторов) влияют на одну зависимую. Мы можем оценить, насколько увеличивается либо уменьшается зависимая переменная при росте независимой/независимых.

Вам даны две переменные: количество друзей ВКонтакте (1, 1, 2, 3, 5, 8, 13, 21, 34, 55) и количество друзей оффлайн (1, 3, 5, 7, 11, 13, 17, 19, 23, 27)

1. Создайте две переменных (два вектора) с любыми удобными вам именами. Они должны быть интервальными. Объедините их в таблицу данных (класс data.frame)

2. Постройте диаграмму рассеивания (scatterplot), где количество друзей ВКонтакте – зависимая переменная (\(y\)) с помощью пакета ggplot2.

3. Добавьте линию тренда (используя geom_smooth(method="lm")). Итоговая картинка должна напоминать эту:

4. Задача: найти уравнение полученной прямой по стандартному образцу \(y = ax+b\)

Чтобы узнать уравнение прямой, создадим регрессионную модель через функцию lm (linear model). Она выдаст нам необходимые коэффициенты:

ols <- lm(y ~ x, data=data1)

Здесь y – зависимая переменная; x – независимая; data1 – название вашей таблицы данных.

Получим описание нашей регрессионной модели, вычисленной методом наименьших квадратов (OLS – Ordinary Least Squares):

summary(ols)  
# обращаем внимание ТОЛЬКО на табличку Сoefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -8.9113     4.5029  -1.979 0.083178 .  
x          1.8422     0.2981   6.180 0.000265 ***

Чтобы составить уравнение прямой, нам необходимо узнать два параметра: точку пересечения (Intercept, в стандартном уравнении прямой – параметр \(b\)) и коэффициент наклона (Slope, \(a\)). Оба показателя находятся в столбце Estimate (крайний левый). \(b\) на строчке (Intercept), коэффициент \(a\) – на строке с наименованием переменной x.

Запишите уравнение прямой!

5. Теперь необходимо вычислить размер ошибки – то, на сколько полученная нами прямая отклоняется от наблюдаемых, т.е. реальных значений наших данных. Таким образом мы поймем, в какой мере регрессионной прямой удается описывать тренды, заключенные в наших данных и насколько можно верить ее предсказаниям. Для этого высчитайте \(E\) (error) для каждого случая: \(E(x)=(y_\text{наблюдаемое} - y_\text{ожидаемое})^2\).

Ожидаемое (предсказываемое) значение можно посчитать простой подстановкой значения x в уравнение прямой.

Вычислите сумму квадратов ошибок (сумму E). На ваш взгляд, ориентируясь на график и данные, какая функция бы наиболее подошла для описания взаимного распределения этих двух переменных?

Задание “Эффекты взаимодействия в регрессионных моделях” (Курс анализа данных)

I. Проинтерпретируйте эффекты взаимодействия

a. \(\text{SelfConcept} = 1.843 + 0.145\cdot\text{MATH} + 0.477 \cdot \text{Best_School} - 0.051 \cdot\text{MATH} \cdot\text{Best_School}\)

Здесь:

  • Self-Concept – уверенность учащегося в своих силах (от «не уверен» до «очень уверен», интервальная)
  • MATH – оценка по математике (от 1 до 10)
  • Best_School – ученик учится в одной из лучших школ города (10% лучших школ)

b. \(\text{Abortion} = 5.502 + 1.420 \cdot \text{Ex_comm} – 0.613 \cdot \text{Chatt} + 0.463 \cdot \text{Ex_comm} \cdot \text{Chatt}\)

Здесь:

  • Abortion – отношение индивида к абортам (от «очень отрицательно» до «положительно»)
  • Ex_comm – входила ли страна раньше в состав СССР (0 – нет, 1 – да)
  • Chatt – насколько часто индивид ходит в церковь (от «очень редко» до «очень часто»)

c. \(\text{Popularity} = 0.232 + 0.027 \cdot\text{BOY} + 0.057 \cdot \text{GPA} – 0.019 \cdot\text{BOY} \cdot \text{GPA}\)

Здесь:

  • Popularity – популярность учащегося по количеству друзей (от 0 до 10)
  • Boy – пол (0 – ж, 1 – м)
  • GPA – средний балл (от 1 до 5, интервальная)

II. Постройте эффект взаимодействия

Используйте набор данных swiss {datasets}.

  1. Сформулируйте гипотезу исследования
  2. Постройте множественную регрессию, объясняющую рождаемость (Fertility). Опишите полученную модель (значимость и направление эффектов)
  3. Добавьте в свою модель эффект взаимодействия (в данном случае между двумя интервальными переменными)
  4. Опишите полученный эффект: направление, значимость
  5. Опишите полученную модель, в том числе эффект взаимодействия. Внимание! Если у вас не получилось значимого эффекта, то все равно опишите содержательно тот эффект, который вы искали, и затем сообщите о его незначимости.
  6. (*) Попробуйте придумать, как можно визуализировать эффект взаимодействия для двух интервальных переменных

Семинар Random Permutations, Means and Hypothesis Testing (Курсы “Программирование для социологов”, “ТВ и МС”)

library(ggplot2)

Acknowledgements first:

  • Example and a part of the code from Ahmed Moustafa. Please, play with the simulation there after you’re done with the task.
  • Data (via the source above) from Śaunak Sen

Data

“”“Systolic blood pressure was measured in progeny from a backcross between two mouse strains. 50 (randomly chosen) mice were genotyped at the D4Mit214 marker. We want to detect association between the D4Mit214 marker genotype and blood pressure. The values show the systolic blood pressure (in mm of Hg) by the marker genotype, BA (heterozygous) or BB (homozygous).”“”

# Heterozygous (BA)
a = c(86, 88, 89, 89, 92, 93, 94, 94, 94, 95, 95, 96, 96, 97, 97, 98, 98, 99, 99, 101, 106, 107, 110, 113, 116, 118)
   
# Homozygous (BB)
b = c(89, 90, 92, 93, 93, 96, 99, 99, 99, 102, 103, 104, 105, 106, 106, 107, 108, 108, 110, 110, 112, 114, 116, 116)

df <- data.frame(
  bpressure=c(a, b),
  type=as.factor(c(rep("BA", length(a)), rep("BB", length(b))))
)

Let’s draw a boxplot to see if there is a difference between two groups.

ggplot(df) + geom_boxplot(aes(type, bpressure, fill=type))

Observed difference

summary(a)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   86.00   94.00   96.50   98.46  100.50  118.00
summary(b)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   89.00   98.25  104.50  103.20  108.50  116.00
diff.observed = mean(b) - mean(a)

Observed difference between two groups: 4.7467949

So, do mice with two genotypes really have different blood pressure (in population), or observed difference (between samples) is due to chance?

Hypothesis testing

How do we test hypothesis? Reminder

Throughout the testing procedure, the null hypothesis is considered to be true.

  • (Ei incumbit probatio qui dicit, non qui negat (the burden of proof is on he who declares, not on he who denies)) (Wikipedia)
  • “This is often expressed in the phrase innocent until proven guilty, coined by the English lawyer Sir William Garrow. Garrow insisted that accusers be robustly tested in court. An objective observer in the position of the juror must reasonably conclude that the defendant almost certainly committed the crime.” (Wikipedia)
  1. Formulate your statements in statistical terms (i.e. \(H_0\) and \(H_A\))
  • For our example: \(H_0: \mu_{BA} = \mu_{BB}\); \(H_A: \mu_{BA} \neq \mu_{BB}\)
  1. Choose a significance level (\(\alpha\))
  • \(\alpha\) is the probability to falsely reject a true \(H_0\), consider its costs
  1. Determine the appropriate test statistic and calculate it. Easy, right?
  2. Check if the test gives us enough evidence to reject the null hypothesis (for chosen \(\alpha\))
  3. Formulate your resulting statement in “normal” language.

Our usual way

t.test(a, b)
## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -2.027, df = 47.926, p-value = 0.04824
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -9.45540466 -0.03818508
## sample estimates:
## mean of x mean of y 
##  98.46154 103.20833

… and we look smart, but have no freaking idea how it works, and what the hell is p-value.

  • What is t-statistics?
  • How it works and why we use it?
  • What assumptions we need to check?
  • When it doesn’t work and how will I know?

Random Permutation Test

Let’s try and use less math we don’t understand and more computing. We will use only sums and means and fractions, I promise :)

Observed difference between two groups: 4.7467949

Throughout the testing procedure, the null hypothesis is considered to be true.

Well, it is easy. If Null hypothesis is true, \(a\) and \(b\) are from the same population, so we can combine our samples.

combined <- c(a, b)

Now we need to check how this difference will behave if \(H_0\) is true in a long run.

So we will make random permutations (with respect to \(a\) and \(b\) lengths), each time count the difference and write it down.

number_of_permutations <- 5000

count_diff <- function(number_of_permutations = 1000) {
  diff.random <- rep(NA_real_, number_of_permutations)
  for (i in 1 : number_of_permutations) {

    # Sample from the combined dataset without replacement
    shuffled = sample(combined, length(combined))
    
    a.random = shuffled[1 : length(a)]
    b.random = shuffled[(length(a) + 1) : length(combined)]
  
    # Null (permuated) difference
    diff.random[i] = mean(b.random) - mean(a.random)
  }
  diff.random 
}

Sampling distribution of the difference between means (provided \(H_0\) is true)

Our observed difference

Critical 5% region

q <- quantile(diff.random, c(0.025, 0.975))

P-values and stuff

P-value is the probability of getting the result you got, or an even more extreme result, if the null hypothesis is true.

Here counting p-value approx. is extremely simple: probability is just a fraction. For each random permutation we check if resulting difference in means is equal or larger than our observed difference and see what share of such cases do we have in the total number of permutations.

pvalue <- function(diff.random, diff.observed, number_of_permutations){
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    (pvalue) 
}

pvalue(diff.random, diff.observed, number_of_permutations)
## [1] 0.0502

Reading (Watching|Playing) List

Tasks:

  • How sample works? What is runif and rnorm? How do they work? Draw them for different \(n\) (i.e. 1e1, 1e2, 1e3, 1e5) and explain.
  • Why making many permutations is neccessary? Test our function for different \(n\), draw histogram for diff.random, count p-values and make conclusions.
  • () Central Limit Theorem in Statistics, actually, proves the result about asymptotic normality of the sampling distribution of means which we observed via simulation. What about sampling distributuion of medians*? Can we use traditional methods here? What about simulations? How our approach will behave in this case for different \(n\)?

Сквозное задание “Организация процесса анализа данных” (Курс “Программирование для социологов”)

Это задание предназначено для систематизации навыков организации процесса анализа данных с использованием системы Data Workflow Drake.

Литература к заданию: главы 4 и 6 из книги ”Data Science at the Command Line”, интернет-источники.

В этом задании вам необходимо организовать процесс анализа и создания отчёта по данным игровых чатов. Исходные данные извлекаются из БД в соответствии с распределением по вариантам. Результат созданного workflow – набор отчётов в формате HTML, сгенерированных из написанных вами шаблнов RMarkdown, и содержащий:

  1. Общий отчёт(‘report.html‘), включающий, как минимум, таблицу с общим количеством сообщений для каждого пользователя и график, показывающий помесячное количество сообщений для каждого из top-5 пользователей в вашем варианте. Приветствуется дополнительная информация, позволяющая составить общее представление об интенсивности коммуникации пользователей. Необходимо помнить, что отчёт будет создаваться динамически, поэтому привязка к конкретным данным недопустима.

  2. Индивидуальный отчёт (‘USERID.html‘) для каждого пользователя, содержащий, как минимум, облако слов для этого пользователя, первые и последние (хронологически) пять сообщений. Не забывайте про стоп-слова. Приветствуется дополнительная информация.

Исходные данные расположены в СУБД и достаточно велики, чтобы вызвать проблемы при непосредственной обработке в R (или сделать её очень медленной). Поэтому необходимо создание workflow, объединяющего разные технологии на разных этапах обработки, в частности, вам нужно:

  • Создать запрос к БД, возвращающий данные для вашего варианта. Для его генерации с учётом конкретных идентификаторов вашего варианта подготовлен скрипт на Python, который нужно добавить в Drakefile с необходимой модификацией.
  • Выполнить его с использованием консольного клиента БД в скрипте оболочки
  • Составить цепочку команд, извлекающую Top-20 самых частых слов в сообщениях вашего варианта(игнорируя все остальные поля, извлеченные из БД), записать Top-20 для исходного и лемматизированного с помощью ‘mystem‘ текста в файл ‘top20.txt‘ (не перезапишите файл при добавлении лемматизированных результатов!)
  • Составить шаблоны RMarkdown для обоих отчётов, с использованием дополнительных источников разобраться с добавлением в ‘Drakefile‘ функции формирования отчётов по шаблону.

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

Контрольные вопросы:

  • Каковы ограничения отчётов, формируемых автоматически?
  • Как автоматизированная часть workflow встраивается в цепочку работы аналитика данных?