Простые статистические модели

результат = (модель) + ошибка

тестовая статистика = дисперсия, объясняемая моделью / дисперсия, не объясняемая моделью

Примеры тестовых статистик: F-тест, t-статистика и другие.

Простая линейная регрессия

Общий вид уравнения, который мы оцениваем \[Y_i = (b_0 + b_1* X_i) + \epsilon \]

\(b_0\) - intercept (свободный член)
\(b_1\) - slope (угол наклона кривой)

Наша цель выбрать такие параметры модели, который минимизируеют разницу между фактическими значениями объясняемой переменной и значениями, предсказанными моделью. В частности, в МНК минимизируется сумма квадратов остатков (могут быть и другие варианты)

\[\sum_{i=1}^{n} (Y_i- \hat{Y_i})^2 \]

Чтобы правильно оценивать значения получаемых коэффициентов, необходимо гарантировать выполнение допущений МНК.

\[ дисперсия = \sum (наблюдения - модель)^2 \]

\[ R^2 = \frac{SS_M}{SS_T} \]

\[ F =\frac{MS_M}{MS_R} \]

\[ t = \frac{b_{наблюд}}{SE_b} \]

Допущения МНК

Допущения - это важно. Если допущения метода не выполняются – необходимо менять метод. Но допущения чаще всего не выполняются на реальных данных. Это ответственность исследователя - определить “критичность” допущений.

Допущения МНК:

  1. Типы переменных на X и Y - количественные(X,Y), факторные (X)
  2. Не-нулевая дисперсия
  3. Нет мультиколлинеарности
  4. Предикторы не коррелированы с факторами, не включенными в модель
  5. Гомокседексичность (стабильность дисперсии ошибок)
  6. Отсутствие автокорреляции ошибок (независимость случаев)
  7. Нормально распределенные ошибки
  8. Независимость Y
  9. Линейная зависимость между предикторами и Y

Проверка нормальности распределения

С начала необходимо просто поосмотреть на данные - к примеру, гистограммы всех переменных

options(digits=2)
dlf <-  read.delim('datasets/Album Sales 1.dat', header = TRUE)

library(ggplot2)
hist.adverts <- ggplot(dlf, aes(adverts))+
   geom_histogram(aes (y = ..density..))+theme_bw()+
   labs(x = 'Затраты на рекламу', y = 'плотность')
hist.adverts
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

hist.adverts + stat_function(fun = dnorm, args = list(mean = mean(dlf$adverts, na.rm = TRUE),
                           sd = sd(dlf$adverts, na.rm = TRUE)), color  = 'red')
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

hist.sales <- ggplot(dlf, aes(sales))+
   geom_histogram(aes (y = ..density..))+theme_bw()+
   labs(x = 'Продажи', y = 'плотность')

hist.sales + stat_function(fun = dnorm, args = list(mean = mean(dlf$sales, na.rm = TRUE),
                           sd = sd(dlf$sales, na.rm = TRUE)), color  = 'red')
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

Посмотреть на данные - Q-Q график (квантиль-квантиль).
Квантиль - доля случаев, меньше определенного значения. Если требования “нормальности выполняется” - данные должны лежать на диагональной линиии

qqplot.day1 <- qplot(sample = dlf$sales, stat = 'qq')
qqplot.day1+theme_bw()

Посмотреть на стастические характеристики данных

Тесты на нормальность - Шапиро-Уилк (Shapiro-Weil) Тест просто проверяет разницу между фактическим и теоретическим нормальным распределением. Если разница значима - тест на нормальность “не проходит”.

Пример расчета регрессии

Пример - влияние рекламы (в $) на продажи музукальных альбомов (в штуках)

album1 <- read.delim('datasets/Album Sales 1.dat', header = TRUE)
head(album1)
##   adverts sales
## 1      10   330
## 2     986   120
## 3    1446   360
## 4    1188   270
## 5     575   220
## 6     569   170
album.sales <- lm(album1$sales ~ album1$adverts)
summary(album.sales)
## 
## Call:
## lm(formula = album1$sales ~ album1$adverts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -152.95  -43.80   -0.39   37.04  211.87 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.34e+02   7.54e+00   17.80   <2e-16 ***
## album1$adverts 9.61e-02   9.63e-03    9.98   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66 on 198 degrees of freedom
## Multiple R-squared:  0.335,  Adjusted R-squared:  0.331 
## F-statistic: 99.6 on 1 and 198 DF,  p-value: <2e-16

результаты модели графически

plot(album1$adverts, album1$sales)
lines(album1$adverts, predict(album.sales), col = 'red')

Пример расчет множественной регрессии

Пример - зависимость уровня убийств по штатам (на примере США) в зависимости от социально-экономических и прочих показателей (численность населения, грамотность, уровень доходов и количество холодных дней).

states<-data.frame(state.x77)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)
##            Murder Population Illiteracy Income Frost
## Murder       1.00       0.34       0.70  -0.23 -0.54
## Population   0.34       1.00       0.11   0.21 -0.33
## Illiteracy   0.70       0.11       1.00  -0.44 -0.67
## Income      -0.23       0.21      -0.44   1.00  0.23
## Frost       -0.54      -0.33      -0.67   0.23  1.00
library(car)
scatterplotMatrix(states,spread = FALSE, lty.smooth = 2, main ='Соотношение переменных')

Рассчитаем модель

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
summary(fit)
## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
##     data = states)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.796 -1.649 -0.081  1.482  7.621 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.23e+00   3.87e+00    0.32    0.751    
## Population  2.24e-04   9.05e-05    2.47    0.017 *  
## Illiteracy  4.14e+00   8.74e-01    4.74  2.2e-05 ***
## Income      6.44e-05   6.84e-04    0.09    0.925    
## Frost       5.81e-04   1.01e-02    0.06    0.954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.5 on 45 degrees of freedom
## Multiple R-squared:  0.567,  Adjusted R-squared:  0.528 
## F-statistic: 14.7 on 4 and 45 DF,  p-value: 9.13e-08
confint(fit)
##                2.5 %  97.5 %
## (Intercept) -6.6e+00 9.02132
## Population   4.1e-05 0.00041
## Illiteracy   2.4e+00 5.90387
## Income      -1.3e-03 0.00144
## Frost       -2.0e-02 0.02083

Значимы только уровень грамотности и численность населения. \(R^2\) довольно высокий (0,57)

Диагностика модели

Нормальность - ошибки имеют нормальное распределение

qqPlot(fit, labels = row.names(states), simulate = TRUE, main = 'График Q-Q')

Результат - похоже, есть выбросы

Независимость ошибок - тест Дарбина-Уотсона

durbinWatsonTest((fit))
##  lag Autocorrelation D-W Statistic p-value
##    1            -0.2           2.3    0.28
##  Alternative hypothesis: rho != 0

Результат - нет автокорреляции

Линейность -

crPlots((fit))

Результат - есть нелинейные связи

Гомокседаксичность

ncvTest((fit))
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.7    Df = 1     p = 0.19
spreadLevelPlot(fit)

## 
## Suggested power transformation:  1.2

Тест показывает на незначимость нулевой гипотезы - остатки постоянны

Глобольная оценка допущений

library(gvlma)
summary(gvlma(fit))
## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
##     data = states)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.796 -1.649 -0.081  1.482  7.621 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.23e+00   3.87e+00    0.32    0.751    
## Population  2.24e-04   9.05e-05    2.47    0.017 *  
## Illiteracy  4.14e+00   8.74e-01    4.74  2.2e-05 ***
## Income      6.44e-05   6.84e-04    0.09    0.925    
## Frost       5.81e-04   1.01e-02    0.06    0.954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.5 on 45 degrees of freedom
## Multiple R-squared:  0.567,  Adjusted R-squared:  0.528 
## F-statistic: 14.7 on 4 and 45 DF,  p-value: 9.13e-08
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = fit) 
## 
##                    Value p-value                Decision
## Global Stat        2.773   0.597 Assumptions acceptable.
## Skewness           1.537   0.215 Assumptions acceptable.
## Kurtosis           0.638   0.425 Assumptions acceptable.
## Link Function      0.115   0.734 Assumptions acceptable.
## Heteroscedasticity 0.482   0.487 Assumptions acceptable.

Стандартные проблемы регрессиий

sqrt(vif(fit))>2
## Population Illiteracy     Income      Frost 
##      FALSE      FALSE      FALSE      FALSE

Что делать?

Выбор лучшей модели

ANOVA - сравнеие вложенных (nested) моделей.

Информационные критерии (AIC, BIC и другие) - меньшие значения параметров означают лучшую модель. Сами абсолютные значения параметров не имеют смысла. Сравнивать по критерями можно не вложенные модели.

Выбор параметров - stepwise regression (но лучше этого не делать)

Cross-validation - разбиение на данных на две части (расчетная и тестовая часть).

Cтатистическая мощность

Величина эффекта (effect size) связана с тремя другими статистическими характеристиками: 1. величина выборки (sample size) 2. уровень вероятности, при котором мы принимаем, что эффект значим (0.05) 3. способность статистического теста обнаружить эффект такого размера (мощность теста)

Если известно любые две из трех характеристик - можно вычислить третью.

Статистическая незначимость - не означает, что эффекта нет! Возможно, просто недостаточно данных.

Мнимая регрессия (spurious regression)

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