Data 606 Lab 7
Wilson Chau
2022-10-23
library(tidyverse)
library(openintro)
library(infer)
(require(stats))## [1] TRUElibrary(stats)The data
Every two years, the Centers for Disease Control and Prevention conduct the Youth Risk Behavior Surveillance System (YRBSS) survey, where it takes data from high schoolers (9th through 12th grade), to analyze health patterns. You will work with a selected group of variables from a random sample of observations during one of the years the YRBSS was conducted.
Load the yrbss data set into your workspace.
data('yrbss', package='openintro')There are observations on 13 different variables, some categorical and some numerical. The meaning of each variable can be found by bringing up the help file:
?yrbss#1. What are the cases in this data set? How many cases are there in our sample?
There are 13,583 cases on this sample. There are 13 variables that define these samples.
glimpse(yrbss)## Rows: 13,583
## Columns: 13
## $ age <int> 14, 14, 15, 15, 15, 15, 15, 14, 15, 15, 15, 1…
## $ gender <chr> "female", "female", "female", "female", "fema…
## $ grade <chr> "9", "9", "9", "9", "9", "9", "9", "9", "9", …
## $ hispanic <chr> "not", "not", "hispanic", "not", "not", "not"…
## $ race <chr> "Black or African American", "Black or Africa…
## $ height <dbl> NA, NA, 1.73, 1.60, 1.50, 1.57, 1.65, 1.88, 1…
## $ weight <dbl> NA, NA, 84.37, 55.79, 46.72, 67.13, 131.54, 7…
## $ helmet_12m <chr> "never", "never", "never", "never", "did not …
## $ text_while_driving_30d <chr> "0", NA, "30", "0", "did not drive", "did not…
## $ physically_active_7d <int> 4, 2, 7, 0, 2, 1, 4, 4, 5, 0, 0, 0, 4, 7, 7, …
## $ hours_tv_per_school_day <chr> "5+", "5+", "5+", "2", "3", "5+", "5+", "5+",…
## $ strength_training_7d <int> 0, 0, 0, 0, 1, 0, 2, 0, 3, 0, 3, 0, 0, 7, 7, …
## $ school_night_hours_sleep <chr> "8", "6", "<5", "6", "9", "8", "9", "6", "<5"…Exploratory data analysis
summary(yrbss$weight)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 29.94 56.25 64.41 67.91 76.20 180.99 1004#2. How many observations are we missing weights from?
sum(is.na(yrbss))## [1] 9476summary(yrbss$weight)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 29.94 56.25 64.41 67.91 76.20 180.99 1004There is 9,476 observations that are missing observations. Under this observation there is 1004 missing weights.
yrbss <- yrbss %>%
mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no"))#3. Make a side-by-side boxplot of physical_3plus and
weight. Is there a relationship between these two
variables? What did you expect and why?
yrbss2 <- yrbss %>%
mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no")) %>%
na.exclude()
ggplot(yrbss2, aes(x=weight, y=physical_3plus)) + geom_boxplot() + theme_bw()
I would say that this relationship is about the physical activity. The
interesting thing I found is that students who said yes to working out,
but the weight gain is similar to the students who said no to working
out.
yrbss %>%
group_by(physical_3plus) %>%
summarise(mean_weight = mean(weight, na.rm = TRUE))## # A tibble: 3 × 2
## physical_3plus mean_weight
## <chr> <dbl>
## 1 no 66.7
## 2 yes 68.4
## 3 <NA> 69.9There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.
Inference
4. Are all conditions necessary for inference satisfied? Comment on
each. You can compute the group sizes with the summarize
command above by defining a new variable with the definition
n().
yrbss %>%
group_by(physical_3plus) %>%
summarise(freq = table(weight)) %>%
summarise(n = sum(freq))## # A tibble: 3 × 2
## physical_3plus n
## <chr> <int>
## 1 no 4022
## 2 yes 8342
## 3 <NA> 215Using the summarize command. I went to see that the respond for yes is still more than double the answer of no.
5. Write the hypotheses for testing if the average weights are different for those who exercise at least times a week and those who don’t.
Null Hypothesis H0 => True unless prove it’s wrong.
H0 = Students who are physically active 3 or more days per week have similar weight as those who are not physically active 3 or more days per week.
Alternative Hypothesis HA => Accept HA, reject H0.
Students who are are physically active 3 or more dyas per week have
Next, we will introduce a new function, hypothesize, that
falls into the infer workflow. You will use this method for
conducting hypothesis tests.
But first, we need to initialize the test, which we will save as
obs_diff.
Notice how you can use the functions specify and
calculate again like you did for calculating confidence
intervals. Here, though, the statistic you are searching for is the
difference in means, with the order being
yes - no != 0.
After you have initialized the test, you need to simulate the test on
the null distribution, which we will save as null.
Here, hypothesize is used to set the null hypothesis as
a test for independence. In one sample cases, the null
argument can be set to “point” to test a hypothesis relative to a point
estimate.
Also, note that the type argument within
generate is set to permute, whichis the
argument when generating a null distribution for a hypothesis test.
We can visualize this null distribution with the following code:
yrbss %>%
group_by(physical_3plus) %>%
summarise(freq = table(weight)) %>%
summarise(n = sum(freq))## # A tibble: 3 × 2
## physical_3plus n
## <chr> <int>
## 1 no 4022
## 2 yes 8342
## 3 <NA> 215devtools::install_github("statswithr/statsr",
dependencies=TRUE,
upgrade_dependencies = TRUE)## pillar (1.8.0 -> 1.8.1 ) [CRAN]
## crayon (1.5.1 -> 1.5.2 ) [CRAN]
## ps (1.7.1 -> 1.7.2 ) [CRAN]
## processx (3.7.0 -> 3.8.0 ) [CRAN]
## digest (0.6.29 -> 0.6.30 ) [CRAN]
## htmltools (0.5.2 -> 0.5.3 ) [CRAN]
## viridisLite (0.4.0 -> 0.4.1 ) [CRAN]
## scales (1.2.0 -> 1.2.1 ) [CRAN]
## isoband (0.2.5 -> 0.2.6 ) [CRAN]
## gtable (0.3.0 -> 0.3.1 ) [CRAN]
## cpp11 (0.4.2 -> 0.4.3 ) [CRAN]
## purrr (0.3.4 -> 0.3.5 ) [CRAN]
## dplyr (1.0.9 -> 1.0.10 ) [CRAN]
## xfun (0.31 -> 0.35 ) [CRAN]
## yaml (2.3.5 -> 2.3.6 ) [CRAN]
## stringr (1.4.0 -> 1.5.0 ) [CRAN]
## evaluate (0.15 -> 0.18 ) [CRAN]
## sass (0.4.1 -> 0.4.4 ) [CRAN]
## jsonlite (1.8.0 -> 1.8.3 ) [CRAN]
## deSolve (NA -> 1.34 ) [CRAN]
## contfrac (NA -> 1.1-12 ) [CRAN]
## elliptic (NA -> 1.4-0 ) [CRAN]
## callr (3.7.2 -> 3.7.3 ) [CRAN]
## knitr (1.39 -> 1.41 ) [CRAN]
## hunspell (NA -> 3.0.2 ) [CRAN]
## commonmark (1.8.0 -> 1.8.1 ) [CRAN]
## bslib (0.3.1 -> 0.4.1 ) [CRAN]
## fontawesome (0.3.0 -> 0.4.0 ) [CRAN]
## httpuv (1.6.5 -> 1.6.6 ) [CRAN]
## tidyr (1.2.0 -> 1.2.1 ) [CRAN]
## tinytex (0.40 -> 0.42 ) [CRAN]
## hypergeo (NA -> 1.2-13 ) [CRAN]
## mvtnorm (NA -> 1.1-3 ) [CRAN]
## pbapply (NA -> 1.6-0 ) [CRAN]
## coda (NA -> 0.19-4 ) [CRAN]
## HistData (NA -> 0.8-7 ) [CRAN]
## spelling (NA -> 2.2 ) [CRAN]
## cubature (NA -> 2.0.4.5 ) [CRAN]
## shiny (1.7.2 -> 1.7.3 ) [CRAN]
## rmarkdown (2.14 -> 2.18 ) [CRAN]
## BayesFactor (NA -> 0.9.12-4.4) [CRAN]
##
## The downloaded binary packages are in
## /var/folders/1m/84tfbh350m597c6cnbw9pnsh0000gn/T//RtmpuVP9gK/downloaded_packages
## ── R CMD build ─────────────────────────────────────────────────────────────────
##
checking for file ‘/private/var/folders/1m/84tfbh350m597c6cnbw9pnsh0000gn/T/RtmpuVP9gK/remotes1f0d507f33d8/StatsWithR-statsr-9cb9eda/DESCRIPTION’ ...
✔ checking for file ‘/private/var/folders/1m/84tfbh350m597c6cnbw9pnsh0000gn/T/RtmpuVP9gK/remotes1f0d507f33d8/StatsWithR-statsr-9cb9eda/DESCRIPTION’
##
─ preparing ‘statsr’:
## checking DESCRIPTION meta-information ...
✔ checking DESCRIPTION meta-information
##
─ checking for LF line-endings in source and make files and shell scripts
##
─ checking for empty or unneeded directories
##
─ building ‘statsr_0.3.0.tar.gz’
##
## #6. How many of these null permutations have a
difference of at least obs_stat?
This the standard workflow for performing hypothesis tests.
#7. Construct and record a confidence interval for the difference between the weights of those who exercise at least three times a week and those who don’t, and interpret this interval in context of the data. CI = {x} z x = sample mean z = confidence level value s = Sample Standard deviation n = sampel size
#Standard deviation
yrbss %>%
group_by(physical_3plus) %>%
summarise(sd_weight = sd(weight, na.rm = TRUE))## # A tibble: 3 × 2
## physical_3plus sd_weight
## <chr> <dbl>
## 1 no 17.6
## 2 yes 16.5
## 3 <NA> 17.6#Mean
yrbss %>%
group_by(physical_3plus) %>%
summarise(mean_weight = mean(weight, na.rm = TRUE))## # A tibble: 3 × 2
## physical_3plus mean_weight
## <chr> <dbl>
## 1 no 66.7
## 2 yes 68.4
## 3 <NA> 69.9#Sample size N
yrbss %>%
group_by(physical_3plus) %>%
summarise(freq = table(weight)) %>%
summarise(n = sum(freq))## # A tibble: 3 × 2
## physical_3plus n
## <chr> <int>
## 1 no 4022
## 2 yes 8342
## 3 <NA> 215# not Active
not_active_mean <- 66.7
not_active_sd <- 17.6
not_active_n <- 4022
# active
active_mean <- 68.4
active_sd <- 16.5
active_n <- 8342
z <- 1.96
# confidence interval for not active
upper_not_active <- not_active_mean + z * (not_active_sd / sqrt(not_active_n))
upper_not_active## [1] 67.24394From this calculation we can state 95% CI, those students who are active at least three times a week have an average weight of 68.05 kg and 68.75 kg. Not so active students have at least 66.16kg to 67.24kg. # 8 Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.
height_table <- as.data.frame(table(yrbss$height))
height_freq <- sum(height_table$Freq)
# mean, standard deviation and sample size
height_mean <- mean(yrbss$height, na.rm = TRUE)
height_mean## [1] 1.691241height_sd <- sd(yrbss$height, na.rm = TRUE)
height_sd## [1] 0.1046973height_n <- yrbss %>%
summarise(freq = table(height)) %>%
summarise(n = sum(freq, na.rm = TRUE))
height_n## # A tibble: 1 × 1
## n
## <int>
## 1 12579z_height <- 1.96Finding the upper/lower height:
# confidence interval for height
upper_height<- height_mean + z_height * (height_sd / sqrt(height_n))
upper_height## n
## 1 1.693071lower_height <- height_mean - z_height * (height_sd / sqrt(height_n))
lower_height## n
## 1 1.689411At 95% CI, the average height in meters for the students is between 1.689411 to 1.693071m
#9 Calculate a new confidence interval for the same parameter at the 90% confidence level. Comment on the width of this interval versus the one obtained in the previous exercis
# set z value to 1.65 for 90% confidence interval
z_90 <- 1.65
#confidence interval for height
upper_height_90 <- height_mean + z_90 * (height_sd / sqrt(height_n))
upper_height_90## n
## 1 1.692781lower_height_90 <- height_mean - z_90 * (height_sd / sqrt(height_n))
lower_height_90## n
## 1 1.689701range_95 <- (upper_height - lower_height)
range_95 ## n
## 1 0.003659302After finding the 95 CI, look for 90
range_90 <- (upper_height_90 - lower_height_90)
range_90## n
## 1 0.003080535# difference between the two ranges
diff_range <- range_95 - range_90
diff_range## n
## 1 0.0005787672Based on calculation at 90% CI vs 95% CI there is a slight difference in range and slight difference in CI.
#10 Conduct a hypothesis test evaluating whether the average height is different for those who exercise at least three times a week and those who don’t.
#11 Now, a non-inference task: Determine the number of different
options there are in the dataset for the hours_tv_per_school_day there
are. There are 7 different options with the 8th being
yrbss %>%
group_by(hours_tv_per_school_day)%>%
summarise(n())## # A tibble: 8 × 2
## hours_tv_per_school_day `n()`
## <chr> <int>
## 1 <1 2168
## 2 1 1750
## 3 2 2705
## 4 3 2139
## 5 4 1048
## 6 5+ 1595
## 7 do not watch 1840
## 8 <NA> 338