###########################
# MS-C1620
# Statistical inference
# Lecture 2
library(scales)
library(EnvStats)
####################################################################################
# Comparison of confidence intervals for expected value
# Highly skew data set
x <- rivers
n <- length(x)
boxplot(x)
# Manual bootstrapping for the expected value
B <- 1000
res <- rep(0, B)
for(b in 1:B){
res[b] <- mean(sample(x, n, replace = TRUE))
}
# Distribution of the bootstrap means
hist(res)
# 95% BS confidence interval
ci_bs <- quantile(res, probs = c(0.025, 0.975))
abline(v = ci_bs, lwd = 2, col = 2)
ci_bs
# Parametric 95% confidence interval assuming normality
ci_ex <- c(mean(x) - 1.96*sd(x)/sqrt(n), mean(x) + 1.96*sd(x)/sqrt(n))
abline(v = ci_ex, lwd = 2, col = 3)
ci_ex
####################################################################################
# BS confidence interval for skewness
library(moments)
# Manual bootstrapping
B <- 1000
res <- rep(0, B)
for(b in 1:B){
res[b] <- skewness(sample(x, n, replace = TRUE))
}
# Distribution of the bootstrap means
hist(res)
# 95% BS confidence interval
ci_bs <- quantile(res, probs = c(0.025, 0.975))
abline(v = ci_bs, lwd = 2, col = 2)
ci_bs
####################################################################
# Illustration of Type I and Type II errors using t-test (optional extra example)
# Single run
n <- 15
theta <- 0
y <- rnorm(n, theta, 3)
# Our test statistic for H0: theta = 0
t_observed_0 <- (mean(y) - 0)/(sd(y)/sqrt(n))
t_observed_0
# Is this too large to happen under H0?
# Find The distribution of t under H0
t_generate <- function(theta){
n <- 15
y <- rnorm(n, theta, 3)
(mean(y) - 0)/(sd(y)/sqrt(n))
}
theta_0 <- replicate(10000, t_generate(0))
hist(theta_0, breaks = 30)
# We choose cut-off points for what is still considered as "normal behavior" under H0 (95% interval)
# (note that this is an approximation and we could also use the true distribution, Student's t with n - 1 degrees of freedom)
quant_0 <- quantile(theta_0, c(0.025, 0.975))
abline(v = quant_0[1], col = "red", lwd = 2)
abline(v = quant_0[2], col = "red", lwd = 2)
# Type I error rate is now roughly 5%
# Where does the observed value fall?
abline(v = t_observed_0, col = "black", lwd = 2)
# p-value measures how far away we are in the tail/how weird the observation is.
mean((theta_0 >= abs(t_observed_0) | theta_0 <= -abs(t_observed_0)))
# Type I error rate is 5% (probability of claiming there is an effect when there is not)
# What about type II errors?
# The distribution of t when H1: theta = 1 is true
theta_1 <- replicate(10000, t_generate(1))
hist(theta_0, breaks = 30, xlim = c(-5, 7))
abline(v = quant_0[1], col = "red", lwd = 2)
abline(v = quant_0[2], col = "red", lwd = 2)
hist(theta_1, breaks = 30, add = TRUE, col = alpha("lightblue", 0.4), lty = 2)
# Type II error rate for theta = 1
mean(theta_1 < quant_0[2] & theta_1 > quant_0[1])
# Conduct a single study where H1: theta = 1 is true
t_observed_1 <- t_generate(1)
abline(v = t_observed_1, col = "black", lwd = 2)
# p-value
mean((theta_0 >= abs(t_observed_1) | theta_0 <= -abs(t_observed_1)))
# Making the "normally behaving" area wider makes true H0 easier to recognize (lower type I)
# but deviations from H0 more difficult to recognize (higher type II) and vice versa
# Recognizing deviations from H0 gets more difficult as we approach theta = 0
theta_01 <- replicate(10000, t_generate(0.1))
hist(theta_0, breaks = 30, xlim = c(-5, 5))
abline(v = quant_0[1], col = "red", lwd = 2)
abline(v = quant_0[2], col = "red", lwd = 2)
hist(theta_01, breaks = 30, add = TRUE, col = alpha("salmon", 0.4), lty = 2)
# Type II error rate for theta = 1
mean(theta_01 < quant_0[2] & theta_01 > quant_0[1])
# Conduct a single study where H0 is false
t_observed_01 <- t_generate(0.1)
abline(v = t_observed_01, col = "black", lwd = 2)
# p-value
mean((theta_0 >= abs(t_observed_01) | theta_0 <= -abs(t_observed_01)))
####################################################################################
# One sample t-test
# Data which show the effect of two soporific drugs (increase in hours of sleep compared to control) on 10 patients.
sleep
# Compare whether the first drugs differs in its effect from placebo at significance level 5%
# H0: mu1 == 0
# H1: mu1 != 0
plot(sleep[sleep$group == 1, 1], rep(1, 10))
t.test(sleep[sleep$group == 1, 1])
# p-value 0.218 - not enough evidence against H0
####################################################################################
# Two-sample t-test
# Data which show the effect of two soporific drugs (increase in hours of sleep compared to control) on 10 patients.
sleep
# Compare whether the drugs differ in their effect at significance level 5%
# H0: mu1 == mu2
# H1: mu1 != mu 2
plot(sleep[sleep$group == 1, 1], rep(1, 10))
plot(sleep[sleep$group == 2, 1], rep(1, 10))
t.test(sleep[sleep$group == 1, 1], sleep[sleep$group == 2, 1])
# p-value 0.079 - not enough evidence against H0
####################################################################################
# Variance test
# This data set provides measurements of the girth, height and volume of timber in 31 felled black cherry trees.
trees
# We test the conjecture that the standard deviation of the height of the felled trees equals 5
# Test the following hypotheses with the significance level 5%
# H0: sigma^2 == 5^2 == 25
# H1: sigma^2 != 25
varTest(trees$Height, sigma.squared = 25)
# p-value 0.034 - H0 not plausible
####################################################################################
# Variance comparison test
# Compare whether the variance of the drug effect differ in the sleep data at significance level 5%
# H0: sigma_1^2 == sigma_2^2
# H1: sigma_1^2 != sigma_2^2
var.test(sleep[sleep$group == 1, 1], sleep[sleep$group == 2, 1])
# p-value 0.743 - no evidence against H0