###########################
# MS-C1620
# Statistical inference
# Lecture 5
# We use the significance level 0.05 throughout the script
################################################
# Normality testing
# Consider three data sets:
# 1. The sepal lengths of 100 irises of the species Versicolor and Virginica
# 2. The January sunspot numbers for the years 1749-1848
# 3. A random sample of size 100 from the standard normal distribution
x1 <- iris[51:150, 1]
x2 <- sunspots[seq(1, by = 12, length.out = 100)]
x3 <- rnorm(100)
# Do the samples look normal?
hist(x1, breaks = 10)
hist(x2, breaks = 10)
hist(x3, breaks = 10)
# Multiple histograms of normal samples
par(mfrow = c(4, 2))
invisible(replicate(8, hist(rnorm(100), breaks = 10)))
par(mfrow = c(1, 1))
# Bowman-Shenton/Jarque-Bera test
library(tseries)
jarque.bera.test(x1)
# No evidence against normality
jarque.bera.test(x2)
# Strong evidence aganist normality
jarque.bera.test(x3)
# QQ-plot
qqnorm(x1)
qqline(x1)
# The sample is slightly right-skew
qqnorm(x2)
qqline(x2)
# The sample is heavily right-skew
qqnorm(x3)
qqline(x3)
# Shapiro-Wilk test
shapiro.test(x1)
# No evidence against normality
shapiro.test(x2)
# Strong evidence against normality
shapiro.test(x3)
#####################################################################
# Multinomial distribution
# Probability of having eye color
# BROWN | GREEN | BLUE
# 0.43 | 0.32 | 0.25
# We observe 292 people
# A sample from Multinomial(292, 0.43, 0.32, 0.25)
rmultinom(1, 292, c(0.43, 0.32, 0.25))
#####################################################################
# Chi-squared tests
##########################
# X^2 goodness-of-fit for discrete data
#Distribution of hair and eye color and sex in 592 statistics students.
HairEyeColor
# H0: the distribution (probabilities) of hair colors is (0.20, 0.50, 0.10, 0.20) for
# (Black, Brown, Red, Blond)?
hair <- margin.table(HairEyeColor, 1)
hair
hair/sum(hair)
# Observed and expected
hair
sum(hair)*c(0.2, 0.5, 0.1, 0.2)
chisq.test(hair, p = c(0.20, 0.50, 0.10, 0.20))
# No evidence against H0
##########################
# X^2 goodness-of-fit for discretized data
# Is the sunspot data we used before from an exponential distribution
x2 <- sunspots[seq(1, by = 12, length.out = 100)]
hist(x2, breaks = 15, probability = TRUE, ylim = c(0, 0.022))
# Estimate the parameters
lambdahat <- 1/mean(x2)
lambdahat
expfit <- function(x){ dexp(x, rate = lambdahat) }
curve(expfit, add = TRUE, col = 3, lwd = 2)
# Break the distribution into 8 categories, (0, 25, 50,...) + leftover category
hist(x2, breaks = 15, xlim = c(-1, 201), ylim = c(0, 23))
abline(v = seq(0, 200, by = 25), col = "red", lwd = 2)
# Observed frequencies
cum_x2_obs <- sapply(seq(0, 200, by = 25), function(i) sum(x2 <= i))[-1]
x2_obs <- c(cum_x2_obs[1], diff(cum_x2_obs), 0)
x2_obs
# Expected frequencies
curve(expfit, col = 3, lwd = 2, from = 0, to = 201)
abline(h = 0, lwd = 1)
segments(seq(0, 200, by = 25), rep(0, 9), seq(0, 200, by = 25), expfit(seq(0, 200, by = 25)), col = 2, lwd = 2)
p_exp <- diff(pexp(seq(0, 200, by = 25), rate = lambdahat))
x2_exp <- 100*p_exp
c(x2_exp, 100 - sum(x2_exp))
# The probabilities don't sum to one so we add one extra category
# Degrees of freedom is "categories" - 1 - "est. parameters" = 9 - 1 - 1 = 7
pchisq(chisq.test(x2_obs, p = c(p_exp, 1 - sum(p_exp)))$statistic, df = 7, lower.tail = FALSE)
# No evidence against the null that the data comes from an exponential distribution
##########################
# X^2 test of homogeneity
# This data set provides information on the fate of passengers on the fatal maiden voyage
# of the ocean liner 'Titanic', summarized according to economic status (class), sex, age and survival.
Titanic
# Let's assume that the 1st, 2nd and 3rd class ticket amounts (and the crew size) were fixed and
# each category had a specific probability to survive
Titanic_2 <- margin.table(Titanic, c(1, 4))
Titanic_2
# Is the probability to survive the same in each category?
chisq.test(Titanic_2)
# Strong evidence against H0
##########################
# X^2 test of independence
# Cross-classification of a sample of British males according to each subject's
# occupational status and his father's occupational status.
occupationalStatus
heatmap(occupationalStatus, Colv = NA, Rowv = NA, revC = TRUE)
# Are the occupational statuses of a father and his son independent?
chisq.test(occupationalStatus)
# Strong evidence against H0
# (the warning has to do with some cell counts being very low,
# and the extreme classes could be combined to give a more reliable result)