Lab Name

#install.packages("infer")

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.0.2

## Warning: package 'tibble' was built under R version 4.0.2

## Warning: package 'tidyr' was built under R version 4.0.2

## Warning: package 'dplyr' was built under R version 4.0.2

library(openintro)

## Warning: package 'openintro' was built under R version 4.0.2

## Warning: package 'airports' was built under R version 4.0.2

## Warning: package 'cherryblossom' was built under R version 4.0.2

## Warning: package 'usdata' was built under R version 4.0.2

library(infer)

## Warning: package 'infer' was built under R version 4.0.2

set.seed(011774)

Introductory Stuff

global_monitor <- tibble(
  scientist_work = c(rep("Benefits", 80000), rep("Doesn't benefit", 20000))
)

write.csv(global_monitor, "gm.csv")
scientist_work <- read.csv("gm.csv")
gm <- read.csv("gm.csv")
rownames(gm) <- NULL

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

global_monitor %>%
  count(scientist_work) %>%
  mutate(p = n/sum(n))

## # A tibble: 2 x 3
##   scientist_work      n     p
##   <chr>           <int> <dbl>
## 1 Benefits        80000   0.8
## 2 Doesn't benefit 20000   0.2

Exercise 1

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

## # A tibble: 2 x 3
##   scientist_work      n     s
##   <chr>           <int> <dbl>
## 1 Benefits           40   0.8
## 2 Doesn't benefit    10   0.2

Describe the distribution of responses from this sample.

This distribution is off by 0.04. But notice, if I multiply the sample size by 10, I get a sampling distribution closer to the population mean. Now it’s less than 0.2.

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

## # A tibble: 2 x 3
##   scientist_work      n    s2
##   <chr>           <int> <dbl>
## 1 Benefits          407 0.814
## 2 Doesn't benefit    93 0.186

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

## # A tibble: 2 x 3
##   scientist_work      n p_hat
##   <chr>           <int> <dbl>
## 1 Benefits           40   0.8
## 2 Doesn't benefit    10   0.2

Exercise 2

According to Lohr, Sampling Design and Analysis, 2nd Edition, p. 29, sampling distributions operate the same as confidence intervals. Take one sampling distribution, it may or may not represent the population parameters. But if multiple sampling distributions are taken, 95% of them will be close to the population parameters. Thus the sampling distribution from another student may or may not be close to the population parameters, but if we took the sampling distributions from all the students in the class, most of them would be at or near the population parameters.

Exercise 3

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

## # A tibble: 2 x 3
##   scientist_work      n    s3
##   <chr>           <int> <dbl>
## 1 Benefits           34  0.68
## 2 Doesn't benefit    16  0.32

In this case, my two samples of 50 were the same.

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" )

Exercise 4

There are 14,999 observations in sample_props50. (Is this what is meant by elements?) The plot of sample_props50 appears to be normal, but see the long-winded explanation below, complete w/ Shapiro-Wilk and a histogram w/ a normal curve superimposed.

g =sample_props50
m <- mean(sample_props50$p_hat)
std <-sd(sample_props50$p_hat)
ggplot (data =sample_props50, aes(x=p_hat)) +
  geom_blank() +
  geom_histogram(aes(y = ..density.. )) +
stat_function(fun = dnorm, args = c(mean = m, sd = std), col = "tomato")

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

write.csv(sample_props50, file = "sample_props50.csv")

  #### I tried running Shapiro-Wilk, but R wouldn't let me do it because the sample size was > 5,000. The NCSS Shapiro-Wilk documentation says "The test was developed by Shapiro and Wilk (1965) for samples up to 20. NCSS uses the approximations

suggested by Royston (1992) and Royston (1995) which allow unlimited sample sizes. Note that Royston only checked the results for sample sizes up to 5000, but indicated that he saw no reason larger sample sizes should not work." Therefore, NCSS allows unlimited sample sizes.

I saved the sample_props50.csv file, imported it into NCSS, and ran the normality tests. The curve looks normal, but the Shapiro-Wilks test said it wasn’t.

  As we discussed in class today, if you have a sufficiently large sample size, you will always have statistical significance. With a sample size of 14,999, the Shapiro-Wilk test found statistical significance, and thus said the sample was not normal. But see below. I ran the Shaprio-Wilk test against a sample of 50 and it showed strong evidence of not being normal. Large or small, samples seem to show the data are not normally distributed, in spite of the normal curve on the histogram.  
  
  This also proves the point of the author of the Stata book I mentioned in class a couple of weeks ago. He  told me we should all know multiple stats packages because one doesn't do everything. R doesn't do Shapiro-Wilk on datasets with > 5,000 observations. NCSS does. Gotta have both. 
   
                Test           Prob   10% Critical  5% Critical Decision

Test Name Value Level Value Value -5% Shapiro-Wilk W 0.9870675 0 Reject normality Anderson-Darling 81.54572 0 Reject normality Martinez-Iglewicz 0.9810345 0.9926747 0.9911315 Can’t reject normality Kolmogorov-Smirnov 0.08221351 0.011 0.012 Reject normality D’Agostino Skewness 10.17995 0 1.645 1.96 Reject normality D’Agostino Kurtosis -0.9085 0.363607 1.645 1.96 Can’t reject normality D’Agostino Omnibus 104.4567 0 4.605 5.991 Reject normality

#write.csv(samp3, file = "samp3.csv")
samp3 <- read.csv("samp3.csv")
shapiro.test(samp3$number)

## 
##  Shapiro-Wilk normality test
## 
## data:  samp3$number
## W = 0.44124, p-value = 1.573e-12

Exercise at the top of page 5, showing the code needed to generate 1 sample of 50. rep_sample_n saves doing this X number of times.

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

## # A tibble: 1 x 3
##   scientist_work      n p_hat
##   <chr>           <int> <dbl>
## 1 Doesn't benefit     8  0.16

Exercise 5 25 Sample proportions from samples of size 10.

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")

There are 23 observations in this dataset. Each observation represents a different sample. p_hat in each row indicates the probability of getting a value of “Doesn’t Benefit” in that sample.

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

### Exercise 6.

Each observation represents a sample As the sample size increases, the histogram bunches together. The mean doesn’t change much – w/ 10 observations, it is 0.22; with the rest (50, 100, 1,000) it stays at 0.2. The SE, however, decreases from 0.11 to 0.06 to 0.04 to 0.01. NB: I ran a sample size of 1,000 just to see what it would do.

Exercise 7. I assumed for this exercise we were to generate only 1 sample.

global_monitor %>%
sample_n(size = 15, replace = TRUE) %>%
count(scientist_work) %>%
mutate(p_hat = n /sum(n)) %>%
filter(scientist_work == "Benefits")

## # A tibble: 1 x 3
##   scientist_work     n p_hat
##   <chr>          <int> <dbl>
## 1 Benefits          11 0.733

I think this one sample would give a higher proportion of “Benefits” than a greater number of samples would give.

Exercise 8

sample_props15 <- global_monitor %>%
                    rep_sample_n(size = 15, reps = 2000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

e8 =sample_props15
m8 <- mean(sample_props15$p_hat)
std8 <-sd(sample_props15$p_hat)
med <- median(sample_props15$p_hat)
ggplot (data =sample_props15, aes(x=p_hat)) +
  geom_blank() +
  geom_histogram(aes(y = ..density.. )) +
stat_function(fun = dnorm, args = c(mean = m8, sd = std8), col = "Blue")

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

In this histogram, it appears the proportion of “Benefits” in the population is ~ 80%. The distribution is skewed to the left, meaning more of the samples have a greater proportion of “Benefits.”

hist(sample_props15$p_hat, main = "Histogram of p_hat, sample size 15", xlab = "p_hat")
abline(v = m8, col = "red")
abline(v = med, col = "blue")

Based on the above histogram w/ the mean and medians overlaying, it seems the proportion of “Benefits” in the population is 80%.

Exercise 9

sample_props150 <- global_monitor %>%
                    rep_sample_n(size = 150, reps = 2000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

e9 =sample_props15
m9 <- mean(sample_props150$p_hat)
std9 <-sd(sample_props150$p_hat)
med9 <- median(sample_props150$p_hat)
ggplot (data =sample_props150, aes(x=p_hat)) +
  geom_blank() +
  geom_histogram(aes(y = ..density.. )) +
stat_function(fun = dnorm, args = c(mean = m9, sd = std9), col = "Green")

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

hist(sample_props150$p_hat, main = "Histogram of p_hat, sample size 150", xlab = "p_hat")
abline(v = m9, col = "red")
abline(v = med9, col = "blue")

Given the 2nd histogram, the population mean seems to always be 0.8, regardless of the size of the sample of the shape of the sampling distribution.

The second sampling distribution is shifted markedly to the left, but the curve seems to show a normal distribution. The 2nd histogram confirms that.

Exercise 10

NB: I assume the authors of this Lab mean 8 & 9, not 2 & 3. The distribution w/ the smaller sample size (15) actually has the larger spread, which is counterintuitive, but in line w/ the results generated by the app. If you want a sampling distribution w/ a more narrow spread, use a larger sample size.