---
title: "Bike Sharing - Data Analysis"
author: "Prof. Radhendushka Srivastava, Department of Mathematics, IIT Bombay"
date: "05/03/2022"
output:
  html_document: default
---

# Loading libraries
```{r,warning=F}
suppressPackageStartupMessages(library(tidyverse))
```

# Reading data
```{r}

# Spoken Tutorials:
# 5. Introduction to Data Frames in R
# 9. Indexing and Slicing Data Frames

# Source: https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset
bike_sharing <- read_csv("day.csv")
bike_sharing <- bike_sharing[,-c(1:8,11,14:15)]
bike_sharing <- bike_sharing %>% filter(weathersit!=3)
```

## 1) Basic statistics of continuous variable
```{r}
# Mean, standard deviation, minimum, maximum, range, quartiles, IQR, skewness and kurtosis
statistics <- function(x)
{
  mean <- mean(x)
  standard_deviation <- sd(x)
  min <- min(x)
  max <- max(x)
  range <- max - min
  first_quartile <- quantile(x)[2]
  third_quartile <- quantile(x)[4]
  IQR <- third_quartile - first_quartile
  skewness <- mean((x-mean(x))^3)/(sd(x)^3)
  kurtosis <- (mean((x-mean(x))^4))/sd(x)^4
  
  return(cat("\t","Mean: ",mean,"\n\t","Standard deviation: ",standard_deviation,"\n\t","Minimum value: ",min,"\n\t","Maximum value: ",max,"\n\t","Range: ",range,"\n\t","First Quartile: ",first_quartile,"\n\t","Third Quartile: ",third_quartile,"\n\t","Interquartile Range: ",IQR,"\n\t","Skewness: ",skewness,"\n\t","Kurtosis: ",kurtosis))
}
```

### 1.1) Basic statistics of the variable "temp" (Normalized temperature in Celsius)
```{r}
statistics(bike_sharing$temp)
```

### 1.2) Basic statistics of the variable "hum" (Normalized humidity)
```{r}
statistics(bike_sharing$hum)
```

### 1.3) Basic statistics of the variable "windspeed" (Normalized wind speed)
```{r}
statistics(bike_sharing$windspeed)
```
 
### 1.4) Basic statistics of the variable "cnt" (Count of total rental bikes)
```{r}
statistics(bike_sharing$cnt)
```

## 2) Histogram (Probability density estimate)
### 2.1) Histogram of the variable "cnt" (Count of total rental bikes)
```{r}

# Spoken Tutorials:
# 15. Plotting Histograms and Pie Chart 

# Distribution of "cnt" variables
hist(bike_sharing$cnt,breaks = 15,col="green",border="black",prob = TRUE,xlab = "",main = "Count of total rental bikes")
x <- seq(min(bike_sharing$cnt),max(bike_sharing$cnt),length=40)
y <- dnorm(x,mean=mean(bike_sharing$cnt),sd=sd(bike_sharing$cnt))
lines(x, y, col="red", lwd=2)
```

## 3) Correlation matrix among continuous variables
```{r}
cor(bike_sharing[,-1])
```

## 4) Scatter plot and linear regression
### 4.1.1) Scatter plot between "temp" (Normalized temperature in Celsius) and "cnt" (Count of total rental bikes) 
```{r}

# Spoken Tutorials:
# 17. Introduction to ggplot2
# 18. Aesthetic Mapping in ggplot2

# "cnt" and "temp" variables
# Scatter plot
ggplot(data = bike_sharing,aes(y=cnt,x=temp)) + geom_point() + labs(x = "Normalized temperature in Celsius", y = "Count of total rental bikes") + theme_bw()
```

### 4.1.2) Linear regression of  "cnt" (Count of total rental bikes) on "temp" (Normalized temperature in Celsius)  
```{r}
cnt_temp_lm <- lm(cnt~temp,data = bike_sharing)
summary(cnt_temp_lm)
```

### 4.1.3) Scatter plot of "temp" (Normalized temperature in Celsius) and "cnt" (Count of total rental bikes) with regression line
```{r}
ggplot(data = bike_sharing,aes(y=cnt,x=temp)) + geom_point() + labs(x = "Normalized temperature in Celsius", y = "Count of total rental bikes") + theme_bw() + geom_smooth(method='lm', formula = y~x)
```

### 4.1.4) Residual Analysis
```{r}
# Plotting squared residuals (Homogeneity Validation)
plot(y=cnt_temp_lm$residuals^2,x=as.numeric(names(cnt_temp_lm$residuals^2)),ylab="Squared Residuals",xlab="Index")
# qq-plot of residuals (Normality Validation)
qqnorm(rstandard(cnt_temp_lm), ylab="Standardized Residuals",xlab="Normal Scores")
qqline(rstandard(cnt_temp_lm))
# Histogram of residuals (Normality Validation)
hist(cnt_temp_lm$residuals,main = "Histogram of Residuals",xlab="",prob=T,ylim=c(0,0.0003))
x <- seq(min(cnt_temp_lm$residuals),max(cnt_temp_lm$residuals),length=40)
y <- dnorm(x,mean=mean(cnt_temp_lm$residuals),sd=sd(cnt_temp_lm$residuals))
lines(x, y, col="red", lwd=2)
```

### 4.2.1) Scatter plot of "hum" (Normalized humidity) and "cnt" (Count of total rental bikes) 
```{r}
ggplot(data = bike_sharing,aes(y=cnt,x=hum)) + geom_point() + labs(x = "Normalized humidity", y = "Count of total rental bikes") + theme_bw()
```

### 4.2.2) Linear regression of  "cnt" (Count of total rental bikes) on "hum" (Normalized humidity)
```{r}
cnt_hum_lm <- lm(cnt~hum,data = bike_sharing)
summary(cnt_hum_lm)
```

### 4.2.3) Scatter plot of "hum" (Normalized humidity) and "cnt" (Count of total rental bikes) with regression line
```{r}
ggplot(data = bike_sharing,aes(y=cnt,x=hum)) + geom_point() + labs(x = "Normalized humidity", y = "Count of total rental bikes") + theme_bw() + geom_smooth(method='lm', formula = y~x)
```

### 4.2.4) Residual Analysis
```{r}
# Plot of squared residuals (Homogeneity Validation)
plot(y=cnt_hum_lm$residuals^2,x=as.numeric(names(cnt_hum_lm$residuals^2)),ylab="Squared Residuals",xlab="Index")
# Normal qq-plot (Normality Validation)
qqnorm(rstandard(cnt_hum_lm), ylab="Standardized Residuals",xlab="Normal Scores")
qqline(rstandard(cnt_hum_lm))
# Histogram of residual (Normality Validation)
hist(cnt_hum_lm$residuals,main = "Histogram of Residuals",xlab="",prob=T)
x <- seq(min(cnt_hum_lm$residuals),max(cnt_hum_lm$residuals),length=40)
y <- dnorm(x,mean=mean(cnt_hum_lm$residuals),sd=sd(cnt_hum_lm$residuals))
lines(x, y, col="red", lwd=2)
```

### 4.3.1) Scatter plot of "windspeed" (Wind speed) and "cnt" (Count of total rental bikes) 
```{r}
ggplot(data = bike_sharing,aes(y=cnt,x=windspeed)) + geom_point() + labs(x = "Normalized wind speed", y = "Count of total rental bikes") + theme_bw()
```

### 4.3.2) Linear regression of "cnt" (Count of total rental bikes) on "windspeed" (Wind speed)  
```{r}
cnt_windspeed_lm <- lm(cnt~windspeed,data = bike_sharing)
summary(cnt_windspeed_lm)
```

### 4.3.3) Scatter plot of "windspeed" (Wind Speed) and "cnt" (Count of total rental bikes) with regression line
```{r}
ggplot(data = bike_sharing,aes(y=cnt,x=windspeed)) + geom_point() + labs(x = "Normalized wind speed", y = "Count of total rental bikes") + theme_bw() + geom_smooth(method='lm', formula = y~x)
```

### 4.3.4) Residual Analysis
```{r}
# Plot of squared residuals (Homogeneity Validation)
plot(y=cnt_windspeed_lm$residuals^2,x=as.numeric(names(cnt_windspeed_lm$residuals^2)),ylab="Squared Residuals",xlab="Index")
# Normal qq-plot (Normality Validation)
qqnorm(rstandard(cnt_windspeed_lm), ylab="Standardized Residuals",xlab="Normal Scores")
qqline(rstandard(cnt_windspeed_lm))
# Histogram of residual (Normality Validation)
hist(cnt_windspeed_lm$residuals,main = "Histogram of Residuals",xlab="",prob=T)
x <- seq(min(cnt_windspeed_lm$residuals),max(cnt_windspeed_lm$residuals),length=40)
y <- dnorm(x,mean=mean(cnt_windspeed_lm$residuals),sd=sd(cnt_windspeed_lm$residuals))
lines(x, y, col="red", lwd=2)
```

## 5) Multiple regression model
### 5.1) Multiple regression model of "cnt" (Count of total rental bikes) on ["temp" (Normalized temperature) + "hum" (Normalized humidity) + "windspeed" (Wind Speed) + "weathersit" (Weather situation)]
```{r}
multi_reg <- lm(cnt~temp+hum+windspeed+weathersit,data = bike_sharing)
summary(multi_reg)
```

### 5.2) Residual Analysis
```{r}
# Plot of squared residual (Homogeneity Validation)
plot(y=multi_reg$residuals^2,x=as.numeric(names(multi_reg$residuals^2)),ylab="Squared Residuals",xlab="Index")
# Normal qq-plot (Normality Validation)
qqnorm(rstandard(multi_reg), ylab="Standardized Residuals",xlab="Normal Scores")
qqline(rstandard(multi_reg))
# Histogram of residual (Normality Validation)
hist(multi_reg$residuals,main = "Histogram of Residuals",xlab="",prob=T,ylim=c(0,0.0003))
x <- seq(min(multi_reg$residuals),max(multi_reg$residuals),length=40)
y <- dnorm(x,mean=mean(multi_reg$residuals),sd=sd(multi_reg$residuals))
lines(x, y, col="red", lwd=2)
```