Example 6.2, Plutonium Radiation in Hanford

# Clear variables
rm(list=ls())

dataHanford <- read.table("hanford.txt", header = FALSE)
## Look at dataHanforda
exposure <- dataHanford[ , 2]
cancerMortality <- dataHanford[ , 3]
hanford <- data.frame(exposure, cancerMortality)
plot(hanford$exposure, hanford$cancerMortality, main = 'Plutonium Radiation in Hanford', xlab = 'Index of Exposure', ylab = 'Cancer Mortality')

y_i = beta_0 + beta_1 * x_i

## design matrix, and obs
y <- hanford$cancerMortality
x <- hanford$exposure
X <- cbind(1, x)  # m = 2
n <- length(y)  # n = 9
## correlation between estimates
cov2cor(t(X) %*% X)

                    x
  1.0000000 0.8143074
x 0.8143074 1.0000000

## remove the correlation ???
X.removeCorrelation <- X
X.removeCorrelation[ ,2] <- X[ ,2] - mean(X[ ,2])
cov2cor(t(X.removeCorrelation) %*% X.removeCorrelation)

                          x
  1.000000e+00 1.199279e-16
x 1.199279e-16 1.000000e+00

## Parameter estimaes
beta <- solve(t(X) %*% X) %*% t(X) %*% y
## 
regressionGLM <- glm(hanford$cancerMortality ~ hanford$exposure)
regressionLM <- lm(hanford$cancerMortality ~ hanford$exposure)

y_i = 157.333333 + 9.232264 * x_i

plot(hanford$exposure, hanford$cancerMortality, main = 'Linear Regression Model of Plutonium Radiation in Hanford', xlab = 'Index of Exposure', ylab = 'Cancer Mortality')
seq_x <- seq(1, 12)
seq_y <- beta[1] + beta[2] * seq_x
seq_y2 <- regressionGLM$coefficients[1] + regressionGLM$coefficients[2] * seq_x
lines(seq_x, seq_y, col = 4, lwd = 2)

## variance parameter
y.hat <- X %*% beta
sigmaSquare.hat <- sum((y - y.hat)^2) / n
## Unbiased Estimate, this is what we use
sigmaSquareUnbiased.hat <- sum((y - y.hat)^2) / (n - 2)
## Standard Error for beta_1 and beta_2
seBeta <- sqrt(diag(sigmaSquareUnbiased.hat * solve(t(X) %*% X)))
# Standard Error line and Confidence Interval Line
lineSE <- function(beta, seBeta, x, y){
    beta_2.se_1 <- beta[2] + seBeta[2]  # Standard error of beta_2
    beta_2.se_2 <- beta[2] - seBeta[2]
    lineSE_1 <- mean(y) - beta_2.se_1 * mean(x)
    lineSE_2 <- mean(y) - beta_2.se_2 * mean(x)
    abline(lineSE_1, beta_2.se_1, lty = 2, col = 'blue')
    abline(lineSE_2, beta_2.se_2, lty = 2, col = 'blue')
}
lineConf <- function(beta, sigmaSquareUnbiased.hat, x, y){
    n <- length(x)
    ssx <- sum(x^2) - sum(x^2) / n
    s.t <- qt(0.975, n-2)
    x.v <- seq(min(x), max(x), length = 100)
    y.v <- beta[1] + beta[2] * x.v
    se <- sqrt(sigmaSquareUnbiased.hat * (1 / n + (x.v - mean(x))^2 / ssx))  # Not standard error ?
    ci <- s.t * se
    y.v.upper <- y.v + ci
    y.v.lower <- y.v - ci
    lines(x.v, y.v.upper, lty = 3, col = 'blue')
    lines(x.v, y.v.lower, lty = 3, col = 'blue')
}
plot(hanford$exposure, hanford$cancerMortality, main = 'Linear Regression Model of Plutonium Radiation in Hanford', xlab = 'Index of Exposure', ylab = 'Cancer Mortality')
seq_x <- seq(1, 12)
seq_y <- beta[1] + beta[2] * seq_x
lines(seq_x, seq_y, col = 4, lwd = 2)

lineSE(beta, seBeta, x, y)
lineConf(beta, sigmaSquareUnbiased.hat, x, y)

## Remember to convert the name of the data, because the name in model must consist with those in the new data.
y <- hanford$cancerMortality
x <- hanford$exposure
regressionLM_2 <- lm(y ~ x)
pred.frame <- data.frame(x = seq(1, 30))
valuePred <- predict(regressionLM_2, pred.frame, se.fit = TRUE)
intervalPred <- predict(regressionLM_2, newdata = pred.frame, interval = "pred")
intervalConf <- predict(regressionLM_2, newdata = pred.frame, interval = "conf")
plot(pred.frame$x, valuePred$fit, pch = 8, ylim = c(min(y, valuePred$fit), max(y, valuePred$fit)), xlim = c(min(x, pred.frame$x), max(x, pred.frame$x)))
matlines(pred.frame$x, intervalPred, lty = c(1, 2, 2), col = c("black", "blue", "blue"), lwd = 2)

points(hanford$exposure, hanford$cancerMortality)
lines(pred.frame$x, intervalConf[, 2], lty = 3, col = "red", lwd = 2)

lines(pred.frame$x, intervalConf[, 3], lty = 3, col = "red", lwd = 2)

## Directly in R
summary(lm(y ~ X[ ,2]))

Call:
lm(formula = y ~ X[, 2])

Residuals:
    Min      1Q  Median      3Q     Max 
-16.283 -12.741   4.020   9.411  18.602 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  114.701      8.039  14.269 1.97e-06 ***
X[, 2]         9.232      1.418   6.513  0.00033 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14 on 7 degrees of freedom
Multiple R-squared:  0.8584,    Adjusted R-squared:  0.8381 
F-statistic: 42.42 on 1 and 7 DF,  p-value: 0.0003301

cbind(beta, seBeta)

               seBeta
  114.700789 8.038525
x   9.232264 1.417528

sqrt(sigmaSquareUnbiased.hat)

[1] 13.9975

## qqplot of residuals
library(car)
qqPlot(y - y.hat)

[1] 1 3

## Note the effect of the small dataHanfordaset...