Анализ аномальных наблюдений

Количественный анализ

Влияние аномальных наблюдений

#сгенерируем выборку из N(1,2)
x=rnorm(100,1,2)
#медиана и среднее не различаются различаются
mean(x);median(x)
## [1] 0.938
## [1] 0.9912
#изменим одно наблюдение на аномально большое
y=x;y[100]=100
#среднее существенно изменилось, а медиана осталась неизменной
#медиана - робастная (устойчивая к исккажениям) статистика
mean(y);median(y)
## [1] 1.95
## [1] 1.009
#изменение дисперсии
var(x);var(y)
## [1] 3.611
## [1] 101.7

Модель Тьюки-Хьюбера

#сгенерируем смесь двух распределений: N(1,2) и Exp(1/10)
n=100
p=0.05 #доля наблюдений из Exp(1/10)
n2=floor(n*p);n1=n-n2
#начальное значение датчика псевдослучайных чисел: для воспроизводимости результатов
set.seed(9995) 
x=c(rnorm(n1,1,2),rexp(n2,1/10))

#выведем "ящик с усами"
#пять аномальных наблюдений, принадлежащих Exp(1/10) отмечены точками вверху
boxplot(x)

plot of chunk unnamed-chunk-2

#просмотрим данные наблюдения
tail(x,5)
## [1] 15.334  4.072 14.418 12.686  7.036
#среднее и медиана сильно отличаются: данные неоднородны
mean(x);median(x)
## [1] 1.478
## [1] 1.133
#используем усеченное среднее
#из вектора x "выбрасываются" по 5% наим. и наиб. наблюдений
#статитика устойчива, но теряется часть информации
mean(x,trim = 0.05)
## [1] 1.186

Робастные статистики

x=rnorm(100,1,2)
y=x;y[100]=100
#медиана - робастная статика для центра распределения значений
median(x);median(y)
## [1] 0.8811
## [1] 0.8811
#робастная статистика разброса значений
#абсолютное медианное отклонение - 
mad(x);mad(y)
## [1] 1.933
## [1] 1.933
#межквартильных размах
IQR(x);IQR(y)
## [1] 2.341
## [1] 2.345

Графический анализ

#сгенерируем n=1000 наблюдений из РВ N(-3,1) и N(3,2)
x=c(rnorm(900,-3,1),rnorm(100,3,2))
#фактор - индикатор РВ
d=as.factor(rep(1:2,c(900,100)))
levels(d)=c("N(-3,1)","N(3,2)")
boxplot(x~d)

plot of chunk unnamed-chunk-4

Проверка гипотез о равенстве средних

Базовые понятия

Тестовая статистика - характеристика данных, которая вычисляется по выборке наблюдений.
Например, статистика t-критерия (t) является нормированной разницей средних и вычисляется по формуле:
t=(mean(v1)-mean(v2))/sqrt(var(v1)/length(v1)+var(v2)/length(v2)),
где v1, v2 - вектора значений двух сравниваемых выборок соответственно.

p-value - вероятность того, что тестовая статистика могла бы принять значение больше фактической (выборочной статистики, оцененной в данном эксперименте) при условии нулевой гипотезы H0.

Возможные решения статистического теста:
1. p-value >= a (уровень значимости, обычно a=0.05): гипотеза H0 не отклоняется.
2. p-value < a : гипотеза H0 отклоняется - принимается альтернативная гипотеза H1.

t-критерий проверяет гипотезу о равенстве средних значений (мат. ожидания) для двух выборок.

#проверка распределения
x=rnorm(100,mean=1)
#проверим, что среднее действительно равно 1 (p-value>0.05 - не отклоняется)
t.test(x,mu=1)
## 
##  One Sample t-test
## 
## data:  x
## t = 0.5967, df = 99, p-value = 0.552
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
##  0.8574 1.2653
## sample estimates:
## mean of x 
##     1.061
#предположение, что среднее равно 0 (p-value<0.05 - отклоняется)
t.test(x,mu=0)
## 
##  One Sample t-test
## 
## data:  x
## t = 10.33, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.8574 1.2653
## sample estimates:
## mean of x 
##     1.061

Модельные данные

#две независимые выборки с различным средним значением
x=rnorm(100,mean=1);y=rnorm(100,mean=2)
#изобразим x и y с помощью "ящиков с усами"
boxplot(x,y)

plot of chunk unnamed-chunk-6

#сравним две выборки на равенство среднего значения
t.test(x,y) #(p-value<0.05 - гипотеза о равенстве средних отклоняется)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -8.316, df = 196.6, p-value = 1.493e-14
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.387 -0.855
## sample estimates:
## mean of x mean of y 
##     1.030     2.151

Реальные данные

Данные о лечении электро-магнитными волнами: “magnets.csv”

Данные о 50 пациентах, испытывающих боли в результате заболевания полиомиелитом. Пациенты получили различное лечение в каждой из двух групп. Первая группа пациентов подверглась облучению электро-магнитными волнами (“активное лечение” при включенном приборе), а остальная часть - получила “пассивное лечение” в виде плацебо (прибор был выключен). Переменные: score1, score2 - интенсивности боли (в баллах от 0 до 10, чем выше балл, тем сильнее боль) для первой и второй частей пациентов соответственно, change - разница между score1 и score2, active - индикатор использования эл.-магн. волн (1 - прибор вкл., 2 - выкл.: плацебо).

#загрузим данные
magnets=read.csv("magnets.csv")
attach(magnets)
## The following object is masked _by_ .GlobalEnv:
## 
##     active
#выведем число пациентов в каждой из групп
table(active)
## active
##  1  2 
## 29 21

Двухвыборочный t-критерий

Для независимых выборок (в данном примере выборки - это показания переменной change в двух разных группах пациентов: ни один из пациентов не находится в двух группах одновременно)

#выведем "ящики с усами"  по переменной change для каждой группы отдельно
boxplot(change~active)

plot of chunk unnamed-chunk-8

mean(change[active=="1"]) #активное лечение
## [1] 5.241
mean(change[active=="2"]) #пассивное лечение
## [1] 1.095

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

#Применим t-критерий для сравнения средних переменной  change
#"~": в двух группах согласно переменной active

t.test(change~active)
## 
##  Welch Two Sample t-test
## 
## data:  change by active
## t = 5.986, df = 42.93, p-value = 3.86e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.749 5.543
## sample estimates:
## mean in group 1 mean in group 2 
##           5.241           1.095

Вывод: p-value близко к нулю, поэтому принимается альтернативная гипотеза о том, что среднее изменение уровня боли в двух группах различается существенно (в том числе на уровне 0.05)

Парный t-критерий

Для зависимых выборок (в данном примере выборки - это показания переменных score1 и score2 для группы одних и тех же пациентов)

#проверим, является ли эффект плацебо значимым, т.е. изменилось ли статистически
#значимо среднее значение уровня боли во второй группе (active=2)
#опция "paired = TRUE" указывает на использование парного критерия
t.test(score1[active=="2"],score2[active=="2"],paired = TRUE)
## 
##  Paired t-test
## 
## data:  score1[active == "2"] and score2[active == "2"]
## t = 3.18, df = 20, p-value = 0.004702
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3769 1.8136
## sample estimates:
## mean of the differences 
##                   1.095

Вывод: p-value близко к нулю, поэтому можно считать, что эффект плацебо также способствует улучшению состояния пациентов, хотя по сравнению с активным лечением незначителен.