R in action -- chapter 8

rm(list = ls())
fit <- lm(weight ~ height, women)
summary(fit)
women$weight
fitted(fit)
confint(fit)
residuals(fit)
with(women,plot(height,weight,
                xlab = "height (in inches)",
                ylab = "weight (in pounds)"))
abline(fit)

fit2 <- lm(weight ~ height + I(height^2), women)
summary(fit2)
with(women,plot(height,weight,
                xlab = "height (in inches)",
                ylab = "weight (in pounds)"))
lines(women$height,fitted(fit2))

library(car)
scatterplot(weight ~ height, data = women,
            spread = FALSE, smoother.agrs=list(lty=2),pch = 20,
            main = "women age 30 - 39",
            xlab = "height (inches)",
            ylab = "weight (lbs.)")

states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
cor(states)
library(car)
scatterplotMatrix(states)

fit <- lm(Murder~Population+Illiteracy+Income+Frost,states)
summary(fit)
coefficients(fit)
confint(fit)

fit <- lm(mpg~hp+wt+hp:wt,mtcars)
summary(fit)

# interpretation about interaction
library(effects)
plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=T)

fit <- lm(weight ~ height, women)
par(mfrow=c(2,2))
plot(fit)


fit <- lm(weight ~ height + I(height^2), women)
par(mfrow=c(2,2))
plot(fit)

newfit  <- lm(weight~ height + I(height^2), data= women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)

states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,states)
par(mfrow=c(2,2))
plot(fit)

par(mfrow=c(1,1))
qqPlot(fit,lables=row.names(states),id.method= "identify",
       simulate=T, main = "Q-Q plot")

states["Nevada",]
fitted(fit)["Nevada"]
residuals(fit)['Nevada']
rstudent(fit)["Nevada"]


residualplot <- function(fit, nbreaks= 10){
  z <- rstudent(fit)
  hist(z, breaks= nbreaks, freq = F, 
       xlab = "Studentized residual",
       main= "Distribution of errors")
  rug(jitter(z),col="brown")
  curve(dnorm(x,mean=mean(z),sd=sd(z)),
        add=T, col="blue",lwd= 2)
  lines(density
        (z)$x,density(z)$y,
        col= "red",lwd= 2,lty= 2)
  legend("topright",
         legend = c("normal curve","kernel density curve"),
         lty= 1:2,col=c("blue","red"),cex=.7)
}
residualplot(fit)

library(car)
durbinWatsonTest(fit)

crPlots(fit)

ncvTest(fit)

spreadLevelPlot(fit)

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
vif(fit)

library(car)
outlierTest(fit)
n <- 50
hat.plot <- function(fit){
  p <- length(coefficients(fit))
  n <- length(fitted(fit))
  plot(hatvalues(fit),main= "Index Plot of hat values")
  abline(h=c(2,3)*p/n,col= "red", lty= 2)
  identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit,which = 4,cook.levels = cutoff)
abline(h=cutoff, lty=2,col="red")

library(car)
avPlots(fit,ask= FALSE, id.method ='identify')
ls("package:car")

influencePlot(fit,id.method="identify", main="Influence plot",
              sub="circle size is proportional to cook's distance")

states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit1 <- lm(Murder~Population+Illiteracy+Income+Frost,states)
fit2 <- lm(Murder~ Population+ Illiteracy +Income , data= states)
anova(fit2,fit1)
AIC(fit1,fit2)

library(MASS)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,states)
stepAIC(fit,direction = "backward")

library(leaps)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
leaps <- regsubsets(Murder ~ Population + Illiteracy +Income + Frost, states,nbest = 4)
plot(leaps, scale = 'adjr2')

library(car)
subsets(leaps,statistic = "cp",
        main = "Cp plot for all subsets regression")
abline(1,1,lty=2,col="red")