Лабораторная работа R
Bredikhina
2022-11-18
Раздел 1. Статистическое оценивание
1.1. Импорт данных для работы
Сначала создадим папку R, куда поместим файл “Данные по вариантам для выполнения типового расчета.xlsx”. При импорте данных из файла надо ответить на два вопроса. Где расположен файл для импорта? Какую команду использовать для импорта? Чтобы ответить на первый вопрос, необходимо определить текущую рабочую директорию (working directory), то есть то место в памяти компьютера, где сохраняются созданные в текущей сессии объекты (векторы, функции, таблицы данных и т.д.). Для этого я воспользовалась командой getwd().
getwd()## [1] "C:/Users/bvict/OneDrive/Рабочий стол/R"setwd('C:/Users/bvict/OneDrive/Рабочий стол/R')Для того, чтобы импортировать данные из Excel, нам необходимо установить программу rio. Устанавливаем пакет через меню RStudio Tools → Install Packages…→ указать имя пакета rio. Запомним, что при каждом сеансе в Rstudio нам необходимо будет заново подключать пакеты
Импортируем файл функцией import пакета rio, аргументом необходимо указать название файла в кавычках.Присваиваем импортированным данным имя dataframe:
library('rio')
dataframe <- import('Данные по вариантам для выполнения типового расчета(3).xlsx')Сохраним теперь отдельным объектом данные по Коэффициенту соотношения заёмных и собственных средств 100 предприятий, которые сохранены в нашем файле в столбце под именем “v5”:
data <- dataframe['v5']Посмотреть на данные целиком можно, вызвав их по имени data. Мы здесь ограничимся лишь просмотром начала таблицы (датафрейма), используя функцию head:
head(data) # показывает первые шесть (по умолчанию) строк каждого столбца…и конца таблицы, используя функцию tail:
tail(data)#показывает последние шесть (по умолчанию) строк каждого столбцаБолее информативен для начального знакомства с данными результат вызова другой функции str, которая позволяет узнать класс объекта, в котором содержатся данные (в нашем случае data.frame), число наблюдений (observations) и столбцов-признаков (variable), а также тип данных каждого признака (указан после имени) и первые четыре наблюдения:
str(data)## 'data.frame': 101 obs. of 1 variable:
## $ v5: chr "5.9" "5.9" "5.9" "8.1" ...Данные успешно импортированы. Теперь необходимо привести их к виду, при котором мы сможем использовать функции статистического анализа R.
1.2. Подготовка данных для статистического анализа
Можно заметить, что к нашим 100 наблюдениям добавилось одно лишнее, которое было внизу столбца, и тип данных не числовой (numeric), а текстовый (character). Исправить все можно двумя строчками кода:
data <- data[1:100,] # отрезает в data первые 100 наблюдений
data <- as.data.frame(as.numeric(data)) # изменяет формат данных на числовойПереименуем столбец данных в соответствии с исследуемой величиной:
names(data) <- 'Коэффициент соотношения заёмных и собственных средств 100 предприятий'Проверим структуру исправленных данных, теперь все верно:
str(data)## 'data.frame': 100 obs. of 1 variable:
## $ Коэффициент соотношения заёмных и собственных средств 100 предприятий: num 5.9 5.9 5.9 8.1 3.9 4.1 4.1 4.1 6.1 11.6 ...Наконец можем перейти к статистическому анализу.
1.3. Дескриптивная статистика
Самое начальное представление об исследуемой величине по имеющейся выборке можно получить, вызвав команду summary базового пакета:
summary(data)## Коэффициент соотношения заёмных и собственных средств 100 предприятий
## Min. : 0.300
## 1st Qu.: 3.900
## Median : 4.800
## Mean : 5.272
## 3rd Qu.: 6.325
## Max. :11.600В результате мы получили минимальное значение в выборке (min), выборочные нижний квартиль (1st qu), медиану (median), среднюю (mean), верхний квартиль (3rd qu), максимальное значение в выборке (max).
Более подробно о распределении изучаемого признака можно узнать, используя функцию Desc (Describe data) пакета DescTools (Tools for Descriptive Statistics):
library('DescTools')
Desc(data)## ------------------------------------------------------------------------------
## Describe data (data.frame):
##
## data frame: 100 obs. of 1 variables
## 100 complete cases (100.0%)
##
## Nr ColName
## 1 Коэффициент соотношения заёмных и собственных средств 100 предприятий
## Class NAs Levels
## numeric . ...
##
##
## ------------------------------------------------------------------------------
## 1 - Коэффициент соотношения заёмных и собственных средств 100 предприятий (numeric)
##
## length n NAs unique 0s mean meanCI'
## 100 100 0 34 0 5.272 4.735
## 100.0% 0.0% 0.0% 5.809
##
## .05 .10 .25 median .75 .90 .95
## 1.200 2.250 3.900 4.800 6.325 9.200 10.600
##
## range sd vcoef mad IQR skew kurt
## 11.300 2.704 0.513 1.631 2.425 0.521 -0.237
##
## lowest : 0.3 (2), 0.84, 0.87, 1.2 (5), 1.8
## highest: 7.9 (6), 8.1 (2), 9.2 (6), 10.6 (5), 11.6 (3)
##
## heap(?): remarkable frequency (13.0%) for the mode(s) (= 4.1)
##
## ' 95%-CI (classic)
Рис. 1. Гистограмма, ящичковая диаграмма и кумулята исследуемой переменной, выдаваемые функцией Desc пакета DescTools
В результате мы получаем таблицу с кратким описанием признаков, имеющихся в датафрейме (у нас один такой признак):
его тип (class), число пропущенных значений (NAs от Not Available),
число категорий,
если признак категориальный (levels).
Также мы получаем расширенную дескриптивную статистику для каждого признака в датафрейме:
объем выборки (length), число и процент наблюденных значений признака, т.е. без учета пропущенных значений (n),
число и процент пропущенных значений признака (NAs),
число уникальных неповторяющихся значений признака (unique),
число и процент нулевых значений признака (0s),
выборочную среднюю,
доверительный интервал для средней, используя оценку среднеквадратического отклонения (meanCI от Confidence Interval),
выборочные процентили, в том числе медиану, размах вариации (range), выборочное среднеквадратическое отклонение (sd от Standard Deviation),
коэффициент вариации (vcoef), медиану абсолютных отклонений значений признака от медианы (mad от Median Absolute Deviation): 𝑚𝑎𝑑(𝑥) = 𝑀𝑒(|𝑥𝑖 − 𝑀𝑒(𝑥)|)
интерквартильный размах (IQR),
выборочные коэффициент асимметрии (skew от Skewness) и эксцесса (kurt от Kurtosis),
наименьшие (lowest) и наибольшие значения признака в выборке (highest) - в скобках указана частота значения, если она больше единицы.
Для удобства дальнейшего анализа запишем данные в вектор под именем volume, а также создадим переменную с объёмом выборки n:
volume <- data[[1]]
n <- length(volume)1.4. Гистограмма эмпирических частот
Рассмотрим, как строить гистограмму частот и теоретическую кривую распределения самостоятельно контролируя необходимые аспекты. Для этого нам будет достаточно базовых функций R.
Построим сначала гистограмму частот. Для этого воспользуемся функцией hist, указав значения следующих аргументов:
x - значения признака (обязательный аргумент, все прочие опциональные),
freq (от frequency) - уточняет построение гистограммы частот или гистограммы плотности: в первом случае необходимо указатьTRUE, во втором - FALSE,
col - цвет гистограммы: например, можно указать blue, или red, или green и т.д.,
breaks - регулирует разбиение наблюдаемой по выборке области значений признака на интервалы: необходимо передать либо желаемое число интервалов, либо вектор из нижних границ интервалов, либо функцию, вычисляющую число интервалов, либо функцию, вычисляющую границы интервалов, либо знакомый вам способ определения ширины интервала – по названию методов - ‘sturges’, ‘scott’, ‘freedman-diaconis’ (отметим, что шаг при этом округляется) (по умолчанию используется формула Стерджеса определения ширины интервала),
xlab - подпись оси абсцисс: необходимо передать сроку текста,
ylab - подпись оси ординат: необходимо передать сроку текста,
ylim - регулирует отрезок оси ординат, который мы видим: необходимо передать вектор с компонентами начала отрезка и его конца,
main - название графика: необходимо передать сроку текста.
hist = hist(volume, freq = TRUE, col = 'blue',
breaks = 'sturges',
xlab = 'Нижняя граница интервала',
ylab = 'Частота',
ylim = c(0, 30),
main = 'Гистограмма эмпирических частот')
Рис. 2. Гистограмма эмпирических частота коэффициента соотношения заёмных средств, построенная с помощью базовой функции hist.
Вектор из абсцисс мы сгенерируем функцией seq, первым аргументом которой необходимо передать значение левого конца отрезка, вторым — правого конца (определим их через минимум и максимум значений по выборке, немного от них отойдя), а последним аргументом — желаемое число точек (чем больше точек мы сгенерируем, тем более гладкой получится соединяющая их линия, наша теоретическая кривая):
x_values <- seq(from = min(volume)-0.1,to = max(volume)+0.1,length = 1000)Значения теоретической кривой нормального распределения определим по формуле:
𝑚𝑖Т = 𝑃(𝑎𝑖≤𝑋≤𝑏𝑖) ∙ 𝑛 = (𝐹(𝑏𝑖, 𝜇^, 𝜎^) − 𝐹(𝑎𝑖, 𝜇^, 𝜎^))∙𝑛 ≈ 𝑓(𝑥𝑖,𝜇^,𝜎^)∙ ℎ ∙ 𝑛
где 𝜇^,𝜎^ — оценки параметров, ℎ — ширина интервала гистограммы, 𝑓(⋅) — функция плотности нормального распределения, 𝐹(⋅) – функция распределения.
Ширину интервалов определим из полученного нами ранее результата вызова функции hist (затем мы и сохранили результат в объекте hist). Объект hist представляет из себя список, содержащий всю основную информацию о построенной гистограмме частот, первый элемент списка — список breaks, который содержит границы интервалов. Шаг мы определим как разницу между вторым и первым элементами этого списка (т.е. верхней и нижней границы первого интервала) (заметим, что в случае иного определения границ интервалов формулу нахождения теоретической кривой, возможно, придется уточнить):
h = unlist(hist['breaks'])[2]-unlist(hist['breaks'])[1]Вектор из значений функции плотности нормального распределения в сгенерированных точках получим, воспользовавшись функцией dnorm, аргументами которой являются вектор из координат точек, значение функции плотности в которых нужно определить, средняя и стандартное отклонение совокупности (мы используем их выборочные оценки) (по умолчанию предполагается, что функция соответствует стандартному нормальному закону):
y_values <- dnorm(x_values, mean = mean(volume), sd = sd(volume))*h*nМожно найти теоретические вероятности через функцию распределения нормального закона (используя функцию в R).
Итак, нанесем теоретическую кривую нормального закона на гистограмму частот (рис. 3) (можно настроить цвет линии аргументом col):
hist = hist(volume, freq = TRUE, col = 'blue',
breaks = 'sturges',
xlab = 'Нижняя граница интервала',
ylab = 'Частота',
ylim = c(0, 30),
main = 'Гистограмма частот и теоретическая кривая Гаусса')
lines(x_values, y_values, col = 'red')
Рис. 3. Гистограмма эмпирических частот коэффициента соотношения заёмных средств и теоретическая кривая нормального распределения, построенные с помощью базовой функции hist
1.5. Интервальные оценки параметров нормального распределения
Для построения доверительных интервалов используем функции пакета DescTools.
Зададим сперва надежность, с которой будем строить все интервальные оценки параметров:
gamma = 0.99Доверительный интервал для средней вычисляется функцией MeanCI (CI от Confidence Interval), аргументами которой необходимо передать вектор значений наблюдений выборки, по которой строится оценка, уровень надежности (conf.level, по умолчанию 0,95), если известно генеральное среднеквадратическое отклонение, то указать его значение, в противном случае, по умолчанию используется его выборочная оценка.
Предположим сначала, что нам известно значение генерального среднеквадратического отклонения 𝜎(𝑋) = 0,2, тогда доверительный интервал для средней строится следующим образом:
MeanCI(volume, sd = 0.2, conf.level = gamma)## mean lwr.ci upr.ci
## 5.272100 5.220583 5.323617Первое значение результата вызова функции соответствует точечной оценке средней (mean) – среднему арифметическому 𝑥 = 5.2721, второе - нижней границе доверительного интервала (lwr.ci от lower bound of the confidence interval), третье - верхней границе интервала (upr.ci от upper bound of the confidence interval).
Таким образом, при известной генеральной дисперсии интервальная оценка генеральной средней: 𝑃(5.221 ≤ 𝜇 ≤ 5.324) = 0.99
В случае, если генеральная дисперсия неизвестна, доверительный интервал определяем следующим образом:
MeanCI(volume, conf.level = gamma)## mean lwr.ci upr.ci
## 5.272100 4.561793 5.982407Доверительный интервал для дисперсии можно найти, воспользовавшись функцией VarCI, указав дополнительным аргументом метод norm:
VarCI(volume, conf.level = gamma)## var lwr.ci upr.ci
## 7.314229 5.209910 10.887198Получаем точечную оценку дисперсии 𝑆̂ 2 = 7.3142, нижнюю и верхнюю границу доверительного интервала для дисперсии с надежностью 0,99. 𝑃(5.2100 ≤ 𝜎2 ≤ 10.8872) = 0.99
Построение асимптотического доверительного интервала для среднеквадратического отклонения не найдено нами в пакетах R, поэтому определим границы доверительного интервала, написав код для их определения самостоятельно (в чем тоже преимущество использования языка R).
Квантиль нормального распределения определим функцией qnorm, аргументом которой передадим 1 − 𝜀/2, где 𝜀 = 1 − 𝛾.
sdCI <- c(sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) + qnorm(1-(1-gamma)/2)),
sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) - qnorm(1-(1-gamma)/2)))
sdCI <- c(sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) + qnorm(1-(1-gamma)/2)),
sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) - qnorm(1-(1-gamma)/2)))
names(sdCI) <- c('sd_min', 'sd_max')
sdCI## sd_min sd_max
## 2.302451 3.337496Таким образом, получаем асимптотический доверительный интервал для генерального среднего квадратического отклонения: 𝑃(2.4663 ≤ 𝜎 ≤ 3.0083) = 0.99
1.6. Интервальные оценки вероятности (генеральной доли) биномиального распределения
Задача из ДКР1: найти с надежностью 0,978 границы, в которых будет лежать доля жителей города,добирающихся на работу общественным транспортом в случае, если из 340 случайноотобранных респондентов 55% добираются на работу общественным транспортом.
Рассмотрим построение доверительного интервала для вероятности.
Поскольку речь идет о выборке большого объема, то мы используем нормальную аппроксимацию (интервал Вальда):
n = 340 # количество случайно отобранных респондентов
m = 187 # количество респондентов, добирающихся на общественном транспорте (55%)
gamma = 0.978 # надежность
BinomCI(m, n, conf.level=gamma, method = 'wald')## est lwr.ci upr.ci
## [1,] 0.55 0.488205 0.611795Итак, получили доверительный интервал, построенный с надежностью 0,978:
𝑃(0,4882 < 𝑝 < 0,6118) = 0,978
Также при большом числе испытаний можно рассмотреть интервал Агрести-Коула (его авторы 3-го источника считают лучшим в случае больших выборок), и в нашем случае он дал интервал меньшей ширины:
BinomCI(m, n, conf.level=gamma, method = 'agresti-coull')## est lwr.ci upr.ci
## [1,] 0.5492403 0.4879072 0.6105734