Связан ли вес целого и вес части
Было исследовано 97 котов
У каждого индивида измеряли:
Пример взят из работы: R. A. Fisher (1947) The analysis of covariance method for the relation between a part and the whole, Biometrics 3, 65–68. Данные представлены в пакете boot
Посмотрим на датасет
library(readxl)
cat <- read_excel("data/catsM.xlsx")
head(cat)
## # A tibble: 6 x 3 ## Sex Bwt Hwt ## <chr> <dbl> <dbl> ## 1 M 2 6.5 ## 2 M 2 6.5 ## 3 M 2.1 10.1 ## 4 M 2.2 7.2 ## 5 M 2.2 7.6 ## 6 M 2.2 7.9
Есть ли пропущенные значения?
sum(!complete.cases(cat))
## [1] 0
Нет пропусков
##Каков объем выборки
nrow(cat)
## [1] 97
##Есть ли отскакивающие значения? {.smaller}
library(ggplot2) gg_dot <- ggplot(cat, aes(y = 1:nrow(cat))) + geom_point() gg_dot + aes(x = Bwt)
gg_dot + aes(x = Hwt)
Цель практически любого исследования - поиск взаимосвязи величин и создание базы для предсказания неизвестного на основе имеющихся данных
| Коэффициент | Функция | Особенности применения |
|---|---|---|
| Коэф. Пирсона | cor(x,y,method="pearson") |
Оценивает связь двух нормально распределенных величин. Выявляет только линейную составляющую взаимосвязи. |
| Ранговые коэффициенты (коэф. Спирмена, Кендалла) | cor(x,y,method="spirman")cor(x,y,method="kendall") |
Не зависят от формы распределения. Могут оценивать связь для любых монотонных зависимостей. |
##Оценка значимости корреляции Для оценки статистической значимости коэффициентов корреляции можно использовать функцию cor.test()
Hint Для построения точечной диаграммы вам понадобится geom_point()
cor.test(cat$Bwt, cat$Hwt)
## ## Pearson's product-moment correlation ## ## data: cat$Bwt and cat$Hwt ## t = 10, df = 100, p-value <2e-16 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.705 0.857 ## sample estimates: ## cor ## 0.793
pl_cat <- ggplot(cat,
aes(x = Bwt, y = Hwt)) +
geom_point() +
xlab("Вес кота") +
ylab("Вес сердца")
pl_cat
##Задание: проанализируйте данные и оцените корреляцию между признаками
Зависимость между диаметром клеток инфузорий и концентрацией клеток в культуре
данные представлены в пакете ISwR
infus <- read_excel("data/hellung.xls")
Опишите связь между концентрацией клеток и их диаметром
В некоторых методах статистики (многомерные методы) вместо корреляции применяется ковариация (согласованное отклонение):
\[ cov(X, Y) = \frac{1}{n - 1}\sum{(x_i - \bar{x})(y_i - \bar{y})} \] В отлиие от корреляции - это не стандартизованная величина, варьирующая в неограниченных пределах
Проведя корреляционный анализ, мы лишь ответили на вопрос “Существует ли статистически значимая связь между величинами?”
Сможем ли мы, используя это знание, предсказть значения одной величины, исходя из знаний другой?
Линейные модели \[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 + \epsilon_i\)
Появляется дополнительный член \(\epsilon_i\) Он вводит в модель влияние неучтенных моделью факторов. Обычно считают, что \(\epsilon \in N(0, \sigma^2)\)
Модель для данного примера имеет такой вид
\(response = 4.6 + 5.3I_{Level2} + 9.9 I_{Level3} + \epsilon_i\)
\(I_{i}\) - dummy variable
В стохастические модели выделяется две части:
Фиксированная часть: \(у_i = 2 + 5x_i\)
Случайная часть: \(\epsilon_i\)
Бывают модели, в которых случайная часть выглядит существенно сложнее (модели со смешаными эффектами). В таких моделях необходимо смоделировать еще и поведение случайной части.
Какая из линий “лучше” описывает облако точек?
“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)
Рост регрессирует (возвращается) к популяционной средней
Угловой коэффициент в зависимости роста потомков от роста родителей- коэффциент регресси
Симулированный пример с использованием geom_violin() 
(из кн. Quinn, Keough, 2002, стр. 85)
Остатки (Residuals):
\[e_i = y_i - \hat{y_i}\]
Линия регрессии (подобранная модель) - это та линия, у которой \(\sum{e_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 |
cat_model <- lm(Hwt ~ Bwt, data = cat) cat_model
## ## Call: ## lm(formula = Hwt ~ Bwt, data = cat) ## ## Coefficients: ## (Intercept) Bwt ## -1.18 4.31
cat_model?coefficients(cat_model) [1]
## (Intercept) ## -1.18
coefficients(cat_model) [2]
## Bwt ## 4.31
as.numeric(coefficients(cat_model) [1] + coefficients(cat_model) [2] * 2.5)
## [1] 9.6
cat$Hwt[10] - fitted(cat_model)[10]
## 10 ## 1.3
residuals(cat_model)[10]
## 10 ## 1.3
summary()summary(cat_model)
## ## Call: ## lm(formula = Hwt ~ Bwt, data = cat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.773 -1.048 -0.298 0.984 4.865 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.184 0.998 -1.19 0.24 ## Bwt 4.313 0.340 12.69 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.56 on 95 degrees of freedom ## Multiple R-squared: 0.629, Adjusted R-squared: 0.625 ## F-statistic: 161 on 1 and 95 DF, p-value: <2e-16
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_cat + geom_smooth(method="lm")
Доверительная зона регрессии с 95% вероятностью накроет истинную регрессионную прямую в генеральной совокупности при повторных выборках.
Возникает из-за неопределенности оценок коэффициентов модели, вследствие выборочного характера оценок.
Линии регрессии, полученные для 100 выборок (по 20 объектов в каждой), взятых из одной и той же генеральной совокупности 
coef(cat_model)
## (Intercept) Bwt ## -1.18 4.31
confint(cat_model)
## 2.5 % 97.5 % ## (Intercept) -3.17 0.798 ## Bwt 3.64 4.987
Если коэффициенты уравнения регрессии - лишь приблизительные оценки параметров, то предсказать значения зависимой переменной можно только с нeкоторой вероятностью.
newdata <- data.frame(Bwt = 2.5)
Predicted <- predict(cat_model, newdata, interval = "prediction", level = 0.95,
se = TRUE)$fit
Predicted
## fit lwr upr ## 1 9.6 6.48 12.7
Подготавливаем данные
cat_predicted <- predict(cat_model, interval="prediction") cat_predicted <- data.frame(cat, cat_predicted) head(cat_predicted)
## Sex Bwt Hwt fit lwr upr ## 1 M 2.0 6.5 7.44 4.28 10.6 ## 2 M 2.0 6.5 7.44 4.28 10.6 ## 3 M 2.1 10.1 7.87 4.72 11.0 ## 4 M 2.2 7.2 8.30 5.16 11.4 ## 5 M 2.2 7.6 8.30 5.16 11.4 ## 6 M 2.2 7.9 8.30 5.16 11.4
EstimateStd. Errort value, Pr(>|t|)Зависимость есть, если \(\beta_1 \ne 0\)
Нулевая гипотеза \(H_0: \beta = 0\)
Тестируем гипотезу
\[t=\frac{b_1-0}{SE_{b_1}}\]
Число степеней свободы: \(df=n-2\)
>- Итак,
>- t value - Значение t-критерия
>- Pr(>|t|) - Уровень значимости
\[Hwt = -1.18 + 4.31 Bwt\]
summary(cat_model)
## ## Call: ## lm(formula = Hwt ~ Bwt, data = cat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.773 -1.048 -0.298 0.984 4.865 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.184 0.998 -1.19 0.24 ## Bwt 4.313 0.340 12.69 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.56 on 95 degrees of freedom ## Multiple R-squared: 0.629, Adjusted R-squared: 0.625 ## F-statistic: 161 on 1 and 95 DF, p-value: <2e-16
| Объясненная дисперсия | Остаточная дисперсия | Полная дисперсия |
|---|---|---|
| \(SS_{Regression}=\sum{(\hat{y}-\bar{y})^2}\) | \(SS_{Residual}=\sum{(\hat{y}-y_i)^2}\) | \(SS_{Total}=\sum{(\bar{y}-y_i)^2}\) |
| \(df_{Regression} = 1\) | \(df_{Residual} = n-2\) | \(df_{Total} = n-1\) |
| \(MS_{Regression} =\frac{SS_{Regression}}{df}\) | \(MS_{Residual} =\frac{SS_{Residual}}{df_{Residual}}\) | \(MS_{Total} =\frac{SS_{Total}}{df_{Total}}\) |
Если зависимости нет, то \(MS _{Regression} = MS_{Residual}\)
\[ F= \frac{MS _{Regression}}{MS_{Residual}}\]
Логика та же, что и с t-критерием
Форма F-распределения зависит от двух параметров: \(df_{Regression} = 1\) и \(df_{Residual} = n-2\)
Коэффициент детерминации описывает какую долю дисперсии зависимой переменной объясняет модель
summary(cat_model)
## ## Call: ## lm(formula = Hwt ~ Bwt, data = cat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.773 -1.048 -0.298 0.984 4.865 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.184 0.998 -1.19 0.24 ## Bwt 4.313 0.340 12.69 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.56 on 95 degrees of freedom ## Multiple R-squared: 0.629, Adjusted R-squared: 0.625 ## F-statistic: 161 on 1 and 95 DF, p-value: <2e-16
Применяется если необходимо сравнить две модели с разным количеством параметров
\[ R^2_{adj} = 1- (1-R^2)\frac{n-1}{n-k}\]
\(k\) - количество параметров в модели
Вводится штраф за каждый новый параметр
Мы показали, что связь между весом сердца и весом тела котов описывается моделью вида
Hwt = -1.18 + 4.31 Bwt (\(F_{1,95}\) = 161, p < 0.001, \(R^2\) = 0.63)