Untitled

rm(list=ls())
library(ggplot2)
library(patchwork)
weather <- read.csv("Data/weather.csv")


(Y.model <- lm(rain ~pressure, data = weather))

p1aY = ggplot(data = weather, aes(x = pressure, y = rain)) + geom_point(size = 1) +
  geom_line(aes(y = Y.model$coefficients[1] + pressure*Y.model$coefficients[2]), linetype = "dashed") +
  labs(title = "Precipitation (mm) against Air Pressure (hPa)")
# Exponential decay

press <- data.frame(weather[4])
Y.pred <-
  cbind(press,
        fit = predict(Y.model, press),
        conf = predict(Y.model, press, interval = "confidence"),
        pred = predict(Y.model, press, interval = "prediction"))
head(Y.pred)
# get rid of the extra fits
Y.pred$conf.fit <- Y.pred$pred.fit <- NULL
head(Y.pred)

Y.pred$e <- Y.model$residuals
head(Y.pred)
(max.e <- max(abs(Y.pred$e)))
(Y.elims <- c(-max.e, max.e))
## Plot against yhat####
# Add a horizontal line at y=0,
# and expand the y-axis to include +/- max residual.

p1aYres = ggplot(data = Y.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Y.elims) +
  xlab(" ") +
  ylab(" ") +
  labs(tag = "identity") +
  labs(title = "Residuals vs Predicted Precipitation") +
  theme(text = element_text(size = 18))


 ggplot(data = weather, aes(x = pressure, y = rain)) + geom_point(size = 1) +
  geom_line(aes(y = Y.pred$fit, col = "fit"), linetype = "dashed") +
  geom_line(aes(y = Y.pred$conf.lwr, col = "conf"), linetype = "dashed") +
  geom_line(aes(y = Y.pred$conf.upr, col = "conf"), linetype = "dashed") +
  geom_line(aes(y = Y.pred$pred.lwr, col = "pred"), linetype = "dashed") +
  geom_line(aes(y = Y.pred$pred.upr, col = "pred"), linetype = "dashed") +
  labs(title = "Precipitation (mm) against Air Pressure (hPa)")



(Ylog.model <- lm(log(rain) ~ pressure, data = weather))

 p1alogY = ggplot(data = weather, aes(x = pressure, y = log(rain))) + geom_point(size = 1) +
  geom_line(aes(y = Ylog.model$coefficients[1] + pressure*Ylog.model$coefficients[2]), linetype = "dashed") +
  labs(title = "Precipitation (log(mm)) against Air Pressure (hPa)")
# Linear to the left, but after 1010 hPa the data seems to have another slope
p1alogY

## resid anal ###
Ylog.pred <-
  cbind(press,
        fit = predict(Ylog.model, press),
        conf = predict(Ylog.model, press, interval = "confidence"),
        pred = predict(Ylog.model, press, interval = "prediction"))
head(Y.pred)
# get rid of the extra fits
Ylog.pred$conf.fit <- Ylog.pred$pred.fit <- NULL
head(Ylog.pred)

Ylog.pred$e <- Ylog.model$residuals
head(Ylog.pred)
(max.e <- max(abs(Ylog.pred$e)))
(Ylog.elims <- c(-max.e, max.e))
## Plot against yhat####
# Add a horizontal line at y=0,
# and expand the y-axis to include +/- max residual.

p1alogYres = ggplot(data = Ylog.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Ylog.elims) +
  xlab(" ") +
  ylab("Residual") +
  labs(tag = "log") +
  #labs(title = "Residuals vs Predicted log-Precipitation") +
  theme(text = element_text(size = 18))


p1alogYints = ggplot(data = weather, aes(x = pressure, y = rain)) + geom_point(size = 1) +
  geom_line(aes(y = exp(Ylog.pred$fit), col = "fit"), linetype = "dashed") +
  geom_line(aes(y = exp(Ylog.pred$conf.lwr), col = "conf"), linetype = "dashed") +
  geom_line(aes(y = exp(Ylog.pred$conf.upr), col = "conf"), linetype = "dashed") +
  geom_line(aes(y = exp(Ylog.pred$pred.lwr), col = "pred"), linetype = "dashed") +
  geom_line(aes(y = exp(Ylog.pred$pred.upr), col = "pred"), linetype = "dashed") +
  labs(title = "Precipitation (mm) against Air Pressure (hPa)")
p1alogYints


(Ycr.model <- lm(rain^(1/3) ~ pressure, data = weather))

p1aYcr = ggplot(data = weather, aes(x = pressure, y = rain^(1/3))) + geom_point(size = 1) +
  geom_line(aes(y = Ycr.model$coefficients[1] +pressure*Ycr.model$coefficients[2]), linetype = "dashed") +
  labs(title = "Precipitation ((mm)^1/3) Against Air Pressure (hPa)")
# Linear, fits our model Y = xB + e
p1aYcr

## resid anal ### 
Ycr.pred <-
  cbind(press,
        fit = predict(Ycr.model, press),
        conf = predict(Ycr.model, press, interval = "confidence"),
        pred = predict(Ycr.model, press, interval = "prediction"))
head(Y.pred)
# get rid of the extra fits
Ycr.pred$conf.fit <- Ycr.pred$pred.fit <- NULL
head(Ycr.pred)

Ycr.pred$e <- Ycr.model$residuals
head(Ycr.pred)
(max.e <- max(abs(Ycr.pred$e)))
(Ycr.elims <- c(-max.e, max.e))
## Plot against yhat####
# Add a horizontal line at y=0,
# and expand the y-axis to include +/- max residual.

p1aYcrres = ggplot(data = Ycr.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Ycr.elims) +
  xlab("Predicted Precipitation") +
  ylab(" ") +
  labs(tag = "cube root") +
  #labs(title = "Residuals vs Predicted Precipitation") +
  theme(text = element_text(size = 18))


  p1aYcrints = ggplot(data = weather, aes(x = pressure, y = rain)) + geom_point(size = 1) +
  geom_line(aes(y = (Ycr.pred$fit)^3, col = "fit"), linetype = "dashed") +
  geom_line(aes(y = (Ycr.pred$conf.lwr)^3, col = "conf"), linetype = "dashed") +
  geom_line(aes(y = (Ycr.pred$conf.upr)^3, col = "conf"), linetype = "dashed") +
  geom_line(aes(y = (Ycr.pred$pred.lwr)^3, col = "pred"), linetype = "dashed") +
  geom_line(aes(y = (Ycr.pred$pred.upr)^3, col = "pred"), linetype = "dashed") +
  labs(title = "Precipitation (mm) against Air Pressure (hPa)")
  
  # FINAL RESIDUAL PLOTS
  p1aYres/p1alogYres/p1aYcrres
  
################## Part 2: #####################