Preparing the model

Let us now prepare the model.

Firstly, we need to deal with the missing data. It is completely pointless to take into account the stations where the satisfaction survey responses aren’t present. So we’ll simply remove these entries. Not that I continue working with both dataframes, one with factors and another one with actual values, in order to see on which one the model performs better.

y <- is.na(df$Global)
y <- !y
df <- df[y,]
dat <- dat[y,]

Next, I examine the rest of the missing values using mice and VIM libraries and complete the entries by generating imputed data using cartesian method for df and pmm for dat (the code for only one of these manipulations is shown):

library(mice)
md.pattern(df)
md.pattern(dat)

library(VIM)
mice_plot <- aggr(df, col=c('navyblue','yellow'),
                  numbers=TRUE, sortVars=TRUE,
                  labels=names(df), cex.axis=.7,
                  gap=3, ylab=c("Missing data","Pattern"))

imputed_Data1 <- mice(df, m=5, maxit = 50, method = 'cart', seed = 500)

df <- complete(imputed_Data1,1)

Let’s check what we’ve got:

summary(df)

##       UIC          Works             P0              P1              P2
##  672006 :  2   Min.   :0.000   Min.   :5.700   Min.   :6.640   Min.   :6.930
##  111849 :  1   1st Qu.:0.000   1st Qu.:8.140   1st Qu.:8.290   1st Qu.:8.134
##  113001 :  1   Median :0.000   Median :8.330   Median :8.545   Median :8.420
##  118000 :  1   Mean   :1.016   Mean   :8.313   Mean   :8.479   Mean   :8.375
##  141002 :  1   3rd Qu.:2.000   3rd Qu.:8.570   3rd Qu.:8.740   3rd Qu.:8.690
##  142109 :  1   Max.   :4.000   Max.   :9.680   Max.   :9.680   Max.   :9.330
##  (Other):119   NA's   :2
##        P3              P4              P5              P6
##  Min.   :6.450   Min.   :5.660   Min.   :5.910   Min.   :6.060
##  1st Qu.:7.830   1st Qu.:7.305   1st Qu.:7.305   1st Qu.:7.475
##  Median :8.168   Median :7.630   Median :7.650   Median :7.812
##  Mean   :8.146   Mean   :7.588   Mean   :7.597   Mean   :7.782
##  3rd Qu.:8.475   3rd Qu.:7.945   3rd Qu.:7.947   3rd Qu.:8.158
##  Max.   :9.600   Max.   :8.680   Max.   :8.600   Max.   :8.840
##
##        P7            Piano           PowStat          BabyFoot
##  Min.   :5.270   Min.   :0.0000   Min.   :0.0000   Min.   :0.000
##  1st Qu.:7.058   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.000
##  Median :7.450   Median :1.0000   Median :0.0000   Median :0.000
##  Mean   :7.396   Mean   :0.7977   Mean   :0.8876   Mean   :0.191
##  3rd Qu.:7.765   3rd Qu.:1.0000   3rd Qu.:2.0000   3rd Qu.:0.000
##  Max.   :8.700   Max.   :3.0000   Max.   :7.0000   Max.   :1.000
##                  NA's   :37       NA's   :37       NA's   :37
##      Distr            People              Clean            Global
##  Min.   :0.0000   Min.   :   538454   Min.   : 82.37   Min.   :6.402
##  1st Qu.:0.0000   1st Qu.:  1899648   1st Qu.: 93.01   1st Qu.:7.834
##  Median :1.0000   Median :  3721942   Median : 95.05   Median :8.086
##  Mean   :0.6517   Mean   : 17384725   Mean   : 94.51   Mean   :8.037
##  3rd Qu.:1.0000   3rd Qu.: 10189076   3rd Qu.: 96.78   3rd Qu.:8.264
##  Max.   :4.0000   Max.   :527278558   Max.   :100.00   Max.   :8.898
##  NA's   :37                           NA's   :4

Note that there are no more ‘NA’ values in df dataframe, but there are still in dat dataframe, as we encoded one of the factor levels to ‘NA’.

Note also that the data in df is not scaled, because the advantage of the linear regression is that it works with unscaled variables, unlike many other methods.

Before proceding to splitting the data into test and training subsets, let’s make sure that it doesn’t contain outliers:

We can’t rely on the software highlighting the data as outliers, but have to use common sense. In our case, for instance, exceptionally high or low satisfaction ratings are most likely real, and the high number of visitors in the train station is, as well. I personally don’t see outliers in this dataset.

I then scale the continuous variables inf df and split the data using the train-and-test regimen of caret library, with 80% of the training data. I do this for both df and dat and check if there’re any problems with partitions. As this is a straghtforward commonly repeated manipulation, I leave the code out.

Model selection and training

I explicitly specify all variables in regression model:

model1 <- {Global ~ Works + Piano + PowStat + BabyFoot + Distr + People + Clean}
# fit linear regression model with scaled and factorized variables
dat_fit <- lm(Global ~ Works + Piano + PowStat + BabyFoot + Distr + People + Clean, data = train_dat)

# fit linear regression model with scaled and non-factor variables
df_fit <- lm(Global ~ Works + Piano + PowStat + BabyFoot + Distr + People + Clean, data = train)

It is fairly obvious that the model with continuous variables performs better in terms of certainty. The residuals are lower and the p-value is smaller. We therefore can choose to continue working with df, leaving dat aside.

The significant variables, at the first glance, are People and Distr. We’ve already seen that all these are correlated, so a number of visitors alone might explain the satisfaction just as well.

Let’s check.

df_simple <- lm(Global ~ People, data = train)

This simple model is performing great. Residual standard error is barely higher than that of the more complex model. Isn’t it actually surprising?

Let’s compare Akaike’s Information Criterion and Bayesian information criterion for both models:

## For df_fit:    AIC =  264.9658 , BIC =  288.5905
## For df_simple: AIC =  258.0035 , BIC =  265.8784

The smaller the values, the better the fit is. However, BIC is punishing the complexity of the model more severely.

Based on Mean Squared Error I would still lean towards the slightly more complex model, even though the difference is not huge:

## For df_fit:    RMSE =  0.8119456
## For df_simple: RMSE =  0.8322459

#let's see if we can build a better model with interactions

full_model <- {Global ~ (Works + Piano + PowStat + BabyFoot + Distr + People + Clean)^2}
chosen_model_fit <- step(lm(full_model, data = train), direction ="backward")
chosen_model_formula <- formula(chosen_model_fit)

Based on the anova comparison of two models, there is significant difference between two.

Checking assumptions of linear regression

Let us first have a look at residuals for three models:

par(mfrow=c(1,3),mar=c(4, 4, 2, 2))
hist(df_fit$residuals,
     main = "df_fit model", xlab = "Residuals")
hist(df_simple$residuals,
     main = "df_simple model", xlab = "Residuals")
hist(chosen_model_fit$residuals,
     main = "chosen_model_fit model", xlab = "Residuals")

All three are normally distributed, as they should be if the fit is good. Let’s have a look at Q-Q plots for all three:

We should also check the model for multicollinearity across explanatory variables and ensure that all values are low:

print(vif(df_fit))

##    Works    Piano  PowStat BabyFoot    Distr   People    Clean
## 1.293768 1.143916 1.230986 1.031874 1.188111 1.478836 1.080747

This obviously won’t make sense for chosen_model_fit, as interactions between variables will have much higher value.