1. Introduction

knitr::include_graphics(paste("~/Desktop/school/pstat_131/pics",
                        "kangaroo.png",sep="/"))
## Warning in knitr::include_graphics(paste("~/Desktop/school/pstat_131/pics", : It
## is highly recommended to use relative paths for images. You had absolute paths:
## "~/Desktop/school/pstat_131/pics/kangaroo.png"

Background

During the time period from 2008 to 2018, Australia experienced a wide range of rainfall patterns and weather events. La Niña and El Niño events occurred which impacted rainfall patterns throughout the country. Several regions in Australia faced prolonged periods of drought which led to significant water shortages, agricultural losses, and overall negative impacts on the environment. In 2010-2011 and 2012-2013, significant flooding occurred in the northeastern part of Australia, which resulted in extensive damages to infrastructure, homes, and agricultural land. Generally, rainfall patterns vary across different regions. The northern regions experience a wet season during the summer months, whereas the southern parts have a Mediterranean climate with winter rainfall. Typically in Australia, the Spring season occurs during the months September, October, and November and the Summer season includes the months December, January, and February. Additionally, the Autumn season is comprised of the months March, April, and May where as the Winter season falls within the months June, July, and August.

Climate change has been a contributing factor to shifts in rainfall patterns all around the world, including Australia. Overall, the time from 2008 to 2018 in Australia had a mix of drought, flooding, and regional variations in rainfall patterns. These events had a significant impact on agriculture, water resources, and the overall environment, emphasizing how vulnerable Australia is to changing weather patterns.

Motive

This project really appealed to me because my parents first moved to Australia right after they got married! Also, I have family in Australia and New Zealand, so I hope to visit both places one day.

Data Description

This dataset contains 10 years of daily weather observations from many locations across Australia. The data collected were drawn from numerous weather stations and the daily observations were drawn from http://www.bom.gov.au/climate/data and http://www.bom.gov.au/climate/dwo/ which was available on Kaggle https://www.kaggle.com/datasets/jsphyg/weather-dataset-rattle-package.

Project Outline

With the background information provided, we can start to look at our goals for this project. Firstly, we will be exploring the data by doing some tidying and getting rid of predictors with the most missing observations and any unnecessary predictors that might not provide much insight. The remaining variables can be used for additional exploration to gain a better understanding of how they can be used to predict the next-day rain. This will mainly include looking at some visualizations of the different predictor variables and the target variable, RainTomorrow. The main goal here is to use the target variable, RainTomorrow, to predict the next-day rain. If it rains the next day for 1mm or more, then the observation for that column is Yes, if not, then the observation for that column is No. Next, we will split the data to create a testing set and a training set, make a recipe, and set folds for the 10-fold cross validation method that we will apply. The different models we will implement to train the data are Logistic Regression, Quadratic Discriminant Analysis, Support Vector Machin (SVM), and Random Forest. After obtaining the results from each model, the one that performed the best will be selected and fit to the test data set to determine how effective the model really is.

2. Exploring the Data

Loading and Exploring the Data

The variable RainTomorrow is the target variable we are trying to predict. If it rains the next day for 1mm or more then the observation for that column is Yes, if not, then the observation for that column is No. There are 145,460 observations and 23 predictors in the original dataset. Also, approximately 10.1% of the observations are missing.

weatherAUS <- read_csv("~/Desktop/school/pstat_131/data/unprocessed/weatherAUS.csv")
## Rows: 145460 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (6): Location, WindGustDir, WindDir9am, WindDir3pm, RainToday, RainTom...
## dbl  (16): MinTemp, MaxTemp, Rainfall, Evaporation, Sunshine, WindGustSpeed,...
## date  (1): Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(weatherAUS)
## # A tibble: 6 × 23
##   Date       Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
##   <date>     <chr>      <dbl>   <dbl>    <dbl>       <dbl>    <dbl> <chr>      
## 1 2008-12-01 Albury      13.4    22.9      0.6          NA       NA W          
## 2 2008-12-02 Albury       7.4    25.1      0            NA       NA WNW        
## 3 2008-12-03 Albury      12.9    25.7      0            NA       NA WSW        
## 4 2008-12-04 Albury       9.2    28        0            NA       NA NE         
## 5 2008-12-05 Albury      17.5    32.3      1            NA       NA W          
## 6 2008-12-06 Albury      14.6    29.7      0.2          NA       NA WNW        
## # … with 15 more variables: WindGustSpeed <dbl>, WindDir9am <chr>,
## #   WindDir3pm <chr>, WindSpeed9am <dbl>, WindSpeed3pm <dbl>,
## #   Humidity9am <dbl>, Humidity3pm <dbl>, Pressure9am <dbl>, Pressure3pm <dbl>,
## #   Cloud9am <dbl>, Cloud3pm <dbl>, Temp9am <dbl>, Temp3pm <dbl>,
## #   RainToday <chr>, RainTomorrow <chr>
library(naniar)
sample_n(weatherAUS, size = 5000) %>% 
  vis_miss()

Data Tidying

Here, we have dropped all the variables with too much missing data. Additionally, we have also changed the binary variable RainTomorrow to be a factor, and reordered it so that Yes is the first level and No is the second level.

drop <- c("Evaporation","Sunshine","Cloud9am","Cloud3pm")
weatherAUS <- weatherAUS[,!(names(weatherAUS)%in%drop)]

# Make RainTomorrow a factor
weatherAUS$RainTomorrow <- factor(weatherAUS$RainTomorrow,levels = c("Yes","No"))

One of the variables in the dataset, date, which notes the month, day, and year that it has rained over these 10 years, from many locations all across Australia, is in the YYYY-MM-DD format. This makes it difficult to use for prediction in a classification model, so here we have extracted the month and year from each value and made them their own predictor. The specific day does not matter as much, but the month will be more beneficial to include since the rain varies across the different seasons over the years. Generally, in the northern part of Australia, majority of the rainfall occurs during the northern wet season from October to April. In the southern regions, rainfall mainly occurs during the southern wet season from April to November. So, the rainfall observed in March 2016 might have a significant difference when compared to the rainfall observed in December 2008.

dates <- as.POSIXct(weatherAUS$Date,format="%Y-%m-%d")

weatherAUS <- weatherAUS %>%
  add_column("Year"=c(format(dates,format="%Y")), .before = 1) %>%
  add_column("Month"=c(format(dates,format="%m")), .before=1)

After adding these two columns, we will now convert them to numeric variables so we can use them in our recipe later on without having to dummy code them.

weatherAUS$Month <- as.numeric(as.character(weatherAUS$Month))
weatherAUS$Year <- as.numeric(as.character(weatherAUS$Year))

So now, we have removed the original Date variable since we have added our new Month and Year columns. Additionally, we have taken out all the missing observations overall. Lastly, we have displayed the current amount of predictors and observations after cleaning our dataset, giving us a little sneak peak of the remaining data we will be working with for the rest of this project.

drop1 <- c("Date")
weatherAUS <- weatherAUS[,!(names(weatherAUS)%in%drop1)]
weatherAUS <- na.omit(weatherAUS)
head(weatherAUS)
## # A tibble: 6 × 20
##   Month  Year Location MinTemp MaxTemp Rainfall WindGustDir WindGustSpeed
##   <dbl> <dbl> <chr>      <dbl>   <dbl>    <dbl> <chr>               <dbl>
## 1    11  2008 Albury      13.4    22.9      0.6 W                      44
## 2    12  2008 Albury       7.4    25.1      0   WNW                    44
## 3    12  2008 Albury      12.9    25.7      0   WSW                    46
## 4    12  2008 Albury       9.2    28        0   NE                     24
## 5    12  2008 Albury      17.5    32.3      1   W                      41
## 6    12  2008 Albury      14.6    29.7      0.2 WNW                    56
## # … with 12 more variables: WindDir9am <chr>, WindDir3pm <chr>,
## #   WindSpeed9am <dbl>, WindSpeed3pm <dbl>, Humidity9am <dbl>,
## #   Humidity3pm <dbl>, Pressure9am <dbl>, Pressure3pm <dbl>, Temp9am <dbl>,
## #   Temp3pm <dbl>, RainToday <chr>, RainTomorrow <fct>
cat("The dimensions for the final set we will be working with are",weatherAUS %>% dim(),"observations and predictors, respectively.")
## The dimensions for the final set we will be working with are 112925 20 observations and predictors, respectively.

Description of Predictors

Now that we have gotten rid of any unnecessary predictors in the data set, we should be able to gain a better understanding of the role that each predictor plays in regards to the response variable, RainTomorrow. The predictors we will be using for the remainder of the project are the following:

  • Month: The month that the rainfall observation was recorded (ranging from Jan-Dec) as represented by its respective number (e.g. 12 for December)
  • Year: The year that the rainfall observation was recorded (ranging from 2007-2017)
  • Location: The common name of the location of the weather station when the observation was recorded
  • MinTemp: The minimum temperature in degrees celsius on the day the observation was recorded
  • MaxTemp: The maximum temperature in degrees celsius on the day the observation was recorded
  • Rainfall: The amount of rainfall recorded for the day in mm
  • WindGustDir: The direction of the strongest wind gust in the 24 hours to midnight
  • WindGustSpeed: The speed in km per hour of the strongest wind gust in the 24 hours to midnight
  • WindDir9am: Direction of the wind at 9am on the day the observation was recorded
  • WindDir3pm: Direction of the wind at 3pm on the day the observation was recorded
  • WindSpeed9am: Wind speed in km per hour averaged over 10 minutes prior to 9am the day of the observation
  • WindSpeed3pm: Wind speed in km per hour averaged over 10 minutes prior to 3pm the day of the observation
  • Humidity9am: Humidity (percent) at 9am on the day of the recorded observation
  • Humidity3pm: Humidity (percent) at 3pm on the day of the recorded observation
  • Pressure9am: Atmospheric pressure (hpa) reduced to mean sea level at 9am on the day the observation was recorded
  • Pressure3pm: Atmospheric pressure (hpa) reduced to mean sea level at 3pm on the day the observation was recorded
  • Temp9am: Temperature in degrees celsius at 9am on the day the observation was recorded
  • Temp3pm: Temperature in degrees celsius at 3pm on the day the observation was recorded
  • RainToday: A boolean variable where Yes represents the precipitation (mm) in the 24 hours to 9am exceeding 1mm, otherwise the observation for that day is No

Visual Exploratory Data Analysis

In order to gain a better understanding of our response variable RainTomorrow, we have created an output variable plot and a correlation matrix to distinguish any possible relationships between the predictor variables. We have also generated additional visualization plots to examine the effect of the predictor variables on the response variable.

RainTomorrow Distribution

Below, we have created a plot of the distribution of whether or not rainfall of 1mm or more was observed over the course of 10 years. To reiterate, a Yes means that rainfall of at least 1mm occurred on a certain day.

weatherAUS %>%
  ggplot(aes(x=RainTomorrow)) +
  geom_bar() +
  labs(x="Rainfall Of 1mm Or More Observed",y="Amount of Rainfall Observations")

From the plot above, we can observe that there are over 75,000 cases of there not being next-day rain of 1mm or more. However, there are approximately 25,000 instances of the next-day rain of 1mm or more being recorded.

Predictor Correlation Plot

In order to see the relationship between our numeric variables, we created a correlation matrix and then generated a heat map of the correlation of these predictors.

weatherAUS %>%
  select_if(is.numeric) %>%
  cor() %>%
  corrplot()

I definitely expected to see significant correlation observed among the predictors because they contain quite a bit of rainfall data that is available for Australia, so they are highly likely to be correlated. One observation that stood out to me was the bigger blue square in the middle. The blue squares represent a strong, positive correlation between the variables. This middle group comprises of the wind speed variables that are averaged over 10 minutes prior to their respective times. The groups of variables that also have a blue square correlation are the minimum and maximum temperatures and the temperatures at 9am and 3pm, as well as the atmospheric pressure reduced to mean sea level at 9am and 3pm. I was also, intrigued by the somewhat stronger, negative correlations between some groups of variables such as humidity and both temperature variables and pressure and both temperature variables as well. I’m surprised that there does not seem to be any correlation between rainfall and any other predictors since the Rainfall variable is one of the most important predictors for RainTomorrow.

Rainfall

The variable Rainfall represents the amount of rainfall, in mm, recorded for the day.

weatherAUS %>%
  dplyr::select("Rainfall","RainTomorrow") %>%
  na.omit(Rainfall) %>%
  filter(Rainfall < 5) %>% 
  ggplot(aes(Rainfall)) + 
  geom_density(aes(fill=RainTomorrow)) +
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

As we can see from the plot above, the more rainfall there is on a certain day, the more likely it is that our response variable, RainTomorrow, will be Yes. This means that the next-day rainfall recorded will be 1mm or more. In order to reach this conclusion, we had filtered out the outliers.

RainToday

The variable RainToday is a boolean variable where 1 represents the precipitation (mm) in the 24 hours to 9am exceeding 1mm, otherwise it’s 0.

ggplot(weatherAUS,aes(x=RainToday,fill=RainTomorrow)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

Based on the percent stacked bar chart displayed above above, we can see that if there is less than 1mm of rainfall on a certain day, then there is a smaller chance that the next-day rain will be Yes. However, if RainToday is Yes, then there is a higher chance (50%) that RainTomorrow will also be Yes.

Temperature

The minimum and maximum temperatures could potentially affect whether or not it rains the next day. As stated before, the northern regions experience a wet season during summer months, and the southern regions have a Mediterranean climate with winter rainfall. The MinTemp and MaxTemp variables represent the minimum and maximum temperature in degrees celsius, respectively, on a certain day.

weatherAUS %>%
  dplyr::select("MaxTemp","RainTomorrow") %>%
  dplyr::mutate(MaxTemp=cut(MaxTemp,breaks=
                               seq(min(MaxTemp),max(MaxTemp),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(MaxTemp) %>%
  na.omit(MaxTemp) %>%
ggplot(aes(y=MaxTemp,fill=RainTomorrow)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

weatherAUS %>%
  dplyr::select("MinTemp","RainTomorrow") %>%
  dplyr::mutate(MinTemp=cut(MinTemp,breaks=
                               seq(min(MinTemp),max(MinTemp),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(MinTemp) %>%
  na.omit(MinTemp) %>%
ggplot(aes(y=MinTemp,fill=RainTomorrow)) +
  geom_bar(position = "fill")+
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

Based on the plot displayed above, we can conclude that when the value of MinTemp is between 23.3 and 27.8 degrees celsius, there is a higher chance that it will rain the next day. For MaxTemp, when the temperature is less than 13.1 degrees celsius, there is a higher chance that RainTomorrow will be Yes.

Humidity

It is known that humidity can have an affect on rainfall patterns. The predictors Humidity9am and Humidity3pm represent the percentage of humidity at 9am and 3pm, respectively, on the particular day an observation was recorded.

weatherAUS %>%
  dplyr::select("Humidity9am","RainTomorrow") %>%
  dplyr::mutate(Humidity9am=cut(Humidity9am,breaks=
                               seq(min(Humidity9am),max(Humidity9am),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(Humidity9am) %>%
  na.omit(Humidity9am) %>%
ggplot(aes(y=Humidity9am,fill=RainTomorrow)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

weatherAUS %>%
  dplyr::select("Humidity3pm","RainTomorrow") %>%
  dplyr::mutate(Humidity3pm=cut(Humidity3pm,breaks=
                               seq(min(Humidity3pm),max(Humidity3pm),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(Humidity3pm) %>%
  na.omit(Humidity3pm) %>%
ggplot(aes(y=Humidity3pm,fill=RainTomorrow)) +
  geom_bar(position = "fill")+
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

We can see that when the humidity was recorded at 9am on a certain day, if the percentage was higher than approximately 75%, then RainTomorrow was more likely to be Yes for that same day. For the humidity recorded at 3pm on a particular day, if the percentage is greater than about 54%, then the next-day rain is more likely to be Yes for that day as well.

Atmospheric Pressure

Atmospheric pressure is an essential factor in predicting rainfall since it influences the movement and behavior of air masses. It is usually considered in conjunction with other meteorological variables such as humidity, temperature, and wind patterns when predicting rainfall.

weatherAUS %>%
  dplyr::select("Pressure9am","RainTomorrow") %>%
  dplyr::mutate(Pressure9am=cut(Pressure9am,breaks=
                               seq(min(Pressure9am),max(Pressure9am),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(Pressure9am) %>%
  na.omit(Pressure9am) %>%
ggplot(aes(y=Pressure9am,fill=RainTomorrow)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

weatherAUS %>%
  dplyr::select("Pressure3pm","RainTomorrow") %>%
  dplyr::mutate(Pressure3pm=cut(Pressure3pm,breaks=
                               seq(min(Pressure3pm),max(Pressure3pm),by=1.5),
                             include.lowest=TRUE)) %>%
  group_by(Pressure3pm) %>%
  na.omit(Pressure3pm) %>%
ggplot(aes(y=Pressure3pm,fill=RainTomorrow)) +
  geom_bar(position = "fill")+
  scale_fill_manual(values=c("#CCFF99","#FFCCFF"))

Based on the plot produced above, if the Pressure recorded at 9am on a certain day is less than 1010 hectopascals, then there is a higher chance of having rainfall of 1mm or more recorded for the next day. It seems that at 986.5 hectopascals, the there is a 50/50 chance that will rain the next day. For the observations recorded at 3pm on a particular day, if the hectopascal value is less than 1022, then RainTomorrow is more likely to be Yes.

3. Setting Up Models

After visually looking at how the predictors effect the response variable RainTomorrow, we started to build our models. We have randomly split up our data into training and testing sets, created a recipe, and established cross-validation within our models.

Splitting into Train/Test

Now that we have started to build the models, we needed to first split our data into separate training and testing sets. Throughout the remainder of the project, we will use the training data up until the very end, where we instead will use our testing data to actually test our models. Firstly, we have set a seed so that the random split can be reproduced each time we train our models. Following this, we have carried out a training/testing split on our data stratified on our response variable RainTomorrow.

set.seed(3123)

weatherAUS_split <- initial_split(weatherAUS, prop=0.75,
                                  strata = "RainTomorrow")
weatherAUS_train <- training(weatherAUS_split)
weatherAUS_test <- testing(weatherAUS_split)

The dimensions of the training set are:

dim(weatherAUS_train)
## [1] 84693    20

The dimensions of the testing set are:

dim(weatherAUS_test)
## [1] 28232    20

Building Our Recipes

Below, we have combined our predictors and response variable to build our recipe which we will use for all our machine learning models. Essentially, we are putting the predictors together in our recipe to create our next-day rainfall prediction models.

We have used 11 out of our 19 predictor variables in our recipe. The predictors Temp9am and Temp3pmas well as WindSpeed9am and WindSpeed3pm have been excluded since these measurements are already represented by MinTemp, MaxTemp, and WindGustSpeed respectively. The fully completed recipe that will be used from now on is included below.

weatherAUS_recipe <- recipe(RainTomorrow~Month+Year+MinTemp+MaxTemp+WindGustSpeed+
                              Humidity9am+Humidity3pm+Pressure9am+Pressure3pm+Rainfall+RainToday,
                            data=weatherAUS_train) %>%
  step_dummy(all_nominal_predictors()) %>% 
  step_scale(all_predictors()) %>%
  step_center(all_predictors())
prep(weatherAUS_recipe) %>% bake(weatherAUS_train)
## # A tibble: 84,693 × 12
##    Month  Year MinTemp MaxTemp WindGustSpeed Humidity9am Humidity3pm Pressure9am
##    <dbl> <dbl>   <dbl>   <dbl>         <dbl>       <dbl>       <dbl>       <dbl>
##  1  1.60 -1.88 -0.847    0.203        0.242       -1.23       -1.23      -0.968 
##  2  1.60 -1.88  0.0327   0.288        0.392       -1.55       -0.994     -1.40  
##  3  1.60 -1.88 -0.559    0.618       -1.26        -1.18       -1.67       0.0274
##  4  1.60 -1.88  0.769    1.23         0.0165       0.771      -0.849     -0.940 
##  5  1.60 -1.88  0.305    0.861        1.14        -0.654      -1.33      -1.17  
##  6  1.60 -1.88  0.257    0.188        0.692       -0.970      -1.52      -1.11  
##  7  1.60 -1.88 -0.799    0.432       -0.434       -1.02       -1.52      -0.570 
##  8  1.60 -1.88  0.0647   0.918       -0.959       -0.496      -1.14      -1.48  
##  9  1.60 -1.88 -0.0153  -0.385        0.242       -0.126      -0.368     -2.31  
## 10  1.60 -1.88 -0.239   -0.170        0.167       -1.08       -0.897     -1.14  
## # … with 84,683 more rows, and 4 more variables: Pressure3pm <dbl>,
## #   Rainfall <dbl>, RainTomorrow <fct>, RainToday_Yes <dbl>

K-Fold Cross Validation

Here, we have stratified our cross validation on our response variable, RainTomorrow. We have also used 10 folds perform the stratified cross validation.

weatherAUS_folds <- vfold_cv(weatherAUS_train,v=10,strata=RainTomorrow)

Since building models have a very long run-time, we have saved the results to an RDA file. This has been done so that we can access it at any point after choosing our desired model.

save(weatherAUS_folds, weatherAUS_recipe, weatherAUS_train, weatherAUS_test, file = "~/Desktop/school/pstat_131/RDA/weatherAUS_model_setup.rda")

4. Building Prediction Models

We have finally been able to build our models! The dataset we are working with is so big that the models required a bit of computing power. So, each model was ran on a separate R file with the loaded data included in the zip file. Each model was saved in an RDA file and will be loaded below when discuss the results. Additionally, all of the models’ files are included in the same zip file that this specific file was in..

As mentioned before, we fit four different models. These include Logistic regression, Quadratic Discriminant Analysis, Support Vector Machine (SVM), and Random Forest. The first two models are relatively simple and will not take much time to run. The second two models are the ones we are more interested in since they are a better fit for a binary classification scenario such as this one.

Performance Metric

In order to evaluate the performance of a model, the two most useful algorithms to do this are accuracy and roc_auc. For our models, we have utilized roc_auc as the metric of performance since it shows the most significant level of efficiency in a binary classification model in which the data is imbalanced. The roc_auc metric, also known as ROC AUC, is a method for calculating the area under the curve (AUC) for the receiver operating characteristic (ROC) curve. The ROC curve is a graphical plot which demonstrates the diagnostic ability of a binary classifier system as its discrimination threshold is varied.

Model Building Process

The Decision Tree and Random Forest models we created both followed a similar process, whereas the Logistic Regression and Quadratic Discriminant Analysis models were a bit different since they are much more simple and have a quicker runtime. The general process for building the model was as follows:

  1. Firstly, we set up the model by specifying the type of model it is, and then we set its engine and its mode. For this project, our mode was always set to ‘classification’.
  2. Next, we set up the workflow for the model, added the new model, and added the weatherAUS recipe from before.

We will skip these next few steps for Logistic Regression and Quadratic Discriminant Analysis since they are simpler models that do not need their hyperparameters to be tuned.

  1. Then, we set up a tuning grid and included the parameters that we need tuned and specified the ranges for how many different levels of tuning we require for each parameter.
  2. We then tuned the model with the chosen hyperparameters.
  3. Following this, we selected the most accurate model from the tuning grid, and then finalized the workflow with the specified tuning parameters.

These last two steps apply to all of our models including the Logistic Regression and Quadratic Discriminant Analysis models.

  1. Afterwards, we fit that most accurate model with our workflow to our weatherAUS training dataset.
  2. Lastly, we saved the results to an RDA file in order to load and use it in our main project file.

5. Model Results

Our dataset is quite large with 112,925 observations, so our models took longer to run. Now, we have all of our completed models saved along with their individual outcomes and scores. From here on, we have loaded in the saved results of each model and analyzed their individual performances.

load("~/Desktop/school/pstat_131/RDA/weatherAUS_model_setup.rda")
load("~/Desktop/school/pstat_131/RDA/weatherAUS_Logistic_Regression.rda")
load("~/Desktop/school/pstat_131/RDA/weatherAUS_Decision_Tree.rda")
load("~/Desktop/school/pstat_131/RDA/weatherAUS_Quadratic_Discriminant.rda")
load("~/Desktop/school/pstat_131/RDA/weatherAUS_Random_Forest.rda")

Visualizing the Results

In order to visualize the results of the models that have been tuned, we have utilized the autoplot function in R. Doing this has visualized the effects that the change in certain parameters had on the metric roc_auc. We have displayed the plots of our models with the three best ROC AUC values, which will be discussed in the next section.

Quadratic Discriminant Analysis

A quadratic discriminant analysis is a statistical classifier model that uses a quadratic decision surface to separate measurements of two or more classes of objects or events. It is a more advanced version of a linear model as it is used to find a non-linear boundary between our classifiers, and assumes that each class follows a Gaussian distribution.

Based on the ROC curve plot below, our Quadratic Discriminant Analysis model predicted the next-day ran relatively well when it was fit back onto our training dataset.

weatherAUS_roc_qda %>%
  roc_curve(RainTomorrow, .pred_Yes) %>%
  autoplot()

Random Forest

A random forest model is a supervised ensemble learning technique that consists of numerous decision trees. Decision tree models tend to struggle with overfitting the training dataset, which is rectified in a random forest model by averaging the prediction results of each decision tree and determining a final output. The random forest algorithm stacks multiple classifiers together to improve its overall performance.

For our random forest model, we have tuned three different hyperparameters, which are: - mtry: amount of predictors that will be sampled randomly during the creation of the models - trees: amount of trees present in the random forest model - min_n: the minimum amount of data values required to be in a tree node in order for it to be split further down the tree

The ROC AUC scores reflected in the visualizations below all relatively follow the same pattern across the different node sizes. Overall, it is apparent that a higher number of trees and a low amount of predictors result in a better ROC AUC value.

autoplot(weatherAUS_rf_tune_auc)

Logistic Regression

Logistic Regression is considered to be a linear method that is used for binary classification problems such as ours. It models the relationship between a set of independent variables and a binary dependent variable using a logistic function. The algorithm estimates the probabilities of the two classes and makes predictions based on a threshold. It uses a maximum likelihood estimation approach to find the optimal coefficients that minimize the error.

Here, we summarized the predicted values by using augment to attach the predicted values to the data, and then generating a confusion matrix to visually understand the summary. We can see the 9252 of the observations that belong to the Yes category were identified correctly. The value of 9512 represents the observations belonging to the Yes class that was incorrectly classified as belonging to the No class. The amount of observations that fall within the No category, but were classified wrongly as belonging to the Yes category is 3501. The value of 62428 refers to the amount of observations belonging to the No class being classified correctly.

weatherAUS_log_fit %>%
  tidy()
## # A tibble: 12 × 5
##    term          estimate std.error statistic   p.value
##    <chr>            <dbl>     <dbl>     <dbl>     <dbl>
##  1 (Intercept)    1.85       0.0128   145.    0        
##  2 Month          0.0115     0.0110     1.05  2.95e-  1
##  3 Year           0.00760    0.0104     0.731 4.64e-  1
##  4 MinTemp       -0.0200     0.0228    -0.875 3.82e-  1
##  5 MaxTemp        0.179      0.0269     6.64  3.09e- 11
##  6 WindGustSpeed -0.564      0.0116   -48.8   0        
##  7 Humidity9am   -0.0882     0.0168    -5.24  1.61e-  7
##  8 Humidity3pm   -1.37       0.0187   -73.2   0        
##  9 Pressure9am   -1.07       0.0401   -26.6   2.28e-156
## 10 Pressure3pm    1.54       0.0404    38.0   0        
## 11 Rainfall      -0.0698     0.0117    -5.99  2.12e-  9
## 12 RainToday_Yes -0.230      0.0115   -20.0   4.90e- 89
predict(weatherAUS_log_fit,new_data=weatherAUS_train,type="prob")
## # A tibble: 84,693 × 2
##    .pred_Yes .pred_No
##        <dbl>    <dbl>
##  1   0.0410     0.959
##  2   0.0317     0.968
##  3   0.00903    0.991
##  4   0.0918     0.908
##  5   0.0758     0.924
##  6   0.0297     0.970
##  7   0.0174     0.983
##  8   0.0368     0.963
##  9   0.219      0.781
## 10   0.0686     0.931
## # … with 84,683 more rows
augment(weatherAUS_log_fit,new_data=weatherAUS_train) %>%
  conf_mat(truth=RainTomorrow,estimate=.pred_class)
##           Truth
## Prediction   Yes    No
##        Yes  9252  3501
##        No   9512 62428
augment(weatherAUS_log_fit,new_data=weatherAUS_train) %>%
  conf_mat(truth=RainTomorrow,estimate=.pred_class) %>%
  autoplot(ype="heatmap")

Model Accuracies

Below, we have summarized the best ROC values from all the models.

weatherAUS_results <- tibble(Model=weatherAUS_mod_names,
                             ROC_AUC=weatherAUS_roc_aucs)

weatherAUS_results <- weatherAUS_results %>%
  dplyr::arrange(-weatherAUS_roc_aucs)

weatherAUS_results
## # A tibble: 4 × 2
##   Model               ROC_AUC
##   <chr>                 <dbl>
## 1 Random Forest         0.999
## 2 Logistic Regression   0.862
## 3 QDA                   0.843
## 4 Decision Tree         0.797

The tibble above shows that the Random Forest model performed the best overall with a ROC AUC value of 0.9989, with the Logistic Regression Model being significantly far behind at the value of 0.8623. Since this is only fitted on the training data, we need to run our models using the testing data. In order to do this, we will be using the Random Forest Model.

6. Results From The Best Models

Since we have established that our best performing model was the Random Forest model, we have analyzed its results and determined the exact model, with the specific hyperparameter values, that had the best performance below.

Random Forest Model

We will now examine how well our Random Forest model performs on unseen data. Above, the high ROC AUC values that were observed was the model’s ability to predict the next-day rain using the data it was originally trained on. This is why the results were very strong.

The Best Model is…

The random forest model #377 is determined to have performed the best out of all of the random forest models, keeping in mind that the Random Forest model performed the best out of all of the four models we ran.

show_best(weatherAUS_rf_tune_auc,metric="roc_auc") %>%
  dplyr::select(-.estimator,.config) %>%
  slice(1)
## # A tibble: 1 × 8
##    mtry trees min_n .metric  mean     n std_err .config               
##   <int> <int> <int> <chr>   <dbl> <int>   <dbl> <chr>                 
## 1     2   100     7 roc_auc 0.880    10 0.00142 Preprocessor1_Model377

After determining the exact best model, we have now fit it to our testing data below to figure out its actual performance in predicting the next-day rain in Australia.

Final ROC AUC Results

The true ROC AUC performance results based on the testing dataset for the Random Forest model #377 has been specified below.

weatherAUS_rf_roc_auc <- augment(weatherAUS_rf_final_fit_auc,new_data=weatherAUS_test,type="prob") %>%
  roc_auc(RainTomorrow,.pred_Yes)%>%
  dplyr::select(.estimate)

weatherAUS_rf_roc_auc
## # A tibble: 1 × 1
##   .estimate
##       <dbl>
## 1     0.879

The ROC AUC score of 0.8789 based off of our testing data means that our model did very well! A value between 0.7 and 0.8 is considered to be a decent result. While there is room for improvement however, our model performed well, so this can be considered a success! Predicting the weather is fairly difficult task, so being able to say that we can predict the next-day rain in Australia very well is a big achievement.

ROC Curve

We plotted our ROC curve un order to visualize the AUC score. The more hugh and to the left the curve lies, the better the model’s AUC will be. The plot below does a fairly good job of following the pattern of the top left right angle of a square and can confirm our computed AUC score above.

augment(weatherAUS_rf_final_fit_auc,new_data=weatherAUS_test,type="prob") %>%
  roc_curve(RainTomorrow,.pred_Yes) %>%
  autoplot()

7. Conclusion

Through the course of this project, we have explored and analyzed our data in order to create and test a model that could predict whether or not it will rain the next day in Australia. Following thorough analysis, testing and computing, we can say that the Random Forest model #377 was the best at predicting the next-day rain in Australia.

Despite the fact that we attained such a high accuracy value for prediction, there is still room for improvement. However, it is difficult when thinking of what specific methods to implement to improve this value. We already started out with a large range of predictors to variables to help anticipate next-day rain. Additionally, majority of the predictors we have used in this project are highly correlated. So, trying to find more predictors to include might not even be beneficial since the minor lack of prediction accuracy of our model could be due to uncontrollable factors of the weather patterns in Australia that is captured by the data itself.

Looking at how our model can be applied realistically, we can say that it performed way better than we could have anticipated. Weather is very hard to predict due to its randomness. Despite having the newest and state-of-the-art technology, weather forecasters still struggle to accurately predict rain, sunshine, or even thunderstorms. We can be very proud of the fact that we were able to achieve a very accurate prediction of next-day rain in Australia.