This is Quantitative Decision Making for Business by using R studio.

Chapter 19: Building Multiple Regression Models

Diamond Prices

The price of jewelry diamond depends on the 4 C’s: Carat weight, Cut, Clarity, and Color. We have collected a sample of diamond prices from a website, together with other related variables. We want to investigate how the price of a diamond is determined by its attributes.

Let’s first read the dataset into R diamond <- read.table('Diamonds_DGK_FE18.txt', sep = '\t', header = TRUE)

and look at its structure str(diamond)

## 'data.frame': 749 obs. of 8 variables: ## $ Price : int 1000 1000 1000 1001 1001 1018 1019 1019 1019 1019 ... ## $ LogPrice : num 3 3 3 3 3 ... ## $ Color : Factor w/ 3 levels "D","G","K": 3 2 1 1 1 2 2 1 1 1 ... ## $ Carat.Weight: num 0.66 0.4 0.36 0.53 0.44 0.41 0.42 0.33 0.32 0.32 ... ## $ Clarity : Factor w/ 7 levels "IF","SI1","SI2",..: 2 4 5 3 2 4 7 7 6 6 ... ## $ Depth : num 62.8 62 61.3 59.4 60.2 61.4 62.3 60.7 62.7 61.5 ... ## $ Table : int 57 59 57 59 60 55 56 57 55 57 ... ## $ Cut : Factor w/ 4 levels "Excellent","Good",..: 1 1 1 4 1 3 1 3 4 1 ...

As we can see, there are 749 observations for each of 8 variables. Some variables are integers (int) (e.g., P rice), some are numerical (num) (e.g., Carat.W eight), and some are categorical (Factor) (e.g., Color).

Let’s first consider modeling the relationship between P rice and Carat W eight. We have both P rice and log10 P rice in the dataset. From the scatterplots below we will use log10 P rice as response variable because the Linearity Assumption is met there. par(mfrow = c(1,2)) plot(diamond$Carat.Weight, diamond$Price, xlab = 'Carat Weight', ylab = 'Price') plot(diamond$Carat.Weight, diamond$LogPrice, xlab = 'Carat Weight', ylab = 'Log10Price')

0.4 0.6 0.8 1.0 1.2 1.4

2 0

0 0

6 0

0 0

1 0

0 0

0

Carat Weight

P ri

ce

0.4 0.6 0.8 1.0 1.2 1.4

3 .0

3 .4

3 .8

Carat Weight

L o

g 1

0 P

ri ce

We now want to add Color to the model as a predictor. To see if it is appropriate to use indicators for Color, we plot log10 P rice against Carat W eight with points labelled by their color levels:

1

library(ggplot2) ggplot(diamond, aes(Carat.Weight, LogPrice, color = Color)) + geom_point()

3.00

3.25

3.50

3.75

4.00

0.25 0.50 0.75 1.00 1.25 Carat.Weight

L o

g P

ri ce

Color

D

G

K

As we can see, the diamonds with different color levels have similar relationships between their prices and carat weights (the slopes are similar). We therefore can use 2 indicators to represent this 3-level categorical predictor. Note that we do not need to make the indicators by ourselves, because R has already recognized Color as a factor variable.

Let’s regress LogP rice on Carat W eight and Color in R as follows: imod1 <- lm(LogPrice ~ Carat.Weight + Color, data = diamond)

and the regression output is given below: summary(imod1)

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + Color, data = diamond) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.27859 -0.07514 -0.00868 0.07294 0.38937 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.733840 0.012396 220.550 < 2e-16 *** ## Carat.Weight 1.028648 0.015887 64.747 < 2e-16 *** ## ColorG -0.075225 0.009224 -8.155 1.48e-15 *** ## ColorK -0.294055 0.014238 -20.653 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1088 on 745 degrees of freedom ## Multiple R-squared: 0.8566, Adjusted R-squared: 0.8561 ## F-statistic: 1484 on 3 and 745 DF, p-value: < 2.2e-16

The model looks good in general. It has R2 = 85.66%. Both F - and t-tests are significant under 5%

2

significance level. So, the model is statistically useful and each predictor contributes significantly given other predictors already in the model.

Note that R takes ColorD as base level and therefore yields slope estimates for ColorG and ColorK. How do we interpret the estimate (-0.294) for ColorK? Given other predictors as constant, the mean log10P rice of diamonds with K-level color is less than that of those with D-level color by 0.294 (unit). But, what does this 0.294 mean on the P rice scale? We can see it from the following equation:

log10 P riceD − log10 P riceK = 0.294

log10 P riceD P riceK

= 0.294

P riceD P riceK

= 100.294 ≈ 1.97

We therefore interpret that the diamonds with D-level color on average cost almost twice as much as those with K-level color, after accounting for other predictors in the model.

It is easy to add interaction terms to the regression model in R. For example, we may add the interaction between Carat W eight and Color to the previous model as follows: summary(lm(LogPrice ~ Carat.Weight + Color + Carat.Weight*Color, data = diamond))

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + Color + Carat.Weight * ## Color, data = diamond) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.27307 -0.07750 -0.00822 0.07186 0.38235 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.70617 0.01929 140.269 < 2e-16 *** ## Carat.Weight 1.07044 0.02742 39.045 < 2e-16 *** ## ColorG -0.03810 0.02694 -1.414 0.15771 ## ColorK -0.16466 0.06302 -2.613 0.00916 ** ## Carat.Weight:ColorG -0.05267 0.03432 -1.535 0.12531 ## Carat.Weight:ColorK -0.14076 0.06356 -2.215 0.02709 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1085 on 743 degrees of freedom ## Multiple R-squared: 0.8577, Adjusted R-squared: 0.8567 ## F-statistic: 895.7 on 5 and 743 DF, p-value: < 2.2e-16

It seems that the interaction terms do not improve the model by any means: the R2 is almost the same, the F -test gets less significant, the t-tests for Color become much less significant, and the t-tests for the interaction terms are not significant. The model does not benefit from its more complexity, so it is not wise to add the interaction terms to the model in this scenario.

Now we want to add more predictors (not interaction terms) to the model in order to improve its R2. We have four variables left in the pool: T able, Depth, Clarity and Cut. Let’s first see whether or not the two continuous variables are linearly related to log10 P rice: par(mfrow = c(1,2)) plot(factor(diamond$Table), diamond$LogPrice, xlab = 'Table', ylab = 'LogPrice') plot(diamond$Depth, diamond$LogPrice, xlab = 'Depth', ylab = 'LogPrice')

3

53 56 59 62 65

3 .0

3 .4

3 .8

Table

L o

g P

ri ce

58 60 62 64

3 .0

3 .4

3 .8

Depth

L o

g P

ri ce

As we can see, both plots show no linear pattern at all, and we therefore do not need those variables in the model. How about the two categorical variables: Clarity and Cut? The new multiple regression model is fitted as follows: imod2 <- lm(LogPrice ~ Carat.Weight + Color + Clarity + Cut, data = diamond)

and the regression output is given below: summary(imod2)

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + Color + Clarity + Cut, ## data = diamond) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.215186 -0.042365 0.005408 0.044710 0.189934 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.780792 0.011558 240.593 < 2e-16 *** ## Carat.Weight 1.200301 0.011365 105.611 < 2e-16 *** ## ColorG -0.088864 0.005866 -15.148 < 2e-16 *** ## ColorK -0.342888 0.009107 -37.649 < 2e-16 *** ## ClaritySI1 -0.237850 0.011736 -20.267 < 2e-16 *** ## ClaritySI2 -0.312076 0.012454 -25.058 < 2e-16 *** ## ClarityVS1 -0.130431 0.012118 -10.764 < 2e-16 *** ## ClarityVS2 -0.178439 0.011684 -15.273 < 2e-16 *** ## ClarityVVS1 -0.064518 0.012445 -5.184 2.81e-07 *** ## ClarityVVS2 -0.076747 0.011878 -6.461 1.89e-10 *** ## CutGood -0.028884 0.011133 -2.594 0.00966 ** ## CutIdeal 0.012127 0.010022 1.210 0.22662 ## CutVery Good -0.015374 0.005503 -2.793 0.00535 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.06806 on 736 degrees of freedom ## Multiple R-squared: 0.9446, Adjusted R-squared: 0.9437 ## F-statistic: 1045 on 12 and 736 DF, p-value: < 2.2e-16

As we can see, the R2 is improved a lot, the overall F -test is highly significant, and all the t-tests are

4

significant except the one for CutIdeal (p-value = 0.23). The insignificant coefficient of CutIdeal indicates that there is no significant difference in the prices between the diamonds with ideal cut and those with excellent cut (Note: the excellent cut is the baseline). We may consider combining those two cut levels.

Don’t forget to check the regression assumptions, although the model “looks wonderful“. Let’s plot the residuals against the fitted values: plot(imod2$fitted.values, imod2$residuals, xlab = 'Fitted values', ylab = 'Residuals') abline(0, 0)

3.0 3.2 3.4 3.6 3.8 4.0

− 0

.2 0

.0 0

.2

Fitted values

R e

si d

u a

ls

and we see an obvious bend that needs to be corrected. Add a squared term of Carat W eight to the model: imod3.tmp <- lm(LogPrice ~ Carat.Weight + I(Carat.Weight^2) +

Color + Clarity + Cut, data = diamond) summary(imod3.tmp)

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + I(Carat.Weight^2) + Color + ## Clarity + Cut, data = diamond) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.152181 -0.031578 -0.002784 0.030240 0.153036 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.388025 0.016190 147.499 < 2e-16 *** ## Carat.Weight 2.365036 0.042400 55.780 < 2e-16 *** ## I(Carat.Weight^2) -0.699846 0.025028 -27.962 < 2e-16 *** ## ColorG -0.095142 0.004092 -23.248 < 2e-16 *** ## ColorK -0.347479 0.006346 -54.755 < 2e-16 *** ## ClaritySI1 -0.286641 0.008359 -34.292 < 2e-16 *** ## ClaritySI2 -0.353371 0.008800 -40.155 < 2e-16 *** ## ClarityVS1 -0.161288 0.008513 -18.947 < 2e-16 *** ## ClarityVS2 -0.215559 0.008246 -26.141 < 2e-16 *** ## ClarityVVS1 -0.077703 0.008682 -8.950 < 2e-16 *** ## ClarityVVS2 -0.103078 0.008327 -12.378 < 2e-16 *** ## CutGood -0.038896 0.007763 -5.010 6.81e-07 *** ## CutIdeal 0.016165 0.006982 2.315 0.020877 * ## CutVery Good -0.014988 0.003834 -3.910 0.000101 *** ## ---

5

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.04741 on 735 degrees of freedom ## Multiple R-squared: 0.9731, Adjusted R-squared: 0.9727 ## F-statistic: 2048 on 13 and 735 DF, p-value: < 2.2e-16

Is there anything wrong with the model? Let’s check if there exists any collinearity in the model: library(car)

## Loading required package: carData vif(imod3.tmp)

## GVIF Df GVIF^(1/(2*Df)) ## Carat.Weight 46.979493 1 6.854159 ## I(Carat.Weight^2) 44.512102 1 6.671739 ## Color 1.354881 2 1.078885 ## Clarity 1.539197 6 1.036592 ## Cut 1.106386 3 1.016993

The VIFs for Carat W eight and its squared term are significantly larger than those for other variables. It’s not surprising because the Carat W eight is highly correlated with its squared term: cor(diamond$Carat.Weight, (diamond$Carat.Weight)^2)

## [1] 0.9879121

To fix this issue, we first center Carat W eight and then square it: diamond$Carat.Weight2 <-(diamond$Carat.Weight - mean(diamond$Carat.Weight))^2

Let’s fit the model with the centered squared term and check its collinearity: imod3 <- lm(LogPrice ~ Carat.Weight + Carat.Weight2 +

Color + Clarity + Cut, data = diamond) vif(imod3)

## GVIF Df GVIF^(1/(2*Df)) ## Carat.Weight 1.660814 1 1.288726 ## Carat.Weight2 1.069692 1 1.034259 ## Color 1.354881 2 1.078885 ## Clarity 1.539197 6 1.036592 ## Cut 1.106386 3 1.016993

Now the predictors are pretty uncorrelated with each other. The regression output of imod3 is given below summary(imod3)

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + Carat.Weight2 + Color + ## Clarity + Cut, data = diamond) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.152181 -0.031578 -0.002784 0.030240 0.153036 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|)

6

## (Intercept) 2.851063 0.008434 338.040 < 2e-16 *** ## Carat.Weight 1.226520 0.007972 153.853 < 2e-16 *** ## Carat.Weight2 -0.699846 0.025028 -27.962 < 2e-16 *** ## ColorG -0.095142 0.004092 -23.248 < 2e-16 *** ## ColorK -0.347479 0.006346 -54.755 < 2e-16 *** ## ClaritySI1 -0.286641 0.008359 -34.292 < 2e-16 *** ## ClaritySI2 -0.353371 0.008800 -40.155 < 2e-16 *** ## ClarityVS1 -0.161288 0.008513 -18.947 < 2e-16 *** ## ClarityVS2 -0.215559 0.008246 -26.141 < 2e-16 *** ## ClarityVVS1 -0.077703 0.008682 -8.950 < 2e-16 *** ## ClarityVVS2 -0.103078 0.008327 -12.378 < 2e-16 *** ## CutGood -0.038896 0.007763 -5.010 6.81e-07 *** ## CutIdeal 0.016165 0.006982 2.315 0.020877 * ## CutVery Good -0.014988 0.003834 -3.910 0.000101 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.04741 on 735 degrees of freedom ## Multiple R-squared: 0.9731, Adjusted R-squared: 0.9727 ## F-statistic: 2048 on 13 and 735 DF, p-value: < 2.2e-16

It is interesting to see that the coefficient of Carat.Weight becomes much more significant in imod3 than in imod3.tmp, because its standard error is reduced a lot after the collinearity is gone. But, the other coefficients and statistics remain the same.

Again, don’t forget to check the regression assumptions, although the model seems to be fantastic. Let’s produce the residual plot and the histogram of residuals: par(mfrow = c(1, 2)) plot(imod3$fitted.values, imod3$residuals, xlab = 'Fitted values', ylab = 'Residuals') abline(0, 0) hist(imod3$residuals, main = "Histogram of residuals", xlab = 'Residuals')

3.0 3.2 3.4 3.6 3.8 4.0

− 0

.1 5

0 .0

0 0

.1 5

Fitted values

R e

si d

u a

ls

Histogram of residuals

Residuals

F re

q u

e n

cy

−0.15 −0.05 0.05 0.15

0 4

0 8

0 1

2 0

The residual plot looks pretty random and the histogram is unimodal and approximately symmetric. Therefore, the regression assumptions are satisfied in this model.

Let’s check if there is any outlier or high-leverage case in the dataset. We compute the standardized residuals and leverages, and produce their histograms as follows:

7

std.res <- rstandard(imod3) lev <- hatvalues(imod3)

par(mfrow = c(1,2)) hist(std.res, xlab = 'Standardized residuals', main = "") hist(lev, xlab = 'Leverages', main = "")

Standardized residuals

F re

q u

e n

cy

−3 −1 0 1 2 3

0 5

0 1

0 0

1 5

0

Leverages F

re q

u e

n cy

0.01 0.03 0.05

0 1

0 0

2 0

0 3

0 0

We can see that there are a few cases that stand out from the rest cases in the histograms. Let’s find out the outliers whose residuals are beyond (-3, 3): diamond[which(abs(std.res) > 3), c(2, 3, 4, 5, 8)]

## LogPrice Color Carat.Weight Clarity Cut ## 193 3.254790 K 0.80 VVS1 Excellent ## 219 3.301464 K 0.77 IF Excellent ## 299 3.400883 D 0.52 VS1 Very Good ## 368 3.500785 K 0.90 VS2 Good ## 375 3.500785 K 0.90 VS2 Good ## 399 3.556423 K 1.06 SI2 Excellent

Regarding the leverages a common rule is to flag any case whose leverage is more than 3 times larger than the mean leverage value, i.e., hii > 3( k+1n ). Let’s use this cutoff to identify the cases with large leverages: lev.cut <- 3*length(imod3$coefficients)/dim(diamond)[1] diamond[which(lev > lev.cut), c(2, 3, 4, 5, 8)]

## LogPrice Color Carat.Weight Clarity Cut ## 574 3.743823 K 1.40 VS1 Good ## 575 3.744058 K 1.16 IF Good ## 615 3.789863 K 1.34 VVS1 Good

It is interesting to note that there is no case that is both outlier and of high leverage. We should investigate those extreme cases further, trying to understand what makes them special.

Are the extreme cases we found out above influential to the regression model? Let’s compute Cook’s distance and DFFITS measure for each case, and plot them in boxplots as follows: cookd <- cooks.distance(imod3) dffit <- dffits(imod3)

par(mfrow = c(1,2)) boxplot(cookd, xlab = "Cook's distance")

8

boxplot(dffit, xlab = 'DFFITS') 0

.0 0

0 0

.0 1

5 0

.0 3

0

Cook's distance

− 0

.6 0

.0 0

.4

DFFITS

From the boxplot of Cook’s distance we see one case stand out a bit from others. Let’s find that case: diamond[which.max(cookd), c(2, 3, 4, 5, 8)]

## LogPrice Color Carat.Weight Clarity Cut ## 368 3.500785 K 0.9 VS2 Good

It is one of the outliers as identified by the standardized residuals. One of the common cutoffs for large DFFITS is 2

√ (k + 1)/(n − k − 1):

dffit.cut <- 2*sqrt(length(imod3$coefficients)/ (dim(diamond)[1] - length(imod3$coefficients)))

This cutoff, however, is a bit too conservative because it yields too many “influential” cases which(abs(dffit) > dffit.cut)

## 166 168 169 193 201 204 219 251 299 309 313 314 329 330 336 350 352 368 ## 166 168 169 193 201 204 219 251 299 309 313 314 329 330 336 350 352 368 ## 375 382 384 392 399 442 453 466 470 536 567 572 573 575 612 623 632 640 ## 375 382 384 392 399 442 453 466 470 536 567 572 573 575 612 623 632 640 ## 644 664 683 707 709 730 733 745 746 ## 644 664 683 707 709 730 733 745 746

We therefore just focus on the case that stands out in the boxplot: diamond[which.max(abs(dffit)), c(2, 3, 4, 5, 8)]

## LogPrice Color Carat.Weight Clarity Cut ## 368 3.500785 K 0.9 VS2 Good

This is the same influential case reflected by Cook’s distance. So, what is so special about that diamond? The residual for that diamond is positive, which means its price is much higher than the estimate given by the model. After a closer look, we find that the price of that diamond is about average, but the color is the worst, the cut is the worst, the clarity is mediocre, and the weight is just above average. A combination of a good price and bad qualities makes this diamond special. We don’t need to remove this case from the data because its influence is mild based on its relatively small Cook’s distance and DFFITS measure.

We are lazy and want to use stepwise regression to help us find a good model, given all the potential predictors. We need to first set up a “null” (simplest possible) model and a “full” (most complicated possible) model:

9

null <- lm(LogPrice ~ 1, data = diamond) full <- lm(LogPrice ~ Carat.Weight + Carat.Weight2 + Color + Clarity +

Cut + Depth + Table , data = diamond)

We should be aware that the stepwise regression implemented in R is a variant from the original stepwise method. At each step it uses AIC, instead of t-test, to select the predictors. It is thus a fully automatic procedure and we do not need to decide the entering or removing significance level as required in the original method. Let’s do it with R: step.reg <- step(null, scope=list(lower=null, upper=full), test = 'F')

## Start: AIC=-1870.33 ## LogPrice ~ 1 ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## + Carat.Weight 1 47.590 13.905 -2981.9 2556.6010 < 2.2e-16 *** ## + Color 2 3.067 58.428 -1904.7 19.5823 5.143e-09 *** ## + Clarity 6 2.918 58.577 -1894.8 6.1614 2.550e-06 *** ## + Carat.Weight2 1 0.585 60.910 -1875.5 7.1752 0.007554 ** ## + Cut 3 0.737 60.758 -1873.4 3.0116 0.029470 * ## <none> 61.495 -1870.3 ## + Depth 1 0.014 61.481 -1868.5 0.1692 0.680926 ## + Table 1 0.001 61.495 -1868.3 0.0069 0.933605 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Step: AIC=-2981.88 ## LogPrice ~ Carat.Weight ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## + Color 2 5.088 8.817 -3319.1 214.9870 < 2.2e-16 *** ## + Clarity 6 3.675 10.230 -3199.8 44.3721 < 2.2e-16 *** ## + Carat.Weight2 1 0.681 13.224 -3017.5 38.4219 9.426e-10 *** ## + Cut 3 0.652 13.253 -3011.8 12.1971 8.482e-08 *** ## + Table 1 0.090 13.815 -2984.7 4.8545 0.02788 * ## + Depth 1 0.079 13.826 -2984.1 4.2495 0.03961 * ## <none> 13.905 -2981.9 ## - Carat.Weight 1 47.590 61.495 -1870.3 2556.6010 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Step: AIC=-3319.13 ## LogPrice ~ Carat.Weight + Color ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## + Clarity 6 5.335 3.482 -4003.0 188.6997 < 2.2e-16 *** ## + Carat.Weight2 1 0.625 8.192 -3372.2 56.7693 1.425e-13 *** ## + Cut 3 0.531 8.286 -3359.6 15.8477 5.366e-10 *** ## + Table 1 0.096 8.721 -3325.3 8.1803 0.004353 ** ## <none> 8.817 -3319.1 ## + Depth 1 0.004 8.812 -3317.5 0.3570 0.550353 ## - Color 2 5.088 13.905 -2981.9 214.9870 < 2.2e-16 *** ## - Carat.Weight 1 49.611 58.428 -1904.7 4192.1171 < 2.2e-16 *** ## ---

10

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Step: AIC=-4002.98 ## LogPrice ~ Carat.Weight + Color + Clarity ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## + Carat.Weight2 1 1.728 1.754 -4514.6 727.0692 < 2.2e-16 *** ## + Cut 3 0.073 3.409 -4012.9 5.2619 0.001351 ** ## + Table 1 0.040 3.442 -4009.7 8.6591 0.003356 ** ## <none> 3.482 -4003.0 ## + Depth 1 0.001 3.481 -4001.1 0.1281 0.720543 ## - Clarity 6 5.335 8.817 -3319.1 188.6997 < 2.2e-16 *** ## - Color 2 6.748 10.230 -3199.8 716.0507 < 2.2e-16 *** ## - Carat.Weight 1 52.411 55.893 -1925.9 11123.3911 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Step: AIC=-4514.58 ## LogPrice ~ Carat.Weight + Color + Clarity + Carat.Weight2 ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## + Cut 3 0.102 1.652 -4553.6 15.1624 1.392e-09 *** ## + Table 1 0.020 1.734 -4521.2 8.5336 0.003593 ** ## <none> 1.754 -4514.6 ## + Depth 1 0.004 1.750 -4514.3 1.6960 0.193223 ## - Carat.Weight2 1 1.728 3.482 -4003.0 727.0692 < 2.2e-16 *** ## - Clarity 6 6.438 8.192 -3372.2 451.4441 < 2.2e-16 *** ## - Color 2 6.889 8.643 -3324.0 1449.3108 < 2.2e-16 *** ## - Carat.Weight 1 53.930 55.684 -1926.7 22691.4952 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Step: AIC=-4553.56 ## LogPrice ~ Carat.Weight + Color + Clarity + Carat.Weight2 + Cut ## ## Df Sum of Sq RSS AIC F value Pr(>F) ## <none> 1.652 -4553.6 ## + Depth 1 0.003 1.649 -4552.9 1.3189 0.2512 ## + Table 1 0.001 1.651 -4552.0 0.3894 0.5328 ## - Cut 3 0.102 1.754 -4514.6 15.1624 1.392e-09 *** ## - Carat.Weight2 1 1.757 3.409 -4012.9 781.8800 < 2.2e-16 *** ## - Clarity 6 5.918 7.569 -3425.4 438.8666 < 2.2e-16 *** ## - Color 2 6.761 8.412 -3338.3 1504.1556 < 2.2e-16 *** ## - Carat.Weight 1 53.195 54.847 -1932.0 23670.6866 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model selected by the stepwise regression is summary(step.reg)

## ## Call: ## lm(formula = LogPrice ~ Carat.Weight + Color + Clarity + Carat.Weight2 + ## Cut, data = diamond)

11

## ## Residuals: ## Min 1Q Median 3Q Max ## -0.152181 -0.031578 -0.002784 0.030240 0.153036 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.851063 0.008434 338.040 < 2e-16 *** ## Carat.Weight 1.226520 0.007972 153.853 < 2e-16 *** ## ColorG -0.095142 0.004092 -23.248 < 2e-16 *** ## ColorK -0.347479 0.006346 -54.755 < 2e-16 *** ## ClaritySI1 -0.286641 0.008359 -34.292 < 2e-16 *** ## ClaritySI2 -0.353371 0.008800 -40.155 < 2e-16 *** ## ClarityVS1 -0.161288 0.008513 -18.947 < 2e-16 *** ## ClarityVS2 -0.215559 0.008246 -26.141 < 2e-16 *** ## ClarityVVS1 -0.077703 0.008682 -8.950 < 2e-16 *** ## ClarityVVS2 -0.103078 0.008327 -12.378 < 2e-16 *** ## Carat.Weight2 -0.699846 0.025028 -27.962 < 2e-16 *** ## CutGood -0.038896 0.007763 -5.010 6.81e-07 *** ## CutIdeal 0.016165 0.006982 2.315 0.020877 * ## CutVery Good -0.014988 0.003834 -3.910 0.000101 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.04741 on 735 degrees of freedom ## Multiple R-squared: 0.9731, Adjusted R-squared: 0.9727 ## F-statistic: 2048 on 13 and 735 DF, p-value: < 2.2e-16

12