Introduction & Description of the Data

This data set is called airq402.dat. It describes the airfares and passengers for U.S. Domestic Routes for 4th Quarter of 2002. That is, in includes information concerning airfare pricing, passenger volume, market share, and distances flown. The data was sourced by the United States Department of Transportation, and presumably comes from a random sample of size n=1,000 of the total number of flights that occured in the fourth quarter of 2002, from the first of October to the 31st of December. The variables and respective descriptions are as follows:

As for practical and analytical questions, I am curious to know which cities had the highest volume of people travelling to and from them, especially considering these data were collected during the fourth fiscal quarter of the year, during the holiday season. I would also like to know why some variables are repeated with different observation values.

The data set does have enough information for me to answer my main statistical inquiry of the relationship between distance and average fare. In the dataset there are elevel variables, and 1,000 total observations. It is stipulated in the guidelines that there are to be at least 15 observations for each variable.


Practical Question

Our objective and primary question is to determine the relationship between airplane ticket price and relevant predictor variables in the dataset.



Exploratory Data Analysis

As shown in last week’s assignment, a pairwise scatterplot of the data show that multiple variables seem highly correlated, and may in fact be repetitive given the nature of the names they were given. Examples include marketShare and marketShare2, which in the original data, we both simply named “market share”.

airfareNumVars <- c(3, 4, 5, 7, 8, 10, 11)
airfareNumeric <- airfare[airfareNumVars]

pairs(airfareNumeric, main ="Pair-wise Association: Scatter Plot")

Because of the extreme correlation between variables like marketShare and marketShare2, I will remove the duplicate variables from the dataset, and create a subset called airfaireMult with only the variables that I am interested in.

airfareMult <- airfare[c(3, 4, 5, 7, 11)]
kable(head(airfareMult))
averageFare distance averageWeeklyPassengers marketShare price
114.47 528 424.56 70.19 111.03
122.47 860 276.84 75.10 118.94
214.42 852 215.76 78.89 167.12
69.40 288 606.84 96.97 68.86
158.13 723 313.04 39.79 145.42
135.17 1204 199.02 40.68 127.69

Multiple Linear Regression Model & Diagnostics

full.model = lm(price ~ ., data = airfareMult)
kable(summary(full.model)$coef, caption ="Statistics of Regression Coefficients")
Statistics of Regression Coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.4524702 4.1537034 9.979641 0.0000000
averageFare 0.6873902 0.0163178 42.125196 0.0000000
distance 0.0042585 0.0016127 2.640632 0.0084048
averageWeeklyPassengers -0.0025610 0.0009608 -2.665499 0.0078119
marketShare -0.2218813 0.0448108 -4.951510 0.0000009

Next I will conduct residual diagnostic tests:

par(mfrow=c(2,2))
plot(full.model)

We can see from the residual plots that there are some violations:

We first perform Box-Cox transformation to correct the non-constant variance and correct the non-normality of the QQ plot.

library(MASS)
par(pty = "s", mfrow = c(2, 2), oma=c(.1,.1,.1,.1), mar=c(4, 0, 2, 0))
##
boxcox(price ~ averageWeeklyPassengers + averageFare + distance +  log(marketShare)  
       , data = airfareMult, lambda = seq(0, 1, length = 10), 
       xlab=expression(paste(lambda, ": log dist2MRT")))
##
boxcox(price ~ averageWeeklyPassengers + averageFare + distance +  marketShare
       , data = airfareMult, lambda = seq(-0.5, 1, length = 10), 
       xlab=expression(paste(lambda, ": dist2MRT")))
##
boxcox(price ~ log(1+averageWeeklyPassengers) + averageFare + distance +  marketShare 
      , data = airfareMult, lambda = seq(-0.5, 1, length = 10), 
       xlab=expression(paste(lambda, ": log-age")))
##
boxcox(price ~ log(1+averageWeeklyPassengers) + averageFare + distance +  log(marketShare)
       , data = airfareMult, lambda = seq(-0.5, 1, length = 10), 
      xlab=expression(paste(lambda, ": log-age, log.dist2MRT")))

sqrt.price.log.share = lm((price)^0.5 ~ averageWeeklyPassengers + averageFare + distance +  log(marketShare)  , data = airfareMult)
kable(summary(sqrt.price.log.share)$coef, caption = "log-transformed model")
log-transformed model
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.0478876 0.4356211 20.770087 0.0000000
averageWeeklyPassengers -0.0001440 0.0000381 -3.777641 0.0001677
averageFare 0.0281128 0.0006526 43.078822 0.0000000
distance 0.0002236 0.0000655 3.413054 0.0006683
log(marketShare) -0.4881688 0.0993120 -4.915508 0.0000010 
par(mfrow = c(2,2))
plot(sqrt.price.log.share)

log.price = lm(log(price) ~ averageWeeklyPassengers + averageFare + distance +  marketShare, data = airfareMult)
kable(summary(log.price)$coef, caption = "log-transformed model")
log-transformed model
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2582865 0.0281403 151.323388 0.00e+00
averageWeeklyPassengers -0.0000323 0.0000065 -4.968216 8.00e-07
averageFare 0.0047022 0.0001105 42.534880 0.00e+00
distance 0.0000362 0.0000109 3.316360 9.45e-04
marketShare -0.0022034 0.0003036 -7.258122 0.00e+00 
par(mfrow = c(2,2))
plot(log.price)

#define plotting area
par(pty = "s", mfrow = c(1, 3))
#Q-Q plot for original model
qqnorm(full.model$residuals, main = "Full-Model")
qqline(full.model$residuals)
#Q-Q plot for Box-Cox transformed model
qqnorm(log.price$residuals, main = "Log-Price")
qqline(log.price$residuals)
#display both Q-Q plots
qqnorm(sqrt.price.log.share$residuals, main = "sqrt price log dist")
qqline(sqrt.price.log.share$residuals)

Goodness of Fit Measures

select=function(m){ # m is an object: model
 e = m$resid                           # residuals
 n0 = length(e)                        # sample size
 SSE=(m$df)*(summary(m)$sigma)^2       # sum of squared error
 R.sq=summary(m)$r.squared             # Coefficient of determination: R square!
 R.adj=summary(m)$adj.r                # Adjusted R square
 MSE=(summary(m)$sigma)^2              # square error
 Cp=(SSE/MSE)-(n0-2*(n0-m$df))         # Mellow's p
 AIC=n0*log(SSE)-n0*log(n0)+2*(n0-m$df)          # Akaike information criterion
 SBC=n0*log(SSE)-n0*log(n0)+(log(n0))*(n0-m$df)  # Schwarz Bayesian Information criterion
 X=model.matrix(m)                     # design matrix of the model
 H=X%*%solve(t(X)%*%X)%*%t(X)          # hat matrix
 d=e/(1-diag(H))                       
 PRESS=t(d)%*%d   # predicted residual error sum of squares (PRESS)- a cross-validation measure
 tbl = as.data.frame(cbind(SSE=SSE, R.sq=R.sq, R.adj = R.adj, Cp = Cp, AIC = AIC, SBC = SBC, PRD = PRESS))
 names(tbl)=c("SSE", "R.sq", "R.adj", "Cp", "AIC", "SBC", "PRESS")
 tbl
 }

output.sum = rbind(select(full.model), select(sqrt.price.log.share), select(log.price))
row.names(output.sum) = c("full.model", "sqrt.price.log.dist", "log.price")
kable(output.sum, caption = "Goodness-of-fit Measures of Candidate Models")
Goodness-of-fit Measures of Candidate Models
SSE R.sq R.adj Cp AIC SBC PRESS
full.model 521669.29147 0.7658960 0.7649548 5 6267.034 6291.5726 529107.91138
sqrt.price.log.dist 822.47761 0.7805446 0.7796624 5 -185.434 -160.8952 833.20936
log.price 23.94319 0.7815141 0.7806357 5 -3722.071 -3697.5327 24.25261

We see that the Goodness of Fit measures of the third model are superior, thus we choose to use the third model as our final model.

kable(summary(log.price)$coef, caption = "Inferential Statistics of Final Model")
Inferential Statistics of Final Model
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2582865 0.0281403 151.323388 0.00e+00
averageWeeklyPassengers -0.0000323 0.0000065 -4.968216 8.00e-07
averageFare 0.0047022 0.0001105 42.534880 0.00e+00
distance 0.0000362 0.0000109 3.316360 9.45e-04
marketShare -0.0022034 0.0003036 -7.258122 0.00e+00

Summary of the Model

\[ \log(price) =4.2583 - 0.0000323\times averageWeeklyPassengers +0.0047022\times averageFare + 0.0000362\times distance - 0.0022034\times marketShare \]

averageWeeklyPassengers and marketShare are both negatively associated with price, while distance and averageFare are positively associated with price.

Discussion

We used regression techniques such as box-cox transformation for the response, and other transformations of the independent variables to identify a suitable final model for this dataset. Due to the nature of some of the variables from the initial data set, I performed variable selection to remove what appeared to be repetitive variables from the data set.

I used common global goodness-of-fit tests for model selection.

I had difficulty addressing and interpreting the true meaning of the final regression coefficients because of the nature of the transformations used.

The violation of the assumption of normality persists.