London Bikes data

We’ll use data from http://www.tfl.gov.uk to analyse usage of the London Bike Sharing scheme. This data has already been downloaded for you and exists in a CSV (Comma Separated Values) file that you have to read in to R.

There is no dropdown menu to read in your data to R, so instead we use functions like read_csv to load data from file into R objects.

#read the CSV file
bike <- read_csv(here::here("data", "london_bikes.csv")) %>% 
  mutate(weekend = if_else((wday == "Sat" | wday == "Sun"), TRUE, FALSE))

Cleaning our data

Sometimes our data needs a bit of ‘cleaning’. For instance, day_of_week is variable type character, or chr. We should, however, treat it as a categorical, or factor variable and relevel it, so Monday is the first level of the factor (or first day of week), etc.

R is fairly sensitive with dates. When you read a CSV file, the date may be in different formats. For instance, Christmas 2017 could be input as 12-25-2017, 25.12.2017, 25 Dec 2017, Dec 25, 2017, etc. To be consistent, we use lubridate’s ymd function, and force variable Day to be a date in the format YYYY-MM-DD

Finally, we can turn season from 1, 2, 3, 4, to words like Winter, Spring, etc.

We will be talking more about data wrangling later, but for now just execute the following piece of code.

# fix dates using lubridate, and generate new variables for year, month, month_name, day, and day_of _week
bike <- bike %>%   
  mutate(
    year=year(date),
    month = month(date),
    month_name=month(date, label = TRUE),
    day_of_week = wday(date, label = TRUE)) 

# generate new variable season_name to turn seasons from numbers to Winter, Spring, etc
bike <- bike %>%  
  mutate(
    season_name = case_when(
      month_name %in%  c("Dec", "Jan", "Feb")  ~ "Winter",
      month_name %in%  c("Mar", "Apr", "May")  ~ "Spring",
      month_name %in%  c("Jun", "Jul", "Aug")  ~ "Summer",
      month_name %in%  c("Sep", "Oct", "Nov")  ~ "Autumn",
    ),
    season_name = factor(season_name, 
                         levels = c("Winter", "Spring", "Summer", "Autumn")))
    
#examine the structure of the datafame
skim(bike)
Data summary
Name bike
Number of rows 5268
Number of columns 41
_______________________
Column type frequency:
character 6
factor 3
logical 2
numeric 27
POSIXct 3
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
preciptype 1411 0.73 4 9 0 3 0
conditions 0 1.00 4 28 0 11 0
description 0 1.00 26 85 0 78 0
icon 0 1.00 4 17 0 6 0
wday 0 1.00 3 3 0 7 0
set 0 1.00 5 5 0 1 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
month_name 0 1 TRUE 12 Aug: 465, Oct: 465, Dec: 464, Sep: 450
day_of_week 0 1 TRUE 7 Sun: 753, Mon: 753, Fri: 753, Sat: 753
season_name 0 1 FALSE 4 Aut: 1365, Sum: 1321, Win: 1294, Spr: 1288

Variable type: logical

skim_variable n_missing complete_rate mean count
severerisk 5268 0 NaN :
weekend 0 1 0.29 FAL: 3762, TRU: 1506

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
bikes_hired 0 1.00 26327.60 9496.92 0.0 19746.75 26070.00 32745.00 73094.00 ▂▇▅▁▁
tempmax 0 1.00 14.73 6.36 -2.6 10.10 14.40 19.50 36.10 ▁▇▇▃▁
tempmin 0 1.00 7.84 5.04 -11.1 4.10 8.00 11.60 22.10 ▁▃▇▇▁
temp 0 1.00 11.23 5.48 -5.9 7.20 11.10 15.50 28.30 ▁▅▇▅▁
feelslikemax 0 1.00 14.09 7.24 -7.2 10.10 14.40 19.50 38.20 ▁▅▇▃▁
feelslikemin 0 1.00 6.18 6.38 -14.6 1.17 6.00 11.60 22.10 ▁▃▇▇▂
feelslike 0 1.00 10.09 6.67 -9.1 4.90 10.60 15.50 28.90 ▁▆▇▇▁
dew 0 1.00 7.41 4.67 -8.0 4.10 7.60 11.00 19.00 ▁▃▇▇▂
humidity 0 1.00 79.28 10.40 39.6 72.20 80.30 87.30 99.90 ▁▂▅▇▅
precip 0 1.00 1.64 4.73 0.0 0.00 0.18 1.34 168.20 ▇▁▁▁▁
precipprob 0 1.00 73.04 44.38 0.0 0.00 100.00 100.00 100.00 ▃▁▁▁▇
precipcover 0 1.00 8.96 11.60 0.0 0.00 8.33 12.50 100.00 ▇▁▁▁▁
snow 0 1.00 0.02 0.25 0.0 0.00 0.00 0.00 8.40 ▇▁▁▁▁
snowdepth 0 1.00 0.05 0.55 0.0 0.00 0.00 0.00 18.10 ▇▁▁▁▁
windgust 1126 0.79 42.08 14.84 9.4 30.72 40.70 51.80 120.90 ▅▇▃▁▁
windspeed 0 1.00 22.95 8.42 5.1 16.88 21.80 27.60 72.20 ▅▇▂▁▁
winddir 0 1.00 196.97 92.33 0.1 127.75 221.60 259.00 359.80 ▃▂▃▇▃
sealevelpressure 0 1.00 1014.90 10.52 966.6 1008.70 1015.70 1022.00 1047.80 ▁▁▇▇▁
cloudcover 0 1.00 60.86 22.45 0.0 46.60 62.70 77.80 100.00 ▁▃▇▇▅
visibility 0 1.00 22.65 10.62 0.2 15.20 22.00 28.90 62.50 ▃▇▆▁▁
solarradiation 0 1.00 112.00 84.89 0.0 40.20 91.50 169.77 358.60 ▇▅▃▂▁
solarenergy 0 1.00 9.66 7.34 0.0 3.50 7.90 14.60 31.00 ▇▅▃▂▁
uvindex 0 1.00 4.26 2.55 0.0 2.00 4.00 6.00 10.00 ▇▆▅▅▁
moonphase 0 1.00 0.48 0.29 0.0 0.25 0.50 0.75 0.98 ▇▇▇▇▇
month 0 1.00 6.62 3.46 1.0 4.00 7.00 10.00 12.00 ▇▅▅▅▇
week 0 1.00 27.01 15.08 1.0 14.00 27.00 40.00 53.00 ▇▇▇▇▇
year 0 1.00 2017.28 4.17 2010.0 2014.00 2017.00 2021.00 2024.00 ▆▇▇▇▇

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
date 0 1 2010-07-30 00:00:00 2024-12-30 00:00:00 2017-10-14 12:00:00 5268
sunrise 0 1 2010-07-30 05:17:39 2024-12-30 08:06:26 2017-10-14 19:21:46 5268
sunset 0 1 2010-07-30 20:50:39 2024-12-30 16:00:26 2017-10-15 06:05:26 5268

Summary Statistics

Besides the number of bikes hired each day, we also have data on the weather for that day

Having loaded and cleaned our data, we can create summary statistics using mosaic’s favstats. This is uni-variate analysis, so in our mosaic syntax we will use favstats(~ bikes_hired, data= bike). We also want to get an overall, time-series plot of bikes over time; for the latter, we just create a scatter plot of bikes_hired vs Day.

favstats(~ bikes_hired, data= bike)
minQ1medianQ3maxmeansdnmissing
01.97e+042.61e+043.27e+047.31e+042.63e+049.5e+0352680

While this analysis shows us overall statistics, what if we wanted to get summary statistics by year, day_of_week, month, or season? Mosaic’s syntax Y ~ X alows us to facet our analysis of a variable Y by another variable X using the syntax favstats( Y ~ Z, data=...)

favstats(bikes_hired ~ year, data=bike)
yearminQ1medianQ3maxmeansdnmissing
20102.76e+039.3e+03 1.4e+04 1.87e+042.8e+04 1.41e+045.62e+031550
20114.56e+031.63e+042.03e+042.37e+042.94e+041.96e+045.5e+03 3650
20123.53e+031.93e+042.62e+043.25e+044.71e+042.6e+04 9.43e+033660
20133.73e+031.76e+042.2e+04 2.74e+043.56e+042.2e+04 7.28e+033650
20144.33e+032.05e+042.77e+043.44e+044.9e+04 2.75e+049.07e+033650
20155.78e+032.21e+042.66e+043.29e+047.31e+042.7e+04 8.55e+033650
20164.89e+032.24e+042.79e+043.51e+044.66e+042.82e+048.85e+033660
20175.14e+032.41e+042.95e+043.45e+044.6e+04 2.86e+048.38e+033650
20185.86e+032.18e+042.92e+043.77e+044.61e+042.9e+04 1.02e+043650
20195.65e+032.42e+042.89e+043.45e+044.47e+042.86e+048.09e+033650
20204.87e+032.01e+042.75e+043.65e+047.02e+042.85e+041.16e+043660
20216.25e+032.17e+043.1e+04 3.82e+045.69e+043e+04       1.1e+04 3650
20220       2.51e+043.14e+044e+04       6.7e+04 3.15e+041.03e+043650
20236.72e+031.99e+042.37e+042.78e+043.53e+042.34e+046.03e+033650
20247.48e+032.04e+042.48e+042.87e+043.52e+042.4e+04 6.14e+033650
 
favstats(bikes_hired ~ day_of_week, data=bike)
day_of_weekminQ1medianQ3maxmeansdnmissing
Sun0       1.41e+042.1e+04 2.97e+046.31e+042.22e+041.05e+047530
Mon3.97e+032.06e+042.57e+043.15e+046.7e+04 2.61e+048.28e+037530
Tue3.76e+032.25e+042.79e+043.44e+046.53e+042.82e+048.62e+037520
Wed4.33e+032.28e+042.8e+04 3.42e+045.44e+042.84e+048.52e+037520
Thu5.65e+032.27e+042.75e+043.44e+047.31e+042.83e+048.79e+037520
Fri5.4e+03 2.14e+042.65e+043.29e+046.7e+04 2.7e+04 8.54e+037530
Sat0       1.6e+04 2.25e+043.08e+047.02e+042.42e+041.11e+047530
 
favstats(bikes_hired ~ month_name, data=bike)
month_nameminQ1medianQ3maxmeansdnmissing
Jan3.73e+031.39e+041.9e+04 2.33e+043.8e+04 1.87e+045.99e+034340
Feb3.53e+031.6e+04 2.06e+042.46e+045.25e+042.03e+046.37e+033960
Mar5.06e+031.76e+042.31e+042.72e+045.66e+042.27e+047.59e+034340
Apr4.87e+032.1e+04 2.59e+043.08e+044.9e+04 2.6e+04 7.61e+034200
May1.07e+042.45e+042.99e+043.6e+04 7.02e+043.04e+048.44e+034340
Jun6.06e+032.77e+043.31e+043.92e+046.53e+043.34e+048.59e+034200
Jul5.56e+032.91e+043.58e+044.13e+047.31e+043.47e+048.35e+034360
Aug4.3e+03 2.56e+043.24e+043.83e+046.7e+04 3.14e+049.6e+03 4650
Sep0       2.49e+043.12e+043.61e+045.18e+043.04e+048.11e+034500
Oct7.07e+032.33e+042.8e+04 3.26e+044.79e+042.75e+046.77e+034650
Nov6.03e+031.88e+042.37e+042.79e+044.47e+042.33e+046.41e+034500
Dec2.76e+031.14e+041.67e+042.29e+043.91e+041.71e+046.91e+034640
 
favstats(bikes_hired ~ season_name, data=bike)
season_nameminQ1medianQ3maxmeansdnmissing
Winter2.76e+031.37e+041.88e+042.35e+045.25e+041.86e+046.58e+0312940
Spring4.87e+032.08e+042.61e+043.13e+047.02e+042.64e+048.51e+0312880
Summer4.3e+03 2.76e+043.37e+043.96e+047.31e+043.31e+048.98e+0313210
Autumn0       2.19e+042.72e+043.26e+045.18e+042.7e+04 7.69e+0313650

Exploratory Data Analysis

Time series plot of bikes rented

While summary statistics allow us to quickly disover key metrics that represent the data, they are often not sufficient by themselves. Often, it is useful to represent the data graphically (or visually) to uncover information that is critical for decision making, such as trends, patterns and outliers.

In this section we will create a time-series scatter-plot that shows the number of bikes hired on each day. We will use the ggplot2 library.

Creating a basic plot is a two-step process. The first step is the specify the data frame and the axes by using the syntax: ggplot(data frame, aes(x=variable, y=variable). The second step is to specify the plot that we want. In this case, we want a scatter (point) plot using geom_point() after a + sign.

bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  
  # Create a ggplot object
  ggplot() +
  
  # Map date to x-axis and bikes_hired to y-axis
  aes(x = date, y = bikes_hired) +
  
  # Add scatter points with 40% transparency to reduce overplotting
  geom_point(alpha = 0.4) +
  
  # Apply a clean black and white theme
  theme_bw() +
  
  # NULL doesn't add anything - often used as a placeholder for easy code modification
  NULL

Further graphs

# Histogram of bikes rented
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_histogram()+
  theme_bw()+
  NULL
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# Histogram faceted by season_name
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_histogram()+
  facet_wrap(~season_name)+
  theme_bw()+
  NULL
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# Histogram faceted by month_name
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_histogram()+
  facet_wrap(~month_name)+
  theme_bw()+
  NULL
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# Histogram faceted by month_name in 4 rows
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_histogram()+
  facet_wrap(~month_name, nrow = 4)+
  theme_bw()+
  NULL
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# Density plot 
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_density()+
  theme_bw()+
  NULL

# Density plot filled by season_name 
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_density(aes(fill=season_name), alpha = 0.3)+
  theme_bw()+
  NULL

# Density plot filled by season_name, and faceted by season_name
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_histogram(aes(fill=season_name), alpha = 0.5)+
  facet_wrap(~season_name, nrow = 4)+
  theme_minimal()+
  
  #remove legend to the right
  theme(legend.position = "none")+
  NULL
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# Density plot filled by season_name, and faceted by month_name
bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x = bikes_hired)+
  geom_density(aes(fill=season_name), alpha = 0.3)+
  facet_wrap(~month_name, nrow = 4)+
  theme_bw()+
  theme(legend.position="none")+
  NULL

#Boxplot of bikes_hired by month_name

bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x=month_name, y= bikes_hired)+
  geom_boxplot()+
  theme_bw()+
  NULL

# bikes_hired vs. weather features

bike %>% 
  
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x=temp, y= bikes_hired)+
  geom_point()+
  geom_smooth(method = "lm", se=FALSE)+
  theme_bw()+
  NULL
## `geom_smooth()` using formula = 'y ~ x'

bike %>% 
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x=temp, y= bikes_hired,
       colour=season_name)+
  geom_point(alpha = 0.1)+
  geom_smooth(method = "lm", se=FALSE)+
  theme_bw()+
  NULL
## `geom_smooth()` using formula = 'y ~ x'

###### bikes on humidity
bike %>% 
  # Filter the bike dataset to include only data from January 1, 2014 onwards
  filter(date >= lubridate::ymd("2014-01-01")) %>%  
  ggplot()+
  aes(x=humidity, y= bikes_hired,
       colour=season_name)+
  geom_point(alpha = 0.1)+
  geom_smooth(method = "lm", se=FALSE)+
  theme_bw()+
  NULL
## `geom_smooth()` using formula = 'y ~ x'

Model building

Correlation - Scatterplot matrix

Besides the number of bikes rented out on a given day, we have also downloaded weather data for London such as mean_temp, humidity, pressure, precipitation, etc. that measure weather conditions on a single day. It may be the case that more bikes are rented out when it’s warmer? Or how can we estimate the effect of rain on rentals?

Your task is to build a regression model that helps you explain the number of rentals per day.

Let us select a few of these numerical variables and create a scatterplot-correlation matrix

bike %>% 
  select(cloudcover, humidity, precip, temp, feelslike, bikes_hired) %>% 
  ggpairs()

Weekend or weekdays? Is there a difference?

t.test(bikes_hired ~ weekend, data= bike)
## 
##  Welch Two Sample t-test
## 
## data:  bikes_hired by weekend
## t = 13.929, df = 2300.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group FALSE and group TRUE is not equal to 0
## 95 percent confidence interval:
##  3739.295 4964.730
## sample estimates:
## mean in group FALSE  mean in group TRUE 
##            27571.74            23219.72

Model 0: using the mean to predict bikes_hired

We start the naive model where we just use the average to predict how many bikes we are going to rent out on a single day

favstats(~bikes_hired, data = bike)
minQ1medianQ3maxmeansdnmissing
01.97e+042.61e+043.27e+047.31e+042.63e+049.5e+0352680
 
# can you create a confidence interval for mean bikes_hired? What is the SE?

model0 <- lm(bikes_hired ~ 1, data= bike)
msummary(model0)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  26327.6      130.8   201.2   <2e-16 ***
## 
## Residual standard error: 9497 on 5267 degrees of freedom
# plot actual data vs predicted data
broom::augment(model0) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

What is the regression’s residual standard error? What is the intercept standard error?

Model 1: bikes_hired on mean_temp

# Define the model
model1 <- lm(bikes_hired ~ feelslikemax, data = bike)

# look at model estimated
msummary(model1)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  14556.64     220.78   65.93   <2e-16 ***
## feelslikemax   835.44      13.94   59.94   <2e-16 ***
## 
## Residual standard error: 7323 on 5266 degrees of freedom
## Multiple R-squared:  0.4056, Adjusted R-squared:  0.4054 
## F-statistic:  3593 on 1 and 5266 DF,  p-value: < 2.2e-16
# plot actual data vs predicted data
broom::augment(model1) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

  • Is the effect of mean_temp significant? Why?
  • What proportion of the overall variability in bike rentals does temperature explain?

Model 2: bikes_hired on mean_temp plus weekend

# Define the model
model2 <- lm(bikes_hired ~ feelslikemax + weekend, data = bike)

# look at model estimated
msummary(model2)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15800.66     221.76   71.25   <2e-16 ***
## feelslikemax   834.53      13.44   62.10   <2e-16 ***
## weekendTRUE  -4306.90     215.28  -20.01   <2e-16 ***
## 
## Residual standard error: 7060 on 5265 degrees of freedom
## Multiple R-squared:  0.4476, Adjusted R-squared:  0.4473 
## F-statistic:  2133 on 2 and 5265 DF,  p-value: < 2.2e-16
# plot actual data vs predicted data
broom::augment(model2) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

  • Fit a regression model called model2 with the following explanatory variables: mean_temp and weekend

    • Is the effect of mean_temp significant? Why?
    • What proportion of the overall variability does this model explain? What is the meaning of the effect (slope) of weekendTRUE? What % of the variability of bikes_hired does your model explain?

Model 3: bikes_hired on mean_temp plus wday

# Define the model
model3 <- lm(bikes_hired ~ temp + wday, data = bike)

# look at model estimated
msummary(model3)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15310.17     335.80  45.594  < 2e-16 ***
## temp         1039.45      18.31  56.760  < 2e-16 ***
## wdayMon      -880.36     375.40  -2.345 0.019057 *  
## wdaySat     -2731.60     375.40  -7.277 3.93e-13 ***
## wdaySun     -4732.83     375.40 -12.607  < 2e-16 ***
## wdayThu      1236.42     375.53   3.293 0.000999 ***
## wdayTue      1205.90     375.52   3.211 0.001330 ** 
## wdayWed      1315.16     375.53   3.502 0.000465 ***
## 
## Residual standard error: 7284 on 5260 degrees of freedom
## Multiple R-squared:  0.4125, Adjusted R-squared:  0.4117 
## F-statistic: 527.6 on 7 and 5260 DF,  p-value: < 2.2e-16
# plot actual data vs predicted data
broom::augment(model3) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

What is the meaning of the effect (slope) of, e.g., wdayMon? What % of the variability of bikes_hired does your model explain?

model4 <- lm(bikes_hired ~ temp + wday + humidity + factor(month), data = bike)

# look at model estimated
msummary(model4)
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     38443.27    1064.54  36.113  < 2e-16 ***
## temp              537.96      31.27  17.205  < 2e-16 ***
## wdayMon          -924.22     344.07  -2.686 0.007250 ** 
## wdaySat         -2661.86     344.07  -7.736 1.22e-14 ***
## wdaySun         -4821.70     344.08 -14.013  < 2e-16 ***
## wdayThu          1308.18     344.18   3.801 0.000146 ***
## wdayTue          1311.24     344.20   3.809 0.000141 ***
## wdayWed          1440.00     344.21   4.184 2.92e-05 ***
## humidity         -256.42      11.00 -23.303  < 2e-16 ***
## factor(month)2    432.82     465.44   0.930 0.352456    
## factor(month)3   1481.66     459.98   3.221 0.001284 ** 
## factor(month)4   2045.02     485.36   4.213 2.56e-05 ***
## factor(month)5   4500.81     517.00   8.706  < 2e-16 ***
## factor(month)6   5800.74     566.18  10.245  < 2e-16 ***
## factor(month)7   5428.74     603.25   8.999  < 2e-16 ***
## factor(month)8   2778.82     588.58   4.721 2.41e-06 ***
## factor(month)9   4373.23     543.67   8.044 1.07e-15 ***
## factor(month)10  4715.45     494.41   9.538  < 2e-16 ***
## factor(month)11  3504.31     459.80   7.621 2.96e-14 ***
## factor(month)12 -1582.60     447.96  -3.533 0.000415 ***
## 
## Residual standard error: 6676 on 5248 degrees of freedom
## Multiple R-squared:  0.5076, Adjusted R-squared:  0.5059 
## F-statistic: 284.8 on 19 and 5248 DF,  p-value: < 2.2e-16
# plot actual data vs predicted data
broom::augment(model4) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

model5 <- lm(bikes_hired ~ feelslikemax + wday + humidity + factor(month) + precip, data = bike)

# look at model estimated
msummary(model5)
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     33695.66    1051.08  32.058  < 2e-16 ***
## feelslikemax      492.05      22.50  21.864  < 2e-16 ***
## wdayMon          -906.85     331.82  -2.733 0.006299 ** 
## wdaySat         -2670.15     331.83  -8.047 1.04e-15 ***
## wdaySun         -4786.87     331.86 -14.424  < 2e-16 ***
## wdayThu          1308.28     331.97   3.941 8.22e-05 ***
## wdayTue          1324.82     331.99   3.991 6.68e-05 ***
## wdayWed          1487.40     331.97   4.481 7.60e-06 ***
## humidity         -196.10      11.04 -17.767  < 2e-16 ***
## factor(month)2    666.31     449.57   1.482 0.138373    
## factor(month)3    963.12     446.34   2.158 0.030988 *  
## factor(month)4   1140.89     476.84   2.393 0.016763 *  
## factor(month)5   3749.68     504.56   7.432 1.25e-13 ***
## factor(month)6   5156.49     543.10   9.495  < 2e-16 ***
## factor(month)7   4962.75     569.75   8.710  < 2e-16 ***
## factor(month)8   2271.48     557.27   4.076 4.65e-05 ***
## factor(month)9   3448.54     524.53   6.574 5.36e-11 ***
## factor(month)10  3938.39     478.80   8.226 2.43e-16 ***
## factor(month)11  3095.74     444.06   6.972 3.52e-12 ***
## factor(month)12 -1649.08     431.65  -3.820 0.000135 ***
## precip           -279.62      19.33 -14.469  < 2e-16 ***
## 
## Residual standard error: 6438 on 5247 degrees of freedom
## Multiple R-squared:  0.5421, Adjusted R-squared:  0.5404 
## F-statistic: 310.7 on 20 and 5247 DF,  p-value: < 2.2e-16
# plot actual data vs predicted data
broom::augment(model5) |> 
  mutate(day=row_number()) |> 
  ggplot()+
  aes(x=day, y = bikes_hired)+
  
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

Further variables/questions to explore on your own

Our dataset has many more variables, so here are some ideas on how you can extend your analysis

  • Are other weather variables useful in explaining bikes_hired?
  • We also have data on days of the week, month of the year, etc. Could those be helpful?
  • What’s the best model you can come up with?
  • Is this a regression model to predict or explain? If we use it to predict, what’s the Residual SE?

Ploting your best model predicitons vs reality

Let us say that your best model is model3. How do the predictions of your model compare to the actual data?

# plot actual data vs predicted data

my_best_model <- broom::augment(model3) %>% 
  mutate(day=row_number())



# Plot fitted line and residuals
ggplot(my_best_model, aes(x=day, y = bikes_hired)) +
  geom_point(alpha = 0.2) +
  geom_point(aes(y = .fitted), 
             shape = 1, 
             colour = "red", alpha = 0.2)+
  theme_bw()

Diagnostics, collinearity, summary tables

As you keep building your models, it makes sense to:

  1. Check the residuals, using autoplot. You will always have some deviation from normality, especially for very high values of n
  2. As you start building models with more explanatory variables, make sure you use car::vif(model_x) to calculate the Variance Inflation Factor (VIF) for your predictors and determine whether you have colinear variables. A general guideline is that a VIF larger than 10 is large, and your model may suffer from collinearity. Remove the variable in question and run your model again without it.
# Residual Analysis
model2 %>% 
  autoplot(which = 1:3) +
  theme_bw()
## Warning: `fortify(<lm>)` was deprecated in ggplot2 3.6.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Check VIF
vif(model2)
## feelslikemax      weekend 
##     1.000011     1.000011
vif(model3)
##          GVIF Df GVIF^(1/(2*Df))
## temp 1.000059  1        1.000029
## wday 1.000059  6        1.000005
vif(model4)
##                   GVIF Df GVIF^(1/(2*Df))
## temp          3.470724  1        1.862988
## wday          1.001354  6        1.000113
## humidity      1.546931  1        1.243757
## factor(month) 3.958549 11        1.064537

Comparison of different models

Create a summary table, using huxtable (https://am01-sep23.netlify.app/example/modelling_side_by_side_tables/) that shows which models you worked on, which predictors are significant, the adjusted \(R^2\), and the Residual Standard Error. If you want to add more models, just make sure you do not forget the comma , after the last model, as shown below

# produce summary table comparing models using huxtable::huxreg()
huxreg(model0, model1, model2, model3, model4,
       statistics = c('#observations' = 'nobs', 
                      'R squared' = 'r.squared', 
                      'Adj. R Squared' = 'adj.r.squared', 
                      'Residual SE' = 'sigma'), 
       bold_signif = 0.05
) %>% 
  set_caption('Comparison of models')
Comparison of models
(1)(2)(3)(4)(5)
(Intercept)26327.597 ***14556.638 ***15800.662 ***15310.174 ***38443.272 ***
(130.846)   (220.781)   (221.757)   (335.796)   (1064.536)   
feelslikemax        835.436 ***834.529 ***                
        (13.938)   (13.438)                   
weekendTRUE                -4306.903 ***                
                (215.285)                   
temp                        1039.450 ***537.960 ***
                        (18.313)   (31.267)   
wdayMon                        -880.364 *  -924.224 ** 
                        (375.400)   (344.065)   
wdaySat                        -2731.599 ***-2661.858 ***
                        (375.400)   (344.074)   
wdaySun                        -4732.832 ***-4821.699 ***
                        (375.400)   (344.076)   
wdayThu                        1236.424 ***1308.183 ***
                        (375.525)   (344.182)   
wdayTue                        1205.904 ** 1311.241 ***
                        (375.524)   (344.204)   
wdayWed                        1315.157 ***1439.998 ***
                        (375.526)   (344.208)   
humidity                                -256.423 ***
                                (11.004)   
factor(month)2                                432.823    
                                (465.442)   
factor(month)3                                1481.660 ** 
                                (459.975)   
factor(month)4                                2045.016 ***
                                (485.364)   
factor(month)5                                4500.815 ***
                                (517.000)   
factor(month)6                                5800.743 ***
                                (566.182)   
factor(month)7                                5428.741 ***
                                (603.254)   
factor(month)8                                2778.824 ***
                                (588.583)   
factor(month)9                                4373.232 ***
                                (543.674)   
factor(month)10                                4715.447 ***
                                (494.406)   
factor(month)11                                3504.314 ***
                                (459.798)   
factor(month)12                                -1582.596 ***
                                (447.961)   
#observations5268        5268        5268        5268        5268        
R squared0.000    0.406    0.448    0.412    0.508    
Adj. R Squared0.000    0.405    0.447    0.412    0.506    
Residual SE9496.916    7322.808    7060.078    7284.098    6675.871    
*** p < 0.001; ** p < 0.01; * p < 0.05.