Конспект занятий по языку R

Часть ноль - импорт

Задача по импортированию файлов решается через вкладку Enviroment и клавишу import. Это позволяет добавить необходимые данные в пути поиска RStudio. Имя задаем по-короче, его надо будет постоянно использовать.

Начнем, основная команда получения синопсиса по загруженным данным. Начало начал, удобно понять как R видит ваши данные. Также можно использовать продвинутый вариант описательной статистики с группирующей переменной.

GuideDATA <- ChickWeight
summary(GuideDATA)
#Группирующая переменная идет через запятую
describeBy (GuideDATA$weight, GuideDATA$Diet)

Данные не всегда определяются правильным образом, по этому иногда ему надо подсказать. Каждый вектор можно назначить переменную фактором, числовой, текстовой, логической, комплексной, необработанной (байты).

#Столбец dose сделать фактором
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
#Столбец color сделать непрерывной переменной
GuideDATA$color <- as.numeric(GuideDATA$color)

Часть первая - линейная статистика

Проверка параметров

Первая задача, проверка распределения на нормальность. Используется: Шапиро-Уилка, гистограмма и Q-Q plot. Далее проверка равенства дисперсий

with(GuideDATA, {
  hist(weight,col="blue",xlab="Вес")
  qqnorm(weight)
  shapiro.test(weight)
})

    Shapiro-Wilk normality test

data:  weight
W = 0.90866, p-value < 2.2e-16

Shapiro-Wilk p-value меньше 0,05 - отвержение H₀, распределение ненормальное, что видно на гистограммах, а Q-Q plot в случае нормального распределения показывает ровную линию от нуля.

Попарное сравнение

Простейший тест Т-тест Стьюдента для попарного сравнения средних или непараметрический Критерий Уилкоксона.

#Первой указывается зависимая переменная, вторым независимая переменная, далее ресурс
t.test(len~supp, data = ToothGrowth)

    Welch Two Sample t-test

data:  len by supp
t = 1.9153, df = 55.309, p-value = 0.06063
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 -0.1710156  7.5710156
sample estimates:
mean in group OJ mean in group VC 
        20.66333         16.96333 

#Непраметрический аналог
wilcox.test(len~supp,data=ToothGrowth)

Предупреждение: не могу подсчитать точное p-значение при наличии повторяющихся наблюдений

    Wilcoxon rank sum test with continuity correction

data:  len by supp
W = 575.5, p-value = 0.06449
alternative hypothesis: true location shift is not equal to 0

#Классическая визуализация ящиками
boxplot(len~supp,data=ToothGrowth) 

Мы видим что t больше 0.05 занчит сохраняем H₀ обе выборки из одной генеральной совокупности. df показывает степени свободы.
Тест применяется с автоматической поправкой Уэма, что позволяет его применять на группы с разной дисперсией, отключается добавлением аргумента ...,var.equal=TRUE

Множественное сравнение

Итак, множественный анализ дисперсий. Главный вопрос - в какую сторону дисперсия в между группами отличается от внутригрупповой дисперсии?
\[ F=\frac{\text{Дисперсия между группами}}{\text{Дисперсия внутри групп}} \]
Если F-критерий Фишера больше критического – Н₀ о равенстве средних в группах отвергаем, внутри есть различия. Требования к выбокам:

1. Надо стремиться к равенству размеров групп. Неравенство размеров сразу делает результаты анализа уязвимыми к отклонениям от требований нормальности и гомогенности.

2. Равенство дисперсий

3. Нормальное распределение в каждой группе
   Если никак, обходим через Краскел-Уолеса, требует только равенства дисперсий kruskal.test(len~dose,data=ToothGrowth)

# Проверка нормальности распределения в каждой группе, нормальность распределения остатков.
shapiro.test(resid(model))

    Shapiro-Wilk normality test

data:  resid(model)
W = 0.96457, p-value = 0.07884

# Проверка равенства дисперсий
leveneTest(len ~ dose, data=ToothGrowth, center="mean")

Levene's Test for Homogeneity of Variance (center = "mean")
      Df F value Pr(>F)
group  2  0.7328  0.485
      57               

# Графичеки
qqPlot(lm(len ~ dose, data = ToothGrowth)) 

[1] 23 32

Если все хорошо - в путь.

model <- aov(len~dose,data=ToothGrowth)
summary(model)

            Df Sum Sq Mean Sq F value   Pr(>F)    
dose         2   2426    1213   67.42 9.53e-16 ***
Residuals   57   1026      18                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Просмотр уровней влияния
TukeyHSD(model)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = len ~ dose, data = ToothGrowth)

$dose
      diff       lwr       upr    p adj
2-1  9.130  5.901805 12.358195 0.00e+00
3-1 15.495 12.266805 18.723195 0.00e+00
3-2  6.365  3.136805  9.593195 4.25e-05

# Простейшая визуализация
plotmeans(len ~ dose, xlab = "Vit. C level", ylab = "Teeth lenght", data = ToothGrowth)

Сначала смотрим на Pr(\>F) - меньше 0.05, значения различаются. F value - насколько сильно различаются (в n раз). В TukeyHSD смотрим на графу diff.

Корреляции и регрессии

Корреляции

Начнем с общих параетров. Коффециент Пирсона, просто и сердито. Оценивает силу линейной связи.

attach(trees)
  # Полный вывод
  cor.test(Height, Volume)

    Pearson's product-moment correlation

data:  Height and Volume
t = 4.0205, df = 29, p-value = 0.0003784
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3095235 0.7859756
sample estimates:
      cor 
0.5982497 

  # Линейный график корреляции, диаграмма рассеивания
  scatterplot(Height, Volume)

detach()

Вывод прос – p-value значимсоть корреляции, cor коэффициент.

Регрессии

Общая формула для построения регрисионной модели:
\[ Y_i=\alpha + \beta X_i + \varepsilon_i \\ \text{Свободный член } + \text{ Регресионный коф. } +\text{ Остатки, ошибки} \] Число параметров сравнения не стоит раздувать. Каждый параметр не только требудет +10 вариантов к выборке, но и увеличивает число моделей к построению в квадрате (все варианты нуждно сравнить со всеми).
Идея проста, сравниваем две изменчивости - с включением остаточной изменчивости или без. Грубо говоря сраниваем сравниваем идеальную модель с реальной. Изменчивости сравниваются:

1. Сравниваем F-критерием, обсуждаемым выше.
   p<0.05 отвергаем H₀ влияние предиктора на переменную значимо

2. Т-статистика попарного сравнения, значимость различий бетта коэффициент от нуля
   бетта отличен от 0 p<0.05 отвергаем H₀ значимое влияние фактора

3. Расчет доверительных интервалов коэффициентов
   

   • Оба 95% CI интервала положительны, не включают ноль – эффект сильный
   • Оба 90% CI интервала положительны, не включают ноль – эффект слабый но есть
   • В противном случае – соси

Теперь в коде, не забываем что распределение должно быть нормальным, а дисперсии гомогенными.

model <- lm(len~dose,data=ToothGrowth)
# Проверка нормальности распределения в каждой группе, нормальность распределения остатков.
shapiro.test(resid(model))

    Shapiro-Wilk normality test

data:  resid(model)
W = 0.96457, p-value = 0.07884

summary(model)

Call:
lm(formula = len ~ dose, data = ToothGrowth)

Residuals:
   Min     1Q Median     3Q    Max 
-8.061 -3.274 -1.061  3.063 10.434 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.3183     1.4542   2.282   0.0262 *  
dose          7.7475     0.6731  11.509   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.257 on 58 degrees of freedom
Multiple R-squared:  0.6955,    Adjusted R-squared:  0.6902 
F-statistic: 132.5 on 1 and 58 DF,  p-value: < 2.2e-16

# Получить доверительные интервалы
confint(model)

                2.5 %   97.5 %
(Intercept) 0.4075019 6.229165
dose        6.4000469 9.094953

F-statistic: 67.42 on 2 and 57 DF, p-value: 9.533e-16 p < 0.05 модель адекватно оценивает данные, если фактор один, то это также означает, что предиктор значительно влияет на выборку.
Multiple R-squared: 0.7029, Adjusted R-squared: 0.6924 - второй коэффициент скоректирован, это коффициент детерминации. Показывает какая доля дисперсии зависимой переменной определяется нашим предиктором. Типа 69,2% всей выбьорки описывается нашей моделью.
(Intercept) – свободный член (длинна зуба при 0 дозе).
Estimate – бетта коэффициент. Std. Error – стандартная ошибка.
Pr(>|t|) – с помощью Стюдента сравниваем бетта коэффциент с нулем.

Изначально доверительные интервалы мы строим для 95%. Если надо построить интервалы для 90% ...,level=0.9
Также хорошо было бы построить диаграмму рассеивания, но ни на уроке, ни у меня сейчас по этой моделе ничего вразумительного получить не удалось. Общая идея такой диаграммы - рассиевание должно быть равным, без каких-либо выплесков.

Теперь попробуем участь две независимые переменные.

model <- lm(Volume~Girth+Height, data=trees)
summary(model)

Call:
lm(formula = Volume ~ Girth + Height, data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4065 -2.6493 -0.2876  2.2003  8.4847 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
Girth         4.7082     0.2643  17.816  < 2e-16 ***
Height        0.3393     0.1302   2.607   0.0145 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.882 on 28 degrees of freedom
Multiple R-squared:  0.948, Adjusted R-squared:  0.9442 
F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

confint(model)

                   2.5 %      97.5 %
(Intercept) -75.68226247 -40.2930554
Girth         4.16683899   5.2494820
Height        0.07264863   0.6058538

Вывод уже понятный. Важно помнить, что полученная модель не может быть “распилена” по коэффициентам. То что мы получили работает только в полном сборе. По этому, если фактора у нас два, то нужно строить 4 моедли: оба фактора вместе влияют, кажждый по отдельности, и ни один фактор не влияет. После того, как мы получили все 4 модели в классическом виде выбирается по F-статистике и по дисперсиям. Но одной золотой модели врят-ли найти.

Есть и ещё одна проблема - автокорреляция. Когда два фактора влияние которых мы сравниваем, сильно кореллируют друг с другом. В этом случае проведем тест корреляции.

cor.test(trees$Girth, trees$Height)

    Pearson's product-moment correlation

data:  trees$Girth and trees$Height
t = 3.2722, df = 29, p-value = 0.002758
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2021327 0.7378538
sample estimates:
      cor 
0.5192801 

p-value < 0.05 корреляция есть, cor = 0.52 да и при том огромная. Следовательно линейную модель по этим двум факторам строить нельяз.
Но можно немного схитрить. Вычесть влияние одного фактора их другого.

model <- lm(Girth~Height, data=trees)
# Часть переменной обхват объясняется высотой дерева, а вот остаточная изменчивость не объясняется высостой дерева. Эти остатки мы и возьмем
model <- lm(Volume~Girth+resid(model), data=trees)
summary(model)

Call:
lm(formula = Volume ~ Girth + resid(model), data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4065 -2.6493 -0.2876  2.2003  8.4847 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -49.7787     5.8039  -8.577 2.54e-09 ***
Girth          6.0347     0.4349  13.876 4.50e-14 ***
resid(model)  -1.3265     0.5089  -2.607   0.0145 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.882 on 28 degrees of freedom
Multiple R-squared:  0.948, Adjusted R-squared:  0.9442 
F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

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

Но если независимая переменная в модели фактор, как анализируется такая переменная? В общем-то из таблицы с данными берется первый фактор и сравнивается с отальными.

ToothGrowth$dose <- as.factor(ToothGrowth$dose)
model <- lm(len~supp+dose, data = ToothGrowth)
summary(model)

Call:
lm(formula = len ~ supp + dose, data = ToothGrowth)

Residuals:
   Min     1Q Median     3Q    Max 
-7.085 -2.751 -0.800  2.446  9.650 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.4550     0.9883  12.603  < 2e-16 ***
suppVC       -3.7000     0.9883  -3.744 0.000429 ***
dose2         9.1300     1.2104   7.543 4.38e-10 ***
dose3        15.4950     1.2104  12.802  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.828 on 56 degrees of freedom
Multiple R-squared:  0.7623,    Adjusted R-squared:  0.7496 
F-statistic: 59.88 on 3 and 56 DF,  p-value: < 2.2e-16

Каждый фактор сравнивается с первым из таблицы и коэффициент показывает насколько все последующие факторы отличаются по влиянию от первого. При этом если вместо + использовать знак * то проявляется эффект взаимодействия. То есть второй вариант для варианта где доза и тип препарата влияют на длинну зуба зависимо. А первый вариант, когда доза влияет независимо от препарата, не важно витамин это или апельсиновый сок.