#Read data set and check structure and names
a3 = read.csv("C:/Users/Ricky/Desktop/STA302/Assignment 3/PPRSA_data_2014.1.1-2014.12.31.csv",header=T) dim(a3)

str(a3)

class(a3)

names(a3)

#Define diagonal and lower panel of the scatter plot matrix
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1))
r <- cor(x, y)
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste("r= ", txt, sep = "") text(0.5, 0.6, txt)
p <- cor.test(x, y)$p.value
txt2 <- format(c(p, 0.123456789), digits = digits)[1] txt2 <- paste("p= ", txt2, sep = "")
if(p<0.01) txt2 <- paste("p= ", "<0.01", sep = "") text(0.5, 0.4, txt2)
}

#Buil Model 0 with only y
m0 <- lm(pm2.5 ~ 1, data = a3) summary(m0)

anova(m1)

pairs(pm2.5 ~ month + hour + DEWP + TEMP + PRES + Iws + Is + Ir, main = "Beijing PM2.5", data = a3)

#Select the best model according to AIC
mboth = step ( m0 ,scope = list ( lower = formula ( m0 ),upper = formula ( m1 )),direction = "both")

#Build Model 2 with log(y) transformation
m2 <- lm(log(pm2.5) ~ month + hour + DEWP + TEMP + PRES + Iws + Is + Ir, data = a3) summary(m2)

anova(m2)

#Select the best model according to AIC
Beijing PM2.5
mboth = step ( m00 ,scope = list ( lower = formula ( m00 ),upper = formula ( m2 )), direction = "both")

formula ( mboth )

summary ( mboth )

#Build Model 3 with sqrt(y)
m3 <- lm(sqrt(pm2.5) ~ month + hour + DEWP + TEMP + PRES + Iws + Is + Ir, data = a3) summary(m3)

#Plot Residual plot, Q-Q plot and Scale-location plot for Model 1, 2 & 3
par(mfrow=c(3,3))
plot(m1,main="Raw data")
plot(m2,which=1,main="Log(pm2.5)")
plot(m3,which=1,main="Sqrt(pm2.5)")
plot(m1,which=2,main="Raw data")
plot(m2,which=2,main="Log(pm2.5)")
plot(m3,which=2,main="Sqrt(pm2.5)")
plot(m1,which=3,main="Raw data")
plot(m2,which=3,main="Log(pm2.5)")
plot(m3,which=3,main="Sqrt(pm2.5)")

#Use box-cox method to find optimal lumbda
library(MASS) bc=boxcox(m1,lambda=seq(-2,2,by=0.01))
bc$x[which.max(bc$y)]

#Perform Log transformation on Iws
m4 <- lm(log(pm2.5) ~ month + hour + DEWP + TEMP + PRES + log(Iws) + Is + Ir, data = a3) summary(m4)
anova(m4)

#Model 5 with Log(Iws), excluding TEMP and PRES
m5 <- lm(log(pm2.5) ~ month + hour + DEWP + log(Iws) + Is + Ir, data = a3) summary(m5)

#compare AICs of all models
c(AIC(m1), AIC(m2), AIC(m3), AIC(m4), AIC(m5))

#Calculating Confidence interval for all models
confint(m1, level = 0.95)
confint(m2, level = 0.95)
confint(m3 ,level = 0.95)
confint(m4, level = 0.95)
confint(m5, level = 0.95)

#Taking prediction interval for all models at new x's
newdata = data.frame(month = 6, hour = 18, DEWP = -10, TEMP = 14, PRES = 1015, Iws = 50, I predict(m1, newdata, interval="predict")
predict(m2, newdata, interval="predict")
predict(m3, newdata, interval="predict")
predict(m4, newdata, interval="predict")
predict(m5, newdata, interval="predict")