###########################
# MS-C1620
# Statistical inference
# Lecture 6


# We use the significance level 0.05 throughout the script


###########################
# Forms of dependency
# EX1

# Compare the following forms of dependency


### Dependency between continuous and continuous variable
### Sepal length vs sepal width of setosa-irises
iris
plot(iris[1:50, 1], iris[1:50, 2], xlab = "Sepal length", ylab = "Sepal width")



### Dependency between continuous and continuous variable
### Ozone content vs. temperature
airquality
plot(airquality[, 1], airquality[, 4], xlab = "Ozone content", ylab = "Temperature")


### Dependency between continuous and continuous variable
### Artificial data with circular shape
phi <- runif(400, 0, 2*pi)
r <- runif(400, 1, 1.5)
x <- r*cos(phi)
y <- r*sin(phi)
plot(x, y, asp = 1)

# Note: no dependency in polar coordinates!
plot(phi, r)


### Dependency between categorial and continuous variable
### Chicken weight vs. feed type
boxplot(weight ~ feed, data = chickwts)

### what if we treat feed as numerical 1...6
chickmatrix = data.matrix(chickwts)
chickx = chickmatrix[,2]
chicky = chickmatrix[,1]
plot(chickx,chicky)


### Dependency within a single sample of continuous variable
### Approval ratings of US presidents
plot(presidents)


### Dependency within and between multiple samples of continuous variables
### Daily closing prices of major European stock indices
plot(EuStockMarkets)





###########################
# Linear dependence

### Perfect linear dependence
x <- rnorm(100)
y <- 10 - 2*x

plot(x, y)


### Zero linear dependence (but some other dependence still)
x <- rnorm(10000)
y <- x^2
plot(x, y)
cor(x, y)


### Correlations for some of the previous examples


# Iris
plot(iris[1:50, 1], iris[1:50, 2], xlab = "Sepal length", ylab = "Sepal width")
cor(iris[1:50, 1], iris[1:50, 2])
# (Note also the correlation matrix:)
cor(iris[1:50, 1:4])


# Air quality
plot(airquality[, 1], airquality[, 4], xlab = "Ozone content", ylab = "Temperature")
cor(airquality[, 1], airquality[, 4], use = "complete.obs")


# Circle data
phi <- runif(400, 0, 2*pi)
r <- runif(400, 1, 1.5)
x <- r*cos(phi)
y <- r*sin(phi)
plot(x, y)
cor(x, y)


# Correlation with the series itself!
plot(presidents)
delay = 1
plot(presidents[1:(120-delay)], presidents[(1+delay):120])
cor(presidents[1:(120-delay)], presidents[(1+delay):120], use = "complete.obs")

delay = 2
cor(presidents[1:(120-delay)], presidents[(1+delay):120], use = "complete.obs")

delay = 4
cor(presidents[1:(120-delay)], presidents[(1+delay):120], use = "complete.obs")

delay = 8
cor(presidents[1:(120-delay)], presidents[(1+delay):120], use = "complete.obs")



###########################
# EX2 Bivariate normal distribution

library(MASS)
library(KernSmooth)

# for rotateable plots (persp3d)
library(rgl)

# Bivariate normal with variances 1, means 0, and correlation 0.5
x <- mvrnorm(10000, c(0, 0), matrix(c(1, 0.5, 0.5, 1), 2, 2))
plot(x)
persp(bkde2D(x, 0.25)$fhat)
persp3d(bkde2D(x, 0.25)$fhat, col="green")   # rotateable plot

# Bivariate normal with variances 1, means 0, and correlation 0.8
x <- mvrnorm(10000, c(0, 0), matrix(c(1, 0.8, 0.8, 1), 2, 2))
plot(x)
persp(bkde2D(x, 0.25)$fhat)
persp3d(bkde2D(x, 0.25)$fhat, col="green")



###########################
# EX3 95% confidence intervals for Pearson correlation


### Sepal length vs sepal width of setosa-irises
x <- iris[1:50, 1]
y <- iris[1:50, 2]
n <- 50

# Is this bivariate normal?
plot(x, y)
persp(bkde2D(data.frame(x, y), 0.25)$fhat, theta = 45)
persp3d(bkde2D(data.frame(x, y), 0.25)$fhat, theta = 45, col="green")

cor(x, y)


### Parametric confidence interval under the assumption of normality
ci_par <- c(tanh(atanh(cor(x, y)) - 1.96/sqrt(n - 3)), tanh(atanh(cor(x, y)) + 1.96/sqrt(n - 3)))
ci_par

# Note: not symmetric around the estimate
plot(c(ci_par, cor(x, y)), rep(1, 3), xlim = c(0, 1))




### Non-parametric confidence interval
B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(cbind(x, y)[sample(1:n, n, replace = TRUE), ])[1, 2]
}

# Distribution of the bootstrap correlations
hist(res, breaks = 20)

# 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
abline(v = ci_par, lwd = 2, col = 3)

# Quite different answers by the parametric and non-parametric approaches -> safer to trust non-parametric







###################################
# Two-sample test for Pearson correlation

### Sepal length/sepal width correlations of setosa vs. versicolor
# Claim the correlations differ for the two species

# Setosa
x_setosa <- iris[1:50, 1]
y_setosa <- iris[1:50, 2]
n <- 50

plot(x_setosa, y_setosa)
cor(x_setosa, y_setosa)


# Versicolor
x_versicolor <- iris[51:100, 1]
y_versicolor <- iris[51:100, 2]
m <- 50

plot(x_versicolor, y_versicolor)
cor(x_versicolor, y_versicolor)


# Test whether the populatation correlation differ
# H0: rho1 == rho2
# H1: rho1 != rho2

z <- (atanh(cor(x_setosa, y_setosa)) - atanh(cor(x_versicolor, y_versicolor)))/sqrt(1/(n - 3) + 1/(m - 3))
2*pnorm(abs(z), lower.tail = FALSE)

# Not quite enough evidence against H0 -> we continue to believe that the correlations could be equally large





##################################
# Significance tests for correlation


### Petal length vs petal width of setosa-irises
# Is the correlation significant?
# H0: rho == 0
# H0: rho != 0

x <- iris[1:50, 3]
y <- iris[1:50, 4]
n <- 50
plot(x, y) # is this bivariate normal?
cor(x, y)


### Parametric test (null = zero correlation)
cor.test(x, y)
# -> enough evidence against null -> the correlation is significant if we can trust the normality of the data






### EX4 Non-parametric permutation test

B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(x, sample(y, n, replace = FALSE))
}

# Distribution of the permutation test replicates
hist(res, breaks = 20)


# Estimated probability of observing a more deviating value for the correlation under H0
abline(v = abs(cor(x, y)), col = 2, lwd = 2)
abline(v = -1*abs(cor(x, y)), col = 2, lwd = 2)
mean(res >= abs(cor(x, y)))
# Very similar p-value as in the parametric test -> same conclusion










##################################
# EX5 Spearman correlation


### Visualization

x1 <- c(1, 2, 3, 4, 5)
y1 <- c(1, 8, 27, 64, 125)
plot(x1, y1, type = "b")
# Spearman correlation ignores the "magnitude" by first "straightening" the data
plot(rank(x1), rank(y1), type = "b")


x2 <- c(1, 2, 3, 4, 5)
y2 <- c(1, 8, 27, 16, 125)
plot(x2, y2, type = "b")
plot(rank(x2), rank(y2), type = "b")



x3 <- c(1, 2, 3, 4, 5)
y3 <- c(1, 1.2, 4, 4.2, 7)
plot(x3, y3, type = "b")
plot(rank(x3), rank(y3), type = "b")






### Ozone content vs. temperature

full_obs <- (!is.na(airquality[, 1])) & (!is.na(airquality[, 4]))
x <- airquality[full_obs, 1]
y <- airquality[full_obs, 4]
n <- length(x)

plot(x, y, xlab = "Ozone content", ylab = "Temperature")





### Pearson correlation, boostrap confidence interval and a permutation test
cor(x, y)


### Bootstrap
B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(cbind(x, y)[sample(1:n, n, replace = TRUE), ])[1, 2]
}

# Distribution of the bootstrap correlations
hist(res, breaks = 20)

# 95% BS confidence interval
ci_bs <- quantile(res, probs = c(0.025, 0.975))
abline(v = ci_bs, lwd = 2, col = 2)
ci_bs



### Permutation test
B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(x, sample(y, n, replace = FALSE))
}

# Distribution of the permutation test replicates
hist(res, breaks = 20)

# Estimated probability of observing a more deviating value for the correlation under H0
mean(res >= abs(cor(x, y)))
# p-value = 0 -> Pearson correlation differs significantly from zero






### Spearman correlation, boostrap confidence interval and a permutation test
rx <- rank(x)
ry <- rank(y)
plot(rx, ry, xlab = "rank(Ozone content)", ylab = "rank(Temperature)")


cor(rx, ry)
# Spearman correlation is larger than Pearson -> there is more monotonic dependency inthe data than
# linear dependency


### Bootstrap
B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(cbind(rx, ry)[sample(1:n, n, replace = TRUE), ])[1, 2]
}

# Distribution of the bootstrap correlations
hist(res, breaks = 20)

# 95% BS confidence interval
ci_bs <- quantile(res, probs = c(0.025, 0.975))
abline(v = ci_bs, lwd = 2, col = 2)
ci_bs



### Permutation test
B <- 1000
res <- rep(0, B)

for(b in 1:B){
  res[b] <- cor(rx, sample(ry, n, replace = FALSE))
}

# Distribution of the permutation test replicates
hist(res, breaks = 20)

# Estimated probability of observing a more deviating value for the correlation under H0
mean(res >= abs(cor(rx, ry)))
# p-value = 0 -> Spearman correlation differs significantly from zero


# Conclusions:
# There is both linear and monotonic dependency in the data
# The monotonic dependency is a somewhat stronger than the linear dependency
# -> When we move to linear models, maybe some transformations should be done to achieve linearity


plot(x, y, xlab = "Ozone content", ylab = "Temperature")
plot(log(x), y, xlab = "log(Ozone content)", ylab = "Temperature")