Time
Victoria
24 09 2021
Testing hypotheses about regression
library(readr)
time <- read_csv("time.csv")## Parsed with column specification:
## cols(
## .default = col_double()
## )## See spec(...) for full column specifications.Ваша задача — ответить на вопрос: правда ли, что дополнительная работа плохо сказывается на сне?
Изучите и опишите загруженные данные: Посмотрите, какие переменные содержатся в данных и какой тип у этих переменных. Выведите таблицу с описательными статистиками количественных переменных. Постройте графики, описывающие количественные переменные (гистограмма, график плотности распределения, бокс-плот и т. п.).
library(broom)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.5 ✓ dplyr 1.0.7
## ✓ tibble 3.1.4 ✓ stringr 1.4.0
## ✓ tidyr 1.1.3 ✓ forcats 0.5.0
## ✓ purrr 0.3.4## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(haven)
library(sandwich)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(dplyr)
library(ggplot2)Посмотрите, какие переменные содержатся в данных и какой тип у этих переменных Выведите таблицу с описательными статистиками количественных переменных.
str(time)## spec_tbl_df [539 × 27] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ age : num [1:539] 32 31 44 30 64 41 47 32 30 23 ...
## $ black : num [1:539] 0 0 0 0 0 0 0 0 0 0 ...
## $ clerical: num [1:539] 0 0 0 0 0 0 0 0 0 0 ...
## $ construc: num [1:539] 0 0 0 0 0 0 0 0 0 0 ...
## $ educ : num [1:539] 12 14 17 12 14 12 13 17 15 16 ...
## $ earns74 : num [1:539] 0 9500 42500 42500 2500 ...
## $ gdhlth : num [1:539] 0 1 1 1 1 1 1 1 1 1 ...
## $ inlf : num [1:539] 1 1 1 1 1 1 1 1 1 1 ...
## $ smsa : num [1:539] 0 0 1 0 0 0 1 0 1 1 ...
## $ lothinc : num [1:539] 10.08 0 0 0 9.33 ...
## $ male : num [1:539] 1 1 1 0 1 1 1 1 1 0 ...
## $ marr : num [1:539] 1 0 1 1 1 1 1 1 1 1 ...
## $ prot : num [1:539] 1 1 0 1 1 1 1 0 0 0 ...
## $ rlxall : num [1:539] 3163 2920 3038 3083 3493 ...
## $ selfe : num [1:539] 0 1 1 1 0 0 1 0 1 0 ...
## $ sleep : num [1:539] 3113 2920 2670 3083 3448 ...
## $ slpnaps : num [1:539] 3163 2920 2760 3083 3493 ...
## $ south : num [1:539] 0 1 0 0 0 0 0 0 0 0 ...
## $ spsepay : num [1:539] 0 0 20000 5000 2400 0 0 0 6000 3000 ...
## $ spwrk75 : num [1:539] 0 0 1 1 1 0 0 0 1 1 ...
## $ totwrk : num [1:539] 3438 5020 2815 3786 2580 ...
## $ union : num [1:539] 0 0 0 0 0 0 0 1 0 0 ...
## $ worknrm : num [1:539] 3438 5020 2815 3786 2580 ...
## $ workscnd: num [1:539] 0 0 0 0 0 ...
## $ yngkid : num [1:539] 0 0 0 0 0 0 0 0 0 0 ...
## $ yrsmarr : num [1:539] 13 0 0 12 33 23 24 11 7 4 ...
## $ hrwage : num [1:539] 7.07 1.43 20.53 9.62 2.75 ...
## - attr(*, "spec")=
## .. cols(
## .. age = col_double(),
## .. black = col_double(),
## .. clerical = col_double(),
## .. construc = col_double(),
## .. educ = col_double(),
## .. earns74 = col_double(),
## .. gdhlth = col_double(),
## .. inlf = col_double(),
## .. smsa = col_double(),
## .. lothinc = col_double(),
## .. male = col_double(),
## .. marr = col_double(),
## .. prot = col_double(),
## .. rlxall = col_double(),
## .. selfe = col_double(),
## .. sleep = col_double(),
## .. slpnaps = col_double(),
## .. south = col_double(),
## .. spsepay = col_double(),
## .. spwrk75 = col_double(),
## .. totwrk = col_double(),
## .. union = col_double(),
## .. worknrm = col_double(),
## .. workscnd = col_double(),
## .. yngkid = col_double(),
## .. yrsmarr = col_double(),
## .. hrwage = col_double()
## .. )summary(time)## age black clerical construc
## Min. :23.00 Min. :0.00000 Min. :0.0000 Min. :0.00000
## 1st Qu.:29.00 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :37.00 Median :0.00000 Median :0.0000 Median :0.00000
## Mean :39.33 Mean :0.04824 Mean :0.1466 Mean :0.02412
## 3rd Qu.:49.00 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.00000
## Max. :65.00 Max. :1.00000 Max. :1.0000 Max. :1.00000
##
## educ earns74 gdhlth inlf
## Min. : 1.00 Min. : 0 Min. :0.0000 Min. :0.0000
## 1st Qu.:12.00 1st Qu.: 2500 1st Qu.:1.0000 1st Qu.:1.0000
## Median :12.00 Median : 8250 Median :1.0000 Median :1.0000
## Mean :12.76 Mean : 9728 Mean :0.8905 Mean :0.7514
## 3rd Qu.:15.00 3rd Qu.:13750 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :17.00 Max. :42500 Max. :1.0000 Max. :1.0000
##
## smsa lothinc male marr
## Min. :0.0000 Min. : 0.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:1.0000
## Median :0.0000 Median : 8.517 Median :1.0000 Median :1.0000
## Mean :0.3989 Mean : 6.251 Mean :0.5733 Mean :0.8126
## 3rd Qu.:1.0000 3rd Qu.: 9.328 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :10.657 Max. :1.0000 Max. :1.0000
##
## prot rlxall selfe sleep slpnaps
## Min. :0.0000 Min. :1380 Min. :0.0000 Min. : 755 Min. :1335
## 1st Qu.:0.0000 1st Qu.:3159 1st Qu.:0.0000 1st Qu.:3029 1st Qu.:3106
## Median :1.0000 Median :3429 Median :0.0000 Median :3278 Median :3375
## Mean :0.6716 Mean :3441 Mean :0.1262 Mean :3266 Mean :3386
## 3rd Qu.:1.0000 3rd Qu.:3698 3rd Qu.:0.0000 3rd Qu.:3532 3rd Qu.:3650
## Max. :1.0000 Max. :6110 Max. :1.0000 Max. :4695 Max. :6110
##
## south spsepay spwrk75 totwrk
## Min. :0.0000 Min. : 0 Min. :0.000 Min. : 0
## 1st Qu.:0.0000 1st Qu.: 0 1st Qu.:0.000 1st Qu.:1594
## Median :0.0000 Median : 0 Median :0.000 Median :2288
## Mean :0.1837 Mean : 5051 Mean :0.475 Mean :2140
## 3rd Qu.:0.0000 3rd Qu.: 8400 3rd Qu.:1.000 3rd Qu.:2700
## Max. :1.0000 Max. :75000 Max. :1.000 Max. :5020
##
## union worknrm workscnd yngkid
## Min. :0.0000 Min. : 0 Min. : 0.00 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:1574 1st Qu.: 0.00 1st Qu.:0.0000
## Median :0.0000 Median :2275 Median : 0.00 Median :0.0000
## Mean :0.2189 Mean :2110 Mean : 30.17 Mean :0.1206
## 3rd Qu.:0.0000 3rd Qu.:2657 3rd Qu.: 0.00 3rd Qu.:0.0000
## Max. :1.0000 Max. :5020 Max. :1337.00 Max. :1.0000
##
## yrsmarr hrwage
## Min. : 0.00 Min. : 0.510
## 1st Qu.: 0.00 1st Qu.: 2.820
## Median : 8.00 Median : 4.240
## Mean :11.88 Mean : 5.045
## 3rd Qu.:21.00 3rd Qu.: 6.210
## Max. :43.00 Max. :35.510
## NA's :134Постройте графики, описывающие количественные переменные (гистограмма, график плотности распределения, бокс-плот и т. п.).
ggplot(data=time,aes(age))+geom_histogram(aes(y = ..density..)) +
geom_density(alpha = 0.1, fill = "red")+ggtitle("Возраст участников")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(time, aes(sleep , after_stat(density),fill= male)) +
geom_freqpoly(binwidth = 500)+ggtitle("Время сна")
time$black=as.factor(time$black)
ggplot(data=time, aes(x=black,age))+geom_boxplot()+ggtitle("Соотношение возраста и цвета кожи участников по полу")+ facet_grid(. ~ male)
ggplot(data = time, aes(x =educ, y= hrwage)) + geom_col( )+ggtitle("Соотношение лет образования и заработной платы")+ylab("Заработная плата")## Warning: Removed 134 rows containing missing values (position_stack).
- Проведите предварительную обработку данных: Проверьте, есть ли в наборе переменные с пропущенными значениями. Если есть, то приведет ли это к каким-то последствиям при последующем анализе?
Проверьте, есть в наборе выбросы. Если есть, то приведет ли это к каким-то последствиям при последующем анализе? Покажите связь между сном и работой на диаграмме рассеяния с линией тренда. Рассчитайте оценку коэффициента корреляции и сделайте вывод о его значимости.
sleep сколько минут тратит на сон в неделю slpnaps то же, но включая дневной короткий сон totwrk сколько минут проводит за работой в неделю worknrm сколько минут тратит на работу по основному месту работы workscnd сколько минут тратит на работу по не основному месту
sum(is.na(time))## [1] 134colSums(is.na(time)) #hrwage почасовая зарплата 134 пропущенных и все там- не будем рассматривать ее в модели## age black clerical construc educ earns74 gdhlth inlf
## 0 0 0 0 0 0 0 0
## smsa lothinc male marr prot rlxall selfe sleep
## 0 0 0 0 0 0 0 0
## slpnaps south spsepay spwrk75 totwrk union worknrm workscnd
## 0 0 0 0 0 0 0 0
## yngkid yrsmarr hrwage
## 0 0 134#выбросы - так как мы анализируем как дополнительная работа влияет на работоспособность , то для анализа выбросов возьмем две переменные - sleep- сон в неделю & workscnd сколько минут тратит на работу по не основному месту - оба из этих данных continious
ggplot(time, aes(x=sleep, y=workscnd)) +
geom_point()
#выбросы как видно есть , кто-то не имеет дополнительной работы, поэтому их показатели на нуле, а у кого-то она есть.Выбросы- те кто работает на дополнительное работе.
#Чтобы проанализировать , как доп. работа влияест на сон, нам нужны те кто не работает на доп работе и кто работает те их workscnd>0
newdata=time %>% filter(time$workscnd!=0) #
times=time %>% mutate(haveextrawork=case_when(time$workscnd!=0 ~ "Yes", TRUE ~ "No")) #создадим новую переменную- разделим на тех , у кого есть доп. работа и у кого ее нет, и посмотрим на данные
ggplot(data=times)+geom_boxplot(aes(y=sleep,x=haveextrawork))
ggplot(newdata, aes(y=sleep, x=workscnd)) +
geom_col()
scatter.smooth(x=newdata$workscnd, y=newdata$sleep, main="totwork~sleep") 
Покажите связь между сном и работой на диаграмме рассеяния с линией тренда. Рассчитайте оценку коэффициента корреляции и сделайте вывод о его значимости.
ggplot(data=time, aes(x=totwrk, y=sleep)) + geom_point(aes(color="red")) + geom_smooth() + labs(title="Scatter",
x = "Xaxis", y = "Y-axis", color="Color") # более точное приближение## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
scatter.smooth(x=time$totwrk, y=time$sleep, main="totwork~sleep")
cor(time$sleep, time$workscnd) #в данном случае у нас почти нудевая корреляция , что означает , что дополнительная работа почти никак не влияет на сон## [1] 0.007587354cor(time$sleep, time$totwrk ) # тогда как при оценке часов на основной работе мы видим отрицательную корреляцию, о есть при увеличении рабочего часа на основной работе на час , то сон в среднем уменьшается на -0.3182219.## [1] -0.3182219- Проведите регрессионный анализ взаимосвязи между переменными сна и работы. Оцените модель парной регрессии. Проинтерпретируйте коэффициент при регрессоре. Что показывают p-value и R2? Какой вывод о качестве регрессии можно сделать?
Linear regression
scatter.smooth(x=time$totwrk, y=time$sleep, main="totwork~sleep") Тогда наша формула линейной регрессии Y=-0.1538 *Sleep + 3595.6109
linearMod <- lm(sleep~totwrk, data=time)
summary(linearMod)##
## Call:
## lm(formula = sleep ~ totwrk, data = time)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2430.95 -228.28 4.21 262.42 1335.22
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3595.61094 45.94460 78.260 < 2e-16 ***
## totwrk -0.15383 0.01978 -7.779 3.78e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 415.4 on 537 degrees of freedom
## Multiple R-squared: 0.1013, Adjusted R-squared: 0.09959
## F-statistic: 60.51 on 1 and 537 DF, p-value: 3.778e-14Модель показывает пересечение и наклон линии регрессии
Что показывают p-value и R2? Какой вывод о качестве регрессии можно сделать?
Так как уровень p-value очень низкий < 2e-16 , что меньше 5 % уровня значимости, что влечет за собой отклонение нулевой гипотезы. Adjusted R-squared: 0.09959 - что говорит нам о том, что данная модель плохо предсказывает нашу объясняемую переменную
h0=0 ho!=h1 Отвергаемм нулеву гипотезу A small p-value <= alpha, which is usually 0.05, indicates that the observed data is sufficiently inconsistent with the null hypothesis, so the null hypothesis may be rejected. The alternate hypothesis is true at the 95% confidence interval
cor.test(time$totwrk, time$sleep,
alternative = c("two.sided"),
method = c("pearson"),
exact = NULL, conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: time$totwrk and time$sleep
## t = -7.7786, df = 537, p-value = 3.778e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3921388 -0.2402221
## sample estimates:
## cor
## -0.3182219