update hyp test

Author: WUIMSM

07-statistical_inference.R

@@ -1,6 +1,8 @@
+## ---- include=FALSE------------------------------------------------------
+knitr::opts_chunk$set(echo = TRUE)
+options(digits = 7)
 
-
-## ------------------------------------------------------------------------
+## ---- message = FALSE, warning=FALSE-------------------------------------
 library(tidyverse)
 library(ggplot2)
 library(latex2exp)
@@ -24,42 +26,19 @@ ggplot(data.frame(hours)) +
   annotate("text", x = mean(hours) + 28, y = 1100, label = "Mean + 2 * SD" )+
   annotate("text", x = mean(hours) -28, y = 1100, label = "Mean - 2 * SD" )
 
-
-## ------------------------------------------------------------------------
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 4, fig.width = 6----
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
-m1 <- mean(hours[student_sample])
-m1
-
-## ------------------------------------------------------------------------
-set.seed(12345)
-samples <- 20000
-means <- matrix(NA, nrow = samples)
-for (i in 1:samples){
-  student_sample <- sample(1:25000, size = 100, replace = FALSE)
-  means[i,] <- mean(hours[student_sample])
-}
-meansdf <- data.frame('true' = mean(hours), 'sample' = mean(means))
-meansdf <- gather(meansdf)
-ggplot(data.frame(means)) +
-  geom_histogram(aes(x = means), bins = 30, fill='white', color='black') +
-  theme_bw() +
-  geom_vline(data = meansdf, aes(xintercept = value, color = key, linetype = key), size=1) +
-  scale_color_discrete(labels = c("Mean of sample means", "Population mean")) +
-  scale_linetype_discrete(labels = c("Mean of sample means", "Population mean")) +
-  theme(legend.title = element_blank(),legend.position = "bottom") +
-  ggtitle('Distribution of sample means')
-
-## ------------------------------------------------------------------------
-head(means, 5)
-min(means)
-max(means)
-mean(means)
+sample_1 <- hours[student_sample]
+ggplot(data.frame(sample_1)) +
+  geom_histogram(aes(x = sample_1), bins = 30, fill='white', color='black') +
+  theme_bw() + xlab("Hours") +
+  geom_vline(aes(xintercept = mean(sample_1)), size=1) +
+  ggtitle(TeX(sprintf("Distribution of listening times ($\\bar{x}$ = %.2f)",round(mean(sample_1),2))))
 
-## ------------------------------------------------------------------------
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 6, fig.width = 8----
+#student_sample <- sample(1:25000, size = 100, replace = FALSE)
+#means <- hours[student_sample]
 library(cowplot)
-library(gridExtra)
-library(grid)
-library(latex2exp)
 set.seed(8830)
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
 means1 <- hours[student_sample]
@@ -102,8 +81,61 @@ title <- ggdraw() + draw_label('Distribution of listening times in four differen
 p <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
 print(p)
 
-## ------------------------------------------------------------------------
+
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 4, fig.width = 6----
 set.seed(12345)
+samples <- 20000
+means <- matrix(NA, nrow = samples)
+for (i in 1:samples){
+  student_sample <- sample(1:25000, size = 100, replace = FALSE)
+  means[i,] <- mean(hours[student_sample])
+}
+
+meansdf <- data.frame('true' = mean(hours), 'sample' = mean(means))
+meansdf <- gather(meansdf)
+ggplot(data.frame(means)) +
+  geom_histogram(aes(x = means), bins = 30, fill='white', color='black') +
+  theme_bw() +
+  geom_vline(data = meansdf, aes(xintercept = value, color = key, linetype = key), size=1) +
+  scale_color_discrete(labels = c("Mean of sample means", "Population mean")) +
+  scale_linetype_discrete(labels = c("Mean of sample means", "Population mean")) +
+  theme(legend.title = element_blank(),
+        legend.position = "bottom") +
+  labs(title = "Histogram of listening times",
+       subtitle = TeX(sprintf("Population mean ($\\mu$) = %.2f; population standard deviation ($\\sigma$) = %.2f",round(mean(hours),2),round(sd(hours),2))),
+       y = 'Number of students', 
+       x = 'Hours') 
+
+## ----message=FALSE, warning=FALSE, eval=TRUE, echo=FALSE, fig.align="center", fig.cap="Relationship between the sample size and the standard error"----
+set.seed(321)
+hours <- rnorm(25000, 50, 10)
+
+R <- 1000
+sems <- numeric()
+replication <- numeric()
+
+for (r in 10:R) {
+  y_sample <- sample(hours, r)
+  sem <- sd(hours)/sqrt(length(y_sample))
+  sems <- rbind(sems, sem)
+  replication <- rbind(replication, r)
+}
+
+df <- as.data.frame(cbind(replication, sems))
+ggplot(data=df, aes(y = sems, x = replication)) + 
+  geom_line() + 
+  ylab("Standard error of the mean") + 
+  xlab("Sample size") + 
+  ggtitle('Relationship between sample size and standard error') +
+  theme_bw()
+
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 6, fig.width = 8----
+library(cowplot)
+library(gridExtra)
+library(grid)
+library(latex2exp)
+set.seed(12345)
+
 sample_size = 10
 samples <- 20000
 means <- matrix(NA, nrow = samples)
@@ -174,13 +206,20 @@ plot4 <- ggplot(data.frame(means)) +
 
 p <- plot_grid(plot1, plot2, plot3, plot4, ncol = 2,
            labels = c("A", "B","C","D"))
-title <- ggdraw() + draw_label('Distribution of sample means', fontface='bold')
+title <- ggdraw() + draw_label('Relationship between sample size and standard error', fontface='bold')
 p <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
 print(p)
+# now add the title
+#title <- ggdraw() + draw_label("", fontface='bold')
+#plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
 
 
-## ------------------------------------------------------------------------
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 3, fig.width = 8----
+library(cowplot)
+library(gridExtra)
+library(grid)
 set.seed(12345)
+
 hours1 <- rnorm(25000, 50, 1)
 sample_size = 100
 samples <- 20000
@@ -218,13 +257,17 @@ plot2 <- ggplot(data.frame(means)) +
   scale_linetype_discrete(labels = c("Mean of sample means", "Population mean")) +
   theme(legend.position = "none") +  ggtitle(TeX(sprintf("n = 100; $\\sigma = 10$; $\\sigma_{\\bar x}$ = %.2f",round(sd(hours2)/sqrt(sample_size),2))))
 
-p <- plot_grid(plot1, plot2, ncol = 1, 
+p <- plot_grid(plot1, plot2, ncol = 2, 
            labels = c("A", "B"))
-title <- ggdraw() + draw_label('Distribution of sample means', fontface='bold')
+title <- ggdraw() + draw_label('Relationship between population SD and standard error', fontface='bold')
 p <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
 print(p)
+# now add the title
+#title <- ggdraw() + draw_label("", fontface='bold')
+#plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
+
 
-## ------------------------------------------------------------------------
+## ----message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE, fig.align="center", fig.height = 4, fig.width = 6----
 set.seed(321)
 hours <- rgamma(25000, shape = 2, scale = 10)
 ggplot(data.frame(hours)) +
@@ -235,14 +278,19 @@ ggplot(data.frame(hours)) +
        y = 'Number of students', 
        x = 'Hours') 
 
+
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 6, fig.width = 8----
+#student_sample <- sample(1:25000, size = 100, replace = FALSE)
+#means <- hours[student_sample]
+
 set.seed(8830)
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
 means1 <- hours[student_sample]
 plot1 <- ggplot(data.frame(means1)) +
   geom_histogram(aes(x = means1), bins = 30, fill='white', color='black') +
   theme_bw() + xlab("Hours") +
   geom_vline(aes(xintercept = mean(means1)), size=1) +
- ggtitle(TeX(sprintf("$\\bar{x}_1$ = %.2f",round(mean(means1),2))))
+ ggtitle(TeX(sprintf("$\\bar{x}$ = %.2f",round(mean(means1),2))))
 
 set.seed(6789)
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
@@ -251,7 +299,7 @@ plot2 <- ggplot(data.frame(means1)) +
   geom_histogram(aes(x = means1), bins = 30, fill='white', color='black') +
   theme_bw() + xlab("Hours") +
   geom_vline(aes(xintercept = mean(means1)), size=1) +
- ggtitle(TeX(sprintf("$\\bar{x}_2$ = %.2f",round(mean(means1),2))))
+ ggtitle(TeX(sprintf("$\\bar{x}$ = %.2f",round(mean(means1),2))))
 
 set.seed(3904)
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
@@ -260,7 +308,7 @@ plot3 <- ggplot(data.frame(means1)) +
   geom_histogram(aes(x = means1), bins = 30, fill='white', color='black') +
   theme_bw() + xlab("Hours") +
   geom_vline(aes(xintercept = mean(means1)), size=1) +
- ggtitle(TeX(sprintf("$\\bar{x}_3$ = %.2f",round(mean(means1),2))))
+ ggtitle(TeX(sprintf("$\\bar{x}$ = %.2f",round(mean(means1),2))))
 
 set.seed(3333)
 student_sample <- sample(1:25000, size = 100, replace = FALSE)
@@ -269,7 +317,7 @@ plot4 <- ggplot(data.frame(means1)) +
   geom_histogram(aes(x = means1), bins = 30, fill='white', color='black') +
   theme_bw() + xlab("Hours") +
   geom_vline(aes(xintercept = mean(means1)), size=1) +
- ggtitle(TeX(sprintf("$\\bar{x}_4$ = %.2f",round(mean(means1),2))))
+ ggtitle(TeX(sprintf("$\\bar{x}$ = %.2f",round(mean(means1),2))))
 
 p <- plot_grid(plot1, plot2, plot3, plot4, ncol = 2,
            labels = c("A", "B","C","D"))
@@ -278,9 +326,14 @@ p <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values con
 print(p)
 
 
-## ------------------------------------------------------------------------
+## ----message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 6, fig.width = 8----
+library(cowplot)
+library(gridExtra)
+library(grid)
 set.seed(321)
+
 hours <- rgamma(25000, shape = 2, scale = 10)
+
 samples <- 10
 means <- matrix(NA, nrow = samples)
 for (i in 1:samples){
@@ -344,7 +397,7 @@ p <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values con
 print(p)
 
 
-## ------------------------------------------------------------------------
+## ---- fig.height = 4, fig.width=6----------------------------------------
 set.seed(321)
 hours <- rgamma(25000, shape = 2, scale = 10)
 
@@ -357,33 +410,21 @@ plot2 <- ggplot(data.frame(hours_s)) +
   geom_histogram(aes(x = hours_s), bins = 30, fill='white', color='black') +
   theme_bw() + xlab("Hours") +
   geom_vline(aes(xintercept = mean(hours_s)), size=1) +
- ggtitle(TeX(sprintf("$n$ = %d; $\\bar{x}$ = %.2f; $s$ = %.2f",sample_size,round(mean(hours_s),2),round(sd(hours_s),2))))
+ ggtitle(TeX(sprintf("Random sample; $n$ = %d; $\\bar{x}$ = %.2f; $s$ = %.2f",sample_size,round(mean(hours_s),2),round(sd(hours_s),2))))
 plot2
 
-
-## ------------------------------------------------------------------------
+## ---- fig.height = 4, fig.width=6----------------------------------------
 qnorm(0.975)
-qnorm(0.025)
 
-## ------------------------------------------------------------------------
+## ---- fig.height = 4, fig.width=6----------------------------------------
 sample_mean <- mean(hours_s)
 se <- sd(hours_s)/sqrt(sample_size)
-ci_upper <- sample_mean + qnorm(0.975)*se
 ci_lower <- sample_mean - qnorm(0.975)*se
-ci_upper
+ci_upper <- sample_mean + qnorm(0.975)*se
 ci_lower
+ci_upper
 
-plot2 <- ggplot(data.frame(hours_s)) +
-  geom_histogram(aes(x = hours_s), bins = 30, fill='white', color='black') +
-  theme_bw() + xlab("Hours") +
-  geom_vline(aes(xintercept = sample_mean), size=1) +
-  geom_vline(aes(xintercept = ci_upper), size=1, color="red") +
-  geom_vline(aes(xintercept = ci_lower), size=1, color="red") +
-  ggtitle(TeX(sprintf("$n$ = %d; $\\bar{x}$ = %.2f; $s$ = %.2f",sample_size,round(mean(hours_s),2),round(sd(hours_s),2))))
-plot2
-
-
-## ------------------------------------------------------------------------
+## ---- fig.height = 15, fig.width=10--------------------------------------
 set.seed(12)
 samples <- 100
 hours <- rgamma(25000, shape = 2, scale = 10)
@@ -404,3 +445,5 @@ ggplot2::ggplot(means_sd, aes(y = y)) +
   scale_color_manual(values = c("red", "black")) +
   guides(color=guide_legend(title="True mean in CI")) +
   theme_bw()
+
+