Анализ данных в R

https://stepik.org/course/129/syllabus

1 Предобработка данных

setwd(“~/R/data_analysis_1”)

1.3 Работа с data frame

Video 1.3.1

https://stepik.org/lesson/11481/step/3?unit=2523

Reading data

?read.table

?read.csv

df_evals <- read.csv("https://stepic.org/media/attachments/lesson/11481/evals.csv", stringsAsFactors=TRUE)

Summaries

head(df_evals, 5) # show 5 first rows (default=6)
tail(df_evals, 5) # show 5 last rows (default=6)
summary(df_evals)  # summary information
names(df_evals)  # show columns names
str(df_evals)  # show structure: columns information

View(df_evals)  # show data frame in table format
# The same table will appear if Data will be clicked in Environment (right upper window)

Video 1.3.2

https://stepik.org/lesson/11481/step/4?unit=2523

Variables, indexes

df_evals$score  # score column
df_evals$score[1:10]  # show first 10 values of 'score' column
df_evals[5, 1]  # show value of 1st column and 5th row in data frame
df_evals[c(10, 15, 64), 1]  # show several values of 1st column with row number 10, 15, 64
df_evals[5,]  # show 5th row
df_evals[,1]  # show 1st column: equivalent df_evals$score
head(df_evals[,2:5]) # show head of columns from second to fifth

mean(df_evals$score) # mean score
summary(df_evals$score) # score statistical information

Creating and deleting columns

df_evals$ten_point_scale <- df_evals$score * 2  # create new column 'ten_point_scale'
df_evals$number <- NULL  # delete column 'number'

nrow(df_evals)  # rows number
ncol(df_evals)  # columns number

df_evals$number <- 1 : nrow(df_evals)  # create column with row number

Video 1.3.3

https://stepik.org/lesson/11481/step/5?unit=2523

Subsetting with condition

head(df_evals[, c(1, 3, 4)]) # show 1, 3, 4 columns
head(df_evals[, c("gender", "age")])  # show gender and age columns
df_evals[df_evals$gender == 'female',]  # show rows with gender == 'female'
head(df_evals[df_evals$gender == 'female', 1:5])  # show head of rows with gender == 'female' and columns from 1st to 5th
head(subset(df_evals, gender == 'female'))  # show only rows with female gender
subset(df_evals, gender == 'female' & score > 4.9)  # two conditions of gender and score
df_evals[c(1, 4, 5), c(1, 4)]  # combined indexation

rbind and cbind

# rbind
df_evals_female <- subset(df_evals, gender == 'female' & score > 4.9)
df_evals_male <- subset(df_evals, gender == 'male' & score > 4.9)

df_rbind = rbind(df_evals_female, df_evals_male)  # concatenate to df_evals_female and df_evals_male by rows
df_rbind

# cbind
df_evals_1 <- head(df_evals[,1:3])
df_evals_2 <- head(df_evals[,5:7])

df_evals_3 <- cbind(df_evals_1, df_evals_2) # concatenate df_evals_1 and df_evals_2 by columns
df_evals_3

Standard R datasets

Когда вы запускаете R, автоматически загружаются основные библиотеки для работы. В их числе библиотека *datasets**, в которой хранятся наборы данных из разных исследований. Чтобы просмотреть весь список доступных датасетов, выполните команду:

library(help = "datasets")  # Standard R datasets

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

data(mtcars)  # добавит датасет в рабочую среду
help(mtcars)  # выведет информацию о датасете
my_data <- mtcars  # запишет датасет в новую переменную

Task 1.3.1

https://stepik.org/lesson/11481/step/7?unit=2523

В этой задаче поработаем со встроенными данными mtcars. В датафрэйме mtcars создайте новую колонку (переменную) под названием even_gear, в которой будут единицы, если значение переменной (gear) четное, и нули если количество нечетное.

data(mtcars)
head(mtcars)

mtcars$even_gear <- (mtcars$gear %% 2 == 0) * 1
mtcars$even_gear[1:5]

[1] 1 1 1 0 0

Task 1.3.2

https://stepik.org/lesson/11481/step/8?unit=2523

Продолжим нашу работу с данными mtcars. Теперь ваша задача создать переменную - вектор mpg_4 и сохранить в нее значения расхода топлива (mpg) для машин с четырьмя цилиндрами (cyl).

mpg_4 <- mtcars$mpg[mtcars$cyl == 4]  # equivalent: mtcars[mtcars$cyl == 4, "mpg"]
mpg_4

 [1] 22.8 24.4 22.8 32.4 30.4 33.9 21.5 27.3 26.0 30.4 21.4

Task 1.3.3

https://stepik.org/lesson/11481/step/9?unit=2523

А теперь научимся отбирать только некоторые строчки из исходных данных.

Ваша задача создать новый dataframe под названием mini_mtcars, в котором будут сохранены только третья, седьмая, десятая, двенадцатая и последняя строчка датафрейма mtcars.

mini_mtcars <- mtcars[c(3, 7, 10, 12, nrow(mtcars)), ]
mini_mtcars

Task 1.3.4

https://stepik.org/lesson/11481/step/10?unit=2523

Укажите какая команда создаст сабсет данных mtcars, только для тех автомобилей, у которых число цилиндров (cyl) не равняется 3, и время разгона автомобиля (qsec) больше среднего по выборке.

s <- subset(mtcars, cyl != 3 & qsec > mean(qsec))
s

1.4 Элементы синтаксиса

Video 1.4.1

https://stepik.org/lesson/11478/step/2?unit=2520

df_evals <- read.csv('https://stepic.org/media/attachments/lesson/11481/evals.csv')  # load dataframe

if

a <- 0

if (a > 0){
  print('positive')
} else if (a < 0) {
  print('negative')
} else {
  print('zero')
}

[1] "zero"

ifelse

ifelse(a > 0, 'positive', 'not positive')  # 'ifelse' is able to work with arrays, see example below

[1] "not positive"

a <- c(1, -1, 2, 3)
ifelse(a > 0, 'positive', 'not positive')

[1] "positive"     "not positive" "positive"     "positive"    

Video 1.4.2

https://stepik.org/lesson/11478/step/3?unit=2520

for

for (i in 1:5) {
  print(i)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

for in dataset

for (i in 460:nrow(df_evals)){
  print(df_evals$score[i])
}

[1] 3.5
[1] 4.4
[1] 4.4
[1] 4.1

for + if

for (i in 1:nrow(df_evals)){
  if (df_evals$score[i] > 4.9){
    print(df_evals$age[i])
  }
}

[1] 47
[1] 47
[1] 47
[1] 46
[1] 46
[1] 50
[1] 50
[1] 57
[1] 58
[1] 58
[1] 58

for + if vs ifelse

for + if

df_evals$quality <- rep(NA, nrow(df_evals))  # create empty array with length of data set df_evals

for (i in 1:nrow(df_evals)){
  if (df_evals$score[i] > 4){
    df_evals$quality[i] <- 'good'
  } else df_evals$quality[i] <- 'bad'
}

df_evals$quality[1:10]

 [1] "good" "good" "bad"  "good" "good" "good" "bad"  "good" "bad"  "good"

for + ifelse

df_evals$quality <- ifelse(df_evals$score < 4, 'good', 'bad')

while

i <- 1

# first four scores
while (i < 5){
  print(df_evals$score[i])
  i <- i + 1
}

[1] 4.7
[1] 4.1
[1] 3.9
[1] 4.8

Task 1.4.1

https://stepik.org/lesson/11478/step/4?unit=2520

Создайте новую числовую переменную new_var в данных mtcars, которая содержит единицы в строчках, если в машине не меньше четырёх карбюраторов (переменная carb) или больше шести цилиндров (переменная cyl). В строчках, в которых условие не выполняется, должны стоять нули.

df <- mtcars
df$new_var <- ifelse(df$carb >= 4 | df$cyl > 6, 1, 0)
head(df)

Task 1.4.2

https://stepik.org/lesson/11478/step/5?unit=2520

В уже существующей переменной my_vector сохранен вектор из 50 чисел. Если среднее значение вектора my_vector больше 20, в переменную result сохраните “My mean is great”, если среднее значение my_vector меньше или равно 20 то в переменную result сохраните строку “My mean is not so great”.

my_vector = 1:50

if (mean(my_vector) > 20){
  result <- "My mean is great"
} else result <- "My mean is not so great"
print(result)

[1] "My mean is great"

Task 1.4.3

https://stepik.org/lesson/11478/step/6?unit=2520

В этой задаче от вас потребуется узнать некоторую информацию о типах данных в R самостоятельно! Встроенные в R данные AirPassengers - это новый для нас формат данных типа Time-Series. Изучите структуру этих данных, прежде чем начать решение задачи! Например напишите команды:

?AirPassengers      # справка о данных
str(AirPassengers)  # структура данных

В встроенных в R данных AirPassengers хранится 144 значения (количество пассажиров в месяц) с 1949 по 1960 год. Данные Time-Series очень похожи на вектор по своей структуре, например мы можем обратиться к любому из 144 элементов используя уже знакомую нам индексацию AirPassengers[1] или AirPassengers[56].

Можно вообще перевести исходные данные в вектор при помощи команды as.vector(AirPassengers) и продолжить с ними работу как с вектором.

И так ваша задача создать переменную good_months и сохранить в нее число пассажиров только в тех месяцах, в которых это число больше, чем показатель в предыдущем месяце.

Важный момент! В R оператор : для создания последовательности имеет приоритет над арифметическими действиями. Таким образом, если у вас есть переменная i, равная 10, и вы хотите создать вектор от 1 до i - 1, воспользуйтесь скобками, чтобы указать последовательность действий.

good_months <- c()

for (i in 2:length(AirPassengers)){
  if (AirPassengers[i] > AirPassengers[i - 1]){
    good_months <- c(good_months, AirPassengers[i])
  }
}

# one row script 1:
# good_months <- AirPassengers[c(FALSE, AirPassengers[2:144] > AirPassengers[1:143])]

# one row script 2:
# good_months <- AirPassengers[-1][AirPassengers[-1] > AirPassengers[-144]] 

good_months

 [1] 118 132 135 148 118 126 141 149 170 140 145 150 178 172 178 199 166 171 180 193 183 218 230 242 194 196 236
[28] 243 264 272 201 204 235 234 264 302 229 242 267 269 270 315 364 278 284 317 318 374 413 306 315 356 355 422
[55] 465 467 336 340 362 363 435 491 505 337 360 406 420 472 548 559 405 417 419 461 472 535 622 432

Task 1.4.4

https://stepik.org/lesson/11478/step/7?unit=2520

Задачка для супер героев, повышенной сложности!

Для встроенных в R данных AirPassengers рассчитайте скользящее среднее с интервалом сглаживания равным 10. Напечатайте получившийся результат (первым значением в выводе должно быть среднее для элементов 1:10, во втором значении - среднее для элементов 2:11 и т.д., в последнем - среднее для элементов 135:144)

Все полученные значения средних сохраните в переменную moving_average.

moving_average <- numeric(135) # empty array with length 135

for (i in 1:(length(AirPassengers) - 9)){
  moving_average[i] <- mean(AirPassengers[i:(i + 9)])
}

moving_average[1:10]

 [1] 129.8 129.0 129.0 127.3 127.0 129.0 129.0 126.7 126.8 130.2

Можно решить и без цикла при помощи разностей кумулятивных сумм!

n <- 10    
d <- AirPassengers    
cx <- c(0, cumsum(d))    
moving_average <- (cx[(n + 1):length(cx)] - cx[1:(length(cx) - n)]) / n
moving_average[1:10]

 [1] 129.8 129.0 129.0 127.3 127.0 129.0 129.0 126.7 126.8 130.2

1.5 Описательные статистики

Video 1.5.2

https://stepik.org/lesson/11479/step/3?unit=2521

Read and prepare data frame

?mtcars

mtcars_df <- mtcars
str(mtcars_df)  # show data frame structure

# change types of vs and am columns
mtcars_df$vs <- factor(mtcars_df$vs, labels = c("V", "S"))  # "vs" from num type to factor type
mtcars_df$am <- factor(mtcars_df$am, labels = c("Auto", "Manual")) # "am" from num type to factor type
str(mtcars_df)

Video 1.5.3

https://stepik.org/lesson/11479/step/4?unit=2521

median, mean, sd, range

median(mtcars_df$mpg)  # median mpg
mean(mtcars_df$disp)   # mean displacement
sd(mtcars_df$hp)       # standard deviation of horse power
range(mtcars_df$cyl)   # range

# mean mpg of 6 cyl cars with V-engines
mean(mtcars_df$mpg[mtcars_df$cyl == 6 & mtcars_df$vs == "V"])

# standard deviation for cars with automatic transmission and number cylinders is not equal 3
sd(mtcars_df$hp[mtcars_df$cyl != 3 & mtcars_df$am == "Auto"])

Task 1.5.1

https://stepik.org/lesson/11479/step/5?unit=2521

Вновь вернемся к данным mtcars. Рассчитайте среднее значение времени разгона (qsec) для автомобилей, число цилиндров (cyl) у которых не равняется 3 и показатель количества миль на галлон топлива (mpg) больше 20.

Получившийся результат (среднее значение) сохраните в переменную result.

result <- mean(mtcars_df$qsec[mtcars_df$cyl != 3 & mtcars_df$mpg > 20])
result

[1] 18.81643

Video 1.5.4

https://stepik.org/lesson/11479/step/6?unit=2521

aggregate

?aggregate

agg <- aggregate(x = mtcars_df$hp, by = list(mtcars_df$vs), FUN = mean)  # mean horse power for V and S engines
agg

mean_hp_vs <- aggregate(x = mtcars_df$hp, by = list(mtcars_df$vs), FUN = mean)
colnames(mean_hp_vs) <- c("VS", "Mean HP")  # rename columns
mean_hp_vs

# equivalent script (shorter way)
aggregate(hp ~ vs, mtcars_df, mean)

aggregate(hp ~ vs + am, mtcars_df, mean) # with two 'by' categories

# equivalent longer script
# aggregate(x = mtcars_df$hp, by = list(mtcars_df$vs, mtcars_df$am), FUN = mean)

# all median values of num type columns
agg <- aggregate(x = mtcars_df[, -c(8, 9)], by = list(mtcars_df$vs, mtcars_df$am), FUN = median)
agg

# sd values of mpg and disp columns
agg <- aggregate(x = mtcars_df[, c("mpg", "disp")], by = list(mtcars_df$vs, mtcars_df$am), FUN = sd)
agg

# equivalent shorter script
#aggregate(cbind(mpg, disp) ~ am + vs, mtcars_df, sd)

Task 1.5.2

https://stepik.org/lesson/11479/step/7?unit=2521

При помощи функции aggregate рассчитайте стандартное отклонение переменной hp (лошадиные силы) и переменной disp (вместимости двигателя) у машин с автоматической и ручной коробкой передач.

Полученные результаты (результаты выполнения функции aggregate) сохраните в переменную descriptions_stat.

descriptions_stat <- aggregate(cbind(hp, disp) ~ am, mtcars, sd)
descriptions_stat

Video 1.5.5

https://stepik.org/lesson/11479/step/8?unit=2521

describe

?describe

d <- describe(x = mtcars_df[, -c(8, 9)])
d

Video 1.5.6

https://stepik.org/lesson/11479/step/9?unit=2521

describeBy

?describeBy

library(psych)
# 'mat' provides output to one matrix
d <- describeBy(x = mtcars_df[, -c(8, 9)], group = mtcars_df$vs, mat =  T, digits = 2)
d

# 'fast' provides less columns of describe function
d <- describeBy(x = mtcars_df[, -c(8, 9)], group = mtcars_df$vs, mat =  T, digits = 2, fast = T)
d

d <- describeBy(mtcars_df$qsec, group = list(mtcars_df$vs, mtcars_df$am), mat = T, digits = 1, fast = T)
d

Video 1.5.7

https://stepik.org/lesson/11479/step/10?unit=2521

check N/A values

sum(is.na(mtcars_df))  # if 0 it means that all values are available

[1] 0

mean(mtcars_df$mpg, na.rm = T)  # na.rm > na remove: remove all na values

[1] 20.09062

# aggregate function removes na values by default

Task 1.5.3

https://stepik.org/lesson/11479/step/11?unit=2521

Воспользуемся встроенными данными airquality. В новую переменную сохраните subset исходных данных, оставив наблюдения только для месяцев 7, 8 и 9.

При помощи функции aggregate рассчитайте количество непропущенных наблюдений по переменной Ozone в 7, 8 и 9 месяце. Для определения количества наблюдений используйте функцию length().

Результат выполнения функции aggregate сохраните в переменную result.

result <- aggregate(Ozone ~ Month, subset(airquality, Month %in% c(7, 8, 9)), length)
result

Task 1.5.4

https://stepik.org/lesson/11479/step/12?unit=2521

Примените функцию describeBy к количественным переменным данных airquality, группируя наблюдения по переменной Month. Чему равен коэффициент асимметрии (skew) переменной Wind в восьмом месяце?

# str(airquality)
d <- describeBy(airquality$Wind, group = list(airquality$Month), mat = T, digits = 3)
d$skew[d$group1 == 8]

[1] 0.035

# library(psych)
d <- skew(subset(airquality, Month == 8)$Wind)
d

[1] 0.03549581

Task 1.5.5

https://stepik.org/lesson/11479/step/13?unit=2521

Обратимся к встроенным данным iris. Соотнесите значения стандартного отклонения переменных.

res <- describe(iris)
res['sd']

Task 1.5.6

https://stepik.org/lesson/11479/step/14?unit=2521

В данных iris расположите по убыванию значения медиан количественных переменных в группе virginica.

#sort(sapply(iris[iris$Species =='virginica',][-5], FUN=median))
d <- describeBy(iris[,1:4], group = iris$Species)$'virginica'['median']
d

Task 1.5.7

https://stepik.org/lesson/11479/step/15?unit=2521

В переменной my_vector сохранен вектор с пропущенными значениями. Вам нужно создать новый вектор fixed_vector, в котором все пропущенные значения вектора my_vector будут заменены на среднее значение по имеющимся наблюдениям.

При этом исходный вектор оставьте без изменений!

?replace

?aggregate

my_vector <- rnorm(30)
my_vector[sample(1:30, 10)] <- NA  # на десять случайных позиций поместим NA
my_vector_mean = mean(my_vector[!is.na(my_vector)])  # mean value of my_vector elements

fixed_vector = replace(my_vector, is.na(my_vector), my_vector_mean)  # replace na elements with mean value

# another solution
# fixed_vector <- replace(my_vector, is.na(my_vector), mean(my_vector, na.rm = T))

fixed_vector

 [1]  0.162019833  0.162019833  0.162019833 -0.994789932 -0.079637674  0.508381304 -0.393351838  0.162019833
 [9]  0.356677781  0.162019833  0.162019833  0.332249004 -0.633532200 -0.673211856  0.355179160  0.382346243
[17]  2.082183352  1.187997779 -0.358954408  0.162019833 -0.773771710 -0.207107907  0.007899274  0.768779799
[25]  0.162019833  0.671175903  0.893696069  0.162019833  0.162019833 -0.191811483

1.6 Описательные статистики. Графики

Video 1.6.1

https://stepik.org/lesson/11787/step/2?unit=2671

# Read data
df <- mtcars
df$vs <- factor(df$vs, labels = c("V", "S"))
df$am <- factor(df$am, labels = c("Auto", "Manual"))

histogram

?hist

diagram <- hist(df$mpg, breaks = 5, main = "Histogram of MPG", xlab = "MPG")

diagram

$breaks
[1] 10 15 20 25 30 35

$counts
[1]  6 12  8  2  4

$density
[1] 0.0375 0.0750 0.0500 0.0125 0.0250

$mids
[1] 12.5 17.5 22.5 27.5 32.5

$xname
[1] "df$mpg"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

boxplot

?boxplot

diagram <- boxplot(mpg ~ am, df, xlab = "Gearbox type", ylab = "MPG")

diagram

$stats
      [,1] [,2]
[1,] 10.40 15.0
[2,] 14.95 21.0
[3,] 17.30 22.8
[4,] 19.20 30.4
[5,] 24.40 33.9

$n
[1] 19 13

$conf
         [,1]    [,2]
[1,] 15.75947 18.6808
[2,] 18.84053 26.9192

$out
numeric(0)

$group
numeric(0)

$names
[1] "Auto"   "Manual"

plot

diagram <- plot(df$mpg, df$hp)

diagram

NULL

Video 1.6.2

https://stepik.org/lesson/11787/step/3?unit=2671

ggplot2

library(ggplot2) # load library

geom_histogram (ggplot2)

ggplot(df, aes(x = mpg))+
  geom_histogram(fill = "white", col = "black", binwidth = 2)

geom_dotplot (ggplot2)

ggplot(df, aes(x = mpg, fill = am))+
  geom_dotplot(binwidth = 1)

geom_density (ggplot2)

ggplot(df, aes(x = mpg, fill = am))+
  geom_density(alpha = 0.3)

Video 1.6.3

https://stepik.org/lesson/11787/step/4?unit=2671

geom_boxplot

ggplot(df, aes(x = am, y = hp, col = vs))+
  geom_boxplot()

geom_point

ggplot(df, aes(x = mpg, y = hp, col = vs, size = qsec))+
  geom_point()

Task 1.6.1

https://stepik.org/lesson/11787/step/5?unit=2671

При помощи функции ggplot() или boxplot() постройте график boxplot, используя встроенные в R данные airquality. По оси x отложите номер месяца, по оси y — значения переменной Ozone.

На графике boxplot отдельными точками отображаются наблюдения, отклоняющиеся от 1 или 3 квартиля больше чем на полтора межквартильных размаха. Сколько таких наблюдений присутствует в сентябре (месяц №9)?

Обратите внимание, что для корректного отображения графика ggplot ожидает факторную переменную по оси x.

df <- airquality
df$Month <- factor(df$Month, labels = 5 : 9)

# library(ggplot2)
ggplot(df, aes(x = Month, y = Ozone))+
  geom_boxplot(na.rm = T)

Task 1.6.2

https://stepik.org/lesson/11787/step/6?unit=2671

Используем знакомые нам данные mtcars.

Нужно построить scatterplot с помощью ggplot из ggplot2, по оси x которого будет mpg, по оси y - disp, а цветом отобразить переменную (hp).

Полученный график нужно сохранить в переменную plot1.

plot1 <- ggplot(mtcars, aes(x = mpg, y = disp, col = hp))+
  geom_point()
plot1

Task 1.6.3

https://stepik.org/lesson/11787/step/7?unit=2671

Укажите, при помощи какого варианта кода мы можем построить следующий график по данным *iris**. Гистограмма распределения переменной Sepal.Length, в которой цвет заполнения столбцов гистограммы зависит от значения переменной Species.

# ggplot(iris, aes(Sepal.Length)) + geom_histogram(aes(fill = Species))
obj <- ggplot(iris, aes(Sepal.Length, fill = Species)) + geom_histogram()
obj

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Task 1.6.4

https://stepik.org/lesson/11787/step/8?unit=2671

Студент Ярослав очень любит строить графики в R. Основываясь на данных *iris** он хочет построить следующий график.

Scatterplot (диаграмма рассеивания), где по оси X будет отложена переменная Sepal.Length, по оси Y переменная Sepal.Width. За цвет точек будет отвечать переменная Species, а за размер точек переменная Petal.Length.

ggplot(iris, aes(Sepal.Length, Sepal.Width, col = Species, size = Petal.Length)) +
  geom_point()

1.7 Сохранение результатов

Video 1.7.1

https://stepik.org/lesson/11480/step/2?unit=2522

getwd()  # show work directory "C:/Users/xxx/Documents"
setwd("~/R/data_analysis_1")  # set work directory

write.csv(mtcars, "mtcars.csv")  # save data frame as csv file

my_mean = mean(10^5 : 10^7)
save(my_mean, "my_mean.RData")  # save value as file RData