## airquality data
Ozone <- airquality[complete.cases(airquality),"Ozone"]
Temp <- airquality[complete.cases(airquality),"Temp"]
plot(Temp, Ozone, pch=19)

## cor
cor(Ozone, Temp)

## cor.test
cor.test(Ozone, Temp)

## By hand
r <- cor(Ozone, Temp)
n <- length(Ozone)
ci.trans <- atanh(r) + qnorm(c(.025, .975))*sqrt(1/(n-3))
tanh(ci.trans)

## lm
fit <- lm(Ozone~Temp)
coef(fit)
confint(fit)
summary(fit)
predict(fit, newdata=data.frame(Temp=c(60, 70, 80, 90)))

## Plot the regression line
plot(Temp, Ozone, pch=19, col="gray")
abline(fit, col="slateblue", lwd=3)

## regression vs correlation
coef(fit)[2]*sd(Temp)/sd(Ozone) ## Same as correlation coefficient

## visreg
require(visreg)
visreg(fit)

## Log transform
fit <- lm(log(Ozone)~Temp)
visreg(fit)
visreg(fit, trans=exp, partial=TRUE, ylab="Ozone") ## Transforming back to original scale

## functions 1
cor.wald <- function(x, y, alpha=.05) {
  r <- cor(x,y)
  n <- length(x)
  SE <- (1-r^2)/sqrt(n)
  z <- qnorm(c(alpha/2, 1-alpha/2))
  r + z*SE
}
cor.score <- function(x, y, alpha=.05) {
  r <- cor(x,y)
  n <- length(x)
  f.u <- function(x) (r-x)/((1-x^2)/sqrt(n)) - qnorm(alpha/2)
  f.l <- function(x) (r-x)/((1-x^2)/sqrt(n)) - qnorm(1-alpha/2)
  c(uniroot(f.l, c(-1,1))$root, uniroot(f.u, c(-1,1))$root)
}

## functions 2
cor.lm <- function(x, y, alpha=0.05) {
  fit <- lm(y~x)
  ci <- confint(fit, 2, level=1-alpha)
  ci*sd(x)/sd(y)
}

## airquality comparison
cor.test(Ozone, Temp)$conf.int
cor.wald(Ozone, Temp)
cor.score(Ozone, Temp)
cor.lm(Ozone, Temp)

## functions 3
genData <- function(n, rho) {
  a <- sqrt(rho/(1-rho))
  U <- rnorm(n)
  V <- matrix(rnorm(n*2), n, 2)
  (V + a*U)/sqrt(1+a^2)
}

## simulation
source("http://myweb.uiowa.edu/pbreheny/571/f14/labs/displayProgressBar.R")
source("http://myweb.uiowa.edu/pbreheny/571/f14/labs/toplegend.R")
N <- 5000
n <- 10
Names <- c("Wald", "Score", "Fisher Z", "LM")
covered <- matrix(NA, N, 4)
rho <- c(0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)
coverage <- matrix(NA, length(rho), 4)
colnames(covered) <- colnames(coverage) <- Names
for (j in 1:length(rho)) {
  for (i in 1:N) {
    X <- genData(n, rho=rho[j])
    x <- X[,1]
    y <- X[,2]
    ci <- cor.wald(x,y)
    covered[i,1] <- (ci[1] < rho[j]) & (ci[2] > rho[j])
    ci <- cor.score(x,y)
    covered[i,2] <- (ci[1] < rho[j]) & (ci[2] > rho[j])
    ci <- cor.test(x,y)$conf.int
    covered[i,3] <- (ci[1] < rho[j]) & (ci[2] > rho[j])
    ci <- cor.lm(x,y)
    covered[i,4] <- (ci[1] < rho[j]) & (ci[2] > rho[j])
  }
  coverage[j,] <- apply(covered, 2, mean)
  displayProgressBar(j, length(rho))
}
col <- hcl(seq(15, 375, len=5), l=60, c=150)[1:4]
matplot(rho, coverage, type="l", lty=1, lwd=3, col=col, xlab=expression(rho), las=1)
toplegend(legend=Names, lwd=3, col=col)