library(car)
Prestige.complete<-na.omit(Prestige)
attach(Prestige.complete)
model<-lm(prestige~income+education+type)
plot(prestige~income)
qqPlot(model$res, main="QQ Plot for Residuals")
hist(model$res, main="Histogram for Residuals")
## Residuals against Regressors
st.res<-rstandard(model)
length(income)
length(st.res)
plot(income,st.res, main="Standardized Residuals and income",ylim=c(-2,2))
plot(education,st.res, main="Standardized Residuals and education",ylim=c(-2,2))
plot(type,st.res, main="Standardized Residuals and type",ylim=c(-2,2))
tapply(st.res,type,sd)
fitted<-model$fitted
plot(fitted,st.res, main="Standardized Residuals and Fitted Values",ylim=c(-2,2))

### calculate all the different residuals
(res<-model$res) # Residuals
(st.res<-rstandard(model)) # Standardized residuals
s<-summary(model)
(hii<-hatvalues(model))   # h_{ii}
(stud.i.res<-res/(s$sigma*sqrt(1-hii))) #internally studentized residuals
(press.res<-res/(1-hii)) # PRESS residuals
(stud.e.res<-rstudent(model))  # function for finding externally studentized residuals

# Create matrix of residuals/diagnostics
res.matrix<-cbind(res, st.res, hii, stud.i.res, press.res, stud.e.res)
res.matrix

# influence.measures(model)
# rstandard(model)
# rstudent(model) 
# hatvalues(model)


library(scatterplot3d)
(p.bc<-Prestige.complete[type=="bc",])
(p.wc<-Prestige.complete[type=="wc",])
(p.prof<-Prestige.complete[type=="prof",])
s3d<-scatterplot3d(p.bc$income,p.bc$education, p.bc$prestige, type="h",highlight.3d=F,angle=60, scale.y=0.7, pch=16,xlim=c(0,15000),ylim=c(6,15),zlim=c(10,100))
s3d$points3d(p.wc$income,p.wc$education, p.wc$prestige, type="h",col="blue",  pch=16)
s3d$points3d(p.prof$income,p.prof$education, p.prof$prestige, type="h",col="red",  pch=16)


##### Partial Residuals Example

model.prestige.income<-lm(prestige~income+type)
model.prestige.income$res
model.education.income<-lm(education~income+type)
model.education.income$res
plot(model.prestige.income$res~model.education.income$res, main="Partial Residual Plot for Prestige and Education")
identify(model.prestige.income$res~model.education.income$res)
abline(lm(model.prestige.income$res~model.education.income$res))

########### These models have the same residuals
# model.res.res<-lm(model.prestige.income$res~model.education.income$res-1)
# summary(model.res.res)
# model.res.res$res
# model.bp<-lm(prestige~income+education)
# model.bp$res
model.prestige.education<-lm(prestige~education)
model.income.education<-lm(income~education)

plot(model.prestige.education$res~model.income.education$res, main="Partial Residual Plot for Prestige and Education")