Зависит ли уровень интеллекта от размера головного мозга? (Willerman et al. 1991)
Scan_03_11 by bucaorg(Paul_Burnett) on Flickr
Было исследовано 20 девушек и 20 молодых людей
У каждого индивида измеряли:
Пример взят из работы: Willerman, L., Schultz, R., Rutledge, J. N., and Bigler, E. (1991), “In Vivo Brain Size and Intelligence,” Intelligence, 15, 223-228.
Данные представлены в библиотеке “The Data and Story Library” http://lib.stat.cmu.edu/DASL/
Посмотрим на датасет
brain <- read.csv("data/IQ_brain.csv", header = TRUE)
head(brain)
## Gender FSIQ VIQ PIQ Weight Height MRINACount ## 1 Female 133 132 124 118 64.5 816932 ## 2 Male 140 150 124 NA 72.5 1001121 ## 3 Male 139 123 150 143 73.3 1038437 ## 4 Male 133 129 128 172 68.8 965353 ## 5 Female 137 132 134 147 65.0 951545 ## 6 Female 99 90 110 146 69.0 928799
Есть ли пропущенные значения?
sum(!complete.cases(brain))
## [1] 2
Где именно?
sapply(brain, function(x) sum(is.na(x)))
## Gender FSIQ VIQ PIQ Weight Height ## 0 0 0 0 2 1 ## MRINACount ## 0
Что это за случаи?
brain[!complete.cases(brain), ]
## Gender FSIQ VIQ PIQ Weight Height MRINACount ## 2 Male 140 150 124 NA 72.5 1001121 ## 21 Male 83 83 86 NA NA 892420
Каков объем выборки
nrow(brain) ## Это без учета пропущенных значений
## [1] 40
Цель практически любого исследования - поиск взаимосвязи величин и создание базы для предсказания неизвестного на основе имеющихся данных
Наличие связи между явлениями не означает, что между ними существует причинно-следственная связь.
Сила связи между явлениями может быть количественно измерена.
Оценивает двух величин от своих средних значений \[cov_{x,y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{n - 1}\]
Коэффициент ковариации варьирует в интервале \(-\infty < cov_{x,y} < +\infty\)
Это стандартизованное значение ковариации
\[r_{x,y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\sum(x_i-\bar{x})^2}\sqrt{\sum(y_i-\bar{y})^2}} = \frac{cov_{x,y}} {\sigma_x \sigma_y}\]
Коэффициент корреляции варьирует в интервале \(-1 \le r_{x,y} \le +1\)
| Коэффициент | Функция | Особенности применения |
|---|---|---|
| Коэф. Пирсона | cor(x,y,method="pearson") |
Оценивает связь двух нормально распределенных величин. Выявляет только линейную составляющую взаимосвязи. |
| Ранговые коэффициенты (коэф. Спирмена, Кендалла) | cor(x,y,method="spearman")cor(x,y,method="kendall") |
Не зависят от формы распределения. Могут оценивать связь для любых монотонных зависимостей. |
Hint 1: Обратите внимание на то, что в датафрейме есть пропущенные значения. Изучите, как работают с NA функции, вычисляющие коэффициенты корреляции.
pl_brain <- ggplot(brain,
aes(x = MRINACount, y = PIQ)) +
geom_point() +
xlab("Brain size") +
ylab("IQ test")
pl_brain
cor.test(brain$PIQ, brain$MRINACount, method = "pearson", alternative = "two.sided")
## ## Pearson's product-moment correlation ## ## data: brain$PIQ and brain$MRINACount ## t = 3, df = 38, p-value = 0.01 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.0856 0.6232 ## sample estimates: ## cor ## 0.387
cor(brain[,2:6], use = "pairwise.complete.obs")
## FSIQ VIQ PIQ Weight Height ## FSIQ 1.0000 0.9466 0.93413 -0.05148 -0.0860 ## VIQ 0.9466 1.0000 0.77814 -0.07609 -0.0711 ## PIQ 0.9341 0.7781 1.00000 0.00251 -0.0767 ## Weight -0.0515 -0.0761 0.00251 1.00000 0.6996 ## Height -0.0860 -0.0711 -0.07672 0.69961 1.0000
Проведя корреляционный анализ, мы лишь ответили на вопрос “Существует ли статистически значимая связь между величинами?”
Сможем ли мы, используя это знание, предсказть значения одной величины, исходя из знаний другой?
Линейные модели \[y = b_0 + b_1x\]
\[y = b_0 + b_1x_1 + b_2x_2\] Нелинейные модели \[y = b_0 + b_1^x\]
\[y = b_0^{b_1x_1+b_2x_2}\]
Простая модель \[y = b_0 + b_1x\]
Множественная модель \[y = b_0 + b_1x_1 + b_2x_2 + b_3x_3 + ... + b_nx_n\]
Модель: \(у_i = 2 + 5x_i\)
Два параметра: угловой коэффициент (slope) \(b_1=5\); свободный член (intercept) \(b_0=2\)
Чему равен \(y\) при \(x=10\)?
Модель: \(у_i = 2 + 5x_i + \varepsilon_i\)
Появляется дополнительный член \(\varepsilon_i\) Он вводит в модель влияние неучтенных моделью факторов. Обычно считают, что \(\epsilon \in N(0, \sigma^2)\)
В стохастические модели выделяется две части:
Фиксированная часть: \(у_i = 2 + 5x_i\)
Случайная часть: \(\varepsilon_i\)
Бывают модели, в которых случайная часть выглядит существенно сложнее (модели со смешанными эффектами). В таких моделях необходимо смоделировать еще и поведение случайной части.
Модель для данного примера имеет такой вид
\(response = 4.6 + 5.3I_{Level2} + 9.9 I_{Level3}\)
\(I_{i}\) - dummy variable
Какая из линий “лучше” описывает облако точек?
## `geom_smooth()` using formula 'y ~ x'
“Essentially, all models are wrong,
but some are useful”
(Georg E. P. Box)
Френсис Галтон (Francis Galton)
“the Stature of the adult offspring … [is] … more mediocre than the stature of their Parents” (цит. по Legendre & Legendre, 1998)
Рост регрессирует (возвращается) к популяционной средней
Угловой коэффициент в зависимости роста потомков от роста родителей- коэффциент регресси
Симулированный пример, в котором
\[ y_i = 10x_i + \varepsilon_i \\ \varepsilon \in N(0, 10) \]
Точки (данные) - это выборки из нескольких “локальных” генеральных совокупностей с нормальным распределением зависимой переменной (каждая совокупность соответствует одному значению предиктора).
Линия регрессии проходит через средние значения нормальных распределений.
Параметры прямой подбираются так, чтобы вероятность существования данных при таких параметрах была бы максимальной.
Аналогичная картинка, но с использованием geom_violin()
## `geom_smooth()` using formula 'y ~ x'
(из кн. Quinn, Keough, 2002, стр. 85)
Остатки (Residuals):
\[\varepsilon_i = y_i - \hat{y_i}\]
Линия регрессии (подобранная модель) - это та линия, у которой \(\sum{\varepsilon_i}^2\) минимальна.
lm()fit <- lm(formula, data)
Модель записывается в виде формулы
| Модель | Формула |
|---|---|
| Простая линейная регрессия \(\hat{y_i}=b_0 + b_1x_i\) |
Y ~ X Y ~ 1 + X Y ~ X + 1 |
| Простая линейная регрессия (без \(b_0\), “no intercept”) \(\hat{y_i}=b_1x_i\) |
Y ~ -1 + X Y ~ X - 1 |
| Уменьшенная простая линейная регрессия \(\hat{y_i}=b_0\) |
Y ~ 1 Y ~ 1 - X |
| Множественная линейная регрессия \(\hat{y_i}=b_0 + b_1x_i +b_2x_2\) |
Y ~ X1 + X2 |
lm()fit <- lm(formula, data)
Элементы формул для записи множественных моделей
| Элемент формулы | Значение |
|---|---|
: |
Взаимодействие предикторов Y ~ X1 + X2 + X1:X2 |
* |
Обозначает полную схему взаимодействий Y ~ X1 * X2 * X3 аналогично Y ~ X1 + X2 + X3+ X1:X2 + X1:X3 + X2:X3 + X1:X2:X3 |
. |
Y ~ . В правой части формулы записываются все переменные из датафрейма, кроме Y |
brain_model <- lm(PIQ ~ MRINACount, data = brain) brain_model
## ## Call: ## lm(formula = PIQ ~ MRINACount, data = brain) ## ## Coefficients: ## (Intercept) MRINACount ## 1.74376 0.00012
brain_model?coefficients(brain_model) [1]
## (Intercept) ## 1.74
coefficients(brain_model) [2]
## MRINACount ## 0.00012
as.numeric(coefficients(brain_model) [1] + coefficients(brain_model) [2] * 900000)
## [1] 110
brain$PIQ[10] - fitted(brain_model)[10]
## 10 ## 30.4
residuals(brain_model)[10]
## 10 ## 30.4
summary()summary(brain_model)
## ## Call: ## lm(formula = PIQ ~ MRINACount, data = brain) ## ## Residuals: ## Min 1Q Median 3Q Max ## -39.6 -17.9 -1.6 17.0 42.3 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.7437570 42.3923825 0.04 0.967 ## MRINACount 0.0001203 0.0000465 2.59 0.014 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 21 on 38 degrees of freedom ## Multiple R-squared: 0.15, Adjusted R-squared: 0.127 ## F-statistic: 6.69 on 1 and 38 DF, p-value: 0.0137
EstimateStd. Errort valuePr(>|t|)
| Параметр | Оценка | Стандартная ошибка |
|---|---|---|
| \(\beta_1\) Slope |
\(b _1 = \frac {\sum _{i=1}^{n} {[(x _i - \bar {x})(y _i - \bar {y})]}}{\sum _{i=1}^{n} {(x _i - \bar x)^2}}\) или проще \(b_0 = r\frac{sd_y}{sd_x}\) |
\(SE _{b _1} = \sqrt{\frac{MS _e}{\sum _{i=1}^{n} {(x _i - \bar {x})^2}}}\) |
| \(\beta_0\) Intercept |
\(b_0 = \bar y - b_1 \bar{x}\) | \(SE _{b _0} = \sqrt{MS _e [\frac{1}{n} + \frac{\bar x}{\sum _{i=1}^{n} {(x _i - \bar x)^2}}]}\) |
| \(\epsilon _i\) | \(e_i = y_i - \hat {y_i}\) | \(\approx \sqrt{MS_e}\) |
pl_brain + geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'
>Доверительная зона регрессии. В ней с 95% вероятностью лежит регрессионная прямая, описывающая связь в генеральной совокупности.
Возникает из-за неопределенности оценок коэффициентов модели, вследствие выборочного характера оценок.
Симулированный пример: Генеральная совокупность, в которой связь между Y и X, описывается следующей зависимостью \[ y_i = 10 + 10x_i + \varepsilon_i \\ \varepsilon \in N(0, 20) \]
pop_x <- rnorm(1000, 10, 3)
pop_y <- 10 + 10*pop_x + rnorm(1000, 0, 20)
population <- data.frame(x=pop_x, y=pop_y)
ggplot(population, aes(x=x, y=y)) +
geom_point(alpha=0.3, color="red") +
geom_abline(aes(intercept=10, slope=10),
color="blue", size=2)
Линии регрессии, полученные для 100 выборок (по 20 объектов в каждой), взятых из одной и той же генеральной совокупности
coef(brain_model)
## (Intercept) MRINACount ## 1.74376 0.00012
confint(brain_model)
## 2.5 % 97.5 % ## (Intercept) -84.0751348 87.562649 ## MRINACount 0.0000261 0.000214
## `geom_smooth()` using formula 'y ~ x' ## `geom_smooth()` using formula 'y ~ x' ## `geom_smooth()` using formula 'y ~ x'
Если коэффициенты уравнения регрессии - лишь приблизительные оценки параметров, то предсказать значения зависимой переменной можно только с нeкоторой вероятностью.
newdata <- data.frame(MRINACount = 900000)
predict(brain_model, newdata,
interval = "prediction",
level = 0.95, se = TRUE)$fit
## fit lwr upr ## 1 110 66.9 153
Подготавливаем данные
brain_predicted <- predict(brain_model, interval="prediction") brain_predicted <- data.frame(brain, brain_predicted) head(brain_predicted)
## Gender FSIQ VIQ PIQ Weight Height MRINACount fit lwr upr ## 1 Female 133 132 124 118 64.5 816932 100 56.1 144 ## 2 Male 140 150 124 NA 72.5 1001121 122 78.2 166 ## 3 Male 139 123 150 143 73.3 1038437 127 81.9 171 ## 4 Male 133 129 128 172 68.8 965353 118 74.5 161 ## 5 Female 137 132 134 147 65.0 951545 116 73.0 159 ## 6 Female 99 90 110 146 69.0 928799 113 70.4 157
## `geom_smooth()` using formula 'y ~ x'
pl_brain +
# 1) Линия регрессии и ее дов. интервал
# Если мы указываем fill внутри aes() и задаем фиксированное значение -
# появится соотв. легенда с названием.
# alpha - задает прозрачность
geom_smooth(method = "lm", aes(fill = "Conf.interval"), alpha = 0.4) +
# 2) Интервал предсказаний создаем при помощи геома ribbon ("лента")
# Данные берем из другого датафрейма - из brain_predicted
# ymin и ymax - эстетики геома ribbon, которые задают нижний и верхний
# край ленты в точках с заданным x (x = MRINACount было задано в ggplot()
# при создании pl_brain, поэтому сейчас его указывать не обязательно)
#
geom_ribbon(data = brain_predicted, aes(ymin = lwr, ymax = upr, fill = "Conf. area for prediction"), alpha = 0.2) +
# 3) Вручную настраиваем цвета заливки при помощи шкалы fill_manual.
# Ее аргумент name - название соотв. легенды, values - вектор цветов
scale_fill_manual(name = "Intervals", values = c("green", "gray")) +
# 4) Название графика
ggtitle("Confidence interval \n and confidence area for prediction")
Модель “работает” только в том диапазоне значений независимой переменной (\(x\)), для которой она построена (интерполяция). Экстраполяцию надо применять с большой осторожностью.
## `geom_smooth()` using formula 'y ~ x'
EstimateStd. Errort value, Pr(>|t|)