#
# File name: lec08.R
#
# data example: patient satisfaction
# KNNL Problem 6.15 and 7.5
# read in data
# study of patient satisfaction by a hospital administrator
# y = patient satisfaction
# x1 = age (in years)
# x2 = severity of illness (an index)
# x3 = anxiety level (an index)
#
patient = read.table("CH06PR15.txt")
names(patient) = c("y","x1","x2","x3")
patient.lm123 = lm(y ~ x1 + x2 + x3, data = patient)
summary(patient.lm123)
# 7.5(a) obtain ANOVA decomposition
# Note: not in the right order!
anova(patient.lm123)
# Note: in the correct order!
patient.lm213 = lm(y ~ x2 + x1 + x3, data = patient)
anova(patient.lm213)
# 7.5(b) Test whether X3 can be dropped from the model
# H_0: beta_3 = 0 versus H_a: beta_3 != 0
# Test Statistic: F* = [SSR(X3|X1,X2)/1]/[SSE(X1,X2,X3)/(n-4)]
# F* = 364.2/101.2 = 3.5997
# Note: could get this from the MS as follows:
b.mses = anova(patient.lm213)$Mean
b.mses
b.mses[3]/b.mses[4]
# Reject H_0 if F* > F(1-0.025,1,42)
qf(0.975,1,42)
# Or use p-value
1-pf(3.599735,1,42)
# Conclusion: fail to reject H_0 at the 0.025 level of significance.
# Not enough evidence to reject H_0.
# 7.6. H_0: beta_2 = beta_3 = 0 versus H_a: at least one beta is not zero
patient.lm123 = lm(y ~ x1 + x2 + x3, data = patient)
patient.lm1 = lm(y ~ x1, data = patient)
anova(patient.lm1, patient.lm123)
# Reject H_0 if F* = 4.1768 > F(1-0.025,2,42)
qf(.975,2,42)
# Or use p-value
1-pf(4.1768,2,42)
# [1] 0.02216124
# Conclusion: reject H_0 at the 0.025 level of significance.
# There is moderate evidence to claim that at least one of the parameters is not zero.
#############################
# VIF and Multicollinearity #
#############################
pairs(patient)
round(cor(patient),2)
# VIF by Prof. Jennifer Hoeting
# The "VIF" function below computes the VIF values for a model fit using the lm function.
# To use this function, copy the entire function into R. If you
# save your workspace, you'll have access to the function in future sessions.
# If you don't save your workspace at the end of your session, you'll need to
# re-evaluate this function code before using it in future R sessions.
VIF = function(lm.object) {
## lm.object is the linear model object
d = dim(model.matrix(lm.object))[2]
X = model.matrix(lm.object)[,2:d]
a = 1/sqrt(dim(X)[1]-1)*scale(X)
b = cbind(diag(solve(t(a)%*%a)))
dimnames(b) = list(dimnames(X)[[2]],"VIF")
return(b)
}
VIF(patient.lm123)
# Or simply, use the vif function in HH (or some other library)
library(HH)
vif(patient.lm123)