HW 3
Milagros Rodriguez Ventura
1. Import the dataset Apartments.xlsx
library(readxl)## Warning: package 'readxl' was built under R version 4.3.3data <- read_xlsx('./Apartments.xlsx')
head(data)## # A tibble: 6 × 5
## Age Distance Price Parking Balcony
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 7 28 1640 0 1
## 2 18 1 2800 1 0
## 3 7 28 1660 0 0
## 4 28 29 1850 0 1
## 5 18 18 1640 1 1
## 6 28 12 1770 0 1Description:
- Age: Age of an apartment in years
- Distance: The distance from city center in km
- Price: Price per m2
- Parking: 0-No, 1-Yes
- Balcony: 0-No, 1-Yes
data$ID <- seq(1, nrow(data)) 2. Change categorical variables into factors. (0.5 p)
data$ParkingF <- factor(data$Parking,
levels = c("0", "1"),
labels = c("NO", "YES"))
data$BalconyF <- factor(data$Balcony,
levels = c("0", "1"),
labels = c("NO", "YES"))
head(data)## # A tibble: 6 × 8
## Age Distance Price Parking Balcony ID ParkingF BalconyF
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <fct> <fct>
## 1 7 28 1640 0 1 1 NO YES
## 2 18 1 2800 1 0 2 YES NO
## 3 7 28 1660 0 0 3 NO NO
## 4 28 29 1850 0 1 4 NO YES
## 5 18 18 1640 1 1 5 YES YES
## 6 28 12 1770 0 1 6 NO YES3. Test the hypothesis H0: Mu_Price = 1900 eur. What can you conclude? (1 p)
Assumptions:
- Variable “Price” is numeric
- Variable “Price” is normally distributed:
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3ggplot(data, aes(x = Price)) +
geom_histogram(binwidth = 500, colour = "black") +
ylab("Frequency") +
xlab("Price")
As the histogram doesn’t show evidently normal distribution of the
variable, I run a Shapiro-Wilk test
shapiro.test(data$Price)##
## Shapiro-Wilk normality test
##
## data: data$Price
## W = 0.94017, p-value = 0.0006513We reject H0 at p = 0.0006513, so the variable is not normally distributed and I need to run a non-parametric test (Wilcoxon Signed Rank Test):
median(data$Price)## [1] 1950- H0: me = 1950
- H1: me != 1950
wilcox.test(data$Price,
mu = 1950,
correct = FALSE)##
## Wilcoxon signed rank test
##
## data: data$Price
## V = 2048, p-value = 0.2408
## alternative hypothesis: true location is not equal to 1950library(effectsize)## Warning: package 'effectsize' was built under R version 4.3.3effectsize(wilcox.test(data$Price,
mu = 1950,
correct = FALSE))## r (rank biserial) | 95% CI
## ---------------------------------
## 0.15 | [-0.10, 0.38]
##
## - Deviation from a difference of 1950.Conclusion: we can’t reject H0 at p = 0.2408, so we conclude that the median price of the apartments is 1950
4. Estimate the simple regression function: Price = f(Age). Save results in object fit1 and explain the estimate of regression coefficient, coefficient of correlation and coefficient of determination. (2 p)
fit1 <- lm(Price ~ Age,
data = data)
summary(fit1)##
## Call:
## lm(formula = Price ~ Age, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -623.9 -278.0 -69.8 243.5 776.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2185.455 87.043 25.108 <2e-16 ***
## Age -8.975 4.164 -2.156 0.034 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 369.9 on 83 degrees of freedom
## Multiple R-squared: 0.05302, Adjusted R-squared: 0.04161
## F-statistic: 4.647 on 1 and 83 DF, p-value: 0.03401cor_coef <- sqrt(summary(fit1)$r.squared)
cor_coef## [1] 0.230255Regression coefficients
- B0: for new departments (age 0), the average price is USD 2185.455 millions (p < 0.001)
- B1: for each additional year that the apartment has, the price per m2 of it decreases, on average, by USD 8.975 millions (p = 0.034)
Coefficient of correlation
There’s a weak correlation between Age and Price. This correlation is negative, which means that when the age increases, the price of the apartment decreases
Coefficient of determination
We can say that 5.3% of the variability of Price is explained by the linear effect of age
5. Show the scateerplot matrix between Price, Age and Distance. Based on the matrix determine if there is potential problem with multicolinearity. (0.5 p)
library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3scatterplotMatrix(data[ , c(-4, -5, -6, -7)],
smooth = FALSE) 
It seems that there are not multicollinearity problems because it looks like there is no correlation between the variables Distance and Age based on the scatterplot matrix
6. Estimate the multiple regression function: Price = f(Age, Distance). Save it in object named fit2.
fit2 <- lm(Price ~ Age + Distance,
data = data)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -603.23 -219.94 -85.68 211.31 689.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2460.101 76.632 32.10 < 2e-16 ***
## Age -7.934 3.225 -2.46 0.016 *
## Distance -20.667 2.748 -7.52 6.18e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 286.3 on 82 degrees of freedom
## Multiple R-squared: 0.4396, Adjusted R-squared: 0.4259
## F-statistic: 32.16 on 2 and 82 DF, p-value: 4.896e-117. Check the multicolinearity with VIF statistics. Explain the findings. (1 p)
vif(fit2) ## Age Distance
## 1.001845 1.001845mean(vif(fit2))## [1] 1.001845As the VIF of the variables is less than 1 and the mean VIF is close to 1, none of the variables should be removed, since the is no multicollinearity
8. Calculate standardized residuals and Cooks Distances for model fit2. Remove any potentially problematic units (outliers or units with high influence). (2 p)
data$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
data$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(data$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(data$StdResid)##
## Shapiro-Wilk normality test
##
## data: data$StdResid
## W = 0.95303, p-value = 0.003645hist(data$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
By looking at the histogram of the standardized residuals we can
conclude that there are no outliers, because there are no values outside
the range (-3, +3)
As regards the histogram of Cooks distances, that gives information about the variables with high influence, we can see that there’s a gap and we need to remove the units that follow this gap
We find the units that make the histogram non-continuous:
head(data[order(-data$CooksD),], 5) ## # A tibble: 5 × 10
## Age Distance Price Parking Balcony ID ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <fct> <fct> <dbl> <dbl>
## 1 5 45 2180 1 1 38 YES YES 2.58 0.32
## 2 43 37 1740 0 0 55 NO NO 1.44 0.104
## 3 2 11 2790 1 0 33 YES NO 2.05 0.069
## 4 7 2 1760 0 1 53 NO YES -2.15 0.066
## 5 37 3 2540 1 1 22 YES YES 1.58 0.061We can see that the first unit is the one that has to be removed
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondata <- data %>%
filter(!ID == 38) We fit the model again:
fit2 <- lm(Price ~ Age + Distance,
data = data)data$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
data$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(data$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(data$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
We can still see a gap in the Cooks distances histogram so we need to
remove the observation that follows that gap
head(data[order(-data$CooksD),], 3) ## # A tibble: 3 × 10
## Age Distance Price Parking Balcony ID ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <fct> <fct> <dbl> <dbl>
## 1 43 37 1740 0 0 55 NO NO 1.60 0.129
## 2 2 11 2790 1 0 33 YES NO 2.22 0.084
## 3 7 2 1760 0 1 53 NO YES -2.24 0.072We remove the first observation
library(dplyr)
data <- data %>%
filter(!ID == 55)And, now, we fit the model again and look if there are still gaps
fit2 <- lm(Price ~ Age + Distance,
data = data)data$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
data$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(data$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(data$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
We can still see gaps in the histograms, so we must remove the units
that follow the gaps
head(data[order(-data$CooksD),], 5) ## # A tibble: 5 × 10
## Age Distance Price Parking Balcony ID ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <fct> <fct> <dbl> <dbl>
## 1 2 11 2790 1 0 33 YES NO 2.17 0.084
## 2 7 2 1760 0 1 53 NO YES -2.35 0.084
## 3 37 3 2540 1 1 22 YES YES 1.57 0.065
## 4 8 26 2300 1 1 25 YES YES 1.85 0.054
## 5 45 21 1910 0 1 31 NO YES 1.07 0.05The first 2 observations are the ones that we must remove so that we can have a continuous histogram
library(dplyr)
data <- data %>%
filter(!ID == 33)
data <- data %>%
filter(!ID == 53) Now, I fit the model again to check if there are no more gaps
fit2 <- lm(Price ~ Age + Distance,
data = data)data$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
data$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(data$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(data$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
Finally, the resulting histogram is continuous, which means that there are no units of high impact. The dataset now contains 81 observations
9. Check for potential heteroskedasticity with scatterplot between standarized residuals and standrdized fitted values. Explain the findings. (1 p)
data$StdFitted <- scale(fit2$fitted.values)
library(car)
scatterplot(y = data$StdResid, x = data$StdFitted,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
By looking at the scatterplot, I can conclude that there’s
heteroskedasticity because the variances of the errors are not constant.
For smaller values of x the variance is lower than for bigger values of
x
10. Are standardized residuals ditributed normally? Show the graph and formally test it. Explain the findings. (1 p)
hist(data$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(data$StdResid)##
## Shapiro-Wilk normality test
##
## data: data$StdResid
## W = 0.94297, p-value = 0.001286By looking and the histogram and also by looking at the results of the Shapiro-Wilk test we can conclude that the standardized residuals are not normally distributed (we reject H0 at p = 0.001)
This violates one of the assumptions of the classical linear regression model, but it is not problematic because the sample size is large enough
11. Estimate the fit2 again without potentially excluded units and show the summary of the model. Explain all coefficients. (2 p)
library(readxl)
data <- read_xlsx('./Apartments.xlsx')
data$ParkingF <- factor(data$Parking,
levels = c("0", "1"),
labels = c("NO", "YES"))
data$BalconyF <- factor(data$Balcony,
levels = c("0", "1"),
labels = c("NO", "YES"))
data$ID <- seq(1, nrow(data)) fit2 <- lm(Price ~ Age + Distance,
data = data)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -603.23 -219.94 -85.68 211.31 689.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2460.101 76.632 32.10 < 2e-16 ***
## Age -7.934 3.225 -2.46 0.016 *
## Distance -20.667 2.748 -7.52 6.18e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 286.3 on 82 degrees of freedom
## Multiple R-squared: 0.4396, Adjusted R-squared: 0.4259
## F-statistic: 32.16 on 2 and 82 DF, p-value: 4.896e-11B0: the price of a new apartment (age 0) and at 0 km from the city center is, on average, USD 2460.101 millions (p < 0.001)
B1: for each additional year that the apartment has, the average price decreases by USD 7.934 millions (assuming that the distance remains constant) (p = 0.016)
B2: for each additional km further away from the city center, the average price of the apartments decreases by USD 20.667 millions (assuming that the age of the apartment remains constant) (p < 0.001)
12. Estimate the linear regression function Price = f(Age, Distance, Parking and Balcony). Be careful to correctly include categorical variables. Save the object named fit3.
fit3 <- lm(Price ~ Age + Distance + ParkingF + BalconyF,
data = data)13. With function anova check if model fit3 fits data better than model fit2. (1 p)
anova(fit2, fit3)## Analysis of Variance Table
##
## Model 1: Price ~ Age + Distance
## Model 2: Price ~ Age + Distance + ParkingF + BalconyF
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 82 6720983
## 2 80 5991088 2 729894 4.8732 0.01007 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1H0: the more complex model doesn’t provide a significantly better fit to the data than the simpler model (no significant difference between fit2 and fit3)
H1: the more complex model provides a significantly better fit to the data than the simpler model (significant difference between fit2 and fit3)
So, looking at the anova test, we reject H0 at p = 0.01, which means that the more complex model (fit 3), fits the data significantly better than the simpler model (fit2)
14. Show the results of fit3 and explain regression coefficient for both categorical variables. Can you write down the hypothesis which is being tested with F-statistics, shown at the bottom of the output? (2 p)
summary(fit3)##
## Call:
## lm(formula = Price ~ Age + Distance + ParkingF + BalconyF, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -459.92 -200.66 -57.48 260.08 594.37
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2301.667 94.271 24.415 < 2e-16 ***
## Age -6.799 3.110 -2.186 0.03172 *
## Distance -18.045 2.758 -6.543 5.28e-09 ***
## ParkingFYES 196.168 62.868 3.120 0.00251 **
## BalconyFYES 1.935 60.014 0.032 0.97436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.7 on 80 degrees of freedom
## Multiple R-squared: 0.5004, Adjusted R-squared: 0.4754
## F-statistic: 20.03 on 4 and 80 DF, p-value: 1.849e-11B3: given the age, distance and the fact that the apartment has or not a balcony, if it has parking the apartment is, on average, USD 196.168 millions more expensive compared to the apartments that don’t have parking (p = 0.003)
B4: given the age, distance and whether the apartment has or not parking, if it has a balcony, the price of the apartment will be, on average, USD 1.935 millions more expensive than the apartments with no balcony (Note: there is no evidence that this coefficient is significant, so it wouldn’t be necessary to explain it)
F-statistic hypotheses:
- H0: \(\rho^2 = 0\) -> we can’t explain anything with the model
- H1: \(\rho^2 \neq 0\) -> we found linear relationships between the explanatory variables and the dependent variable
In this case, as p < 0.001, we reject H0, which means that the price can be explained with a linear relationship between the variables Age, Distance, BalconyF and ParkingF
15. Save fitted values and calculate the residual for apartment ID2. (1 p)
data$Fitted <- fitted.values(fit3)
residualID2 <- data[data$ID == 2, 'Price'] - data[data$ID == 2, 'Fitted']
residualID2## Price
## 1 442.5889