library(car)

save_figs = FALSE

source('../../Data/data_loaders.R')
DF = load_appendix_cement_data()

pairs(DF)
scatterplotMatrix(DF)

# Form the standardized predictors:
#
DFs = scale(DF)
DFs = data.frame(DFs)

m = lm( Y ~ ., data=DF )
print(summary(m))

# Part (a):
#
print('Variance inflation factors:')
print(vif(m))

XTX = cor(DF)
print('Correlation matrix:')
print(XTX)

ES = eigen(XTX[ 1:4, 1:4 ])
print('Eigenvalues:')
print(ES$values)
ES$vectors[, 1] = -ES$vectors[, 1] # change the sign of all components of the first eigenvector
print('Eigenvectors:')
print(ES$vectors)

# Part (b):
#
# Compute the primcipal components (and a scatter plot of them):
#
Z = as.matrix(DFs[ , 1:4 ]) %*% ES$vectors
colnames(Z) = c('Z1', 'Z2', 'Z3', 'Z4')
print(summary(Z))

if( save_figs ){ postscript("../../WriteUp/Graphics/Chapter5/ex_5_6_scatter_plot_of_Z.eps", onefile=FALSE, horizontal=FALSE) }
pairs(Z, xlim=c(-3, +2), ylim=c(-3, +2))
if( save_figs ){ dev.off() }


# Compute the coefficients of the linear constraint:
#
ss = scale(DF)
smeans = attr(ss, "scaled:center")
sstd = attr(ss, "scaled:scale")
lev = ES$vectors[, 4] # coefficients of the last eigenvector

coefficients_of_x_i = lev / sstd[1:4]
rhs = sum( smeans[1:4] / sstd[1:4] )
print('coefficients_of_x_i:')
print(coefficients_of_x_i)
print('rhs:')
print(rhs)
lhs_mean = mean(coefficients_of_x_i)
print(lhs_mean)
print('new coefficients_of_x_i:')
print(round(coefficients_of_x_i/lhs_mean, 3))
print('new rhs:')
print(round(rhs/lhs_mean, 3))


# Part (c):
#
DF$X5 = 100 - rowSums(DF[, 1:4 ])
DF = DF[, c( 1:4, 6, 5 ) ] # make the response Y the last column

m = lm( Y ~ . - 1, data=DF )
print('Model with X5 and no intercept:')
print(summary(m))

m = lm( Y ~ ., data=DF )
print('Model with X5 and an intercept:')
print(summary(m))


# Part (d):
#
DF = DF[, c(1, 2, 3, 5, 6) ]

m = lm( Y ~ ., data=DF )
print(summary(m))

print('Variance Inflation Factors:')
print(vif(m))

# Form the standardized predictors:
#
DFs = scale(DF)
DFs = data.frame(DFs)

XTX = cor(DF)
print('Correlation matrix:')
print(XTX)

ES = eigen(XTX[ 1:4, 1:4 ])
print('Eigenvalues:')
print(ES$values)