library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(lubridate, warn.conflicts = FALSE)
library(jtools)
library(corrplot)
## corrplot 0.90 loaded
The data for this analysis is downloaded from https://archive.ics.uci.edu/ml/machine-learning-databases/00275/
The data is first read via the following code chunk
bsd <- read.csv('Bike-Sharing-Dataset/day.csv')
The summary for the bike sharing data is given below.
summary(bsd)
## instant dteday season yr
## Min. : 1.0 Length:731 Min. :1.000 Min. :0.0000
## 1st Qu.:183.5 Class :character 1st Qu.:2.000 1st Qu.:0.0000
## Median :366.0 Mode :character Median :3.000 Median :1.0000
## Mean :366.0 Mean :2.497 Mean :0.5007
## 3rd Qu.:548.5 3rd Qu.:3.000 3rd Qu.:1.0000
## Max. :731.0 Max. :4.000 Max. :1.0000
## mnth holiday weekday workingday
## Min. : 1.00 Min. :0.00000 Min. :0.000 Min. :0.000
## 1st Qu.: 4.00 1st Qu.:0.00000 1st Qu.:1.000 1st Qu.:0.000
## Median : 7.00 Median :0.00000 Median :3.000 Median :1.000
## Mean : 6.52 Mean :0.02873 Mean :2.997 Mean :0.684
## 3rd Qu.:10.00 3rd Qu.:0.00000 3rd Qu.:5.000 3rd Qu.:1.000
## Max. :12.00 Max. :1.00000 Max. :6.000 Max. :1.000
## weathersit temp atemp hum
## Min. :1.000 Min. :0.05913 Min. :0.07907 Min. :0.0000
## 1st Qu.:1.000 1st Qu.:0.33708 1st Qu.:0.33784 1st Qu.:0.5200
## Median :1.000 Median :0.49833 Median :0.48673 Median :0.6267
## Mean :1.395 Mean :0.49538 Mean :0.47435 Mean :0.6279
## 3rd Qu.:2.000 3rd Qu.:0.65542 3rd Qu.:0.60860 3rd Qu.:0.7302
## Max. :3.000 Max. :0.86167 Max. :0.84090 Max. :0.9725
## windspeed casual registered cnt
## Min. :0.02239 Min. : 2.0 Min. : 20 Min. : 22
## 1st Qu.:0.13495 1st Qu.: 315.5 1st Qu.:2497 1st Qu.:3152
## Median :0.18097 Median : 713.0 Median :3662 Median :4548
## Mean :0.19049 Mean : 848.2 Mean :3656 Mean :4504
## 3rd Qu.:0.23321 3rd Qu.:1096.0 3rd Qu.:4776 3rd Qu.:5956
## Max. :0.50746 Max. :3410.0 Max. :6946 Max. :8714
sum(is.na(bsd))
## [1] 0
As we can see there are no missing values in the dataset.
str(bsd)
## 'data.frame': 731 obs. of 16 variables:
## $ instant : int 1 2 3 4 5 6 7 8 9 10 ...
## $ dteday : chr "2011-01-01" "2011-01-02" "2011-01-03" "2011-01-04" ...
## $ season : int 1 1 1 1 1 1 1 1 1 1 ...
## $ yr : int 0 0 0 0 0 0 0 0 0 0 ...
## $ mnth : int 1 1 1 1 1 1 1 1 1 1 ...
## $ holiday : int 0 0 0 0 0 0 0 0 0 0 ...
## $ weekday : int 6 0 1 2 3 4 5 6 0 1 ...
## $ workingday: int 0 0 1 1 1 1 1 0 0 1 ...
## $ weathersit: int 2 2 1 1 1 1 2 2 1 1 ...
## $ temp : num 0.344 0.363 0.196 0.2 0.227 ...
## $ atemp : num 0.364 0.354 0.189 0.212 0.229 ...
## $ hum : num 0.806 0.696 0.437 0.59 0.437 ...
## $ windspeed : num 0.16 0.249 0.248 0.16 0.187 ...
## $ casual : int 331 131 120 108 82 88 148 68 54 41 ...
## $ registered: int 654 670 1229 1454 1518 1518 1362 891 768 1280 ...
## $ cnt : int 985 801 1349 1562 1600 1606 1510 959 822 1321 ...
Using the lubridate package to first parse the column into year,month,day(ymd) and then using the month function of lubridate package along with the label argument to return month names and finally storing these month names into a new column named month_name of the dataset.
Finally, the column month_name is converted from factor to character type using as.character command.
bsd$month_name <- month( ymd(bsd[,c("dteday")]), label=TRUE )
bsd$month_name <- as.character(bsd[,c("month_name")])
# Below code can be uncommented to verify whether the column has been converted to character type
# str(bsd)
model1 <- lm(cnt ~ month_name, data = bsd)
summary(model1)
##
## Call:
## lm(formula = cnt ~ month_name, data = bsd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5177.2 -1095.2 -249.3 1290.0 4669.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4484.9 196.7 22.799 < 2e-16 ***
## month_nameAug 1179.5 275.9 4.275 2.17e-05 ***
## month_nameDec -1081.1 275.9 -3.918 9.79e-05 ***
## month_nameFeb -1829.6 281.8 -6.492 1.58e-10 ***
## month_nameJan -2308.6 275.9 -8.366 3.09e-16 ***
## month_nameJul 1078.8 275.9 3.909 0.000101 ***
## month_nameJun 1287.5 278.2 4.628 4.38e-06 ***
## month_nameMar -792.6 275.9 -2.873 0.004192 **
## month_nameMay 864.9 275.9 3.134 0.001793 **
## month_nameNov -237.7 278.2 -0.854 0.393113
## month_nameOct 714.3 275.9 2.589 0.009829 **
## month_nameSep 1281.6 278.2 4.607 4.83e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1524 on 719 degrees of freedom
## Multiple R-squared: 0.3906, Adjusted R-squared: 0.3813
## F-statistic: 41.9 on 11 and 719 DF, p-value: < 2.2e-16
As we can observe the R-squared in our case is 0.3906 which means that ~39% of the variation in cnt is explained by the month name variable. R-squared basically is the square of the correlation coefficient and informs us how much percent of the variation of dependent variable is explained by the indepenedent variable.
It can be observed that the month names are appearing in lexicographic order. So it is easy to identify that month of April is set as the reference. This can be confirmed also by going through the list and finding the missing month from the 12 months.
In our case, the intercept represents the estimate of cnt for the month of April, and the coefficient estimates for the other variables represent the difference between that month from the month of April.
Hence, 4484.9 is the reference month’s prediction.
2176.339 for the month of January and 5772.367 for the month of June
newObservations = data.frame(month_name=c('Jan', 'Jun'))
newObservations$predictedValues = predict(model1, newObservations)
# model1$fitted.values
There can be multiple ways of seeing the predicted values for the month of January and june, all of which arrive at the same answer.
The first one is theoretical where we can see that our model is cnt = 4484.9 + (-2308.6)(1 if month if june, 0 otherwise) + (coefficient of Feb)(1 if month if june, 0 otherwise) + ……
So basically the above linear model uses the coefficients that we find out by fitting our regression model. So the prediction for the month of January would be cnt = 4484 - 2308.6 = 2176.339 Similarly, for the month of June, it would be cnt = 4484 + 1287.5 = 5772.367
This can also be achieved with the predict() function as used in the code chunk above or the fitted.values column of the model1.
model2 = lm(cnt ~ month_name + temp, data = bsd)
summary(model2)
##
## Call:
## lm(formula = cnt ~ month_name + temp, data = bsd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4896.6 -1080.0 -228.4 1245.2 3372.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1554.39 390.76 3.978 7.66e-05 ***
## month_nameAug -308.08 315.42 -0.977 0.3290
## month_nameDec -170.96 283.80 -0.602 0.5471
## month_nameFeb -764.81 296.15 -2.582 0.0100 *
## month_nameJan -852.31 313.41 -2.719 0.0067 **
## month_nameJul -701.18 335.50 -2.090 0.0370 *
## month_nameJun -47.47 307.78 -0.154 0.8775
## month_nameMar -297.20 269.38 -1.103 0.2703
## month_nameMay 86.73 278.37 0.312 0.7555
## month_nameNov 390.66 275.22 1.419 0.1562
## month_nameOct 620.72 263.30 2.357 0.0187 *
## month_nameSep 368.25 285.93 1.288 0.1982
## temp 6235.14 729.40 8.548 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1453 on 718 degrees of freedom
## Multiple R-squared: 0.4469, Adjusted R-squared: 0.4377
## F-statistic: 48.35 on 12 and 718 DF, p-value: < 2.2e-16
First, we notice that the R-squared for model2 is .4469, which is higher than the R-squared of 0.3906 in model1. This shows that a lot more variance of cnt variable is explained by our current model since it combines two variables month_name and temp to perform the prediction rather than using just a single variable to predict cnt. This means that we can have a little more confidence in the accuracy of these predictions than in the predictions from model1.
tempmodel = lm(cnt ~ temp, data = bsd)
summary(tempmodel)
##
## Call:
## lm(formula = cnt ~ temp, data = bsd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4615.3 -1134.9 -104.4 1044.3 3737.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1214.6 161.2 7.537 1.43e-13 ***
## temp 6640.7 305.2 21.759 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1509 on 729 degrees of freedom
## Multiple R-squared: 0.3937, Adjusted R-squared: 0.3929
## F-statistic: 473.5 on 1 and 729 DF, p-value: < 2.2e-16
Also, we can observe with the help of the above code chunk that the R-squared is not simply a sum of the R-squared from the simple regression models regressed using month_name or temp. This could be because month_name and temp are correlated to each other.
export_summs(model1, model2)
| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 4484.90 *** | 1554.39 *** |
| (196.71) | (390.76) | |
| month_nameAug | 1179.52 *** | -308.08 |
| (275.94) | (315.42) | |
| month_nameDec | -1081.09 *** | -170.96 |
| (275.94) | (283.80) | |
| month_nameFeb | -1829.60 *** | -764.81 * |
| (281.83) | (296.15) | |
| month_nameJan | -2308.56 *** | -852.31 ** |
| (275.94) | (313.41) | |
| month_nameJul | 1078.78 *** | -701.18 * |
| (275.94) | (335.50) | |
| month_nameJun | 1287.47 *** | -47.47 |
| (278.19) | (307.78) | |
| month_nameMar | -792.64 ** | -297.20 |
| (275.94) | (269.38) | |
| month_nameMay | 864.87 ** | 86.73 |
| (275.94) | (278.37) | |
| month_nameNov | -237.72 | 390.66 |
| (278.19) | (275.22) | |
| month_nameOct | 714.33 ** | 620.72 * |
| (275.94) | (263.30) | |
| month_nameSep | 1281.62 *** | 368.25 |
| (278.19) | (285.93) | |
| temp | 6235.14 *** | |
| (729.40) | ||
| N | 731 | 731 |
| R2 | 0.39 | 0.45 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||
It can be seen that the coefficient estimated for the month_nameJan is -2308.56 for model1 whereas it increase to -852.31 in the case of model2.
This could be because month_name and temp are correlated with each other, and the coefficient estimates represent the unique effect of each predictor variable after considering the effect(s) of the other predictor variables.
newObservations2 = data.frame(month_name=c('Jan'), temp=0.25)
newObservations2$predictedValues = predict(model2, newObservations2)
The predicted value for count is 2260.863 when the month is January and the temp is 0.25