STAT270 - Project 3 - Modeling with Regression

# load the packages for graphing and data wrangling
library(ggplot2)
library(PASWR2)
library(car) # install the car package first, if not already installed
library(dplyr) 
library(lattice)
library(boot)
library(MASS)

Note: If you Rmd file submission knits you will receive total of (5 points)

Case Study: Real Estate

Data and ideas for this case study come from (Militino et al., 2004).

Problem 16, page 899 (not entire problem, only specified parts below)

The goal of this case study is to walk the user through the creation of a parsimonious multiple linear regression model that can be used to predict the total price (totalprice) of apartments by their hedonic (structural) characteristics. The data frame VIT2005 contains several variables, and further description of the data can be found in the help file (listed below).

A data frame with 218 observations on the following 15 variables:

• totalprice (the market total price (in Euros) of the apartment including garage(s) and storage room(s))

• area (the total living area of the apartment in square meters)

• zone (a factor indicating the neighborhood where the apartment is located with levels Z11, Z21, Z31, Z32, Z34, Z35, Z36, Z37, Z38, Z41, Z42, Z43, Z44, Z45, Z46, Z47, Z48, Z49, Z52, Z53, Z56, Z61, and Z62)

• category (a factor indicating the condition of the apartment with levels 2A, 2B, 3A, 3B, 4A, 4B, and 5A ordered so that 2A is the best and 5A is the worst)

• age (age of the apartment in years)

• floor (floor on which the apartment is located)

• rooms (total number of rooms including bedrooms, dining room, and kitchen)

• out (a factor indicating the percent of the apartment exposed to the elements: The levels E100, E75, E50, and E25, correspond to complete exposure, 75% exposure, 50% exposure, and 25% exposure, respectively.)

• conservation (is an ordered factor indicating the state of conservation of the apartment. The levels 1A, 2A, 2B, and 3A are ordered from best to worst conservation.)

• toilets (the number of bathrooms)

• garage (the number of garages)

• elevator (indicates the absence (0) or presence (1) of elevators.)

• streetcategory (an ordered factor from best to worst indicating the category of the street with levels S2, S3, S4, and S5)

• heating (a factor indicating the type of heating with levels 1A, 3A, 3B, and 4A which correspond to: no heating, low-standard private heating, high-standard private heating, and central heating, respectively.)

• storage (the number of storage rooms outside of the apartment)

Complete the parts below

(10 pts) Quiz-Project 3 Pr.1

1. Characterize the shape, center, and spread of the variable totalprice.

Hint: Use ggplot function from ggplot2 package to graph the totalprice density function. Use median and IQR to find the median and IQR for the totalprice. Fill in the blank spaces in the observation below:

Solution: Hint: Use the template for the graph:

ggplot(data = YOUR DATA, aes(x = variable to plot)) + geom_density(fill = "your favorite color") + theme_bw()

YOUR CODE HERE:

library(ggplot2) # Plot the density function of ‘totalprice’ ggplot(data = VIT2005, aes(x = totalprice)) + geom_density(fill = “lightblue”) + theme_bw() + ggtitle(“Density Plot of Total Price”)

Calculate the median and IQR for ‘totalprice’

median_totalprice <- median(VIT2005\(totalprice) iqr_totalprice <- IQR(VIT2005\)totalprice)

median_totalprice iqr_totalprice

# MD <- median(VIT2005$totalprice)
# 
# iqr <- IQR(VIT2005$totalprice)
# 
# c(MD, iqr)

Observation: The distribution of totalprice is … with a median of ... and an IQR of ...

(10 pts) Quiz-Project 3 Pr.2

2. Use scatterplotMatrix() from car package or pairs() to explore the relationships between totalprice and the numerical explanatory variables area, age, floor, rooms, toilets, garage, elevator, and storage.

Hint: To use scatterplotMatrix type scatterplotMatrix( ~ totalprice + var1 + var2 + ... + var n, data = VIT2005), do not use more than 5 variables to produce input that fits the screen and can be reviewed. Use the command as many times as you need to review how totalprice correlates with other variables in the data.

Solution:

YOUR CODE HERE: library(car)

scatterplotMatrix( ~ totalprice + area + age + floor + rooms, data = VIT2005, smooth = FALSE)

scatterplotMatrix( ~ totalprice + toilets + garage + elevator + storage, data = VIT2005, smooth = FALSE)

# e.g. matrix can be produced with
scatterplotMatrix( ~ totalprice + area + age + floor + rooms, data = VIT2005)

scatterplotMatrix( ~ totalprice + elevator + storage + toilets, data = VIT2005 )

# or use 
# pairs(~ totalprice + area + age + floor + rooms, data = VIT2005)
# 
# pairs(~ totalprice + toilets + garage + elevator + storage, data = VIT2005)

Observation: The variable totalprice appears to have a moderate linear relationship with area.

3. Total of (55 pts) Quiz-Project 3 Pr.3 to 8, 10 pts for each correctly removed variable, 5 pts to find the correlation.

Compute the correlation between totalprice and all of the other numerical variables. List the three variables in order along with their correlation coefficients that have the highest correlation with totalprice.

Model (A): Use backward elimination to develop a model that predicts totalprice using the data frame VIT2005. Use a “P-value-to remove” of 5%. Store the final model in the object modelA.

(5 pts) Quiz-Project 3 Pr.3

The correlation coefficients are:

NUM <- c("area", "age", "floor", "rooms", "toilets", "garage","elevator", "storage")
COR <- cor(VIT2005[, "totalprice"], VIT2005[, NUM])
COR

Observation: The highest three correlations with totalprice occur with area (0.8092), toilets (0.6876), and rooms(0.5256).

Model (A) The functions drop1() and update() are used to create a model using backward elimination.

Hint: drop1(model_.be_name, test = "F") test for significance of all individual predictors given all others are already in the model.

# use a model where totalprice is regressed on all other variables
model.be <- lm(totalprice ~ ., data = VIT2005)

# 
drop1(model.be, test = "F")

(10 pts) (part c.1) Quiz-Project 3 Pr.4

Which one appears most insignificant (biggest P-value)? Drop it from the model.

E.g. If age is most insignificant, use the update function to update the model with the following code:

model.be <- update(model.be, .~. - age) - the first dot means we use the same response variable, the second all previously included predictors minus age. The again use drop1 and so on, until all remaining variables are of significance 0.05 as specified in the directions of the problem

YOUR CODE HERE: model.be <- lm(totalprice ~ area + age + floor + rooms + toilets + garage + elevator + storage, data = VIT2005)

drop1(model.be, test = “F”) model.be <- update(model.be, . ~ . - age) drop1(model.be, test = “F”) model.be <- update(model.be, . ~ . - floor) drop1(model.be, test = “F”)

model.be <- update(model.be, . ~ . - rooms)

drop1(model.be, test = “F”)

model.be <- update(model.be, . ~ . - storage)

drop1(model.be, test = “F”)

model.be <- update(model.be, . ~ . - elevator)

drop1(model.be, test = “F”)

totalprice ~ area + toilets + garage

model.be <- update(model.be, .~. - age)
drop1(model.be, test = "F")

Single term deletions

Model:
totalprice ~ area + rooms + toilets + garage + elevator + storage
         Df  Sum of Sq        RSS    AIC  F value    Pr(>F)
<none>                 2.3894e+11 4551.7                   
area      1 1.4652e+11 3.8546e+11 4653.9 129.3900 < 2.2e-16
rooms     1 4.4164e+05 2.3894e+11 4549.7   0.0004 0.9842626
toilets   1 1.2991e+10 2.5193e+11 4561.2  11.4725 0.0008425
garage    1 3.2312e+10 2.7125e+11 4577.3  28.5338 2.376e-07
elevator  1 1.5731e+10 2.5467e+11 4563.6  13.8914 0.0002487
storage   1 7.0112e+09 2.4595e+11 4556.0   6.1915 0.0136104
            
<none>      
area     ***
rooms       
toilets  ***
garage   ***
elevator ***
storage  *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(10 pts) (part c.2) Quiz-Project 3 Pr.5

Which one is to be dropped next?

Answer: model.be <- update(model.be, . ~ . - floor) drop1(model.be, test = “F”)

(10 pts) (part c.3) Quiz-Project 3 Pr.6 Which one is to be dropped next?

Answer: model.be <- update(model.be, . ~ . - rooms)

drop1(model.be, test = “F”)

(10 pts) (part c.4) Quiz-Project 3 Pr.7 Which one is to be dropped next?

Answer:

(10 pts) (part c.5) Quiz-Project 3 Pr.8 Which one is to be dropped next?

Answer: model.be <- update(model.be, . ~ . - storage)

If all variable are significant create the object holding the optimal model as per the description of the problem using the code:

formula(model.be)

modelA <- lm(formula(model.be), data = VIT2005)

Question 1: (5 pts) Quiz-Project 3 Pr.9 Which variables are left in the model, list them?

Observation: Backward elimination suggests using the variables area, zone, category, out, toilets, garage, elevator, streetcorner, and heating to best predict totalprice. area

zone

category

out

toilets

garage

elevator

streetcorner

heating

1. (5 pts) Quiz-Project 3 Pr.10

Set the seed for reproducibility

set.seed(5)

Assuming modelA is already fitted, e.g.:

modelA <- lm(totalprice ~ area + toilets + garage + elevator, data = VIT2005)

Load the necessary library

library(boot)

Compute leave-one-out cross-validation error (CV_n)

cv_error_n <- cv.glm(data = VIT2005, glmfit = modelA) cv_error_n$delta[1]

Compute five-fold cross-validation error (CV_5)

cv_error_5 <- cv.glm(data = VIT2005, glmfit = modelA, K = 5) cv_error_5$delta[1]

(ii) **(5 pts)** **Quiz-Project 3 Pr.11**

Compute $R^2$, $R^2_a$, the `AIC`, and the `BIC` for **Model (A)**. What is the proportion of total variability explained by **Model (A)**?

**Your Solution:**
# Fit the model (if not already fitted)
# modelA <- lm(totalprice ~ area + toilets + garage + elevator, data = VIT2005)

# Compute R-squared (R^2)
R2 <- summary(modelA)$r.squared

# Compute Adjusted R-squared (R^2_a)
R2_a <- summary(modelA)$adj.r.squared

# Compute AIC
AIC_value <- AIC(modelA)

# Compute BIC
BIC_value <- BIC(modelA)

# Proportion of total variability explained by Model (A)
proportion_variability_explained <- R2

# Print results
cat("R-squared (R^2):", R2, "\n")
cat("Adjusted R-squared (R^2_a):", R2_a, "\n")
cat("AIC:", AIC_value, "\n")
cat("BIC:", BIC_value, "\n")
cat("Proportion of total variability explained by Model (A):", proportion_variability_explained, "\n")

(i) 


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-output-begin eyJkYXRhIjoiXG48IS0tIHJuYi1zb3VyY2UtYmVnaW4gZXlKa1lYUmhJam9pWUdCZ2NseHViVzlrWld4Qlp5QThMU0JuYkcwb1ptOXliWFZzWVNodGIyUmxiQzVpWlNrc0lHUmhkR0VnUFNCV1NWUXlNREExS1Z4dVhHNGpJQ01nUm05eUlHTjJMbWRzYlNBdElIUm9aU0JrWldaaGRXeDBJR2x6SUhSdklITmxkQ0JMSUdWeGRXRnNJSFJ2SUhSb1pTQnVkVzFpWlhJZ2IyWWdiMkp6WlhKMllYUnBiMjV6SUdsdUlHUmhkR0VnZDJocFkyZ2daMmwyWlhNZ2RHaGxJSFZ6ZFdGc0lHeGxZWFpsTFc5dVpTMXZkWFFnWTNKdmMzTXRkbUZzYVdSaGRHbHZiaTVjYmlNZ1hHNXpaWFF1YzJWbFpDZzFLU0FqSUhWelpTQm1iM0lnY21Wd2JHbGpZWFJwYjI0Z2NIVnljRzl6WlhOY2JtTjJMbVZ5Y205eVRpQThMU0JqZGk1bmJHMG9aR0YwWVNBOUlGWkpWREl3TURVc0lHZHNiV1pwZENBOUlDQnRiMlJsYkVGbktWeHVYRzVEVms1aElEd3RJR04yTG1WeWNtOXlUaVJrWld4MFlWc3hYVnh1UTFaT1lWeHVZR0JnSW4wPSAtLT5cblxuYGBgclxubW9kZWxBZyA8LSBnbG0oZm9ybXVsYShtb2RlbC5iZSksIGRhdGEgPSBWSVQyMDA1KVxuXG4jICMgRm9yIGN2LmdsbSAtIHRoZSBkZWZhdWx0IGlzIHRvIHNldCBLIGVxdWFsIHRvIHRoZSBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zIGluIGRhdGEgd2hpY2ggZ2l2ZXMgdGhlIHVzdWFsIGxlYXZlLW9uZS1vdXQgY3Jvc3MtdmFsaWRhdGlvbi5cbiMgXG5zZXQuc2VlZCg1KSAjIHVzZSBmb3IgcmVwbGljYXRpb24gcHVycG9zZXNcbmN2LmVycm9yTiA8LSBjdi5nbG0oZGF0YSA9IFZJVDIwMDUsIGdsbWZpdCA9ICBtb2RlbEFnKVxuXG5DVk5hIDwtIGN2LmVycm9yTiRkZWx0YVsxXVxuQ1ZOYVxuYGBgXG5cbjwhLS0gcm5iLXNvdXJjZS1lbmQgLS0+XG4ifQ== -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubW9kZWxBZyA8LSBnbG0oZm9ybXVsYShtb2RlbC5iZSksIGRhdGEgPSBWSVQyMDA1KVxuXG4jICMgRm9yIGN2LmdsbSAtIHRoZSBkZWZhdWx0IGlzIHRvIHNldCBLIGVxdWFsIHRvIHRoZSBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zIGluIGRhdGEgd2hpY2ggZ2l2ZXMgdGhlIHVzdWFsIGxlYXZlLW9uZS1vdXQgY3Jvc3MtdmFsaWRhdGlvbi5cbiMgXG5zZXQuc2VlZCg1KSAjIHVzZSBmb3IgcmVwbGljYXRpb24gcHVycG9zZXNcbmN2LmVycm9yTiA8LSBjdi5nbG0oZGF0YSA9IFZJVDIwMDUsIGdsbWZpdCA9ICBtb2RlbEFnKVxuXG5DVk5hIDwtIGN2LmVycm9yTiRkZWx0YVsxXVxuQ1ZOYVxuYGBgIn0= -->

```r
modelAg <- glm(formula(model.be), data = VIT2005)

# # For cv.glm - the default is to set K equal to the number of observations in data which gives the usual leave-one-out cross-validation.
# 
set.seed(5) # use for replication purposes
cv.errorN <- cv.glm(data = VIT2005, glmfit =  modelAg)

CVNa <- cv.errorN$delta[1]
CVNa

modelAg <- glm(formula(model.be), data = VIT2005)

# UNCOMMENT the code below and run it

# set.seed(5) # use for replication purposes
# cv.error5 <- cv.glm(data = VIT2005, glmfit = modelAg, K = 5)
# CV5a <- cv.error5$delta[1]
# CV5a

Observation: The \(CV_n = ...\) for Model (A), and \(CV_5 =...\) for Model (A).

2. Since this problem and a few more will request many goodness of fit statistics, a function called mgof() is written to compute the requested values. Use it as shown below.

mgof <- function(model = model, data = DF, ...){
  R2a <- summary(model)$adj.r.squared
  R2 <- summary(model)$r.squared
  aic <- AIC(model)
  bic <- AIC(model, k = log(nrow(data)))
  se <- summary(model)$sigma
  form <- formula(model)
  ANS <- c(R2 = R2, R2.adj = R2a, AIC = aic, BIC = bic, SE = se)
  ANS
}

MGOF <- mgof(model = modelA, data = VIT2005)
MGOF

Observation: The total proportion of variability (\(R^2\)) explained by modelA is 0.9138.

4. (10 pts) Quiz-Project 3 Pr.12

Explore the residuals of the Models (A) using the function residualPlot() or residualPlots() from the package car. Comment on the results. (Diagnostics: Checking the model assumptions)

# uncomment to run the residual plot for the modelA
# residualPlot(modelA, main = "Model A")

Load the necessary library

library(car)

Plot the residuals for Model A

residualPlot(modelA, main = “Residual Plot for Model A”)

(d) Question 2: Is there curvature of the residuals on the above plot?

Observation: The residuals versus the fitted values for Model (A) have a definite curvature indicating the model is not quite adequate.

Extra Credit (10 pts) Quiz-Project 3 Pr.13

(e) Use the function boxCox() from car to find a suitable transformation for totalprice. Build Model (E) Use backward elimination to develop a model that predicts log(totalprice) using the data frame VIT2005. Use a “P-value-to remove” of 5%. Store the final model in the object modelE.

For details on boxCox() see page 856 in text

boxCox(modelA, lambda = seq(-0.5, 0.5, length = 200))

Observation:

A log transformation is suggested for the response totalprice in Model (A). # Load necessary libraries library(car)

Step 1: Find the Suitable Transformation using boxCox()

Applying the Box-Cox transformation for modelA to determine lambda

boxCox(modelA, lambda = seq(-0.5, 0.5, length = 200))

Step 2: Create Model (E) by transforming the response variable using log

Log transformation of totalprice

modelE <- lm(log(totalprice) ~ area + age + floor + rooms + toilets + garage + elevator + storage, data = VIT2005)

Step 3: Apply Backward Elimination using p-value of 5% for removal

Start with the full model and use step() to perform backward elimination

modelE_final <- step(modelE, direction = “backward”, k = log(nrow(VIT2005)))

Display the final model after backward elimination

summary(modelE_final)

Extra Credit (20 pts) Quiz-Project 3 Pr.14 - 4 pts for excluding the correct variable at each step.

Model (E):

Use backward elimination to develop a model that predicts log(totalprice) using the data frame VIT2005. Use a “P-value-to remove” of 5%. Store the final model in the object modelE.

Model (E) The functions drop1() and update() are used to create a model using backward elimination. (as shown at the beginning)

Load necessary libraries

library(car)

Step 1: Log transformation of totalprice

modelE <- lm(log(totalprice) ~ area + age + floor + rooms + toilets + garage + elevator + storage, data = VIT2005)

Step 2: Backward elimination using drop1() and update()

Initial model

summary(modelE)

Step 3: Perform backward elimination

Drop variable with highest p-value > 0.05

drop1(modelE, test = “F”) # Check for p-values

Suppose ‘storage’ has the highest p-value (greater than 0.05) in the initial model

modelE <- update(modelE, . ~ . - storage)

Check the model again

drop1(modelE, test = “F”)

Suppose ‘age’ has the next highest p-value > 0.05

modelE <- update(modelE, . ~ . - age)

Check the model again

drop1(modelE, test = “F”)

Suppose ‘floor’ has the next highest p-value > 0.05

modelE <- update(modelE, . ~ . - floor)

Check the model again

drop1(modelE, test = “F”)

Continue this process until all variables have p-values < 0.05

The final model (modelE) should have predictors with significant p-values

Step 4: Final Model

summary(modelE)

# transform the response variable using natural log transformation

VITNEW <- VIT2005 %>% mutate(logtotalprice = log(totalprice))

# build the fitted full model
model.be <- lm(logtotalprice ~ ., data = VITNEW[ ,-1]) # exclude the totalprice variable since it is no longer response variable
drop1(model.be, test = "F")

and so on.