library(lattice)
library(MASS)

# Save pictures in the current working directory
myDir <- './'

# Load data Fig 7-1
data1<-read.table('Fig7-1.txt')

# Run regression
z1<-lm(V2~V1,data = data1)
summary(z1)

# Plot X vs Y with regression line and 'residuals', and save the picture.
png(file = paste(myDir,'Fig7-1_R.png', sep ='/'), width = 500, height = 500)
xyplot(V2 ~ V1, data = data1, xlab = 'Stretch', ylab = 'Strength',
       panel = function(x, y, ...) {
         panel.grid(h = -1, v = -1)
         panel.xyplot(x, y, pch = 16, cex = 1.25)
         panel.lmline(x, y)
         z1 <- lm(V2 ~ V1, data = data1)
         ind = order(x)[1:3]
         panel.segments(x, y, x, z1$fitted.values,
                        col = 'cyan4')
       }
       )
graphics.off()


# Load data Fig 7-2
data2<-read.table('Fig7-2.txt')

# Run regression
z2<-lm(V2~V1,data = data2)
summary(z2)

# Plot X vs Y with regression line and 'residuals'.
png(file = paste(myDir,'Fig7-2_R.png', sep ='/'), width = 500, height = 500)
xyplot(V2 ~ V1, data = data2, xlab = 'X', ylab = 'Y',
       panel = function(x, y, ...) {
         panel.grid(h = -1, v = -1)
         panel.xyplot(x, y, pch = 16, cex = 1.25)
         panel.lmline(x, y)
         z1 <- lm(V2 ~ V1, data = data2)
         ind = order(x)[1:3]
         panel.segments(x, y, x, z1$fitted.values,
                        col = 'cyan4')
       }
       )
graphics.off()

# Regress 1/Y on 1/X
#z21 <- lm(1/V2 ~ 1/V1, data = data2)
#summary(z21)
Y<-1/data2$V2;
X<-1/data2$V1;
z21<-lm(Y~X);
summary(z21)

# Plot X vs Y with regression line and 'residuals'.
png(file = paste(myDir,'Fig7-2_R+.png', sep ='/'), width = 500, height = 500)
xyplot(Y ~ X, xlab = '1/X', ylab = '1/Y',
       panel = function(x, y, ...) {
         panel.grid(h = -1, v = -1)
         panel.xyplot(x, y, pch = 16, cex = 1.25)
         panel.lmline(x, y)
         z1 <- lm(Y ~ X)
         ind = order(x)[1:3]
         panel.segments(x, y, x, z1$fitted.values,
                        col = 'cyan4')
       }
       )
graphics.off()


# Load data Fig 7-4
data3<-read.table('Fig7-4.txt')

# Run regression
z3<-lm(V2~V1,data = data3)
summary(z3)


# Plot Time vs Carbon (Fig7-4) with regression line and 'residuals'.
png(file = paste(myDir,'Fig7-4_R.png', sep ='/'), width = 500, height = 500)
xyplot(V2 ~ V1, data = data3, xlab = 'Carbon', ylab = 'Time',
       panel = function(x, y, ...) {
         panel.grid(h = -1, v = -1)
         panel.xyplot(x, y, pch = 16, cex = 1.25)
         panel.lmline(x, y)
         z1 <- lm(V2 ~ V1, data = data3)
         ind = order(x)[1:3]
         panel.segments(x, y, x, z1$fitted.values,
                        col = 'cyan4')
       }
       )
graphics.off()

# Plot residuals vs fitted.
png(file = paste(myDir,'Fig7-4_residues.png', sep ='/'),width = 500, height = 500)
xyplot(residuals ~ fitted.values, data = z3, xlab = 'FITTED', ylab = 'RESIDUALS',
       pch = 16, cex = 1.25,
       panel = function(x, y, ...) {
         panel.grid(h = -1, v = -1)
         panel.xyplot(x, y, ...)
         panel.abline(h = 0)
       }
       )
graphics.off()

# Check normality of residuals.
png(file = paste(myDir,'Fig7-4_normality.png', sep ='/'),
    width = 350, height = 650)
p1 <- histogram(~residuals, data = z3,
                xlab = 'RESIDUALS', ylab = 'PERCENT OF TOTAL')
p2 <- qqmath(~residuals, data = z3, distribution = qnorm,
       pch = 16, cex = 1.25, xlab = 'NORMAL QUANTILES', ylab = 'RESIDUALS',
       panel = function(x, ...) {
         panel.qqmathline(x, ...)
         panel.qqmath(x, ...)
         panel.grid()
       }
       )
print(p1, position = c(0, .49, 1, 1), more = TRUE)
print(p2, position = c(0, 0, 1, .51))
graphics.off()

# Predict the mean Y for X=145 and confidence intervals.
ciPred <- predict(z3,interval = 'prediction',
	newdata=data.frame(V1=c(145)))

# Get covariance matrix of regression coefficients.
z3Summary <- summary(z3)
z3Correlation <- z3Summary$sigma^2 * z3Summary$cov.unscaled
# Simulate 30 regression lines from the approximate distrib of the coefficients.
simEstimates <- mvrnorm(30, z3$coefficients, z3Correlation)

# Plot the simulated regression lines.
png(file = paste(myDir,'Fig7-4_Rlines.png', sep ='/'),
    width = 350, height = 350)
xyplot(V2 ~ V1, data = data3, xlab = 'Carbon', ylab = 'Time',
       panel = function(x, y, ...) {
         panel.grid()
         panel.xyplot(x, y, cex = 1.25, pch = 16)
          for (i in 1:nrow(simEstimates)) {
           panel.abline(simEstimates[i, ], col = gray(.9))
         }
         panel.abline(z3$coefficients, col = 'red', lwd = 1.5)
 
       }
       )
graphics.off()