Housing
Anoosha Qaisar
2026-07-09
housing <- read.csv("C:/Users/anoos/OneDrive/Desktop/Housing.csv")
str(housing)## 'data.frame': 545 obs. of 13 variables:
## $ price : int 13300000 12250000 12250000 12215000 11410000 10850000 10150000 10150000 9870000 9800000 ...
## $ area : int 7420 8960 9960 7500 7420 7500 8580 16200 8100 5750 ...
## $ bedrooms : int 4 4 3 4 4 3 4 5 4 3 ...
## $ bathrooms : int 2 4 2 2 1 3 3 3 1 2 ...
## $ stories : int 3 4 2 2 2 1 4 2 2 4 ...
## $ mainroad : chr "yes" "yes" "yes" "yes" ...
## $ guestroom : chr "no" "no" "no" "no" ...
## $ basement : chr "no" "no" "yes" "yes" ...
## $ hotwaterheating : chr "no" "no" "no" "no" ...
## $ airconditioning : chr "yes" "yes" "no" "yes" ...
## $ parking : int 2 3 2 3 2 2 2 0 2 1 ...
## $ prefarea : chr "yes" "no" "yes" "yes" ...
## $ furnishingstatus: chr "furnished" "furnished" "semi-furnished" "furnished" ...head(housing)## price area bedrooms bathrooms stories mainroad guestroom basement
## 1 13300000 7420 4 2 3 yes no no
## 2 12250000 8960 4 4 4 yes no no
## 3 12250000 9960 3 2 2 yes no yes
## 4 12215000 7500 4 2 2 yes no yes
## 5 11410000 7420 4 1 2 yes yes yes
## 6 10850000 7500 3 3 1 yes no yes
## hotwaterheating airconditioning parking prefarea furnishingstatus
## 1 no yes 2 yes furnished
## 2 no yes 3 no furnished
## 3 no no 2 yes semi-furnished
## 4 no yes 3 yes furnished
## 5 no yes 2 no furnished
## 6 no yes 2 yes semi-furnished# 1. Convert all categoricals to factor
housing$mainroad <- as.factor(housing$mainroad)
housing$guestroom <- as.factor(housing$guestroom)
housing$basement <- as.factor(housing$basement)
housing$hotwaterheating <- as.factor(housing$hotwaterheating)
housing$airconditioning <- as.factor(housing$airconditioning)
housing$prefarea <- as.factor(housing$prefarea)
housing$furnishingstatus <- as.factor(housing$furnishingstatus)
str(housing)## 'data.frame': 545 obs. of 13 variables:
## $ price : int 13300000 12250000 12250000 12215000 11410000 10850000 10150000 10150000 9870000 9800000 ...
## $ area : int 7420 8960 9960 7500 7420 7500 8580 16200 8100 5750 ...
## $ bedrooms : int 4 4 3 4 4 3 4 5 4 3 ...
## $ bathrooms : int 2 4 2 2 1 3 3 3 1 2 ...
## $ stories : int 3 4 2 2 2 1 4 2 2 4 ...
## $ mainroad : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ guestroom : Factor w/ 2 levels "no","yes": 1 1 1 1 2 1 1 1 2 2 ...
## $ basement : Factor w/ 2 levels "no","yes": 1 1 2 2 2 2 1 1 2 1 ...
## $ hotwaterheating : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ airconditioning : Factor w/ 2 levels "no","yes": 2 2 1 2 2 2 2 1 2 2 ...
## $ parking : int 2 3 2 3 2 2 2 0 2 1 ...
## $ prefarea : Factor w/ 2 levels "no","yes": 2 1 2 2 1 2 2 1 2 2 ...
## $ furnishingstatus: Factor w/ 3 levels "furnished","semi-furnished",..: 1 1 2 1 1 2 2 3 1 3 ...n <- nrow(housing)
names(housing)## [1] "price" "area" "bedrooms" "bathrooms"
## [5] "stories" "mainroad" "guestroom" "basement"
## [9] "hotwaterheating" "airconditioning" "parking" "prefarea"
## [13] "furnishingstatus"library(corrplot)## corrplot 0.95 loadedlibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.6.1##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnumeric_data <- housing %>% select_if(is.numeric)
categorical_data <- housing %>% select_if(is.factor)
col_col <- c("price" ,"area","bedrooms" ,"bathrooms", "stories" , "parking" )
corrplot(cor(numeric_data), method = "number")
corrplot(cor(housing[,col_col]), method = "number")
covariates1 <- names(categorical_data)
for (var in covariates1) {
boxplot(housing$price ~ categorical_data[[var]] , xlab = var )
}
pairs(numeric_data) 
Based on the corrplot, the response variable price shows a strong positive correlation with area and bathrooms. However, the correlation between numerical predictors suggests the presence of multicollinearity.
According to the pairs plot, there is a strong positive correlation between area and price. The observed ‘funnel shape’ indicates heteroscedasticity and right-skewness in the data. There is weak correlation among the other predictors, as evidenced by the vertical clusters. Finally, the boxplots reveal heteroscedasticity, as indicated by the differing medians and unequal variability between ‘yes’ and ‘no’ categories.
model.fit <- lm(price ~ . , data = housing)
summary(model.fit) ##
## Call:
## lm(formula = price ~ ., data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2619718 -657322 -68409 507176 5166695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42771.69 264313.31 0.162 0.871508
## area 244.14 24.29 10.052 < 2e-16 ***
## bedrooms 114787.56 72598.66 1.581 0.114445
## bathrooms 987668.11 103361.98 9.555 < 2e-16 ***
## stories 450848.00 64168.93 7.026 6.55e-12 ***
## mainroadyes 421272.59 142224.13 2.962 0.003193 **
## guestroomyes 300525.86 131710.22 2.282 0.022901 *
## basementyes 350106.90 110284.06 3.175 0.001587 **
## hotwaterheatingyes 855447.15 223152.69 3.833 0.000141 ***
## airconditioningyes 864958.31 108354.51 7.983 8.91e-15 ***
## parking 277107.10 58525.89 4.735 2.82e-06 ***
## prefareayes 651543.80 115682.34 5.632 2.89e-08 ***
## furnishingstatussemi-furnished -46344.62 116574.09 -0.398 0.691118
## furnishingstatusunfurnished -411234.39 126210.56 -3.258 0.001192 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1068000 on 531 degrees of freedom
## Multiple R-squared: 0.6818, Adjusted R-squared: 0.674
## F-statistic: 87.52 on 13 and 531 DF, p-value: < 2.2e-16library(car)## Warning: package 'car' was built under R version 4.6.1## Loading required package: carData## Warning: package 'carData' was built under R version 4.6.1##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(model.fit)## GVIF Df GVIF^(1/(2*Df))
## area 1.325250 1 1.151195
## bedrooms 1.369477 1 1.170246
## bathrooms 1.286621 1 1.134293
## stories 1.478055 1 1.215753
## mainroad 1.172728 1 1.082926
## guestroom 1.212838 1 1.101289
## basement 1.323050 1 1.150239
## hotwaterheating 1.041506 1 1.020542
## airconditioning 1.211840 1 1.100836
## parking 1.212837 1 1.101289
## prefarea 1.149196 1 1.072006
## furnishingstatus 1.109639 2 1.026350The linear regression summary confirms that area, bathrooms, and airconditioning are the most influential predictors of price (\(p < 0.001\)). Notably, bedrooms and furnishingstatussemi-furnished were not statistically significant in this model, suggesting they may provide redundant information. The adjusted \(R^2\) of 0.674 indicates a moderate model fit, though the wide range in residuals further corroborates the presence of heteroscedasticity.”To formally assess multicollinearity, Generalized Variance Inflation Factors (GVIF) were calculated. All \(GVIF^{\frac{1}{2 \cdot Df}}\) values were found to be below 1.25, indicating that multicollinearity is not present to a problematic degree. This confirms that the model’s standard errors are not inflated and that the individual coefficient estimates are stable and reliable for interpretation.
# Create the new log-transformed model
model.log <- lm(log(price) ~ ., data = housing)
# Compare the summary
summary(model.log)##
## Call:
## lm(formula = log(price) ~ ., data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60723 -0.12531 0.00631 0.13088 0.63018
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.435e+01 5.115e-02 280.577 < 2e-16 ***
## area 4.949e-05 4.701e-06 10.528 < 2e-16 ***
## bedrooms 2.947e-02 1.405e-02 2.098 0.036411 *
## bathrooms 1.632e-01 2.000e-02 8.159 2.47e-15 ***
## stories 9.013e-02 1.242e-02 7.258 1.41e-12 ***
## mainroadyes 1.175e-01 2.752e-02 4.268 2.34e-05 ***
## guestroomyes 6.979e-02 2.549e-02 2.738 0.006391 **
## basementyes 8.931e-02 2.134e-02 4.185 3.34e-05 ***
## hotwaterheatingyes 1.634e-01 4.319e-02 3.784 0.000172 ***
## airconditioningyes 1.748e-01 2.097e-02 8.336 6.62e-16 ***
## parking 4.477e-02 1.133e-02 3.953 8.77e-05 ***
## prefareayes 1.268e-01 2.239e-02 5.661 2.46e-08 ***
## furnishingstatussemi-furnished 1.797e-02 2.256e-02 0.796 0.426115
## furnishingstatusunfurnished -1.099e-01 2.443e-02 -4.500 8.36e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2067 on 531 degrees of freedom
## Multiple R-squared: 0.699, Adjusted R-squared: 0.6916
## F-statistic: 94.84 on 13 and 531 DF, p-value: < 2.2e-16# Plot diagnostics to see if the funnel shape is gone
plot(model.log)
The log-transformed regression model provides a more accurate and
interpretable fit, successfully addressing the heteroscedasticity found
in the initial analysis. The model demonstrates that key amenities such
as air conditioning, the number of bathrooms, and the preference area
have the most substantial positive impacts on property value, while the
log-transformation allows these effects to be expressed as intuitive
percentage changes.”“The ‘Residuals vs. Leverage’ plot indicates that
the model is stable. There are no observations exceeding the Cook’s
distance boundaries, which confirms that no single data point is
disproportionately influencing the regression results. While points 24,
140, and 379 represent houses with slightly unusual characteristics or
larger prediction errors, they do not constitute influential outliers
that necessitate removal from the dataset.
mod.red <- update(model.log, ~ . -furnishingstatus)
summary(mod.red)##
## Call:
## lm(formula = log(price) ~ area + bedrooms + bathrooms + stories +
## mainroad + guestroom + basement + hotwaterheating + airconditioning +
## parking + prefarea, data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.67787 -0.12197 0.01574 0.12792 0.67737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.427e+01 4.726e-02 302.057 < 2e-16 ***
## area 5.055e-05 4.855e-06 10.412 < 2e-16 ***
## bedrooms 3.409e-02 1.451e-02 2.349 0.01920 *
## bathrooms 1.684e-01 2.066e-02 8.150 2.60e-15 ***
## stories 9.177e-02 1.284e-02 7.148 2.91e-12 ***
## mainroadyes 1.302e-01 2.835e-02 4.593 5.47e-06 ***
## guestroomyes 7.393e-02 2.634e-02 2.807 0.00519 **
## basementyes 9.948e-02 2.201e-02 4.520 7.63e-06 ***
## hotwaterheatingyes 1.783e-01 4.460e-02 3.997 7.31e-05 ***
## airconditioningyes 1.769e-01 2.159e-02 8.195 1.88e-15 ***
## parking 5.113e-02 1.166e-02 4.386 1.39e-05 ***
## prefareayes 1.276e-01 2.314e-02 5.515 5.45e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2138 on 533 degrees of freedom
## Multiple R-squared: 0.6767, Adjusted R-squared: 0.6701
## F-statistic: 101.4 on 11 and 533 DF, p-value: < 2.2e-16anova(mod.red, model.log)## Analysis of Variance Table
##
## Model 1: log(price) ~ area + bedrooms + bathrooms + stories + mainroad +
## guestroom + basement + hotwaterheating + airconditioning +
## parking + prefarea
## Model 2: log(price) ~ area + bedrooms + bathrooms + stories + mainroad +
## guestroom + basement + hotwaterheating + airconditioning +
## parking + prefarea + furnishingstatus
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 533 24.356
## 2 531 22.682 2 1.6738 19.592 6.174e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To evaluate the importance of the furnishingstatus variable, an Analysis of Variance (ANOVA) was conducted comparing a model excluding this feature against the full model. The resulting F-statistic (\(F(2, 531) = 19.59, p < 0.001\)) indicates that furnishingstatus significantly improves the model’s explanatory power. Therefore, retaining this variable is necessary to minimize residual error and maintain the model’s predictive integrity.
# Instead of using just 'area'
model.poly <- lm(log(price) ~ poly(area, 2) + bedrooms + bathrooms + stories +
mainroad + guestroom + basement + hotwaterheating +
airconditioning + parking + prefarea + furnishingstatus,
data = housing)
# View the results
summary(model.poly)##
## Call:
## lm(formula = log(price) ~ poly(area, 2) + bedrooms + bathrooms +
## stories + mainroad + guestroom + basement + hotwaterheating +
## airconditioning + parking + prefarea + furnishingstatus,
## data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60764 -0.12468 0.00489 0.12529 0.61363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.63185 0.05097 287.067 < 2e-16 ***
## poly(area, 2)1 2.60834 0.23754 10.980 < 2e-16 ***
## poly(area, 2)2 -0.73949 0.21732 -3.403 0.000718 ***
## bedrooms 0.03081 0.01392 2.214 0.027267 *
## bathrooms 0.16205 0.01981 8.180 2.12e-15 ***
## stories 0.08756 0.01232 7.107 3.85e-12 ***
## mainroadyes 0.10340 0.02757 3.751 0.000196 ***
## guestroomyes 0.05440 0.02564 2.122 0.034337 *
## basementyes 0.09172 0.02115 4.337 1.73e-05 ***
## hotwaterheatingyes 0.16485 0.04277 3.855 0.000130 ***
## airconditioningyes 0.16415 0.02100 7.817 2.93e-14 ***
## parking 0.04140 0.01126 3.677 0.000260 ***
## prefareayes 0.12565 0.02217 5.667 2.38e-08 ***
## furnishingstatussemi-furnished 0.01376 0.02237 0.615 0.538931
## furnishingstatusunfurnished -0.11231 0.02420 -4.642 4.36e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2047 on 530 degrees of freedom
## Multiple R-squared: 0.7054, Adjusted R-squared: 0.6976
## F-statistic: 90.65 on 14 and 530 DF, p-value: < 2.2e-16“To refine the model’s accuracy, a quadratic term for area was introduced using the poly() function. This adjustment proved statistically significant (\(p < 0.001\)), confirming that the relationship between house area and price is non-linear, likely reflecting diminishing marginal returns for larger properties. This model achieves an Adjusted \(R^2\) of 0.6976, representing the most effective and statistically valid fit for the data.”
# Create new data with all variables used in the model
new.data <- data.frame(
area = 17000,
bedrooms = 3, # Typical value
bathrooms = 2, # Typical value
stories = 2, # Typical value
mainroad = factor("yes", levels = levels(housing$mainroad)),
guestroom = factor("no", levels = levels(housing$guestroom)),
basement = factor("no", levels = levels(housing$basement)),
hotwaterheating = factor("no", levels = levels(housing$hotwaterheating)),
airconditioning = factor("yes", levels = levels(housing$airconditioning)),
parking = 1, # Typical value
prefarea = factor("no", levels = levels(housing$prefarea)),
furnishingstatus = factor("furnished", levels = levels(housing$furnishingstatus))
)
# Now the prediction will work
predicted_log_price <- predict(model.poly, newdata = new.data)
# To get the price in actual dollars, use exp()
predicted_price <- exp(predicted_log_price)
predicted_price## 1
## 6888963confint <- predict(model.poly, newdata = new.data, interval = "confidence")
confint## fit lwr upr
## 1 15.74543 15.49901 15.99185exp(confint)## fit lwr upr
## 1 6888963 5384358 8814016pred <- predict(model.poly, newdata = new.data, interval = "prediction")
pred## fit lwr upr
## 1 15.74543 15.27389 16.21697exp(pred)## fit lwr upr
## 1 6888963 4298995 11039282# Basic scatter plot
plot(housing$area, log(housing$price),
main = "Model Fit: Area vs. Log(Price)",
xlab = "Area", ylab = "Log(Price)",
col = "gray", pch = 16)
# Predict across the range of your data
x_range <- seq(min(housing$area), max(housing$area), length.out = 100)
# We use the median of other variables for a "typical" house
line_pred <- predict(model.poly, newdata = data.frame(
area = x_range, bedrooms = median(housing$bedrooms),
bathrooms = median(housing$bathrooms), stories = median(housing$stories),
mainroad = "yes", guestroom = "no", basement = "no",
hotwaterheating = "no", airconditioning = "yes",
parking = median(housing$parking), prefarea = "no",
furnishingstatus = "furnished"))
# Draw the blue line
lines(x_range, line_pred, col = "blue", lwd = 3)
plot(model.poly)
prediction interval is always wider then the confidence interval . This is because a confidence interval only reflects the uncertainty in estimating the average price (reflects the uncertainty in the parameter estimators) for a given housing charateristics. A prediction interval, however, must also account for the random variability (the error term) of a single, specific new observation around that average.
“The ‘Residuals vs. Leverage’ diagnostic plot for the refined quadratic model confirms overall model stability. While observation #8 exhibits high leverage due to its unique characteristics,You still do not see the faint, curved dashed lines (Cook’s distance boundaries) encompassing any of your points. it remains well within the Cook’s distance boundaries. This indicates that no single observation is exerting undue influence on the model’s coefficients, and the results remain robust.” In your housing model, the fact that you did not see these curved lines (or rather, that all your points remained within the boundaries) is a strong indicator of model stability. It confirms that your findings are based on the overall trend of the housing data, rather than being skewed by a few extreme or erroneous observations.If you ever see a plot
where points are outside these lines, it means those specific observations (e.g., a specific house in your dataset) have a disproportionate impact on your model results compared to the rest of the data.
an <- aov(price ~ furnishingstatus, data = housing)
summary(an)## Df Sum Sq Mean Sq F value Pr(>F)
## furnishingstatus 2 1.798e+14 8.99e+13 28.27 2.09e-12 ***
## Residuals 542 1.723e+15 3.18e+12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukeyf <- TukeyHSD(an)
plot(tukeyf)
“An ANOVA test confirmed that furnishingstatus is a significant
predictor of housing prices (\(F(2, 542) =
28.27, p < 0.001\)). Subsequent Tukey post-hoc analysis
revealed significant differences between all pairs of furnishing
categories. Specifically, prices decrease progressively from fully
furnished to semi-furnished, and finally to unfurnished, with all
pairwise differences being statistically significant.”
The Goal: This test determines if the average price of houses is significantly different across the three categories of furnishingstatus (furnished, semi-furnished, unfurnished).
F-value (28.27): This is high, indicating that the variation between the groups is much larger than the variation within the groups.
Pr(>F) (2.09e-12): Since this is significantly less than 0.05, we conclude that furnishing status has a statistically significant impact on house price.
model.glm <- glm(price ~ poly(area, 2) + bedrooms + bathrooms + stories +
mainroad + guestroom + basement + hotwaterheating +
airconditioning + parking + prefarea + furnishingstatus,
data = housing, family = Gamma(link = "log"))
summary(model.glm)##
## Call:
## glm(formula = price ~ poly(area, 2) + bedrooms + bathrooms +
## stories + mainroad + guestroom + basement + hotwaterheating +
## airconditioning + parking + prefarea + furnishingstatus,
## family = Gamma(link = "log"), data = housing)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.63419 0.05115 286.102 < 2e-16 ***
## poly(area, 2)1 2.63119 0.23838 11.038 < 2e-16 ***
## poly(area, 2)2 -0.70867 0.21809 -3.249 0.001230 **
## bedrooms 0.03558 0.01397 2.547 0.011133 *
## bathrooms 0.15954 0.01988 8.025 6.59e-15 ***
## stories 0.08830 0.01236 7.142 3.05e-12 ***
## mainroadyes 0.10634 0.02766 3.844 0.000136 ***
## guestroomyes 0.05622 0.02573 2.185 0.029346 *
## basementyes 0.09583 0.02122 4.516 7.78e-06 ***
## hotwaterheatingyes 0.17266 0.04292 4.023 6.58e-05 ***
## airconditioningyes 0.15978 0.02107 7.582 1.53e-13 ***
## parking 0.04610 0.01130 4.080 5.21e-05 ***
## prefareayes 0.11528 0.02225 5.181 3.14e-07 ***
## furnishingstatussemi-furnished 0.01254 0.02245 0.559 0.576637
## furnishingstatusunfurnished -0.10591 0.02428 -4.362 1.55e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.04217857)
##
## Null deviance: 76.501 on 544 degrees of freedom
## Residual deviance: 22.091 on 530 degrees of freedom
## AIC: 16520
##
## Number of Fisher Scoring iterations: 5“The coefficients represent the effect of each predictor on the natural log of the house price. For example, the airconditioning variable has a coefficient of 0.16, indicating a positive contribution to the log-price, while furnishingstatusunfurnished has a coefficient of -0.10, indicating a negative impact on the log-price, holding other factors constant.”
par(mfrow = c(2,2))
plot(model.glm)
par(mfrow = c(1,1))“Diagnostic analysis of the Gamma GLM confirms the model’s validity and robustness. The Residuals vs. Fitted and Scale-Location plots demonstrate consistent variance, confirming the absence of heteroscedasticity. The Q-Q Residuals plot indicates that the model errors align well with the theoretical Gamma distribution, and the Residuals vs. Leverage plot confirms that no single observation exerts undue influence on the model’s coefficients (all Cook’s distances < 0.5). Consequently, the model is considered highly reliable for predictive purposes
# 1. Create a sequence of area values (as you did)
xp <- seq(min(housing$area), max(housing$area), length = 100)
# 2. Create the data frame with ALL variables
new <- data.frame(
area = xp,
bedrooms = median(housing$bedrooms),
bathrooms = median(housing$bathrooms),
stories = median(housing$stories),
mainroad = factor("yes", levels = levels(housing$mainroad)),
guestroom = factor("no", levels = levels(housing$guestroom)),
basement = factor("no", levels = levels(housing$basement)),
hotwaterheating = factor("no", levels = levels(housing$hotwaterheating)),
airconditioning = factor("yes", levels = levels(housing$airconditioning)),
parking = median(housing$parking),
prefarea = factor("no", levels = levels(housing$prefarea)),
furnishingstatus = factor("furnished", levels = levels(housing$furnishingstatus))
)
# 3. Predict using 'type = "response"'
# This gives you the predicted price directly in currency (not log scale)
predictions <- predict(model.glm, newdata = new, type = "response")
# 4. Plot the result
plot(housing$area, housing$price, col = "gray", pch = 16,
xlab = "Area", ylab = "Price", main = "Gamma GLM Predicted Price Curve")
lines(xp, predictions, col = "red", lwd = 3)
# 1. Create a sequence for the variable you want to explore
# Since bathrooms are discrete numbers, we use unique values instead of seq()
bath_levels <- sort(unique(housing$bathrooms))
# 2. Create the data frame, keeping area and others at their median
new_bath <- data.frame(
area = median(housing$area),
bedrooms = median(housing$bedrooms),
bathrooms = bath_levels, # This is the variable we are exploring
stories = median(housing$stories),
mainroad = factor("yes", levels = levels(housing$mainroad)),
guestroom = factor("no", levels = levels(housing$guestroom)),
basement = factor("no", levels = levels(housing$basement)),
hotwaterheating = factor("no", levels = levels(housing$hotwaterheating)),
airconditioning = factor("yes", levels = levels(housing$airconditioning)),
parking = median(housing$parking),
prefarea = factor("no", levels = levels(housing$prefarea)),
furnishingstatus = factor("furnished", levels = levels(housing$furnishingstatus))
)
# 3. Predict and Plot
predictions_bath <- predict(model.glm, newdata = new_bath, type = "response")
plot(new_bath$bathrooms, predictions_bath, type = "b", pch = 19, col = "blue",
xlab = "Number of Bathrooms", ylab = "Predicted Price",
main = "Effect of Bathrooms on Predicted Price")
# Poisson regression for count data
model.poisson <- glm(bedrooms ~ area + stories,
data = housing,
family = poisson(link = "log"))
# Check for overdispersion
summary(model.poisson)##
## Call:
## glm(formula = bedrooms ~ area + stories, family = poisson(link = "log"),
## data = housing)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.198e-01 7.876e-02 10.409 < 2e-16 ***
## area 1.301e-05 1.140e-05 1.142 0.254
## stories 1.079e-01 2.735e-02 3.947 7.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 99.167 on 544 degrees of freedom
## Residual deviance: 81.770 on 542 degrees of freedom
## AIC: 1697.2
##
## Number of Fisher Scoring iterations: 4“The model goodness-of-fit was assessed by comparing the residual deviance to the degrees of freedom. With a residual deviance of 81.77 on 542 degrees of freedom, the ratio is approximately 0.15. This suggests the model captures the data structure well, though the low ratio indicates potential underdispersion, which confirms the model is not over-estimating the variability of the bedroom counts.”
# Poisson regression for count data
ad.poisson <- glm(bedrooms ~ area + stories + bathrooms,
data = housing,
family = poisson(link = "log"))
# Check for overdispersion
summary(ad.poisson)##
## Call:
## glm(formula = bedrooms ~ area + stories + bathrooms, family = poisson(link = "log"),
## data = housing)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.331e-01 8.700e-02 8.426 < 2e-16 ***
## area 8.017e-06 1.164e-05 0.689 0.49092
## stories 8.623e-02 2.900e-02 2.973 0.00295 **
## bathrooms 1.168e-01 5.039e-02 2.317 0.02049 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 99.167 on 544 degrees of freedom
## Residual deviance: 76.531 on 541 degrees of freedom
## AIC: 1693.9
##
## Number of Fisher Scoring iterations: 4peaser_chi <- residuals(ad.poisson, type = "pearson")
phi <- sum((peaser_chi)^2)/ad.poisson$df.residual
phi## [1] 0.1444885lambdahat <- fitted(ad.poisson)
plot(lambdahat, peaser_chi,
xlab = expression(hat(lambda)), ylab = "Pearson Residuals",
main = "Mean vs. Pearson Residuals")
abline(h = 0, lty = 2)
A dispersion check was conducted on the Poisson regression model, yielding a dispersion parameter (\(\phi\)) of 0.144. This indicates the presence of underdispersion, likely due to the restricted range of the bedroom count variable. To ensure robust statistical inference, a Quasipoisson model was subsequently employed to adjust standard errors. The results confirm that the significance of stories and bathrooms persists under this conservative estimation, validating the findings.
“Diagnostic analysis of Pearson residuals vs. predicted means revealed distinct diagonal banding, a characteristic pattern indicating that the Poisson model is failing to account for the discrete, ordinal nature of the bedrooms variable. Given that the count range is tightly constrained and follows an ordered hierarchy, an Ordinal Logistic Regression approach is more appropriate. This method respects the ranking of bedroom counts and avoids the distributional assumptions of the Poisson model, which were not met by the observed data.”
validation set approach
n <- nrow(housing)
Train <- sample(1:n, 0.5 * n )
Train_house <- housing[Train,]
Test <- housing[-Train, ]
lm_fit <- lm(price ~ area, data = Train_house)
summary(lm_fit)##
## Call:
## lm(formula = price ~ area, data = Train_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4973015 -1029597 -305765 739996 6323038
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.336e+06 2.636e+05 8.863 <2e-16 ***
## area 4.741e+02 4.693e+01 10.104 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1639000 on 270 degrees of freedom
## Multiple R-squared: 0.2744, Adjusted R-squared: 0.2717
## F-statistic: 102.1 on 1 and 270 DF, p-value: < 2.2e-16pred_val <- predict(lm_fit , newdata = Test)
mse <- mean((Test$price - pred_val )^2)
mse## [1] 2.312563e+12po_fit <- lm(price ~ poly(area,degree = 2), data = Train_house)
summary(po_fit)##
## Call:
## lm(formula = price ~ poly(area, degree = 2), data = Train_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3736278 -1057860 -260657 785636 6089390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4802617 96734 49.648 < 2e-16 ***
## poly(area, degree = 2)1 16562472 1595371 10.382 < 2e-16 ***
## poly(area, degree = 2)2 -6394179 1595371 -4.008 7.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1595000 on 269 degrees of freedom
## Multiple R-squared: 0.3152, Adjusted R-squared: 0.3102
## F-statistic: 61.92 on 2 and 269 DF, p-value: < 2.2e-16pred_val_po <- predict(po_fit , newdata = Test)
mse_po <- mean((Test$price - pred_val_po)^2)
mse_po## [1] 2.262949e+12po3_fit <- lm(price ~ poly(area,degree = 3), data = Train_house)
summary(po3_fit)##
## Call:
## lm(formula = price ~ poly(area, degree = 3), data = Train_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3901932 -995963 -230368 832135 5948878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4802617 96369 49.836 < 2e-16 ***
## poly(area, degree = 3)1 16562472 1589358 10.421 < 2e-16 ***
## poly(area, degree = 3)2 -6394179 1589358 -4.023 7.48e-05 ***
## poly(area, degree = 3)3 -2770883 1589358 -1.743 0.0824 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1589000 on 268 degrees of freedom
## Multiple R-squared: 0.3229, Adjusted R-squared: 0.3153
## F-statistic: 42.61 on 3 and 268 DF, p-value: < 2.2e-16pred_val_po3 <- predict(po3_fit , newdata = Test)
mse_po3 <- mean((Test$price - pred_val_po3)^2)
mse_po3## [1] 2.722392e+12In housing data, while larger homes generally cost more, the relationship is rarely perfectly linear. However, higher-order polynomials (like degree 3 or 4) are notorious for becoming unstable at the edges of your data range (the “tails” of your area distribution). By choosing the degree 2 model, you have successfully identified the model that best balances capturing the trend without introducing unnecessary complexity.
LOOCV
library(boot)##
## Attaching package: 'boot'## The following object is masked from 'package:car':
##
## logitglm_fit <- glm(price ~ area, data = housing)
cv.error <- cv.glm(housing, glm_fit)
cv.error$delta## [1] 2.513734e+12 2.513711e+12cv.error1 <- rep(0,10)
for ( i in 1:10){
glm_fit_i <- glm(price ~ poly(area, degree = i ) , data = housing)
cv.error1[i] <- cv.glm(housing, glm_fit_i)$delta[1]
}
cv.error1## [1] 2.513734e+12 2.414752e+12 2.441341e+12 2.342341e+12 2.379141e+12
## [6] 2.459825e+12 2.416350e+12 2.727776e+12 4.053577e+12 5.473217e+13K FOLD CROSS VALIDATION
library(boot)
cv.error2 <- rep(0,10)
for ( i in 1:10){
glm_fit_2 <- glm(price ~ poly(area, degree = i ) , data = housing)
cv.error2[i] <- cv.glm(housing, glm_fit_2, K =10 )$delta[1]
}
cv.error2## [1] 2.518716e+12 2.420287e+12 2.483398e+12 2.333778e+12 2.385050e+12
## [6] 2.496823e+12 2.409573e+12 2.708670e+12 3.848721e+12 1.105725e+13Looking at both your LOOCV and 10-Fold results, Degree 4 consistently provides the lowest error.Degrees 1–3: The model is “learning” the trend, and error drops as the model gains the flexibility to follow the curve of the data.Degree 4 (The Sweet Spot): This is where the model is complex enough to capture the true underlying pattern but hasn’t yet started to “over-memorize” the random noise.Degrees 5–10: You notice the error begins to increase (especially drastically by degree 8, 9, and 10). This is a classic indicator of overfitting. The model has become so flexible that it is now “fitting the noise” rather than the trend, causing its performance on unseen data to degrade significantly.3. Comparing Validation MethodsLOOCV vs. 10-Fold: You’ll notice the values are quite similar, which gives you high confidence in the result.Stability: 10-Fold CV is generally preferred over LOOCV because it is computationally faster and often provides a more stable estimate of the test error by reducing the correlation between the training folds.The “Explosion” of Error: The massive jump in error at Degree 10 (\(5.47 \times 10^{13}\)) is a classic example of why you should never use high-degree polynomials without regularization. They become highly erratic at the boundaries of your area variable, leading to massive prediction errors for houses that are very small or very large.SummaryYour systematic testing proves that while area is a strong predictor, the relationship is non-linear. You should choose the Degree 4 model for the best predictive performance on new data.
dim(housing)## [1] 545 13sum(is.na(housing))## [1] 0housing <- na.omit(housing)
library(leaps)## Warning: package 'leaps' was built under R version 4.6.1subset_sel <- regsubsets( price ~ . , housing)
reg_sub <- summary(subset_sel)
reg_sub## Subset selection object
## Call: regsubsets.formula(price ~ ., housing)
## 13 Variables (and intercept)
## Forced in Forced out
## area FALSE FALSE
## bedrooms FALSE FALSE
## bathrooms FALSE FALSE
## stories FALSE FALSE
## mainroadyes FALSE FALSE
## guestroomyes FALSE FALSE
## basementyes FALSE FALSE
## hotwaterheatingyes FALSE FALSE
## airconditioningyes FALSE FALSE
## parking FALSE FALSE
## prefareayes FALSE FALSE
## furnishingstatussemi-furnished FALSE FALSE
## furnishingstatusunfurnished FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
## area bedrooms bathrooms stories mainroadyes guestroomyes basementyes
## 1 ( 1 ) "*" " " " " " " " " " " " "
## 2 ( 1 ) "*" " " "*" " " " " " " " "
## 3 ( 1 ) "*" " " "*" " " " " " " " "
## 4 ( 1 ) "*" " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " "*" "*" " " " " " "
## 7 ( 1 ) "*" " " "*" "*" " " " " "*"
## 8 ( 1 ) "*" " " "*" "*" " " " " "*"
## hotwaterheatingyes airconditioningyes parking prefareayes
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " "*" " " " "
## 4 ( 1 ) " " "*" " " " "
## 5 ( 1 ) " " "*" " " "*"
## 6 ( 1 ) " " "*" "*" "*"
## 7 ( 1 ) " " "*" "*" "*"
## 8 ( 1 ) " " "*" "*" "*"
## furnishingstatussemi-furnished furnishingstatusunfurnished
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " " "
## 3 ( 1 ) " " " "
## 4 ( 1 ) " " " "
## 5 ( 1 ) " " " "
## 6 ( 1 ) " " " "
## 7 ( 1 ) " " " "
## 8 ( 1 ) " " "*"subset_full <- regsubsets( price ~ . , housing, nvmax = 13)
reg_sub_full <- summary(subset_full)
reg_sub_full$adjr2## [1] 0.2859806 0.4631117 0.5421602 0.5784222 0.6113350 0.6325451 0.6475058
## [8] 0.6584091 0.6666483 0.6711105 0.6736184 0.6745276 0.6740117reg_sub_full$cp## [1] 648.34486 353.65001 222.81663 163.34417 109.63175 75.43509 51.66305
## [8] 34.65430 22.08445 15.75246 12.64309 12.15805 14.00000reg_sub_full$bic## [1] -171.9818 -322.0760 -403.5843 -443.2624 -482.2731 -507.5682 -524.9351
## [8] -536.7742 -544.7978 -546.8612 -545.7537 -541.9967 -535.8581plot(reg_sub_full$adjr2 , xlab = "Number of variables", ylab = "adj R^2", type = "b")
adj_max <- which.max(reg_sub_full$adjr2)
points( adj_max ,reg_sub_full$adjr2[adj_max], pch = 20, col = "red")
plot(reg_sub_full$bic , xlab = "Number of variables", ylab = "bic", type = "b")
bic_min <- which.min(reg_sub_full$bic)
points( bic_min ,reg_sub_full$bic[ bic_min], pch = 20, col = "red" )
plot(reg_sub_full$cp , xlab = "Number of variables", ylab = "cp", type = "b")
cp_min <- which.min(reg_sub_full$cp)
points( cp_min ,reg_sub_full$cp[ cp_min], pch = 20, col = "red" )
plot(subset_full , scale = "bic")
coef(subset_full , 10)## (Intercept) area
## 196289.1411 252.1058
## bathrooms stories
## 1034173.9440 494356.1424
## mainroadyes basementyes
## 406815.3534 459481.4981
## hotwaterheatingyes airconditioningyes
## 859533.5760 885120.3631
## parking prefareayes
## 277731.5878 667628.2774
## furnishingstatusunfurnished
## -396086.1661subset_frd <- regsubsets( price ~ . , housing, method = "forward", nvmax = 13)
reg_frd <- summary(subset_frd)
reg_frd## Subset selection object
## Call: regsubsets.formula(price ~ ., housing, method = "forward", nvmax = 13)
## 13 Variables (and intercept)
## Forced in Forced out
## area FALSE FALSE
## bedrooms FALSE FALSE
## bathrooms FALSE FALSE
## stories FALSE FALSE
## mainroadyes FALSE FALSE
## guestroomyes FALSE FALSE
## basementyes FALSE FALSE
## hotwaterheatingyes FALSE FALSE
## airconditioningyes FALSE FALSE
## parking FALSE FALSE
## prefareayes FALSE FALSE
## furnishingstatussemi-furnished FALSE FALSE
## furnishingstatusunfurnished FALSE FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: forward
## area bedrooms bathrooms stories mainroadyes guestroomyes basementyes
## 1 ( 1 ) "*" " " " " " " " " " " " "
## 2 ( 1 ) "*" " " "*" " " " " " " " "
## 3 ( 1 ) "*" " " "*" " " " " " " " "
## 4 ( 1 ) "*" " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " "*" "*" " " " " " "
## 7 ( 1 ) "*" " " "*" "*" " " " " "*"
## 8 ( 1 ) "*" " " "*" "*" " " " " "*"
## 9 ( 1 ) "*" " " "*" "*" " " " " "*"
## 10 ( 1 ) "*" " " "*" "*" "*" " " "*"
## 11 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## hotwaterheatingyes airconditioningyes parking prefareayes
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " "*" " " " "
## 4 ( 1 ) " " "*" " " " "
## 5 ( 1 ) " " "*" " " "*"
## 6 ( 1 ) " " "*" "*" "*"
## 7 ( 1 ) " " "*" "*" "*"
## 8 ( 1 ) " " "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*"
## furnishingstatussemi-furnished furnishingstatusunfurnished
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " " "
## 3 ( 1 ) " " " "
## 4 ( 1 ) " " " "
## 5 ( 1 ) " " " "
## 6 ( 1 ) " " " "
## 7 ( 1 ) " " " "
## 8 ( 1 ) " " "*"
## 9 ( 1 ) " " "*"
## 10 ( 1 ) " " "*"
## 11 ( 1 ) " " "*"
## 12 ( 1 ) " " "*"
## 13 ( 1 ) "*" "*"coef(subset_frd , 10)## (Intercept) area
## 196289.1411 252.1058
## bathrooms stories
## 1034173.9440 494356.1424
## mainroadyes basementyes
## 406815.3534 459481.4981
## hotwaterheatingyes airconditioningyes
## 859533.5760 885120.3631
## parking prefareayes
## 277731.5878 667628.2774
## furnishingstatusunfurnished
## -396086.1661subset_back <- regsubsets( price ~ . , housing, method = "backward", nvmax = 13)
reg_back <- summary(subset_back)
reg_back## Subset selection object
## Call: regsubsets.formula(price ~ ., housing, method = "backward", nvmax = 13)
## 13 Variables (and intercept)
## Forced in Forced out
## area FALSE FALSE
## bedrooms FALSE FALSE
## bathrooms FALSE FALSE
## stories FALSE FALSE
## mainroadyes FALSE FALSE
## guestroomyes FALSE FALSE
## basementyes FALSE FALSE
## hotwaterheatingyes FALSE FALSE
## airconditioningyes FALSE FALSE
## parking FALSE FALSE
## prefareayes FALSE FALSE
## furnishingstatussemi-furnished FALSE FALSE
## furnishingstatusunfurnished FALSE FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: backward
## area bedrooms bathrooms stories mainroadyes guestroomyes basementyes
## 1 ( 1 ) "*" " " " " " " " " " " " "
## 2 ( 1 ) "*" " " "*" " " " " " " " "
## 3 ( 1 ) "*" " " "*" " " " " " " " "
## 4 ( 1 ) "*" " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " "*" "*" " " " " " "
## 7 ( 1 ) "*" " " "*" "*" " " " " "*"
## 8 ( 1 ) "*" " " "*" "*" " " " " "*"
## 9 ( 1 ) "*" " " "*" "*" " " " " "*"
## 10 ( 1 ) "*" " " "*" "*" "*" " " "*"
## 11 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## hotwaterheatingyes airconditioningyes parking prefareayes
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " "*" " " " "
## 4 ( 1 ) " " "*" " " " "
## 5 ( 1 ) " " "*" " " "*"
## 6 ( 1 ) " " "*" "*" "*"
## 7 ( 1 ) " " "*" "*" "*"
## 8 ( 1 ) " " "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*"
## furnishingstatussemi-furnished furnishingstatusunfurnished
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " " "
## 3 ( 1 ) " " " "
## 4 ( 1 ) " " " "
## 5 ( 1 ) " " " "
## 6 ( 1 ) " " " "
## 7 ( 1 ) " " " "
## 8 ( 1 ) " " "*"
## 9 ( 1 ) " " "*"
## 10 ( 1 ) " " "*"
## 11 ( 1 ) " " "*"
## 12 ( 1 ) " " "*"
## 13 ( 1 ) "*" "*"coef(subset_back, 10)## (Intercept) area
## 196289.1411 252.1058
## bathrooms stories
## 1034173.9440 494356.1424
## mainroadyes basementyes
## 406815.3534 459481.4981
## hotwaterheatingyes airconditioningyes
## 859533.5760 885120.3631
## parking prefareayes
## 277731.5878 667628.2774
## furnishingstatusunfurnished
## -396086.1661validation set approach
set.seed(1)
Train_vali <- sample(c(TRUE,FALSE),nrow(housing), replace = TRUE)
Test <- (!Train_vali)
subset_sel <- regsubsets( price ~ . , housing[Train_vali, ], nvmax = 13)
reg_sub <- summary(subset_sel)
reg_sub## Subset selection object
## Call: regsubsets.formula(price ~ ., housing[Train_vali, ], nvmax = 13)
## 13 Variables (and intercept)
## Forced in Forced out
## area FALSE FALSE
## bedrooms FALSE FALSE
## bathrooms FALSE FALSE
## stories FALSE FALSE
## mainroadyes FALSE FALSE
## guestroomyes FALSE FALSE
## basementyes FALSE FALSE
## hotwaterheatingyes FALSE FALSE
## airconditioningyes FALSE FALSE
## parking FALSE FALSE
## prefareayes FALSE FALSE
## furnishingstatussemi-furnished FALSE FALSE
## furnishingstatusunfurnished FALSE FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: exhaustive
## area bedrooms bathrooms stories mainroadyes guestroomyes basementyes
## 1 ( 1 ) "*" " " " " " " " " " " " "
## 2 ( 1 ) "*" " " "*" " " " " " " " "
## 3 ( 1 ) "*" " " "*" " " " " " " " "
## 4 ( 1 ) "*" " " "*" " " " " " " " "
## 5 ( 1 ) "*" " " "*" " " " " " " " "
## 6 ( 1 ) "*" " " "*" "*" " " " " " "
## 7 ( 1 ) "*" " " "*" "*" " " " " " "
## 8 ( 1 ) "*" " " "*" "*" " " " " "*"
## 9 ( 1 ) "*" " " "*" "*" " " " " "*"
## 10 ( 1 ) "*" " " "*" "*" "*" " " "*"
## 11 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## hotwaterheatingyes airconditioningyes parking prefareayes
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " "*" " " " "
## 4 ( 1 ) " " "*" " " "*"
## 5 ( 1 ) "*" "*" " " "*"
## 6 ( 1 ) "*" "*" " " "*"
## 7 ( 1 ) "*" "*" "*" "*"
## 8 ( 1 ) "*" "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*"
## furnishingstatussemi-furnished furnishingstatusunfurnished
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " " "
## 3 ( 1 ) " " " "
## 4 ( 1 ) " " " "
## 5 ( 1 ) " " " "
## 6 ( 1 ) " " " "
## 7 ( 1 ) " " " "
## 8 ( 1 ) " " " "
## 9 ( 1 ) " " "*"
## 10 ( 1 ) " " "*"
## 11 ( 1 ) " " "*"
## 12 ( 1 ) " " "*"
## 13 ( 1 ) "*" "*"test_mat <- model.matrix( price ~ . , data = housing[Test, ])
val.error <- rep(NA,13)
for (i in 1:13) {
coefi <- coef(subset_sel, id = i)
pred <- test_mat[, names(coefi)] %*% coefi
val.error[i] <- mean((housing$price[Test] - pred)^2)
}
val.error## [1] 2.413849e+12 1.847891e+12 1.621165e+12 1.491646e+12 1.509419e+12
## [6] 1.360067e+12 1.280644e+12 1.234845e+12 1.192712e+12 1.165999e+12
## [11] 1.154399e+12 1.144285e+12 1.144419e+12which.min(val.error)## [1] 12coef(subset_sel, 12)## (Intercept) area
## -116551.1318 288.7282
## bedrooms bathrooms
## 66600.0462 1162193.3081
## stories mainroadyes
## 398279.9102 334143.9382
## guestroomyes basementyes
## 294706.3624 403769.7360
## hotwaterheatingyes airconditioningyes
## 1056406.7197 1003416.7058
## parking prefareayes
## 312384.1860 622521.3793
## furnishingstatusunfurnished
## -335493.6186predict.regsubsets <- function(object, newdata, id, ...) {
# Fix: Extract formula from the 'call' attribute
form <- as.formula(object$call[[2]])
# Create the model matrix
mat <- model.matrix(form, newdata)
# Get coefficients for the chosen subset size (id)
coefi <- coef(object, id = id)
# Get the names of the coefficients
xvar <- names(coefi)
# Perform matrix multiplication to get predictions
mat[, xvar] %*% coefi
}
set.seed(2)
k = 10
n <- nrow(housing)
folds <- sample(rep(1:k , length = n))
cv.error3 <- matrix(NA, k, 13, dimnames = list(NULL, paste(1:13)))
for (j in 1:k){
best.fit <- regsubsets( price ~ . ,housing[folds != j, ], nvmax = 13)
for (i in 1:13 ) {
pred <- predict(best.fit, housing[folds == j, ] , id = i)
cv.error3[j, i] <- mean((housing$price[folds == j] - pred)^2)
}
}
mean.cv.error3 <- apply(cv.error3,2, mean)
mean.cv.error3## 1 2 3 4 5 6
## 2.593896e+12 1.894128e+12 1.616311e+12 1.535269e+12 1.375187e+12 1.319692e+12
## 7 8 9 10 11 12
## 1.297249e+12 1.268441e+12 1.215801e+12 1.193985e+12 1.183862e+12 1.175012e+12
## 13
## 1.182018e+12which.min(mean.cv.error3)## 12
## 12