title: “Lab 5: Statistical Inference”
author: “Evan McLaughlin”
date: “10.4.2020
knitr::opts_chunk$set(eval = TRUE, results = FALSE, fig.show = "show", message = FALSE)

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.0.4
library(openintro)
library(ggplot2)
library(dplyr)
library(infer)
## Warning: package 'infer' was built under R version 4.0.3
library(trelliscopejs)
## Warning: package 'trelliscopejs' was built under R version 4.0.3
global_monitor <- tibble(scientist_work = c(rep("Benefits", 80000), rep("Doesn't benefit", 20000)))

ggplot(global_monitor, aes(x = scientist_work)) + geom_bar() +  labs(x = "", y = "", title = "Do you believe that the work scientists do benefit people like you?") +  coord_flip()

global_monitor %>%
  count(scientist_work) %>%
  mutate(p = n /sum(n))
samp1 <- global_monitor %>%
  sample_n(50)

samp1 %>%
  count(scientist_work) %>%
  mutate(p_samp = n /sum(n))

Exercise A-1

Describe the distribution of responses in this sample. How does it compare to the distribution of responses in the population. Hint: Although the sample_n function takes a random sample of observations (i.e. rows) from the dataset, you can still refer to the variables in the dataset with the same names. Code you presented earlier for visualizing and summarising the population data will still be useful for the sample, however be careful to not label your proportion p since you’re now calculating a sample statistic, not a population parameters. You can customize the label of the statistics to indicate that it comes from the sample.

The sample is largely comparable to the larger population, as we would assume.

###samp1 %>%
  ###count(scientist_work) %>%
  ###mutate(p_hat = n /sum(n))

Exercise A-2

Would you expect the sample proportion to match the sample proportion of another student’s sample? Why, or why not? If the answer is no, would you expect the proportions to be somewhat different or very different? Ask a student team to confirm your answer.

We would expect the samples to align fairly closely with one another as we are all drawing from the same sample with the same proportion of respondents who believe and don’t believe that scientists help them.

Exercise A-3

Take a second sample, also of size 50, and call it samp2. How does the sample proportion of samp2 compare with that of samp1? Suppose we took two more samples, one of size 100 and one of size 1000. Which would you think would provide a more accurate estimate of the population proportion?

The new sample of 50 strays a bit more from the first sample and the population as a whole, but not overwhelmingly. We would expect the larger sample proportion to hew more closely to the full population since its larger sample size makes it more similar to the population it seeks to describe.

samp2 <- global_monitor %>%
  sample_n(50)
samp2 %>%
  count(scientist_work) %>%
  mutate(p_samp2 = n /sum(n))

samp_h <- global_monitor %>%
  sample_n(100)
samp_h %>%
  count(scientist_work) %>%
  mutate(p_samp_h = n /sum(n))

samp_m <- global_monitor %>%
  sample_n(1000)
samp_m %>%
  count(scientist_work) %>%
  mutate(p_samp_m = n /sum(n))
sample_props50 <- global_monitor %>%
    rep_sample_n(size = 50, reps = 15000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Doesn't benefit")

ggplot(data = sample_props50, aes(x = p_hat)) +
geom_histogram(binwidth = 0.02) + labs(x = "p_hat (Doesn't benefit)", title = "Sampling distribution of p_hat", subtitle = "Sample size = 50, Number of samples = 15000"
  )

sd <- sqrt((mean(sample_props50$p_hat)) * (1-mean(sample_props50$p_hat)) / 50)
### Not sure why this didn't work - Evan to check

sd

sample_props50 %>% 
  ggplot() +  geom_histogram(aes(x = p_hat), binwidth = 0.075) +  xlab("sample props")

### Exercise A-4

How many elements are there in sample_props50? Describe the sampling distribution, and be sure to specifically note its center. Make sure to include a plot of the distribution in your answer.

“Replicate”, “scientist_work”, “p_hat”, and “n” constitute four elements. Sample_props50 is normally distributed around a mean of .2, with a standard deviation of ~.06.

global_monitor %>%
  sample_n(size = 50, replace = TRUE) %>%
  count(scientist_work) %>%
  mutate(p_hat = n /sum(n)) %>%
  filter(scientist_work == "Doesn't benefit")

Exercise A-5

To make sure you understand how sampling distributions are built, and exactly what the rep_sample_n function does, try modifying the code to create a sampling distribution of 25 sample proportions from samples of size 10, and put them in a data frame named sample_props_small. Print the output. How many observations are there in this object called sample_props_small? What does each observation represent?

25 observations, each of which represents an instance that a “Doesn’t Benefit” response is recorded.

sample_props_small <- global_monitor %>%
  rep_sample_n(size = 10, reps = 25, replace = TRUE) %>%
  count(scientist_work) %>%
  mutate(p_hat = n /sum(n)) %>%
  filter(scientist_work == "Doesn't benefit")
sample_props_small

 ggplot(data = sample_props50, aes(x = p_hat)) +
    geom_histogram(binwidth = 0.02)

### Exercise A-6

Use the app below to create sampling distributions of proportions of Doesn’t benefit from samples of size 10, 50, and 100. Use 5,000 simulations. What does each observation in the sampling distribution represent? How does the mean, standar error, and shape of the sampling distribution change as the sample size increases? How (if at all) do these values change if you increase the number of simulations? (You do not need to include plots in your answer.)

The distribution with 10 is right skewed, while the 50 and 100 samples are more normally distributed.

Exercise B

us_adults <- tibble(climate_change_affects = c(rep("Yes", 62000), rep("No", 38000)))
ggplot(us_adults, aes(x = climate_change_affects)) +
  geom_bar() +  labs(x = "", y = "", title = "Do you think climate change is affecting your local community?") + coord_flip()

us_adults %>%
  count(climate_change_affects) %>%
  mutate(p = n /sum(n))