aa989190f363e46d/pharm-simpson-para-pub.Rmd

## pharm-simpson-para-pub.Rmd
---
title: "Иллюстрация парадокса Симпсона"
author: "mz"
date: '`r Sys.time()`'
output:
  html_document: default
  pdf_document:
    keep_tex: yes
header-includes: \usepackage[T1,T2A]{fontenc}
---

```{r,echo=FALSE}
suppressMessages(library(ggplot2))
suppressMessages(library(magrittr))
suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(knitr))
library(parallel)
```

#От чего это

Весенний семестр, фарм. фак, Витебск, ОЭФ. Тема про основы экономического анализа деятельности аптечной организации, не вызывает никакого слюноотделения у преподавателя, а у студента тем более. В этом блюде не хватает специй, чего-нибудь жгучего, вроде смутного дискомфорта от того, что интуиция подвела тебя на простом примере. И, коль скоро на первом занятии в семестре речь идет о показателях относительных, а еще и о разнообразных группировках, то идея добавить в качестве перца парадокс Симпсона, приходит сама собой.

Чтобы создать тематический пример, нужно достаточно глубоко разобрать с сущности моделируемого явления, а разбираться лучше всего с определения. Парадоксом Симпсона принято называть такое явление, когда средние значения в группировках по, по меньшей мере, двум признакам для определенной категории одного признака меньше, вне зависимости от категории по второму признаку, но, в то же время, при группировке только по первому признаку средние обращаются, то есть категория наименьшим средним перестает быть таковой. Легче понять на примерах, а их исторически набралось уже довольно много, самые известные из них: дискриминация по половому признаку при процента поступивших в университет Беркли [[5], [6]], дискриминация по полу в отношении часовых ставок оплаты труда в различных регионах США [[4]], выживаемость пассажирова Титаника [[3]], соревнования улиток по восхождению на гору Фудзи [[2]], набор в хор [[1]]. Есть даже прямой пример из фармакологии про метаанализ исследований росиглитазона [[12]], но, это не про экономику.

[1]: http://shvarz.livejournal.com/433213.html
[2]: http://avva.livejournal.com/1919943.html
[3]: http://statistics.about.com/od/HelpandTutorials/a/What-Is-Simpsons-Paradox.htm
[4]: onlinelibrary.wiley.com/doi/10.1111/j.1468-0084.1985.mp47004005.x/full
[5]: http://vudlab.com/simpsons/
[6]: https://habrahabr.ru/post/279665/
[12]: http://bmcmedresmethodol.biomedcentral.com/articles/10.1186/1471-2288-8-34

Примеров не мало, и я бы мог прибегнуть к магии преобразования мяса в шоколад [[7], [8], [9]], но колебания в силе чувствую я, думая как стану заниматься таким, тем более, в статье «PRACTITIONERS CORNER: On Simpson's Paradox in Economic Statistics»[[4]], отлично описаны и проанализированы условия возникновения парадоксальной ситуации, и условий этих ровно два:

  1. Максимальное значение в группе с меньшими значениями ($L_{max}$) должно быть больше, чем минимальное значение в группе с бóльшими значениями ($B_{min}$)
  2. Разница между $L_{max}$ и $B_{min}$ ($\Delta$) должна быть больше, чем сумма произведений долей групп по первому признаку и соответствующих значений средних, более формально это удобно вразить как: $L_{max} - B_{min} > \sum_{i=1}^{n} ({B_i}-B_{min}) \times \omega_{B_i} + \sum_{i=1}^{n} (L_{max}-{L_i}) \times \omega_{L_i}$ , где $i$ — номер группы по второму признаку.


[7]: https://www.dissernet.org/publications/plagiat-geit.htm
[8]: https://ru.wikipedia.org/wiki/%D0%98%D0%B3%D0%BE%D1%88%D0%B8%D0%BD,_%D0%98%D0%B3%D0%BE%D1%80%D1%8C_%D0%9D%D0%B8%D0%BA%D0%BE%D0%BB%D0%B0%D0%B5%D0%B2%D0%B8%D1%87
[9]: https://republic.ru/fast/russia/kak-deputatu-napisat-dissteratsiyu-menyaem-shokolad-na-myaso-gotovo-917855.xhtml

#Симулякры и симуляции

Перед тем как сказку сделать былью, нужно её сочинить. Итак, преамбула пафоса мистического эпоса про анализ данных в аптеке будет следующая: начитавшись про естественный отбор, руководство небольшой аптечной сети «Фхтагн р'льех» приняло решение что держаться нету больше сил и пора сокращать ассортимент, и штаты. И если со штатми было просто: планировалось применить старейшую и самую надежную стратегию «е-отбор», то по второму пункту сложно определиться, какой именно ассортимент сокращать. Делить решено было по возрастным группам: олдыри и адалты (первый признак для группировки). Сказано — сделано! После тщательнейшего проектирования и сбора данных было установлено, что каждой отдельной аптеке (второй признак для группировки) одно посещение покупателем из одной из двух групп (средний чек) приносит в среднем рублей:

```{r make data,cache=FALSE, echo=FALSE,results='hide'}

# Calculate the number of cores
no_cores <- detectCores() - 1

#set.seed(15)
set.seed(199)
#set.seed(14)

base <- 30
count <- 5
#adult   <- c(32.9,36.2,38.5,31.1)
#eldery  <- c(30.8,35.1,37.7,30.2)
#adult <- eldery * 1.1
adult   <- base + sample.int(base,count)/5
eldery  <- adult - sample.int(base,count)/10

eldery.d <- abs(eldery - max(eldery))
adult.d  <- abs(adult  - min(adult))

aim.d <- max(eldery) - min(adult)

weighter.fabr <- function(sample.size = 30){

  #adult   <- c(32.9,36.2,38.5,31.1)
  #eldery  <- c(30.8,35.1,37.7,30.2)

  function(stub = NA){
    list(eldery.weights = as.numeric(c(1:length(eldery),
                                       sample(1:length(eldery),
                                              sample.size ,replace = T,
                                              pnorm(rnorm(length(eldery))))) %>% table()/(sample.size +length(eldery))),
          adult.weights = as.numeric(c(1:length(adult),
                                       sample.int(length(adult),
                                              sample.size ,replace = T,
                                              pnorm(rnorm(length(eldery))))) %>% table()/(sample.size +length(adult))))
  }
}

check.candidate <- function(candidates){
  all(sapply(candidates,
             function(candidate){aim.d < sum(candidate[["eldery.weights"]]*eldery.d + candidate[["adult.weights"]]*adult.d)}))
}

# Initiate cluster
cl <- makeCluster(no_cores)
clusterExport(cl, "eldery")
clusterExport(cl, "adult")

candidates <- list(list(eldery.weights = rep(1,length(eldery)),
                  adult.weights = rep(1,length(adult))))

i <- 0
while (check.candidate(candidates)){
  #candidate <- get.weights()
  candidates <- parLapply(cl,1:no_cores,weighter.fabr(3000))
  i = i+7
}

candidates[[which.max(sapply(candidates,
                             function(candidate){aim.d -
                                                  sum(candidate[["eldery.weights"]]*eldery.d +
                                                      candidate[["adult.weights"]]*adult.d)}))]] -> good.example

adult.weights <- good.example[["adult.weights"]]
eldery.weights <- good.example[["eldery.weights"]]

stopCluster(cl)
```

```{r insert table, echo=FALSE, results='asis', error=FALSE,warning=FALSE}

kable(data.frame(подразделение = paste0("Аптека №",1:count),
                  олдыри = eldery,
                  адалты = adult))

```

У самых эффективных среди менеджеров возникло однозначное решение, посколько только его можно принять глядя на данные. Сразу понятно что олдыри тут попали в группу $L$. Но смотреть и видеть не одно и то же, даже Джек «ретард» Салли знал это от синемордых. Некоторые из менеджеров постарше, больше, чем просто решать, любят складывать. Короче, решили они посчитать, сколько всей сети в среднем приносит один олдовый ходок за один заход, результат обескуражил:

```{r, echo=FALSE, results='asis'}
kable(data.frame(показатель = c('адалты','олдыри','разность'),
      величина = c(sum(adult*adult.weights),
      sum(eldery*eldery.weights),
      sum(eldery*eldery.weights)-sum(adult*adult.weights))))
```

[Закипел тимбилдинг,… полетел фалафель…](http://joyreactor.cc/post/2952256) Вес! Ответ на загадку. Готовь вазелин и заплатку.

#Что к чему

Реально, всему причиной вес, но не в том смысле, что олдыри жирнее. Просто то ли в аптеках, куда их ходит больше, сотрудники научились лучше их окучивать, то ли олдырям больше нравится ходить туда, где сотрудники лучше вынимают деньги, точно ответить нельзя по имеющимся данным, но факт в том, что нежиданный результат перегруппировки связан с тем, что среди олдовых в «дорогую» аптеку ходит большя доля олдовых, чем доля адалтов. У меня и график есть, чтоб было понятнее:

```{r graphing, echo=FALSE, error=FALSE,warning=FALSE}
categories <- factor(c("олдыри","адалты"))
data.frame(cat = rep(categories,c(15000,15000)),
           apt = factor(c(sample(1:count,15000,replace = T,prob=eldery.weights),
                          sample(1:count,15000,replace = T,prob=adult.weights)),
                        levels=1:count,labels=paste0('Аптека №',1:count),ordered=T)) -> df

df %>% ggplot(aes(cat,fill=apt)) +
  geom_histogram(stat="count",alpha=.55,position = 'fill') +
  xlab('Категория покупателя') + ylab('Вес в категории') +
  scale_fill_discrete('Аптека')

```

Или, если кому ближе дух Луки Пачоли, табличка:

```{r, echo=FALSE,results='asis'}

kable(data.frame("Аптека" = paste0('Аптека №',1:count),
                 "Олдыри" = eldery.weights*100,
                 "Адалты" = adult.weights*100),digits = 2)

```


И для закрепления инсайта, приведу диаграмму на которой видно сразу все. Бирюзовые элементы относятся к олдырям, коралловые — к адалтам, диаметры кругов — доля в групы (по первому признаку), посетившая соответствующую аптеку, линиям показаны уровни средних чеков без группировки по аптекам:

```{r, echo=FALSE}
apts <- factor(1:count,labels=paste0('Аптека №',1:count),ordered = T,levels = 1:count)
data.frame(apt = rep(apts,2),
           cat = rep(c('eldery','adult'),c(count,count)),
           weights = c(eldery.weights, adult.weights),
           price = c(eldery,adult)) -> sdf

sdf %>%
  ggplot(aes(apt,price,group=cat,color=cat)) +
  geom_point(aes(size=weights),alpha = .55) +
  geom_hline(yintercept = sum(adult*adult.weights),color='coral') +
  geom_hline(yintercept = sum(eldery*eldery.weights),color='cyan') +
  xlab('Аптека')+ylab('Средняя цена')+
  scale_color_hue(c=45, l=80,"Категоря посетителей", labels=c('юные','пожилые')) +
  scale_size_area('Доля, (%)',breaks = c(.025,.050,.1,.25,.50),labels=c('2.5','5','10','25','50'),max_size = 15)
```

Или то же, но в профиль [[10]]. Коралловый и бирюзовый использованы аналогично предыдущей диаграмме.

[10]: http://www.tandfonline.com/doi/abs/10.1080/00031305.1985.10479387

```{r, echo=FALSE}

sdf %>% mutate(wp=weights*price) %>% group_by(cat) %>% summarise(wp=sum(wp),weights=sum(weights)) -> sdfg
sdf %>%
  ggplot(aes(cat,price,group=apt,color=apt)) +
  geom_line() +
  geom_point(aes(size=weights),alpha = .55) +
  xlab('Категория посетителей')+ylab('Средняя цена')+
  scale_color_hue(c=45, l=80,"Аптека", labels=paste0('Аптека №',1:count)) +
  scale_size_area('Доля, (%)',breaks = c(.025,.050,.1,.25,.50),labels=c('2.5','5','10','25','50'),max_size = 15) +
  scale_x_discrete(labels=c('юные','пожилые')) +
  geom_line(data = sdfg,
            aes(cat,wp,group='1'),
            alpha = rep(.75),
            colour='gray',
            inherit.aes = F) +
  geom_point(data = sdfg,
             aes(cat,wp),
             size=15,
             colour=c('coral','cyan'),
             alpha = .55,
             inherit.aes = F)

```

К сожалению я не смог освоить хард-кор из «Visualizing Statistical Mix Effects and Simpson's Paradox»[[11]].

[11]: https://www.ncbi.nlm.nih.gov/pubmed/26356927

Тем, кто захочет воcпроизвести мои построения поможет исходный [Rmd-файл](https://gist.github.com/aa989190f363e46d/b299f5fc3a80a4d67aebca3f499752a4) для запуска в R через knit (Rstudio делает это особенно удобно).

Осталось решить только еще один вопрос — что придумать про планирование рецептуры, на следующее занятие.
	---
	title: "Иллюстрация парадокса Симпсона"
	author: "mz"
	date: '`r Sys.time()`'
	output:
	html_document: default
	pdf_document:
	keep_tex: yes
	header-includes: \usepackage[T1,T2A]{fontenc}
	---

	```{r,echo=FALSE}
	suppressMessages(library(ggplot2))
	suppressMessages(library(magrittr))
	suppressMessages(library(dplyr))
	suppressMessages(library(tidyr))
	suppressMessages(library(knitr))
	library(parallel)
	```

	#От чего это

	Весенний семестр, фарм. фак, Витебск, ОЭФ. Тема про основы экономического анализа деятельности аптечной организации, не вызывает никакого слюноотделения у преподавателя, а у студента тем более. В этом блюде не хватает специй, чего-нибудь жгучего, вроде смутного дискомфорта от того, что интуиция подвела тебя на простом примере. И, коль скоро на первом занятии в семестре речь идет о показателях относительных, а еще и о разнообразных группировках, то идея добавить в качестве перца парадокс Симпсона, приходит сама собой.

	Чтобы создать тематический пример, нужно достаточно глубоко разобрать с сущности моделируемого явления, а разбираться лучше всего с определения. Парадоксом Симпсона принято называть такое явление, когда средние значения в группировках по, по меньшей мере, двум признакам для определенной категории одного признака меньше, вне зависимости от категории по второму признаку, но, в то же время, при группировке только по первому признаку средние обращаются, то есть категория наименьшим средним перестает быть таковой. Легче понять на примерах, а их исторически набралось уже довольно много, самые известные из них: дискриминация по половому признаку при процента поступивших в университет Беркли [[5], [6]], дискриминация по полу в отношении часовых ставок оплаты труда в различных регионах США [[4]], выживаемость пассажирова Титаника [[3]], соревнования улиток по восхождению на гору Фудзи [[2]], набор в хор [[1]]. Есть даже прямой пример из фармакологии про метаанализ исследований росиглитазона [[12]], но, это не про экономику.

	[1]: http://shvarz.livejournal.com/433213.html
	[2]: http://avva.livejournal.com/1919943.html
	[3]: http://statistics.about.com/od/HelpandTutorials/a/What-Is-Simpsons-Paradox.htm
	[4]: onlinelibrary.wiley.com/doi/10.1111/j.1468-0084.1985.mp47004005.x/full
	[5]: http://vudlab.com/simpsons/
	[6]: https://habrahabr.ru/post/279665/
	[12]: http://bmcmedresmethodol.biomedcentral.com/articles/10.1186/1471-2288-8-34

	Примеров не мало, и я бы мог прибегнуть к магии преобразования мяса в шоколад [[7], [8], [9]], но колебания в силе чувствую я, думая как стану заниматься таким, тем более, в статье «PRACTITIONERS CORNER: On Simpson's Paradox in Economic Statistics»[[4]], отлично описаны и проанализированы условия возникновения парадоксальной ситуации, и условий этих ровно два:

	1. Максимальное значение в группе с меньшими значениями ($L_{max}$) должно быть больше, чем минимальное значение в группе с бóльшими значениями ($B_{min}$)
	2. Разница между $L_{max}$ и $B_{min}$ ($\Delta$) должна быть больше, чем сумма произведений долей групп по первому признаку и соответствующих значений средних, более формально это удобно вразить как: $L_{max} - B_{min} > \sum_{i=1}^{n} ({B_i}-B_{min}) \times \omega_{B_i} + \sum_{i=1}^{n} (L_{max}-{L_i}) \times \omega_{L_i}$ , где $i$ — номер группы по второму признаку.


	[7]: https://www.dissernet.org/publications/plagiat-geit.htm
	[8]: https://ru.wikipedia.org/wiki/%D0%98%D0%B3%D0%BE%D1%88%D0%B8%D0%BD,_%D0%98%D0%B3%D0%BE%D1%80%D1%8C_%D0%9D%D0%B8%D0%BA%D0%BE%D0%BB%D0%B0%D0%B5%D0%B2%D0%B8%D1%87
	[9]: https://republic.ru/fast/russia/kak-deputatu-napisat-dissteratsiyu-menyaem-shokolad-na-myaso-gotovo-917855.xhtml

	#Симулякры и симуляции

	Перед тем как сказку сделать былью, нужно её сочинить. Итак, преамбула пафоса мистического эпоса про анализ данных в аптеке будет следующая: начитавшись про естественный отбор, руководство небольшой аптечной сети «Фхтагн р'льех» приняло решение что держаться нету больше сил и пора сокращать ассортимент, и штаты. И если со штатми было просто: планировалось применить старейшую и самую надежную стратегию «е-отбор», то по второму пункту сложно определиться, какой именно ассортимент сокращать. Делить решено было по возрастным группам: олдыри и адалты (первый признак для группировки). Сказано — сделано! После тщательнейшего проектирования и сбора данных было установлено, что каждой отдельной аптеке (второй признак для группировки) одно посещение покупателем из одной из двух групп (средний чек) приносит в среднем рублей:

	```{r make data,cache=FALSE, echo=FALSE,results='hide'}

	# Calculate the number of cores
	no_cores <- detectCores() - 1

	#set.seed(15)
	set.seed(199)
	#set.seed(14)

	base <- 30
	count <- 5
	#adult <- c(32.9,36.2,38.5,31.1)
	#eldery <- c(30.8,35.1,37.7,30.2)
	#adult <- eldery * 1.1
	adult <- base + sample.int(base,count)/5
	eldery <- adult - sample.int(base,count)/10

	eldery.d <- abs(eldery - max(eldery))
	adult.d <- abs(adult - min(adult))

	aim.d <- max(eldery) - min(adult)

	weighter.fabr <- function(sample.size = 30){

	#adult <- c(32.9,36.2,38.5,31.1)
	#eldery <- c(30.8,35.1,37.7,30.2)

	function(stub = NA){
	list(eldery.weights = as.numeric(c(1:length(eldery),
	sample(1:length(eldery),
	sample.size ,replace = T,
	pnorm(rnorm(length(eldery))))) %>% table()/(sample.size +length(eldery))),
	adult.weights = as.numeric(c(1:length(adult),
	sample.int(length(adult),
	sample.size ,replace = T,
	pnorm(rnorm(length(eldery))))) %>% table()/(sample.size +length(adult))))
	}
	}

	check.candidate <- function(candidates){
	all(sapply(candidates,
	function(candidate){aim.d < sum(candidate[["eldery.weights"]]eldery.d + candidate[["adult.weights"]]adult.d)}))
	}

	# Initiate cluster
	cl <- makeCluster(no_cores)
	clusterExport(cl, "eldery")
	clusterExport(cl, "adult")

	candidates <- list(list(eldery.weights = rep(1,length(eldery)),
	adult.weights = rep(1,length(adult))))

	i <- 0
	while (check.candidate(candidates)){
	#candidate <- get.weights()
	candidates <- parLapply(cl,1:no_cores,weighter.fabr(3000))
	i = i+7
	}

	candidates[[which.max(sapply(candidates,
	function(candidate){aim.d -
	sum(candidate[["eldery.weights"]]*eldery.d +
	candidate[["adult.weights"]]*adult.d)}))]] -> good.example

	adult.weights <- good.example[["adult.weights"]]
	eldery.weights <- good.example[["eldery.weights"]]

	stopCluster(cl)
	```

	```{r insert table, echo=FALSE, results='asis', error=FALSE,warning=FALSE}

	kable(data.frame(подразделение = paste0("Аптека №",1:count),
	олдыри = eldery,
	адалты = adult))

	```

	У самых эффективных среди менеджеров возникло однозначное решение, посколько только его можно принять глядя на данные. Сразу понятно что олдыри тут попали в группу $L$. Но смотреть и видеть не одно и то же, даже Джек «ретард» Салли знал это от синемордых. Некоторые из менеджеров постарше, больше, чем просто решать, любят складывать. Короче, решили они посчитать, сколько всей сети в среднем приносит один олдовый ходок за один заход, результат обескуражил:

	```{r, echo=FALSE, results='asis'}
	kable(data.frame(показатель = c('адалты','олдыри','разность'),
	величина = c(sum(adult*adult.weights),
	sum(eldery*eldery.weights),
	sum(elderyeldery.weights)-sum(adultadult.weights))))
	```

	[Закипел тимбилдинг,… полетел фалафель…](http://joyreactor.cc/post/2952256) Вес! Ответ на загадку. Готовь вазелин и заплатку.

	#Что к чему

	Реально, всему причиной вес, но не в том смысле, что олдыри жирнее. Просто то ли в аптеках, куда их ходит больше, сотрудники научились лучше их окучивать, то ли олдырям больше нравится ходить туда, где сотрудники лучше вынимают деньги, точно ответить нельзя по имеющимся данным, но факт в том, что нежиданный результат перегруппировки связан с тем, что среди олдовых в «дорогую» аптеку ходит большя доля олдовых, чем доля адалтов. У меня и график есть, чтоб было понятнее:

	```{r graphing, echo=FALSE, error=FALSE,warning=FALSE}
	categories <- factor(c("олдыри","адалты"))
	data.frame(cat = rep(categories,c(15000,15000)),
	apt = factor(c(sample(1:count,15000,replace = T,prob=eldery.weights),
	sample(1:count,15000,replace = T,prob=adult.weights)),
	levels=1:count,labels=paste0('Аптека №',1:count),ordered=T)) -> df

	df %>% ggplot(aes(cat,fill=apt)) +
	geom_histogram(stat="count",alpha=.55,position = 'fill') +
	xlab('Категория покупателя') + ylab('Вес в категории') +
	scale_fill_discrete('Аптека')

	```

	Или, если кому ближе дух Луки Пачоли, табличка:

	```{r, echo=FALSE,results='asis'}

	kable(data.frame("Аптека" = paste0('Аптека №',1:count),
	"Олдыри" = eldery.weights*100,
	"Адалты" = adult.weights*100),digits = 2)

	```


	И для закрепления инсайта, приведу диаграмму на которой видно сразу все. Бирюзовые элементы относятся к олдырям, коралловые — к адалтам, диаметры кругов — доля в групы (по первому признаку), посетившая соответствующую аптеку, линиям показаны уровни средних чеков без группировки по аптекам:

	```{r, echo=FALSE}
	apts <- factor(1:count,labels=paste0('Аптека №',1:count),ordered = T,levels = 1:count)
	data.frame(apt = rep(apts,2),
	cat = rep(c('eldery','adult'),c(count,count)),
	weights = c(eldery.weights, adult.weights),
	price = c(eldery,adult)) -> sdf

	sdf %>%
	ggplot(aes(apt,price,group=cat,color=cat)) +
	geom_point(aes(size=weights),alpha = .55) +
	geom_hline(yintercept = sum(adult*adult.weights),color='coral') +
	geom_hline(yintercept = sum(eldery*eldery.weights),color='cyan') +
	xlab('Аптека')+ylab('Средняя цена')+
	scale_color_hue(c=45, l=80,"Категоря посетителей", labels=c('юные','пожилые')) +
	scale_size_area('Доля, (%)',breaks = c(.025,.050,.1,.25,.50),labels=c('2.5','5','10','25','50'),max_size = 15)
	```

	Или то же, но в профиль [[10]]. Коралловый и бирюзовый использованы аналогично предыдущей диаграмме.

	[10]: http://www.tandfonline.com/doi/abs/10.1080/00031305.1985.10479387

	```{r, echo=FALSE}

	sdf %>% mutate(wp=weights*price) %>% group_by(cat) %>% summarise(wp=sum(wp),weights=sum(weights)) -> sdfg
	sdf %>%
	ggplot(aes(cat,price,group=apt,color=apt)) +
	geom_line() +
	geom_point(aes(size=weights),alpha = .55) +
	xlab('Категория посетителей')+ylab('Средняя цена')+
	scale_color_hue(c=45, l=80,"Аптека", labels=paste0('Аптека №',1:count)) +
	scale_size_area('Доля, (%)',breaks = c(.025,.050,.1,.25,.50),labels=c('2.5','5','10','25','50'),max_size = 15) +
	scale_x_discrete(labels=c('юные','пожилые')) +
	geom_line(data = sdfg,
	aes(cat,wp,group='1'),
	alpha = rep(.75),
	colour='gray',
	inherit.aes = F) +
	geom_point(data = sdfg,
	aes(cat,wp),
	size=15,
	colour=c('coral','cyan'),
	alpha = .55,
	inherit.aes = F)

	```

	К сожалению я не смог освоить хард-кор из «Visualizing Statistical Mix Effects and Simpson's Paradox»[[11]].

	[11]: https://www.ncbi.nlm.nih.gov/pubmed/26356927

	Тем, кто захочет воcпроизвести мои построения поможет исходный [Rmd-файл](https://gist.github.com/aa989190f363e46d/b299f5fc3a80a4d67aebca3f499752a4) для запуска в R через knit (Rstudio делает это особенно удобно).

	Осталось решить только еще один вопрос — что придумать про планирование рецептуры, на следующее занятие.