Цели, задачи и процесс анализа данных

За последние несколько десятилетий арсенал прикладной статистики значительно пополнился за счет развития новых методов, краткий список которых включает в себя: линейную регрессию, обобщённые линейные и смешанные модели (generalized linear models), обобщённые аддитивные модели (generalized additive models), модели со смешанными эффектами (mixed-effects models), регрессионные и классификационные деревья (regression and classification trees), нейронные сети (neural networks), анализ выживаемости (survival analysis), многомерный анализ (multivariate analysis) с такими его методами, как: анализ главных компонент (principal component analysis), канонический анализ соответствий (canonical correspondence analysis), различные временные ряды (time series) и т.д. Часть из них появилась достаточно давно, но появление быстрых современных компьютеров и свободного программного обеспечения, такого как GNU R, сделало все эти требующие вычислительных ресурсов методы доступными практически для каждого исследователя. Однако такая доступность ещё больше обострила хорошо известную проблему всех статистических методов, которую на английском языке часто описывают как “rubbish in, rubbish out”“, т.е.”мусор на входе - мусор на выходе“.

Речь здесь идет о следующем: чудес не бывает, и если мы не будем уделять должного внимания тому, как тот или иной метод работает и какие требования предъявляет к анализируемым данным, то получаемые с его помощью результаты нельзя будет воспринимать всерьез. Поэтому каждый раз исследователю следует начинать свою работу с тщательного ознакомления со свойствами полученных данных и проверки необходимых условий применимости соответствующих статистических методов. Этот начальный этап анализа называют разведочным (Exploratory Data Analysis).

Выявление значений-промахов

При обработке реальных экономических, финансовых, технических и иных данных, полученных в процессе наблюдений, измерений, расчётов, иногда один или несколько результатов наблюдений резко выделяются, т.е. далеко отстоят от основной массы данных. В рассматриваемых случаях возникает естественная мысль о том, что подобные наблюдения не относятся к изучаемой совокупности, поскольку содержат грубую погрешность, а получены они в результате ошибки, промаха. В справочнике по метрологии об этом явлении говорится так: “Грубая погрешность (промах) - это случайная погрешность результата отдельного наблюдения, входящего в ряд измерений, которая для данных условий резко отличается от остальных результатов этого ряда. Они, как правило, возникают из-за ошибок или неправильных действий оператора (его психофизиологического состояния, неверного отсчета, ошибок в записях или вычислениях, неправильного включения приборов или сбоев в их работе и др.). Возможной причиной возникновения промахов также могут быть кратковременные резкие изменения условий проведения измерений. Если промахи обнаруживаются в процессе измерений, то результаты, их содержащие, отбрасывают. Однако чаще всего промахи выявляют только при окончательной обработке результатов измерений с помощью специальных критериев” [1, C. 88].

Чувствительность разных статистических методов к наличию аномальных наблюдений в данных неодинакова. Так, при использовании обобщённой линейной модели для анализа зависимой переменной, распределённой по закону Пуассона, наличие промахов может вызвать избыточную дисперсию, что сделает модель неприменимой [2]. В то же время при использовании непараметрического многомерного шкалирования, основанного на индексе Жаккара, все исходные данные переводятся в номинальную шкалу с двумя значениями (1/0), и наличие промахов никак не сказывается на результат анализа [3]. Исследователь должен чётко понимать эти различия между разными методами и при необходимости выполнять проверку на наличие аномальных наблюдений в данных.

Обычно для выявления промахов используют диаграммы размахов [4], на которых отображаются медианы и разбросы значений анализируемых совокупностей. Диаграммы размахов, или “ящики с усами”" (box-whisker plots), получили свое название за характерный вид: точку или линию, соответствующую медиане или средней арифметической, окружает прямоугольник (“ящик”“), длина которого соответствует одному из показателей разброса или точности оценки генерального параметра. Дополнительно от этого прямоугольника отходят”усы“”, также соответствующие по длине одному из показателей разброса или точности. Графики этого типа очень популярны, поскольку позволяют дать очень полную статистическую характеристику анализируемой совокупности. Кроме того, диаграммы размаха можно использовать для визуальной экспресс-оценки разницы между двумя и более группами (например, между датами отбора проб, экспериментальными группами, участками пространства, и т.п.).Обычно медиана в этом графике представлена в виде горизонтальной линии. Эта линия окаймляется прямоугольником, длина которого соответствует интерквартильному размаху (ИКР). От прямоугольника отходят линии, ограничивающие интервал длиной 1,5ИКР. Наблюдения, находящиеся за пределами интервала 1,5ИКР, рассматриваются как потенциальные промахи. Однако всегда следует внимательно относиться к такого рода нестандартным наблюдениям. поскольку поскольку эти “промахи”" вполне могут оказаться в порядке вещей для исследуемой совокупности. Так в нормальном распределении шансы случайной величины превысить два стандартных отклонения составляют 1 к 21, а шансы выйти за три стандартных отклонения составляют 1 к 369.

В R для построения диаграмм размаха используется функция boxplot.

Подключим стандартные данные, имеющиеся в R:

data("cars")
head(cars)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

В данном случае мы имеем дело с данными зависимости тормозного пути автомобиля от скорости. Подробно сведени можно получить с помощью команды ?cars

cars {datasets} R Documentation
Speed and Stopping Distances of Cars

Description

The data give the speed of cars and the distances taken to stop. Note that the data were recorded in the 1920s.

Usage

cars
Format

A data frame with 50 observations on 2 variables.

[,1]     speed   numeric     Speed (mph)
[,2]     dist    numeric     Stopping distance (ft)
Source

Ezekiel, M. (1930) Methods of Correlation Analysis. Wiley.

References

McNeil, D. R. (1977) Interactive Data Analysis. Wiley.

Построим диаграммы для обоих колонок значений

boxplot(cars$speed)

boxplot(cars$dist)

Обратите внимание, что на второй диаграмме явно виден потенциальный выброс.

Другим очень полезным, но, к сожалению, недостаточно используемым графическим средством выявления промахов является точечная диаграмма Кливленда. На таком графике по оси ординат откладывают порядковые номера отдельных наблюдений, а по оси абсцисс - значения этих наблюдений. Наблюдения, “значительно”" выделяющиеся из основного облака точек, потенциально могут быть промахами. Для построения точечных диаграмм Кливленда в R служит функция dotchart.

dotchart(cars$dist)
points(cars[49,2],49,pch=19,col="red")

Обратите внимание на точку в правом верхнем углу, специально покрашенную второй командой в красный цвет. Возможно, что это выброс.

Другим не менее полезным инструментом визуализации яляется диаграмма типа «скрипка» (violin plot), совмещенная с классической диаграммой размахов (box-and-whisker plot). Эта комбинированная графическая форма представления дает больше информации о характере распределения, чем диаграмма размахов или гистограмма, т.к. помимо данных о медиане и квартилях, отражает еще и показатели ядерной плотности распределения [5]. К сожалению, выбросы она не показывает.

library(vioplot)
## Loading required package: sm
## Package 'sm', version 2.2-5.4: type help(sm) for summary information
vioplot(cars$dist)

Поэтому можно использовать другой график, beanplot из пакета beanplot который следует предварительно установить из меню RStudio Tools.

library(beanplot)
beanplot(cars$dist,cars$speed)

Подробнее об этом графике смотри https://cran.r-project.org/web/packages/beanplot/vignettes/beanplot.pdf

До этого момента мы называли “промахом”" или “выбросом” наблюдение, которое существенно отличается от большинства других наблюдений в исследуемой совокупности. Однако более строгий подход к определению промахов состоит в оценке того, какое влияние эти необычные наблюдения оказывают на результаты анализа. При этом следует делать различие между необычными наблюдениями для зависимых и независимых переменных (предикторов). Решение отбросить какие-либо данные в конечном счете всегда субъективно. Наиболее строгий подход к проблеме анализа подобных аномальных данных - повторить измерения много раз. Если подобная аномалия повторится снова, вероятно, можно будет выяснить её причину. Если же аномальное показание не повторится, то даже оставляя его, оно мало скажется на дальнейших вычислениях из-за большого объёма измерений. Однако при небольшом объёме опытных данных и отсутствии явных причин появления грубых ошибок целесообразно использовать статистические методы выявления промахов [6, c.245]. Рассмотрим несколько статистических тестов на наличие аномальных наблюдений, имеющихся в R.

Q тест Диксона

(Dixon’s Q test) применим для малых выборок (меньше 30 наблюдений). Он применим только для данных, распределённых нормально. Если тест выявляет выбросы, то после их устранения, повторно этот тест использовать нельзя [7]. Для его проведения используют функцию dixon.test из пакета outliers [8]. Нулевая гипотеза для этого и последующего тестов на наличие аномальных значений заключается в том, что наименьшее или наибольшее отклонение от среднего в выборке не является промахом.

Поскольку у нас больше 30 значений, возьмём для примера другую случайную выборку и проведем этот тест.

library(outliers)
sample<-c(28.28,29.10,30.95,31.49,32.49,32.79,32.79,31.99,31.99)
dixon.test(sample)
## 
##  Dixon test for outliers
## 
## data:  sample
## Q = 0.18182, p-value = 0.966
## alternative hypothesis: lowest value 28.28 is an outlier

Как мы видим, наименьшее значение 28.28 не является выбросом (p-value > 0.05).

Тест хи-квадрат для промахов

Другой тест, тест хи-квадрат для промахов, требует знаний о дисперсии генеральной совокупности. В противном случае его использовать нельзя. Тест применим для тестирования остатков моделей, входит в пакет outliers и вызывается командой chisq.out.test.

chisq.out.test(cars$dist)
## 
##  chi-squared test for outlier
## 
## data:  cars$dist
## X-squared = 8.933, p-value = 0.002801
## alternative hypothesis: highest value 120 is an outlier

Тест Граббса

Для проверки наличия промахов в нормально распределённых данных также используется тест Граббса (Grubbs test). Тест также входит в пакет outliers и вызывается командой grubbs.test. Вызов команды имеет вид grubbs.test( sample, type = 10, opposite = FALSE, two.sided = TRUE). Параметр sample должен содержать исследуемые данные, параметр type принимает одно из трех значений: 10 - для одного промаха, 11 - для двух промахов с противоположных концов распределения, 20 - для двух промахов с одной стороны распределения.

grubbs.test(cars$dist)
## 
##  Grubbs test for one outlier
## 
## data:  cars$dist
## G = 2.98880, U = 0.81397, p-value = 0.04413
## alternative hypothesis: highest value 120 is an outlier

Если удаление необычных значений предиктора по тем или иным причинам не представляется возможным, помочь может определенное преобразование этого предиктора (например, логарифмирование). С необычными значениями зависимой переменной бороться сложнее, особенно при построении регрессионных моделей. Преобразование путём, например, логарифмирования, может помочь, но поскольку зависимая переменная представляет особый интерес при построении регрессионных моделей, лучше попытаться подобрать метод анализа, который основан на распределении вероятностей, допускающем больший разброс значений для больших средних величин (например, гамма-распределение для непрерывных переменных или распределение Пуассона для дискретных количественных переменных). Такой подход позволит работать с исходными значениями зависимой переменной.

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

Экспресс-обзор регрессионный анализ

Поскольку windows XP продолжло радовать мы перешли к регрессии, немного забежав вперед.

Выделим следующие этапы построения эконометрической модели для срезов данных:

  1. Спецификация модели:
  1. Оценивание структурных параметров модели.
  2. Верификация эконометрической модели
  1. Формулирование выводов и их интерпретация.

В R линейная регрессионная модель классическим методом наименьшиъ квадратов строится ледующим образом:

lm1<-lm(dist~speed,data=cars)
summary(lm1)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Интерпретация результатов, полученных выше с помощью команды summary, осуществляется по следующему алгоритму. Рассмотрим полученные результаты по порядку сверху вниз, за исключением наиболее важной статистики, F-статистики, которая находится в самом низу.

F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

F-статистика вычисляется по формуле \(F=\frac{\frac{Q_R}{k-1}}{\frac{Q_e}{n-k}}\), далее результат сравнивается с квантилью распределения Фишера с данными (1 и 48) степенями свободы и вычисляется вероятность того, что F превышает соответсвующий квантиль распределения Фишера. При необходимость, тот же результат можно получить используя команду pf c параметром lower.tail=F. Команда pf возвращает значение функции распределения Фишера, т.е. вероятность того, что значения функции не превосходят заданную величину. Для решения противоположной задачи используется параметр lower.tail=T.

pf(89.57,1,48,lower.tail=F)
## [1] 1.489077e-12

Полученное p-значение p-value показывает вероятность того, что модель является незначимой. Иными словами, если оно получилось меньше чем 0,05, то это означает, что наиболее вероятен факт, что модель значима. Если же, напротив, она превышает 0,05, то модель, скорее всего, назначаемая. В приведённом выше примере p-value: < 1.49e-12, откуда следует, что модель, скорее всего, является значимой. Следует отметить, что большинство людей в первую очередь смотрят на \(R^2\), считая его самой важной статистикой, однако, если модель является незначимой, то все остальные её показатели, какими бы хорошими они не были, не имеют никакого значения.

Самые первые строки, которые выдаёт команда summary,содержат в себе описание модели, которая была построена с помощью команды lm.

Call:
lm(formula = dist ~ speed, data = cars)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-29.069  -9.525  -2.272   9.215  43.201 

это условие не соблюдается. Медиана (Median) смещена влево относительно нуля, первый квартиль (1Q) не совсем симметричен относительно третьего (3Q), с незначительным смещением влево, в то время как для нормального распределения характерна симметрия. Максимумы (Max) и минимумы (Min) же позволяют достаточно быстро оценить наличие потенциальных промахов в исходных данных, поскольку промахи в исходных данных порождают промахи в остатках модели.

Следующий блок содержит информацию о коэффициентах модели и их значимости.

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791     6.7584  -2.601   0.0123 *  
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Колонка Estimate содержит подобранные коэффициенты регрессии \(\tilde{\beta}_j\), колонка Std. Error содержит их среднеквадратичные отклонения \(s_j\). Колонка t value рассчитывается как их отношение (без учёта модуля, который будет учтен далее). В колонке Pr(>|t|) показана вероятность того, что данный коэффициент является незначимым. Для удобства интерпретации значимость коэффициентов также отображается символами справа в соответствии с легендой под чертой.

Дальше идут знакомые вам Эр-квардрат и скорретированный Эр-квадрат.

Поговорим о статистике R-квадрат или, как её иногда называют, коэффициента детерминации. Она показывает, насколько условная дисперсия модели отличается от дисперсии реальных значений Y. Более точно — это единица минус доля необъяснённой дисперсии (дисперсии случайной ошибки модели, или условной по факторам дисперсии зависимой переменной) в дисперсии зависимой переменной. Его рассматривают как универсальную меру зависимости одной случайной величины от множества других. В частном случае линейной зависимости \(R^2\) является квадратом так называемого множественного коэффициента корреляции между зависимой переменной и объясняющими переменными. В частности, для модели парной линейной регрессии коэффициент детерминации равен квадрату обычного коэффициента корреляции между y и x. Если этот коэффициент близок к 1, то условная дисперсия модели достаточно мала и весьма вероятно, что модель неплохо описывает данные.Для приемлемых моделей предполагается, что коэффициент детерминации должен быть хотя бы не меньше 50% (в этом случае коэффициент множественной корреляции превышает по модулю 70%). Модели с коэффициентом детерминации выше 80% можно признать достаточно хорошими (коэффициент корреляции превышает 90%). Значение коэффициента детерминации 1 означает функциональную зависимость между переменными.

Однако, у статистики R-квадрат есть один серьезный недостаток: при увеличении числа предикторов эта статистика может только возрастать. Поэтому, может показаться, что модель с большим количеством предикторов лучше, чем модель с меньшим, даже если все новые предикторы никак не влияют на зависимую переменную. Тут можно вспомнить про принцип бритвы Оккама. Следуя ему, по возможности, стоит избавляться от лишних предикторов в модели, поскольку она становится более простой и понятной. Для этих целей была придумана статистика скорректированный R-квадрат (скорректированный коэффициент детерминации). Она представляет собой обычный R-квадрат, но со штрафом за большое количество предикторов. Основная идея: если новые независимые переменные дают большой вклад в качество модели, значение этой статистики растет, если нет — то наоборот уменьшается.

Для построения диагностических графиков остатков используется команда plot.

plot(lm1)

Для построения доверительного интервала коэффициентов команда confint

confint(lm1)
##                  2.5 %    97.5 %
## (Intercept) -31.167850 -3.990340
## speed         3.096964  4.767853

В качестве иллюстрации команд построим модель без свободного члена и выведем параметры модели, причем сделаем это одной командой

summary(lm2<-lm(dist~speed-1,data=cars))
## 
## Call:
## lm(formula = dist ~ speed - 1, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.183 -12.637  -5.455   4.590  50.181 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## speed   2.9091     0.1414   20.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.26 on 49 degrees of freedom
## Multiple R-squared:  0.8963, Adjusted R-squared:  0.8942 
## F-statistic: 423.5 on 1 and 49 DF,  p-value: < 2.2e-16

И регрессионную зависимость от квадрата скорости

summary(lm3<-lm(dist~I(speed^2)-1,data=cars))
## 
## Call:
## lm(formula = dist ~ I(speed^2) - 1, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.350  -7.988   1.325   8.080  49.939 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## I(speed^2) 0.153374   0.007122   21.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.61 on 49 degrees of freedom
## Multiple R-squared:  0.9044, Adjusted R-squared:  0.9025 
## F-statistic: 463.8 on 1 and 49 DF,  p-value: < 2.2e-16
plot(lm3)

Устраним 49 измерение,содержащее выброс и построим регрессию заново.

cars2<-cars[-49,]
summary(lm4<-lm(dist~I(speed^2)-1,data=cars2))
## 
## Call:
## lm(formula = dist ~ I(speed^2) - 1, data = cars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.718  -6.372   1.780   8.769  50.738 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## I(speed^2) 0.149295   0.007114   20.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.04 on 48 degrees of freedom
## Multiple R-squared:  0.9017, Adjusted R-squared:  0.8997 
## F-statistic: 440.4 on 1 and 48 DF,  p-value: < 2.2e-16
plot(lm4)

Список источников

  1. Сергеев А. Г., Крохин В. В. Метрология: Учеб. пособие для вузов.— М.: Логос, 2001. — 408 с. — ISBN: 5-94010-039-2.
  2. Hilbe J. M. Negative Binomial Regression. – 1st. edition. – Cambridge Universit Press, 2007. – 251 p. – ISBN: 978-0521857727. 3.Legendre P., Legendre L. F. J. Numerical Ecology (Developments in EnvironmentalbModelling) [Kindle Edition]. – Elsevier Science, 1998. – 870 p. – ASIN : B0092L8KFA.
  3. Graphical Methods for Data Analysis (Wadsworth & Brooks/Cole Statistics/Probability Series) / John M. Chambers, William S. Cleveland, Paul A. Tukey, Beat Kleiner. – Duxbury Press, 1983. – 395 p. – ISBN: 978-0534980528.
  4. Hintze J.L., Nelson R.D. Violin Plots: A Box Plot-Density Trace Synergism / Jerry L. Hintze, Ray D. Nelson // The American Statistician. – 1998. – Vol. 52. – P. 181-184
  5. Куликов Е. И. Прикладной статитический анализ. Учебное пособие для вузов. — 2-е. изд. — М. : Горячая линия–Телеком, 2008. — 464 с. — ISBN: 978-5-99120021-9. 7.Lewis N. D. 100 Statistical Tests in R [Kindle Edition]. – Heather Hills Press, 2013. – 496 p. – ASIN : B00C143F14. 8.Komsta L. – outliers: Tests for outliers, 2011. – R package version 0.14. URL: http://CRAN.R-project.org/package=outliers