Task1

1. Выбор данных

https://vincentarelbundock.github.io/Rdatasets/datasets.html - здесь есть множество разных датасетов с описанием. При этом можно предварительно посмотреть на количество наблюдений, количество признаков, сколько из них являются факторами и т.д.

2. Считывание и просмотр

df <- read.csv(file = "Angell.csv", header = TRUE, as.is = FALSE) #Скачиваем csv, сохраняем в папку, где лежит r markdown и считываем данные. 
head(df) #Head() Показывает первые шесть строк данных.

3. Описание данных

Рассматриваемые данные содержат различную информацию о 43 городах в Америке около 1950г.

moral - Moral integration - Composite of crime rate and welfare expenditures

hetero - Ethnic Heterogenity - From percentages of nonwhite and foreign-born white residents

mobility - Geographic Mobility - From percentages of residents moving into and out of the city

region - A factor with levels: E Northeast; MW Midwest; S Southeast; W West.

4. Типы признаков

moral - количественный hetero - количественный mobility - количественный region - качественный

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.0.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

summarize(df, across(moral:region, function(x) max(table(x)))) #считает сколько раз встречается мода у каждого признака

# df %>% summarize(across(moral:region, function(x) max(table(x)))) #считает сколько раз встречается мода у каждого признака

Из таблицы выше можем заключить, что hetero - абсолютно непрерывный признак. Кроме того, можем считать, что moral и mobility ближе к непрерывным признакам.

5. Порядковые признаки

Здесь необходимо проверить соответствие текстовых меток порядкового признака (если такие использованы) их естественному порядку.
В данном датасете таких признаков нет.

6. Matrix plot, outliers, etc.

library(lattice)
panel.hist.splom <- function(x, ...) {
  yrng <- current.panel.limits()$ylim
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks
  nB <- length(breaks)
  y <- h$counts
  y <- yrng[1] + 0.95 * diff(yrng) * y / max(y)
  panel.rect(breaks[-nB], yrng[1], breaks[-1], y, col="cyan", ...)
}
splom(~subset(df, select=c(moral:mobility)), 
      data = df, groups = region,
      diag.panel = panel.hist.splom, pch = 21, auto.key = list(space="top", columns=4))

7. Симметричность распределений.

Из матрикс плота видно, что распределение hetero - сильно несимметричное с хвостом вправо, поэтому прологарифмируем его и построим заново матрикс плот.

dfl <- transform(df, hetero=log(hetero))

splom(~subset(dfl, select=c(moral:mobility)), 
      data = dfl, groups = region,
      diag.panel = panel.hist.splom, pch = 21, auto.key = list(space="right", columns=1, title = "Region"))

8. Аутлаеры

В данном случае на одном скаттерплоте (hetero & moral) видна линейная зависимость, не имеющая аутлаеров относительно нее. На скаттерплоте (mobility & hetero) отсутствует зависимость и, соответственно, аутлаеры. На (moral & mobility) также линейная зависимость, но вместо аутлаеров наблюдается скорее неоднородность.

10. Матрикс плот и аутлаеры для отдельной группы.

splom(~subset(dfl, select=c(moral:mobility)), 
      data = dfl, subset = region == 'S')

dflo <- dfl #сохраняем датасет
dflo[38,] <- NA #удаляем аутлаер относительно линейной зависимости между hetero и mobility
splom(~dflo[,2:4], 
      data = dflo, subset = region == 'S')

11. Descriptive statistics

summary(df)

##           X          moral           hetero         mobility     region 
##  Akron     : 1   Min.   : 4.20   Min.   :10.60   Min.   :12.10   E : 9  
##  Atlanta   : 1   1st Qu.: 8.70   1st Qu.:16.90   1st Qu.:19.45   MW:14  
##  Baltimore : 1   Median :11.10   Median :23.70   Median :25.90   S :14  
##  Birmingham: 1   Mean   :11.20   Mean   :31.37   Mean   :27.60   W : 6  
##  Bridgeport: 1   3rd Qu.:13.95   3rd Qu.:39.00   3rd Qu.:34.80          
##  Buffalo   : 1   Max.   :19.00   Max.   :84.50   Max.   :49.80          
##  (Other)   :37

summary(dfl)

##           X          moral           hetero         mobility     region 
##  Akron     : 1   Min.   : 4.20   Min.   :2.361   Min.   :12.10   E : 9  
##  Atlanta   : 1   1st Qu.: 8.70   1st Qu.:2.827   1st Qu.:19.45   MW:14  
##  Baltimore : 1   Median :11.10   Median :3.165   Median :25.90   S :14  
##  Birmingham: 1   Mean   :11.20   Mean   :3.267   Mean   :27.60   W : 6  
##  Bridgeport: 1   3rd Qu.:13.95   3rd Qu.:3.663   3rd Qu.:34.80          
##  Buffalo   : 1   Max.   :19.00   Max.   :4.437   Max.   :49.80          
##  (Other)   :37

library(moments)
summarize(dfl, across(moral:mobility, list(kurtosis = kurtosis, skewness = skewness)))

2.1 Выбор категоризующей переменной

В качестве категоризующего признака возьмем region. Сравнивать будем центральные штаты и юго-восточный регион, как имеющие наибольшее количество наблюдений (14 и 13, соответственно).

2.2 Boxplot

dfcomp <- dflo %>% filter(region == "S" | region == "MW") #Выбираем из дата фрейма только наблюдения из юго-восточного и центрального регионов. Filter - вытаскивает подвыборку дата фрейма. %>% - передает dflo как первый параметр функции, идущей после этого знака.
dfcomp$region <- droplevels(dfcomp$region) #Изменяет параметр-фактор, отвечающий за принадлежность городу региона, оставляя только присутствующие регионы (в нашем случае юго-восточный и центральный)
bwplot(moral ~ region, data = dfcomp, col = c("forestgreen", "gold"), main = "Moral Integration", xlab = "Region") #строим боксплоты моральной интеграции по региону. moral ~ region - моральная интеграция как главный параметр, ~ region - разбивает боксплоты на группы в зависимости от региона. data - передаем наш дата фрейм. col - раскрашивает боксплоты в зависимости от группы. main - название над картинкой. xlab - подпись оси х.

#Далее, аналогичные картинки относительно двух других признаков.

bwplot(hetero ~ region, data = dfcomp, col = c("forestgreen", "gold"), main = "Ethnic Heterogenity", xlab = "Region")

bwplot(mobility ~ region, data = dfcomp, col = c("forestgreen", "gold"), main = "Geographic Mobility", xlab = "Region")

2.3 Нормальность признаков

library(tidyr)

## Warning: package 'tidyr' was built under R version 4.0.2

qqmath(~value | name, data = pivot_longer(dfcomp, moral:mobility), subset = region == "S")

qqmath(~value | name, data = pivot_longer(dfcomp, moral:mobility), groups = region, auto.key = TRUE,
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, ...)
          panel.qqmath(x, ...)
       })

shapiro.test(subset(dfcomp, region == "S")$moral)  #Критерий шапиро-уилка позволяет проверить гипотезу о том, что случайная величина имеет нормальное распределение.

## 
##  Shapiro-Wilk normality test
## 
## data:  subset(dfcomp, region == "S")$moral
## W = 0.88893, p-value = 0.09432

shapiro.test(subset(dfcomp, region == "S")$mobility)

## 
##  Shapiro-Wilk normality test
## 
## data:  subset(dfcomp, region == "S")$mobility
## W = 0.8902, p-value = 0.09816

#Далее, аналогично для другого региона

qqmath(~value | name, data = pivot_longer(dfcomp, moral:mobility), subset = region == "MW")

shapiro.test(subset(dfcomp, region == "MW")$moral)  #Критерий шапиро-уилка позволяет проверить гипотезу о том, что случайная величина имеет нормальное распределение.

## 
##  Shapiro-Wilk normality test
## 
## data:  subset(dfcomp, region == "MW")$moral
## W = 0.95037, p-value = 0.5665

shapiro.test(subset(dfcomp, region == "MW")$mobility)

## 
##  Shapiro-Wilk normality test
## 
## data:  subset(dfcomp, region == "MW")$mobility
## W = 0.90973, p-value = 0.1562

2.4 t-test, критерий Манна-Уитни

t.test(moral ~ region, data = dfcomp)$p.value #Т-тест критерий позволяет проверить гипотезу о равенстве средних по одному признаку у двух групп. Moral - признак, по которому сравниваем, region - признак, по которому разбиваем на две группы. $p.value - печатает только p-value. Иначе печатает все результаты, предоставляемые тестом.

## [1] 2.062762e-05

t.test(hetero ~ region, data = dfcomp)

## 
##  Welch Two Sample t-test
## 
## data:  hetero by region
## t = -5.3758, df = 23.917, p-value = 1.628e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.249586 -0.556182
## sample estimates:
## mean in group MW  mean in group S 
##         2.999347         3.902231

t.test(mobility ~ region, data = dfcomp)

## 
##  Welch Two Sample t-test
## 
## data:  mobility by region
## t = -2.5619, df = 24.245, p-value = 0.01704
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.366225  -1.442566
## sample estimates:
## mean in group MW  mean in group S 
##         26.05714         33.46154

wilcox.test(mobility ~ region, data = dfcomp) #непараметрический аналог, работает аналогичным образом.

## 
##  Wilcoxon rank sum exact test
## 
## data:  mobility by region
## W = 36, p-value = 0.006618
## alternative hypothesis: true location shift is not equal to 0

2.5 Критерий Колмогорова-Смирнова

ks.test(dfcomp[dfcomp$region == "S",4], dfcomp[dfcomp$region == "MW", 4]) #Критерий Колмогорова-Смирнова для проверки гипотезы о равенстве функций распределения двух случайных величин.

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  dfcomp[dfcomp$region == "S", 4] and dfcomp[dfcomp$region == "MW", 4]
## D = 0.57143, p-value = 0.009595
## alternative hypothesis: two-sided

3.1 Анализ зависимостей по группе.

Посмотрим на матрикс плот данных по юго-восточному региону.

splom(~dflo[,2:4], 
      data = dflo, subset = region == 'S') #Пользуемся ранее построенным матрикс плотом в пункте 10.

3.2 Коэффициент корреляции Пирсона.

cor.test(~ hetero + mobility, data = dflo, subset = region == 'S', method = "pearson") #Считает коэффициент корреляции Пирсона между двумя столбцами матрицы.

## 
##  Pearson's product-moment correlation
## 
## data:  hetero and mobility
## t = -5.7287, df = 11, p-value = 0.0001325
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9590838 -0.6010331
## sample estimates:
##        cor 
## -0.8654242

library("Hmisc") #Библиотека, позволяющая строить корреляционные м атрицы вместе с матрицами соответствующих p-value теста о значимости коэффициента корреляции.

## Warning: package 'Hmisc' was built under R version 4.0.2

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.0.2

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

rcorr(as.matrix(subset(dflo, region == 'S', select=moral:mobility))) #Выдает корреляционную матрицу, а также матрицу подсчитанных p-value для теста о значимости полученного коэффициента корреляции (Нулевая гипотеза: коэффициент корреляции не значим).

##          moral hetero mobility
## moral     1.00  -0.25    -0.10
## hetero   -0.25   1.00    -0.87
## mobility -0.10  -0.87     1.00
## 
## n= 13 
## 
## 
## P
##          moral  hetero mobility
## moral           0.4116 0.7443  
## hetero   0.4116        0.0001  
## mobility 0.7443 0.0001

Таким образом, гипотеза о том, что коэффициент корреляции между hetero и mobility равняется нулю, отвергается.

3.3 Коэффициент корреляции Спирмена

cor.test( ~ hetero + mobility, data = dflo, subset = region == "S", method = "spearman") #Считает коэффициент корреляции Спирмена между двумя столбцами матрицы.

## 
##  Spearman's rank correlation rho
## 
## data:  hetero and mobility
## S = 644, p-value = 0.00325
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.7692308

cor.test( ~ moral + mobility, data = dflo, subset = region == "S", method = "spearman")

## Warning in cor.test.default(x = c(10.4, 10.2, 10.2, 9.7, 9.3, 8.6, 8.2, : Cannot
## compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  moral and mobility
## S = 455.13, p-value = 0.4094
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2503441

rcorr(as.matrix(subset(dflo, region == 'S', select=moral:mobility)), type = "spearman")

##          moral hetero mobility
## moral     1.00  -0.26    -0.25
## hetero   -0.26   1.00    -0.77
## mobility -0.25  -0.77     1.00
## 
## n= 13 
## 
## 
## P
##          moral  hetero mobility
## moral           0.3936 0.4094  
## hetero   0.3936        0.0021  
## mobility 0.4094 0.0021

3.4 Причины, следствия и частные корреляции.

cor.test(~ hetero + mobility, data = dfl, method = "pearson") #смотрим на проверку гипотезы о значимости коэффициента корреляции между hetero и mobility для всех групп.

## 
##  Pearson's product-moment correlation
## 
## data:  hetero and mobility
## t = -0.7233, df = 41, p-value = 0.4736
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3991340  0.1946604
## sample estimates:
##        cor 
## -0.1122459

library(ppcor) #Библиотека, позволяющая проверять гипотезы о значимости коэффициента частной корреляции.

## Warning: package 'ppcor' was built under R version 4.0.2

## Loading required package: MASS

## Warning: package 'MASS' was built under R version 4.0.2

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

dft <- subset(dfl, region != "W")
dft$region <- as.numeric(dft$region) 
pcor.test(dft$hetero, dft$mobility, dft$region, method = "pearson") #Смотрим на проверку гипотезы о значимости коэффициента частной корреляции между hetero и mobility, с убранным влиянием region.

dfpart <- na.omit(subset(dflo, region == "S"))
pcor.test(dfpart$hetero, dfpart$mobility, dfpart$moral, method = "pearson") #Аналогично для юго-восточного региона, но убираем влияние moral