Untitled

mail@pastecode.io avatar
unknown
r
2 years ago
5.2 kB
1
Indexable
Never
library(ggplot2)
weather <- read.csv("Data/weather.csv")


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

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

#INTE FÄRDIGT, cbind(BLANK) ska va ngt annat, körde Y.data, men existerar ej.
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.

ggplot(data = Y.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Y.elims) +
  xlab("Predicted weight loss (g)") +
  ylab("Residual") +
  labs(tag = "B") +
  labs(title = "Residuals vs predicted values Y-hat") +
  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), linetype = "dashed") +
  geom_line(aes(y = Y.pred$conf.lwr), linetype = "dashed") +
  geom_line(aes(y = Y.pred$conf.upr), linetype = "dashed") +
  geom_line(aes(y = Y.pred$pred.lwr), linetype = "dashed") +
  geom_line(aes(y = Y.pred$pred.upr), linetype = "dashed") +
  labs(title = "Precipitation (mm) against Air Pressure (hPa)")



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

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



## 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.

ggplot(data = Ylog.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Ylog.elims) +
  xlab("Predicted weight loss (g)") +
  ylab("Residual") +
  labs(tag = "B") +
  labs(title = "Residuals vs predicted values Y-hat") +
  theme(text = element_text(size = 18))


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



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

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


## 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.

ggplot(data = Ycr.pred, aes(x = fit, y = e)) +
  geom_point(size = 1) +
  geom_hline(yintercept = 0) +
  expand_limits(y = Ycr.elims) +
  xlab("Predicted weight loss (g)") +
  ylab("Residual") +
  labs(tag = "B") +
  labs(title = "Residuals vs predicted values Y-hat") +
  theme(text = element_text(size = 18))

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