Question 11.6: Consider the pizza delivery data described in Appendix A.4. ......

The output shows that lower temperature, higher bills, and more ordered pizzas increase the delivery times. The branch in the East is the fastest, and so is the driver Domenico. While there are differences with respect to day, discount customers, and the operator, they are not significant at the 5 % level.

(b) The confidence intervals are calculated as: \hat{\beta }_{i} ± t_{n−p−1;1−α/2} ·\hat{\sigma }_{\hat{\beta }} . We know from the model output from (a) that there are 1248 degrees of freedom (1266 observations − 18 estimated coefficients). The respective quantile from the t-distribution is obtained with the qt() function. The coefficients are accessed via the coefficients command (alternatively: mp$coefficients); the variances of the coefficients are either accessed via the diagonal elements of the covariance matrix (diag(vcov(mp))) or the model summary (summary(mp)[[4]][,2]—both of which are laborious. The summary of coefficients, lower confidence limit (lcl), and upper confidence limit (ucl) may be summarized in a matrix, e.g. via merging the individual columns with the cbind command.

(c) The variance is estimated as the residual sum of squares divided by the degrees of freedom, see also (11.27). Applying the residuals command to the model and using other basic operations yields an estimated variance of 28.86936.

\hat{\sigma}^{2}=\frac{\left(y-X\hat{\beta } \right)^{\prime } \left(y-X\hat{\beta } \right)}{n-\left(p+1\right) } =\frac{\hat{e}^{\prime } \hat{e}}{n-\left(p+1\right)} =\frac{1}{n-\left(p+1\right)}\sum\limits_{i=1}^{n}{\hat{e}^{2}_{i}}.  (11.27)

Taking the square root of the result yields \sqrt{ 28.86936} = 5.37 which is also reported in the model output from (a) under “Residual standard error”.

(d) The sum of squares error is defined as \Sigma ^{n}_{i=1} \left(y_{i}-\hat{y}_{i} \right)^{2} . The total sum of squares is \Sigma ^{n}_{i=1} \left(y_{i}-\bar{y} \right)^{2}. This can be easily calculated in R. The goodness of fit is then obtained as R² = 1 − {S Q_{Error}}/{S Q_{Total}} = 0.3178. Dividing SQ_{Error} by n − p − 1 = 1266 − 17 − 1 = 1248 and SQ_{Total} by n − 1 = 1265 yields R^{2}_{ adj} = 0.3085.
This corresponds to the model output from (a).

(e) Applying stepAIC to the fitted model (with option “back” for backward selection) executes the model selection by means of AIC.

The output shows that the full model has an AIC of 4275.15. The smallest AIC is achieved by removing the operator variable from the model.

The reduced model has an AIC of 4273.37. Removing the discount customer variable from the model yields an improved AIC (4272.0 < 4273.37).

The model selection procedure stops here as removing any variable would only increase the AIC, not decrease it.

The final model, based on backward selection with AIC, includes day, driver, branch, number of pizzas ordered, temperature, and bill.

(f) Fitting the linear model with the variables obtained from (e) and obtaining the summary of it yields an R^{2}_{ adj} of 0.3092.

This is only marginally higher than the goodness of fit from the full model (0.3092 > 0.3085). While the selected model is better than the model with all variables, both, with respect to AIC and R^{2}_{ adj}, the results are very close and remind us of the possible instability of applying automated model selection procedures.

(g) Both the normality assumption and heteroscedasticity can be checked by applying plot() to the model. From the many graphs provided we concentrate on the second and third of them:

Figure B.23a shows that the residuals are approximately normally distributed because the dots lie approximately on the bisecting line. There are some smaller deviations at the tails but they are not severe. The plot of the fitted values versus the square root of the absolute values of the residuals shows no pattern; it is a random plot (Fig. B.23b). There seems to be no heteroscedasticity.

(h) Not all variables identified in (e) represent necessarily a “cause” for delayed or improved delivery time. It makes sense to speculate that because many pizzas are being delivered (and need to be made!) the delivery time increases. There might also be reasons why a certain driver is improving the delivery time: maybe he does not care about red lights. This could be investigated further given the results of the model above. However, high temperature does not cause the delivery time to be shorter; likely it is the other way around: the temperature is hotter because the delivery time is shorter. However, all of these considerations remain speculation. A regression model only exhibits associations. If there is a significant association, we know that given an accepted error (e.g. 5 %), values of x are higher when values of y are higher. This is useful but it does not say whether x caused y or vice versa.

(i) To check whether it is worth to add a polynomial, we simply add the squared temperature to the model. To make R understand that we apply a transformation we need to use I().

It can be seen that the null hypothesis H_{0} : β_{temp^{2}} = 0 is rejected. This indicates that it is worthwhile to assume a quadratic relationship between temperature and delivery time.

(j) The prediction can be obtained by the predict command as follows:

The prediction is 36.5 min and therefore 0.8 min higher than the real delivery time.