---
title: "Solar Radiation Prediction"
author: "Alireza Mohammadi (250406221)"
email: "amoham69@uwo.ca"
date: "November 4, 2017"
output:
  pdf_document:
    number_sections: true  
  html_document: default
bibliography: MyReferences.bib
header-includes:
- \usepackage{color}
- \usepackage{titling}
- \usepackage{gensymb}
- \setlength{\parindent}{2em}
- \usepackage{amsmath}
- \usepackage{amssymb}
- \usepackage{wasysym}
---

<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: { 
      equationNumbers: {
 
            autoNumber: "all",
            formatNumber: function (n) {return '9.'+n}
      } 
  }
});
</script>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(GGally)
library(knitr)
library(stargazer)
library(dplyr)
library(lubridate)
library(bestglm)
library(car)
library(tseries)
library(glmnet)
vifx <- function (X) 
{#Author: A.I. McLeod
    Xs <- scale(X)
    xx <- t(Xs) %*% Xs/(nrow(X) - 1)
    diag(solve(xx))
}
rmse <- function(y, yh) sqrt(mean( (y - yh)^2 ))
```

# Introduction

With the rapid depletion of natural resources and at the same time the air pollution at all-time high, more and more emphasis is shifting towards methods that can efficiently produce energy from environmentally friendly renewable resources such as sun, wind, ocean waves, etc. Aside from that, the main source of energy on the surface of Mars is the solar energy and therefore identifying and understanding effects of key variables on solar radiation are of paramount values as they can help with building a model that can accurately predict the amount of solar radiation as a function of some measurable predictors. Such model can equip the astronauts or colonists with crucial information about when and where to deploy the solar energy harvesting equipment for optimal performance. Building models like this is also very useful on Earth for the same reason. 

In this project, we are interested to build a suitable model using multiple linear regression for prediction of the average solar irradiance received in a day as a function of basic meteorological variables of that day. Solar irradiance is power per unit area received from the Sun. The solar power is in the form of the electromagnetic radiation which is emitted from the Sun at various range of wavelengths [@Incropera02 pp. 700-787]. As a result, the amount of solar power measured by a measuring instrument depends on the range of wavelengths that the device can actually detect. For convenience, we shall refer to solar irradiance as solar radiation or just radiation in the rest of this report. 

![The [HI-SEAS](http://hi-seas.org) Habitat \label{fig:HiSEAS} (image from @HiSEAS).](habitat.jpg)

To build a simple model that can predict the solar radiation, we use the dataset that was collected by [NASA HI-SEAS](http://hi-seas.org/) mission (Hawai'i Space Exploration Analog and Simulation, see @HiSEAS) and was provided on [Kaggle](https://www.kaggle.com/dronio/SolarEnergy). HI-SEAS is a Habitat located on the Big Island of Hawaii at approximately 2500 meters above the sea level (see Figure \ref{fig:HiSEAS}). This isolated and unique location was carefully chosen by NASA research team so that it resembles environment on a Mars site. In particular, this location has little variation in weather (with usually cool and dry climate) which makes it suitable for solar radiation study. This report has been prepared in **R** [@Rprogram].

# Description of the Dataset

The dataset contains the measurements of the solar radiation along with meteorological measurements recorded daily at several time intervals from September 1^st^, 2016 until December 31^st^, 2016. In total there are 32,686 observations in the dataset. The raw data consists of UNIX time (indicating computer time in seconds measured from January 1^st^, 1970), date and time of observation, radiation ($W/m^2$), temperature (\degree F), pressure (in-Hg), wind direction (degree), wind speed (mph), time of sunrise and sunset. The UNIX time column is very helpful especially in sorting the data. The rest of columns were produced in MS Excel from these columns. In particular, *daylight* denotes the duration of daylight, that is the time between sunrise and sunset in seconds. The column with heading *dt* is the time interval between two consecutive measurements in seconds and the rest of columns is basically the numerical integration of each variables in each time intervals which will be used later. The goal of this project is to develop a model that can predict the radiation as a function of explanatory variables. A quick check of the data reveals that no observations were taken on these four days: September 30^th^, November 30^th^ December 6^th^ and 7^th^. The total number of days that data were collected is $n=118$ days. Majority of data was recorded at approximately 5 minute time intervals. However, we notice that the time intervals are not the same for all observations. There are time intervals of approximately 10, 15, 20 minutes and even for a few cases of several hours. This is important especially in computation of daily average values. This means that we cannot simply use the sample mean (R function *mean()*) and we need to do a numerical integration in order to compute the correct daily average values. Before computing the integrals, we group the data conveniently by date. The daily average value then is defined by $\bar f = (\int_{t_1}^{t_2} f dt) / (\int_{t_1}^{t_2} dt)$ where $f$ denotes any of the input/output variables and $t_1$ and $t_2$ denote the time of the first and last observations in each day, respectively. The time integrals are numerically computed using simple trapezoidal rule.

```{r DATAINPUT, echo=FALSE, results="asis", warning=FALSE}
rawdata <- read.csv("SolarPrediction.csv", header=TRUE)
#head(rawdata) # there are no results on Dec 6 and Dec 7. The total number of days is 118
names(rawdata)[names(rawdata)=="�..UNIXTime"] <- "UNIXTime" # importing the data from excel for some reason has changed the heading of the first column, I am using this line to correct this
rawdata$Date <- as.Date(rawdata$Date, "%m/%d/%Y") # It is important to use correct format for Date
#head(rawdata)
#summary(rawdata)
data_by_day <- rawdata %>%  group_by(Date) %>% summarise(R = sum(R_dt)/sum(dt), # average daily radiation
            T = sum(T_dt)/sum(dt), # average daily temperature
            P = sum(P_dt)/sum(dt), # average daily pressure
            H = sum(H_dt)/sum(dt), # average daily humidity
            WD = sum(WD_dt)/sum(dt), # average daily wind direction
            WS = sum(WS_dt)/sum(dt), # average daily wind speed
            DayRatio = mean(Daylight)/(24*3600), # ratio of daylight over total day length
            TimeObsRatio = sum(dt)/(24*3600)) # ratio of length of time data was recorded over total day length
#head(data_by_day)
#summary(data_by_day)
mydata <- data.frame(data_by_day)
stargazer(mydata, type="latex", title="Dataframe Summary.", header=FALSE, label="tab:main_stat")
```
In total, we can define up to six linear predictors in our model ($p = 6$). The predictors are daily average temperature (T), pressure (P), humidity (H), wind direction (WD), wind speed (WS) and day ratio (defined as ${\rm DayRatio} = \frac{{\rm daylight}}{24 \times 3600}$). The response variable is the daily average solar radiation (R). All the input/output variables in this analysis are the daily average values and, for convenience, we may not explicitly mention **\textcolor{blue}{daily average}** everywhere we want to refer to them throughout the rest of this report. The statistical summary of the data is given in Table \ref{tab:main_stat}. The last row (TimeObsRatio) indicates the ratio of each day that data was recorded defined by $(\int_{t_1}^{t_2} dt) / (24 \times 3600)$. This variable is not an explanatory variable but rather is defined to help with interpretation of other variables. The maximum value of this ratio is slightly above 1 and the reason is that for some cases the last point of integration was about a couple of seconds into the next day. Moreover, we notice that the duration of data collection in a day can be as low as `r round(100*min(mydata$TimeObsRatio),2)`% of a day (look at the minimum value of TimeObsRatio). It is worthwhile to note that the daily average and instantaneous pressure fluctuations are very small and essentially negligible. One may quickly expect that pressure should not play a significant role in the model.  

We will perform a systematic analysis in a quest to find the best model for prediction of the solar radiation using available explanatory variables.

# Preliminary Analysis of Dataset

## Scatterplot Matrix

As the first step of exploration of the dataset, we use the scatterplot matrix shown in
Figure \ref{fig:scatterplot}. 

```{r SCATTERPLOTMATRIX, echo=FALSE, fig.height=7, fig.width=7, fig.cap="\\label{fig:scatterplot}Scatterplot Matrix."}
pairs(mydata[,c(2,3,4,5,6,7,8)], panel=panel.smooth)
```

The main conclusions from the above scatterplot matrix can be summarized as:

* Reading across the top row for **radiation**, it seems that **humidity** is the most important predictor as the points are more tightly clustered around the loess. We further notice that **humidity** is negatively correlated with **radiation** and thus the lower the humidity, the higher is the amount of radiation received on the surface of the Earth. This physically makes sense as in this case higher portion of solar radiation is either absorbed by the water molecules in the air or by the clouds in the case of cloudy sky or is just simply reflected back into the sky and away from the Earth. **Pressure** and **wind direction** show the weakest influence on **radiation**. There is a moderate positive association between **wind speed** and **radiation**. **Temperature** and **day ratio** are both positively correlated with **radiation**.

* Reading across the second row from the top for **temperature**, we note that except the wind variables (wind speed and wind direction), **temperature** is associated with all other variables. **Temperature** is positively associated with **radiation**, **pressure** and **day ratio**, and is negatively associated with **humidity**. Positive association of **temperature** with **pressure** may remind us of the [ideal gas law](https://en.wikipedia.org/wiki/Ideal_gas_law) $PV = n\bar{R}T$, where $P$ is pressure, $V$ is volume of gas, $n$ is the number of kmol of gas, $\bar{R}$ is the universal gas constant and $T$ is the absolute temperature [@Sonntag03 pp. 61-66].

* **Pressure** (third row from the top), except for the positive association with **temperature** mentioned above, does not illustrate any noticeable association with any other variables.

* Reading across the fourth row from the top for **humidity**, except for the negative association with **radiation** and **temperature** discussed above, we do not see any other clear association. 

* Wind variables (**wind direction** and **wind speed**, rows five and six, respectively) demonstrate the weakest associations with any other variables.

* **Day ratio** is clearly positively associated with **radiation** and **temperature** (see the last row). These associations are not surprising since the length of a day is longer during warmer seasons and thus both average **radiation** and **temperature** are higher during warmer seasons. 

In summary, the scatterplot matrix (Fig. \ref{fig:scatterplot}) suggests that: (i) radiation is most closely predicted by humidity, and (ii) pressure and wind direction are the least important explanatory variables.

## Variance Inflation Factor (VIF)

There is multicollinearity in the explanatory variables when one of them can be linearly predicted with high accuracy using other explanatory variables. The variance inflation factor (VIF) for the design matrix indicates which variables contribute to multicollinearity. The VIF is shown in Figure \ref{fig:vif} using bar chart. As can be seen, there is no significant multicollinearity in the dataset as all VIF values are much smaller than the empirical threshold ${\rm VIF} = 10$ (marked by red dotted line in Fig.\ref{fig:vif}). 

```{r VIF-barchartfig, echo=FALSE, fig.height=3, fig.width=5.25, fig.cap="\\label{fig:vif}VIF for solar radiation predictors."}
ans_vif <- vifx(mydata[,c("T", "P", "H", "WD", "WS", "DayRatio")])
barplot(ans_vif, col="blue", ylab="VIF", ylim=c(0, 12))
abline(h=10, lty=2, col="red", lwd=2)
```

\newpage

Correlation matrix is shown in Table \ref{tab:cormat} which provides another useful tool for detecting possible multicollinearity. A close inspection of the correlation matrix reveals that **temperature** and **day ratio** have a correlation `r round(cor(mydata[,"T"], mydata[,"DayRatio"]),4)` which is considered a moderate correlation since the coefficient of determination for the regression of **temperature** on **day ratio** is only `r round(100*cor(mydata[,"T"], mydata[,"DayRatio"])^2,2)`%.

```{r CORMAT, echo=FALSE, results="asis"}
X <- mydata[,c("T", "P", "H", "WD", "WS", "DayRatio")]
stargazer(cor(X), type="latex", title="Correlation Matrix.", header=FALSE, label="tab:cormat")
```

It may be helpful to also visualize the correlation matrix as illustrated in Figure \ref{fig:ggcor}. The moderate positive correlations between **temperature** and **pressure**, and also **temperature** and **day ratio** are easily noticeable in this figure.

```{r GGCOR, echo=FALSE, fig.height=3.5, fig.width=4.5, fig.cap="\\label{fig:ggcor}Visualization of the Correlation Matrix."}
ggcorr(X)
```

\newpage

# Model Building

## Initial Model

```{r TEST-TRAIN, echo=FALSE, results="asis"}
R <- mydata[,2] # Y = R is the response variable
Xy <- as.data.frame(cbind(X,R))
set.seed(922157)
nTe <- round(nrow(Xy)/3) #39
nTr <- nrow(Xy)-nTe #79
indTr <- sample(1:nrow(Xy), nTr) # sample the indeces
indTe <- !'%in%'(1:nrow(Xy), indTr)
XyTr <- Xy[indTr,]
XyTe <- Xy[indTe,]
p <- ncol(XyTr)-1
XTr <- as.matrix.data.frame(XyTr[,1:p])
yTr <- XyTr[,p+1]
XTe <- as.matrix.data.frame(XyTe[,1:p])
yTe <- XyTe[,p+1]
```

```{r BEST-SUBSET, echo=FALSE, results="asis"}
BEST_glm<-bestglm(XyTr, IC="BIC")
ans0<-lm(R ~ H + WD + WS + DayRatio, data=XyTr)
```

As the first step towards constructing our model, we divide the dataset into two parts: train and test datasets. The training data is used for model building. The test data is used for model validation. We randomly pick $\frac{2}{3}$ of the data as training set, i.e. number of observations in the training set is $n_{tr}=$ `r nTr`. The remaining observations are stored in the test set, $n_{te}=$ `r nTe`. It is suggested that when $p>5$, the best subset method can be used for model selection [see @note:McLeod17, pp. 9]. For this purpose, we use function *bestglm::bestglm()* developed by @rpack:McLeod17 with BIC as criterion. The best subset method suggests the following model: 

\begin{equation}
  R = \beta_0 + \beta_1 {\rm H} + \beta_2 ({\rm WD}) + \beta_3 ({\rm WS}) + \beta_4 ({\rm DayRatio}) + error.
  \label{eq:model0}
\end{equation}

Following this suggestion, we fit a linear model with predictors as described in Eq.\eqref{eq:model0}. Table \ref{tab:model_0} summarizes the fitted linear regression model `r as.character(ans0$call)[2]`. It shows that all predictors are significant at less than 1%. However $R^2 =$ `r round(100*summary(ans0)$r.squared, 1)`% is not that impressive and hence the proportion of variability explained by the model is low. We use basic model diagnostic checks to see how we may improve this model.

```{r ans0Table, echo=FALSE, results="asis"}
stargazer(ans0, type="latex", title="Summary of the fitted linear regression model ${\\rm R} \\sim {\\rm H} + {\\rm WD} + {\\rm WS} + {\\rm DayRatio}$.", header=FALSE, label = "tab:model_0")
```

\newpage

### Basic Diagnostic Checks

The basic diagnostic checks for the model `r as.character(ans0$call)[2]` are shown
in Figure \ref{fig:basic_diag0}.

```{r BASIC-DIAG0, echo=FALSE, fig.height=6, fig.width=6, fig.cap="\\label{fig:basic_diag0}Basic regression diagnostic plots for model ${\\rm R} \\sim {\\rm H} + {\\rm WD} + {\\rm WS} + {\\rm DayRatio}$."}
layout(matrix(1:4, ncol=2, nrow=2))
plot(ans0)
layout(1)
```

A problem is indicated in the *Residuals vs Fitted* diagnostic plot since the loess
trend is not flat. The curve suggests possible nonlinearity due to interaction and/or some nonlinearity
present in the inputs. Also, the *Normal Q-Q* plot suggests that the distribution of the residuals does not follow a normal distribution. However, the *Residuals vs leverage* plot shows that none of the observations have a strong influence on the fit.

\newpage

### Residual Dependency Checks

The residual dependency plots shown in Figure \ref{fig:resdep0} indicates that
indeed **humidity** as well as **day ratio** and possibly **wind direction** exhibit nonlinear effects.

```{r RESDEP-ans0, echo=FALSE, fig.height=5, fig.width=5, fig.cap="\\label{fig:resdep0}Residual dependency plots for model ${\\rm R} \\sim {\\rm H} + {\\rm WD} + {\\rm WS} + {\\rm DayRatio}$."}
ANS <- ans0 #set lm output
XyTr %>%
 select(H, WD, WS, DayRatio, T, P) %>%
 mutate(res=resid(ANS)) %>%
 gather(value, key=var, -res) %>%
 ggplot(mapping=aes(x=value, y=res)) +
 geom_point() +
 stat_smooth(method=loess, se=FALSE, method.args=list(span=2)) +
 facet_wrap(~var, scales="free") +
 xlab("Explanatory variables") +
 ylab("Residual")
 # ggtitle("Residual dependency plots", subtitle = toString(ANS$call))
```

\newpage

## Revised Model 1

The diagnostic checks suggest that the initial model may be improved by including
a quadratic term with humidity so the regression model can be represented as

\begin{equation}
  R = \beta_0 + \beta_1 {\rm H} + \beta_2 {\rm H}^2 + \beta_3 ({\rm WD}) + \beta_4 ({\rm WS}) + \beta_5         ({\rm DayRatio}) + error.
  \label{eq:model1}
\end{equation}

```{r FIT-ans1, echo=FALSE, results="asis"}
ans1 <- lm(R ~ poly(H,2) + WD + WS + DayRatio, data=XyTr)
anova0_1 <- anova(ans0,ans1)
pval0_1 <- anova0_1$`Pr(>F)`[2]
```

Table \ref{tab:model_1} summarizes the fitted linear regression model 1 (see Eq.\eqref{eq:model1}). We see that the quadratic term is significant at less than 1%. The $R^2$ has increased by `r round(100*(summary(ans1)$r.squared-summary(ans0)$r.squared), 1)`% from the previous model to `r round(100*summary(ans1)$r.squared, 1)`%  and the residual sum squares (RSS) has dropped to ${\rm RSS}=$ `r toString(round(sum(resid(ans1)^2), 1))` from ${\rm RSS}=$ `r toString(round(sum(resid(ans0)^2), 1))` for the previous model. The analysis of variance (ANOVA) lack-of-fit test comparing this model with the initial model gives F-statistic $F=$ `r round(anova0_1$F[2], 2)` on (`r ans0$df.residual - ans1$df.residual`, `r ans1$df.residual`) DF and has a two-sided p-value `r format(100*pval0_1, scientific = FALSE, digits = 4)`% which shows that the null-hypothesis  (H~0~: the initial model is true) is rejected at level less than 0.1%. Thus far this model has shown some improvement compared with the previous model. Next, we will look at the basic diagnostic plots for this model.

```{r ans1Table, echo=FALSE, results="asis"}
stargazer(ans1, type="latex", title="Summary of the fitted linear regression for model R $\\sim$ poly(H,2) + WD + WS + DayRatio.", header=FALSE, label = "tab:model_1")
```

\newpage 

### Basic Diagnostic Checks for Revised Model 1

The basic regression plots for model `r as.character(ans1$call)[2]` clearly show patters in *Residuals vs Fitted* as well as *Scale-Location* diagnostic plots (see Figure \ref{fig:basic_diag1}). It seems that there is still some nonlinearity that has not been adequately accounted for by the model. The *Normal Q-Q* plot has improved and suggests that the distribution of the residuals for this model is closer to a normal distribution compared with the initial model. The *Residuals vs leverage* plot shows that none of the observations have a strong influence on the fit.

```{r BASIC-DIAG1, echo=FALSE, fig.height=6, fig.width=6, fig.cap="\\label{fig:basic_diag1}Basic regression diagnostic plots for model R $\\sim$ poly(H,2) + WD + WS + DayRatio."}
layout(matrix(1:4, ncol=2, nrow=2))
plot(ans1)
layout(1)
```

\newpage 

### Residual Dependency Checks for Revised Model 1

To further analyze which variables may still have nonlinear behaviour, we employ the residual dependency plots shown in Figure \ref{fig:resdep1}. We can see from the residual dependency plots for **day ratio** and **wind direction** that they both display possible nonlinear behaviour. We will add a quadratic term for day ratio as the next step which is discussed in the next section.

```{r RESDEP-ans1, echo=FALSE, fig.height=5, fig.width=5, fig.cap="\\label{fig:resdep1}Residual dependency plots for model R $\\sim$ poly(H,2) + WD + WS + DayRatio"}
ANS <- ans1 #set lm output
XyTr %>%
 select(H, WD, WS, DayRatio, T, P) %>%
 mutate(res=resid(ANS)) %>%
 gather(value, key=var, -res) %>%
 ggplot(mapping=aes(x=value, y=res)) +
 geom_point() +
 stat_smooth(method=loess, se=FALSE, method.args=list(span=2)) +
 facet_wrap(~var, scales="free") +
 xlab("Explanatory variables") +
 ylab("Residual")
 # ggtitle("Residual dependency plots", subtitle = toString(ANS$call))
```

\newpage 

## Revised Model 2

```{r FIT-ans2, echo=FALSE, results="asis"}
ans2 <- lm(R ~ poly(H,2) + WD + WS + poly(DayRatio,2), data=XyTr)
anova1_2 <- anova(ans1,ans2)
pval1_2 <- anova1_2$`Pr(>F)`[2]
```

As the next step in our model building procedure, we add a quadratic term for day ratio. Thus, the revised model takes the form

\begin{equation}
  R = \beta_0 + \beta_1 {\rm H} + \beta_2 {\rm H}^2 + \beta_3 ({\rm WD}) + \beta_4 ({\rm WS}) + \beta_5         ({\rm DayRatio}) + \beta_6 ({\rm DayRatio})^2 + error.
  \label{eq:model2}
\end{equation}

The summary of the fitted model based on Eq.\eqref{eq:model2} is provided in Table \ref{tab:model_2}. Note that all regression coefficients are significant at less than 1%. The $R^2$ has increased by `r round(100*(summary(ans2)$r.squared-summary(ans1)$r.squared), 1)`% from the previous model to `r round(100*summary(ans2)$r.squared, 1)`%  and the residual sum squares has dropped to ${\rm RSS}=$ `r toString(round(sum(resid(ans2)^2), 1))` from ${\rm RSS}=$ `r toString(round(sum(resid(ans1)^2), 1))` for the previous model. The ANOVA lack-of-fit test comparing this model with the revised model 1 gives F-statistic $F=$ `r round(anova1_2$F[2], 2)` on (`r ans1$df.residual - ans2$df.residual`, `r ans2$df.residual`) DF and has a two-sided p-value `r format(100*pval1_2, scientific = FALSE, digits = 4)`% which shows that the null-hypothesis (H~0~: the revised model 1 is true) is rejected at level less than 1%. Thus, this model is preferred compared with the previous model. The next logical step is to analyze the basic diagnostic plots for this model.

```{r ans2Table, echo=FALSE, results="asis"}
stargazer(ans2, type="latex", title="Summary of the fitted linear regression for model R $\\sim$ poly(H,2) + WD + WS + poly(DayRatio,2).", header=FALSE, label = "tab:model_2")
```

\newpage

### Basic Diagnostic Checks for Revised Model 2

The basic diagnostic plots for model `r as.character(ans2$call)[2]` are shown in Figure \ref{fig:basic_diag2}. Patterns are still easily detectable especially in *Scale-Location* plot (see Figure \ref{fig:basic_diag1} top right plot). This signals the possibility of presence of some other form of nonlinearity and/or interaction effects that have not been  handled well by the model. The *Normal Q-Q* plot has improved from the previous model 1 and suggests that the distribution of the residuals of this model is close enough to a normal distribution. The *Residuals vs leverage* plot shows that none of the observations have a strong influence on the fit.

```{r BASIC-DIAG2, echo=FALSE, fig.height=6, fig.width=6, fig.cap="\\label{fig:basic_diag2}Basic regression diagnostic plots for model  R $\\sim$ poly(H,2) + WD + WS + poly(DayRatio,2)."}
layout(matrix(1:4, ncol=2, nrow=2))
plot(ans2)
layout(1)
```

\newpage

### Residual Dependency Checks for Revised Model 2

Residual dependence plots are shown for revised model 2 in Figure \ref{fig:resdep2}. The only predictor that still demonstrates some nonlinear behaviour is wind direction. In the next section we will add a quadratic term for wind direction.

```{r RESDEP-ans2, echo=FALSE, fig.height=5, fig.width=5, fig.cap="\\label{fig:resdep2}Residual dependency plots for model R $\\sim$ poly(H,2) + WD + WS + poly(DayRatio,2)"}
ANS <- ans2 #set lm output
XyTr %>%
 select(H, WD, WS, DayRatio, T, P) %>%
 mutate(res=resid(ANS)) %>%
 gather(value, key=var, -res) %>%
 ggplot(mapping=aes(x=value, y=res)) +
 geom_point() +
 stat_smooth(method=loess, se=FALSE, method.args=list(span=2)) +
 facet_wrap(~var, scales="free") +
 xlab("Explanatory variables") +
 ylab("Residual")
 # ggtitle("Residual dependency plots", subtitle = toString(ANS$call))
```

\newpage 

## Revised Model 3

```{r FIT-ans3, echo=FALSE, results="asis"}
ans3 <- lm(R ~ poly(H,2) + poly(WD,2) + WS + poly(DayRatio,2), data=XyTr)
anova2_3 <- anova(ans2,ans3)
pval2_3 <- anova2_3$`Pr(>F)`[2]
# fit the same model on the whole dataset
ans_complete <- lm(R ~ poly(H,2) + poly(WD,2) + WS + poly(DayRatio,2), data=Xy)
```

Now we add a quadratic term for wind direction. The revised model 3 takes the form

\begin{equation}
  R = \beta_0 + \beta_1 {\rm H} + \beta_2 {\rm H}^2 + \beta_3 ({\rm WD}) + \beta_4 ({\rm WD})^2 + \beta_5 ({\rm WS}) + \beta_6 ({\rm DayRatio}) + \beta_7 ({\rm DayRatio})^2 + error.
  \label{eq:model3}
\end{equation}

```{r ans3Table, echo=FALSE, results="asis"}
stargazer(ans3, ans_complete, type="latex", title="Summary of the fitted linear regression for model R $\\sim$ poly(H,2) + poly(WD,2) + WS + poly(DayRatio,2).", header=FALSE, label = "tab:model_3", column.labels=c("train dataset","complete dataset"))
```

The summary of the fitted model 3 described in Eq.\eqref{eq:model3} is provided in Table \ref{tab:model_3} (see the train dataset column). Note that all regression coefficients are significant at less than 1%. The coefficient of determination for this model is $R^2=$ `r round(100*summary(ans3)$r.squared, 1)`% which shows an increase of `r round(100*(summary(ans3)$r.squared-summary(ans2)$r.squared), 1)`% from the previous model. The residual sum squares has dropped to ${\rm RSS}=$ `r toString(round(sum(resid(ans3)^2), 1))` from ${\rm RSS}=$ `r toString(round(sum(resid(ans2)^2), 1))` for the previous model. The ANOVA lack-of-fit test comparing this model with the revised model 2 gives F-statistic $F=$ `r round(anova2_3$F[2], 2)` on (`r ans2$df.residual - ans3$df.residual`, `r ans3$df.residual`) DF and has a two-sided p-value `r format(100*pval2_3, scientific = FALSE, digits = 3)`% which indicates that the null-hypothesis (H~0~: the revised model 2 is true) is rejected at level less than 0.1%. As will be discussed later, this model provides the best performance compared with all other linear regression models we tested in this project. The results of the same model (Eq.\eqref{eq:model3}) fitted on the complete dataset is also provided in the second column of Table \ref{tab:model_3} for comparison purposes and as can be seen, for instance, $R^2=$ `r round(100*summary(ans_complete)$r.squared, 1)`% for the complete dataset which shows that the proportion of variability explained by this model is very significant and therefore we may have found a useful model provided that all other diagnostic tests are OK. Next we will look at the basic diagnostic plots for this model.

### Basic Diagnostic Checks for Revised Model 3

Basic diagnostic plots for model `r as.character(ans3$call)[2]` are shown in Figure \ref{fig:basic_diag3}. No noticeable patters or fan-shape behaviour are detected in *Residuals vs Fitted* and *Scale-Location* plots. Hence, there is no monotonic variance change. The *Normal Q-Q* plot demonstrates a satisfactory shape. The *Residuals vs leverage* plot shows that none of the observations have a strong influence on the fit.

```{r BASIC-DIAG3, echo=FALSE, fig.height=6, fig.width=6, fig.cap="\\label{fig:basic_diag3}Basic regression diagnostic plots for model R $\\sim$ poly(H,2) + poly(WD,2) + WS + poly(DayRatio,2)."}
layout(matrix(1:4, ncol=2, nrow=2))
plot(ans3)
layout(1)
```

\newpage

### Residual Dependency Checks for Revised Model 3

Residual dependence plots are shown for the revised model 3 in Figure \ref{fig:resdep3}. The loess curves are approximately flat for the explanatory variables of this model (day ratio, humidity and wind direction) and hence the model has sufficiently handled the nonlinearity.

```{r RESDEP-ans3, echo=FALSE, fig.height=5, fig.width=5, fig.cap="\\label{fig:resdep3}Residual dependency plots for model R $\\sim$ poly(H,2) + poly(WD,2) + WS + poly(DayRatio,2)"}
ANS <- ans3 #set lm output
XyTr %>%
 select(H, WD, WS, DayRatio, T, P) %>%
 mutate(res=resid(ANS)) %>%
 gather(value, key=var, -res) %>%
 ggplot(mapping=aes(x=value, y=res)) +
 geom_point() +
 stat_smooth(method=loess, se=FALSE, method.args=list(span=2)) +
 facet_wrap(~var, scales="free") +
 xlab("Explanatory variables") +
 ylab("Residual")
 # ggtitle("Residual dependency plots", subtitle = toString(ANS$call))
```

Next, we compare the AIC and BIC values of all the above four models. Table 7 summarizes these results and shows that the revised model 3 has the lowest AIC and BIC values and are preferred based on these criteria. Backward stepwise/stagewise methods based on AIC and BIC criteria also suggests that the revised model 3 has the best performance (not shown).

```{r AIC-BIC, echo=FALSE}
AIC_results <- AIC(ans0, ans1, ans2, ans3)
BIC_results <- BIC(ans0, ans1, ans2, ans3)
AIC_BIC <- matrix(c(AIC_results$AIC, BIC_results$BIC), ncol=4, byrow=TRUE)
colnames(AIC_BIC) <- c("initial model", "revised model 1", "revised model 2", "revised model 3 (best model)")
rownames(AIC_BIC) <- c("AIC", "BIC")
kable(AIC_BIC, type="latex", digits=2, caption="Summary of the AIC and BIC values for different models.", header=FALSE, label = "tab:AIC_BIC")
#stargazer(AIC_BIC, type="latex", title="Summary of the AIC and BIC values for different models.", header=FALSE, label = "tab:AIC_BIC")
```

```{r ADDITIONAL-STEPWISE, echo=FALSE, results="asis", eval=FALSE}
ansAIC <- step(ans3)
ansBIC <- step(ans3, k=log(nrow(XyTr)))
```

```{r ADDITIONAL-NORM, echo=FALSE, results="asis"}
# let us test the normality assumption again using shapiro-wilk test
shapiro_ans <- shapiro.test(resid(ans3))
# Since p-value is large and therefore the null-hypothesis cannot be rejected. There is no evidence that that the data are not from a normally distributed population.
jarque_bera_ans <- jarque.bera.test(resid(ans3))
```

Shapiro-Wilk and Jarque-Bera normality tests give p-values of `r round(100*shapiro_ans$p.value, 2)`% and `r round(100*jarque_bera_ans$p.value, 2)`%, respectively, and therefore the null-hypothesis (H~0~: data are from a normally distributed population) cannot be rejected. In other words, there is no evidence that the data are not from a normally distributed population as far as these tests are concerned.

```{r ADDITIONAL-SERIAL, echo=FALSE, results="asis"}
# let us check whether there is serial correlation
dw_ans_normal <- durbinWatsonTest(ans3, method = "normal") # use normal
dw_ans_boot <- durbinWatsonTest(ans3, method = "resample") # use bootstrap
# Since the p-value is greater than 10% so the null-hypothesis cannot be rejected which means there is no evidence to support possible serial dependence.
```

Durbin-Watson test is useful to detect serial correlation. The null-hypothesis for this test is H~0~: there is no correlation among residuals (no serial time dependence). We performed this test using both normal and bootstrap methods and the p-values are `r round(100*dw_ans_normal$p, 2)`% and  `r round(100*dw_ans_boot$p, 2)`%, respectively. The D-W statistic is $d=$ `r round(dw_ans_normal$dw, 2)` and the sample autocorrelation of the residuals is $r=$ `r round(dw_ans_normal$r, 2)`. Thereby, the null-hypothesis cannot be rejected which means there is no evidence to support possible serial dependence in the residuals based on this test.

## Overfitting Lack-of-Fit Test

```{r FIT-ans4, echo=FALSE, results="asis"}
ans4 <- lm(R ~ poly(H,2) + poly(WD,2) + poly(WS,2) + poly(DayRatio,2), data=XyTr)
anova3_4 <- anova(ans3,ans4)
pval3_4 <- anova3_4$`Pr(>F)`[2]
```

In this section we overfit the best model (revised model 3) and perform various tests to judge whether the overfitted model is useful or not. To make the model overfit, let us add a quadratic term for wind speed. So the model can be represented as

\begin{equation}
  R = \beta_0 + \beta_1 {\rm H} + \beta_2 {\rm H}^2 + \beta_3 ({\rm WD}) + \beta_4 ({\rm WD})^2 + \beta_5 ({\rm WS}) + \beta_6 ({\rm WS})^2 + \beta_7 ({\rm DayRatio}) + \beta_8 ({\rm DayRatio})^2 + error.
  \label{eq:model4}
\end{equation}

The summary of the fitted model described by Eq.\eqref{eq:model4} is provided in Table \ref{tab:model_4}. We notice that based on the two-sided p-value for the t-test which is `r round(100*summary(ans4)$coefficients[7,4], 2)`%, the quadratic term for wind speed is not significant at 10% and therefore should be dropped. The $R^2$ has increased by only `r round(100*(summary(ans4)$r.squared-summary(ans3)$r.squared), 1)`% from the previous model to `r round(100*summary(ans4)$r.squared, 1)`%  and the residual sum squares has dropped to ${\rm RSS}=$ `r toString(round(sum(resid(ans4)^2), 1))` from ${\rm RSS}=$ `r toString(round(sum(resid(ans3)^2), 1))` for the revised model 3. 

```{r ans4Table, echo=FALSE, results="asis"}
stargazer(ans4, type="latex", title="Summary of the fitted linear regression for model R $\\sim$ poly(H,2) + poly(WD,2) + poly(WS,2) + poly(DayRatio,2).", header=FALSE, label = "tab:model_4")
```

The ANOVA lack-of-fit test comparing this model with the revised model 3 gives F-statistic $F=$ `r round(anova3_4$F[2], 2)` on (`r ans3$df.residual - ans4$df.residual`, `r ans4$df.residual`) DF and has a two-sided p-value `r format(100*pval3_4, scientific = FALSE, digits = 4)`% (the same as p-value for the t-test discussed above) which implies that the null-hypothesis (H~0~: the revised model 3 is true) cannot be rejected at 10%. The BIC values are ${\rm BIC}_{3}=$ `r round(BIC(ans3), 2)` and ${\rm BIC}_{\rm overfit}=$ `r round(BIC(ans4), 2)` for the revised model 3 and the overfitted model, respectively, which  further reinforces the conclusion obtained by the ANOVA test and suggests that the simpler model (revised model 3) is preferred. We have also tried other possibilities for the overfit model such as adding temperature or pressure and found out that adding these predictors substantially deteriorate the performance as compared with the revised model 3 and these overfitted models can be easily rejected using ANOVA or AIC/BIC tests (not shown).

\newpage

# Comparison with Other Methods

We showed in the previous sections, following a step by step approach, that the revised model 3 (best model) provides superior performance compared with all other linear regression models considered in this project. In this section, we will compare our best model with two penalized regression models as well as a Random Forest model. 

```{r PENALIZED-LASSO, echo=FALSE, fig.height=4, fig.width=5, fig.cap="\\label{fig:PENALIZED_lasso}Variations of the mean-squared error as a function of $\\log(\\lambda)$ for L1-penalized (LASSO) regression. The selected tuning parameter is $\\lambda = 0.38$."}
# Penalized regression
#
# The data frame defined below is used for penalized regression. We have added three more columns with quadratic values to accomodate proper comparison with our best linear regression model
data_by_day2 <- rawdata %>%  group_by(Date) %>% summarise(T = sum(T_dt)/sum(dt), # average daily temp
            P   = sum(P_dt)/sum(dt), # average daily pressure
            H   = sum(H_dt)/sum(dt), # average daily humidity
            H2  = (sum(H_dt)/sum(dt))^2, # average daily humidity^2
            WD  = sum(WD_dt)/sum(dt), # average daily wind direction
            WD2 = (sum(WD_dt)/sum(dt))^2, # average daily wind direction^2
            WS  = sum(WS_dt)/sum(dt), # average daily wind speed
            DayRatio  = mean(Daylight)/(24*3600), # ratio of daylight over total day length
            DayRatio2 = (mean(Daylight)/(24*3600))^2, # ratio of (daylight over total day length)^2
            R = sum(R_dt)/sum(dt)) # average daily radiation
mydata2 <- data.frame(data_by_day2)
X2 <- mydata2[,-c(1,11)]
X2y <- as.data.frame(cbind(X2,R))
X2yTr <- X2y[indTr,]
X2yTe <- X2y[indTe,]
X2Tr <- as.matrix.data.frame(X2yTr[,1:(p+3)])
X2Te <- as.matrix.data.frame(X2yTe[,1:(p+3)])
#
#LASSO, alpha=1
ans_lasso <- glmnet(X2Tr, yTr, alpha=1) # fit a GLM with lasso
#plot(ans_lasso,label=TRUE, lwd=2) #coefficient path
ansCV_lasso <- cv.glmnet(X2Tr, yTr, alpha=1) # Cross-validation for glmnet and returns a value for lambda
plot(ansCV_lasso)
#title(sub=paste0("LASSO, lambdaHat = ", round(ansCV_lasso$lambda.1se, 2)))
#coef(ansCV_lasso)
yH_LASSO_Te <- as.vector(predict(ans_lasso, newx=X2Te, type="response", s=ansCV_lasso$lambda.1se))
RMSELasso_Te <- rmse(yTe, yH_LASSO_Te)
yH_LASSO_Tr <- as.vector(predict(ans_lasso, newx=X2Tr, type="response", s=ansCV_lasso$lambda.1se))
RMSELasso_Tr <- rmse(yTr, yH_LASSO_Tr)
```

In order to be able to make proper comparison between penalized regression models and the best model, we needed to redefine the design matrix to include the three quadratic variables for humidity, wind direction and day ratio. For penalized regressions, we used *glmnet* package which, by default, implements regularized 10-fold cross-validation to select the tuning parameter $\lambda$. Figure \ref{fig:PENALIZED_lasso} shows the variation of mean-squared error as a function of $\lambda$ for L1-penalized (LASSO) regression. The tuning parameter $\lambda=$ `r round(ansCV_lasso$lambda.1se, 2)` is chosen based on the one-standard-deviation rule, that is the simplest model with average mean-squared error within one standard deviation of the lowest mean-squared error is selected. Based on the one-standard-deviation rule, a model with eight predictors is chosen (quadratic term in day ratio is removed).

```{r PENALIZED-Ridge, echo=FALSE, fig.height=4, fig.width=5, fig.cap="\\label{fig:PENALIZED_ridge}Variations of mean-squared errors as a function of $\\log(\\lambda)$ for L2-penalized (Ridge) regression. The selected tuning parameter is $\\lambda = 70.46$."}
#
# Ridge, alpha=0
ans_ridge <- glmnet(X2Tr, yTr, alpha=0) # fit a GLM with ridge
#plot(ans_ridge,label=TRUE, lwd=2)#coefficient path 
ansCV_ridge <- cv.glmnet(X2Tr, yTr, alpha=0) # Cross-validation for glmnet and returns a value for lambda
plot(ansCV_ridge)
#title(sub=paste0("RReg, $\\lambda$ = ", round(ansCV_ridge$lambda.1se, 2)))
#coef(ansCV_ridge)
yH_ridge_Te <- as.vector(predict(ans_ridge, newx=X2Te, type="response", s=ansCV_ridge$lambda.1se))
rmseRidge_Te <- rmse(yTe, yH_ridge_Te)
yH_ridge_Tr <- as.vector(predict(ans_ridge, newx=X2Tr, type="response", s=ansCV_ridge$lambda.1se))
rmseRidge_Tr <- rmse(yTr, yH_ridge_Tr)
```

The results of L2-penalized (Ridge) regression is shown in Figure \ref{fig:PENALIZED_ridge}. The tuning parameter chosen based on the one-standard-deviation rule is $\lambda=$ `r round(ansCV_ridge$lambda.1se, 2)` which corresponds to a model with all the nine predictors (six linear terms plus three quadratic terms) retained in the model.

```{r RANDOM-FOREST, echo=FALSE}
ansRF <- randomForest::randomForest(R ~ ., data=XyTr, ntree=1000)
yH_RF_Te <- predict(ansRF, newdata=XyTe)
rmseRF_Te <- rmse(yTe, yH_RF_Te)
yH_RF_Tr <- predict(ansRF, newdata=XyTr)
rmseRF_Tr <- rmse(yTr, yH_RF_Tr)
```

\newpage

As the last step, we use one of the machine learning technique, the Random Forest method. The results of the root-mean-square errors (RMSE) for all the above methods are summarized in Table \ref{tab:rmses}. As can be seen, the best model (revised model 3) has the best performance on prediction of test dataset. However, as expected, Random Forest shows better performance on the train dataset. The best model also performs much better than both L1/L2 penalized regression models on both test and train datasets. The model selected by the backward stagewise method is identical to the best model and thus their performances are also identical.


```{r COMPARISON, echo=FALSE, results="asis"}
# OLS (best model 3)
yH_OLS_Te <- predict(ans3, newdata=XyTe)
rmseOLS_Te <- rmse(yTe, yH_OLS_Te) # RMSE for the test data
#RMSE for the training data
yH_OLS_Tr <- ans3$fitted.values
rmseOLS_Tr <- rmse(yTr, yH_OLS_Tr) # RMSE for the train data
#backward stagewise
ansStep <- step(ans3, trace=0)
yH <- predict(ansStep, newdata=XyTe)
rmseStep_Te <- rmse(yTe, yH)
yH <- predict(ansStep, newdata=XyTr)
rmseStep_Tr <- rmse(yTr, yH)
#
rmses <- matrix(c(rmseOLS_Te, rmseStep_Te, rmseRidge_Te, RMSELasso_Te, rmseRF_Te, rmseOLS_Tr, rmseStep_Tr, rmseRidge_Tr, RMSELasso_Tr, rmseRF_Tr), ncol=5, byrow=TRUE)
colnames(rmses) <- c("Best Model", "BackStage", "RidgeReg", "LASSO", "Random Forest")
rownames(rmses) <- c("test", "train")
#rmses
stargazer(rmses, type="latex", title="Summary of the root-mean-square errors.", header=FALSE, label = "tab:rmses")
```

\newpage

The last plot (Figure \ref{fig:prediction}) shows the values predicted for the radiation by the best model and their corresponding actual values from the test dataset. The line $y=x$ is added to assist with analyzing the performance of the model. As can be seen, the points are clustered around the line $y=x$ which further shows that the best multiple linear regression model we built in this project is actually useful for prediction purposes on a new dataset (test dataset).

```{r COMPARISON-test-train, echo=FALSE, fig.height=4.6, fig.width=4, fig.cap="\\label{fig:prediction}Actual response variable values of the test dataset versus their corresponding values predicted by the revised model 3 (the best model)."}
plot(yH_OLS_Te, yTe, xlab = "Predicted values", ylab = "Actual values", xlim = c(0, 400), ylim = c(0, 400))
abline(0,1,col="red")
```



# Conclusions

The main goal of this project was to build a model using multiple linear regression for prediction of solar radiation as a function of some basic meteorological variables as well as information about the length of a day. We used the dataset that was collected  recently (less than a year ago) by NASA HI-SEAS mission in a period of four months on an isolated site located on the Big  Island of Hawaii. The data was made available recently and thus is considered as a live dataset. 

Through a rigorous analysis, we constructed a model using multiple linear regression. We discovered that temperature is a confounding variable as it is strongly correlated with the output variable (radiation) and all the predictors excluding wind variables. This causes a spurious association between temperature and radiation. This is in agreement with the findings of @Sun15 who suggested that day ratio is a much more important variable than temperature as one may intuitively consider it as an important variable. We found that the best model has the form `r as.character(ans3$call)[2]` with coefficient of determination $R^2=$ `r round(100*summary(ans3)$r.squared, 1)`% and $R^2=$ `r round(100*summary(ans_complete)$r.squared, 1)`% on train and complete datasets, respectively. The $R^2$ of the best model is quite significant which shows that this model explains large proportion of variability of the data. 

We performed various diagnostic tests to check the validity of the model as well as to verify whether the model obeys the standard assumptions of linear regression such as normality, constant variance and no serial correlation. In particular, we implemented Wilk-Shapiro and Jarque-Bera tests for normality as well as Durbin-Watson test for serial correlation. All of these tests showed that the dataset and the best model complied really well with the linear regression assumptions. 

In order to arrive at the best model, we used many variable selection methods such as, iterative model building, stepwise/stagewise and best subset regression using AIC/BIC or repeated hold cross-validation, and L1/L2 penalized regression. In addition, we used Random Forest method as a benchmark for comparison purposes. We found that the best model has the superior prediction performance (lowest RMSE) on the test dataset against all other methods. However, as expected, the Random Forest has the best prediction performance on the train dataset. It was also shown that the best model performs better on both test and train datasets when compared with L1-penalized (LASSO) and L2-penalized (Ridge) regression models.  

We had no reason to drop wind variables (wind direction and wind speed) in the best model from statistical point of view, however we suspect that these variables may be just manifestation of some lurking variables. It is hard to imagine that radiation and wind direction or wind speed are causally related. For instance, a high wind speed could be an indication of a sunny sky. If the dataset can be extended to at least one full year and/or more input variables such as cloud cover, precipitation and solar zenith angle of the sun can be considered, more accurate model may be constructed. It may be also useful to build a model that can provide accurate hourly average solar radiation prediction based on the explanatory variables using time series regression.  

\section*{References}

\bibliography{MyReferences}
