HW 3
Santina Bardengo
1. Import the dataset Apartments.xlsx
library(readxl)## Warning: package 'readxl' was built under R version 4.3.3mydata <- read_xlsx("./Apartments.xlsx")
head(mydata,10)## # A tibble: 10 × 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 1
## 7 14 20 1850 0 1
## 8 18 6 1970 1 1
## 9 22 7 2270 1 0
## 10 25 2 2570 1 0Description:
- 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
mydata$ID <- seq(1, nrow(mydata))summary(mydata)## Age Distance Price Parking Balcony ID
## Min. : 1.00 Min. : 1.00 Min. :1400 Min. :0.0000 Min. :0.0000 Min. : 1
## 1st Qu.:12.00 1st Qu.: 4.00 1st Qu.:1710 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:22
## Median :18.00 Median :12.00 Median :1950 Median :1.0000 Median :0.0000 Median :43
## Mean :18.55 Mean :14.22 Mean :2019 Mean :0.5059 Mean :0.4353 Mean :43
## 3rd Qu.:24.00 3rd Qu.:20.00 3rd Qu.:2290 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:64
## Max. :45.00 Max. :45.00 Max. :2820 Max. :1.0000 Max. :1.0000 Max. :852. Change categorical variables into factors. (0.5 p)
mydata$ParkingF <- factor(mydata$Parking,
levels = c("0", "1"),
labels = c("NO", "YES"))
mydata$BalconyF <- factor(mydata$Balcony,
levels = c("0", "1"),
labels = c("NO", "YES"))
head(mydata)## # 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)
Check the assumptions: - Variable price is numerical - Check that the variable price has a normal distribution
library(ggplot2)
ggplot(mydata, aes(x = Price)) +
geom_histogram(binwidth = 500, colour = "black") +
ylab("Frequency") +
xlab("Price")
As the graph doesn’t show significant evidence of normality I’ll do the Shapiro test.
shapiro.test(mydata$Price)##
## Shapiro-Wilk normality test
##
## data: mydata$Price
## W = 0.94017, p-value = 0.0006513We reject H0 (price is normally distributed) at p < 0.001. Then, price is not normally distributed so we need to make the non-parametric test: Wilcoxon Signed Rank Test.
- H0: me = 1950
- H1: me != 1950
wilcox.test(mydata$Price,
mu = 1950,
correct = FALSE)##
## Wilcoxon signed rank test
##
## data: mydata$Price
## V = 2048, p-value = 0.2408
## alternative hypothesis: true location is not equal to 1950We can’t accept H0 at p = 0.2408. Based on the sample data I found that the median of price of the Apartments is 1950. (p=0.2408)
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 = mydata)
summary(fit1)##
## Call:
## lm(formula = Price ~ Age, data = mydata)
##
## 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.230255B0: If the Apartment is new, it means that there is no Age, the average price is 2185.455 (p < 0.001)
B1: When we increase in one year the age of the apartment, the price will decree on average by 8.974. (p = 0.04)
Coefficient of correlation: The sign of the correlation is determine with the sign of the regression coefficient (B1) so we infer that the correlation is negative and weak.
Coefficient of determination: 5.302% of the variability of Price can be explained by 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(mydata[ , c(-4, -5, -6, -7)], #Scatterplot matrix
smooth = FALSE) 
It seems that there is no problem of multicolinearity as the graph doesn’t show significance evidence of correlation between Age and Distance.
6. Estimate the multiple regression function: Price = f(Age, Distance). Save it in object named fit2.
fit2 <- lm(Price ~ Age + Distance,
data = mydata)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata)
##
## 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. Chech the multicolinearity with VIF statistics. Explain the findings. (1 p)
vif(fit2)## Age Distance
## 1.001845 1.001845mean(vif(fit2))## [1] 1.001845Since the VIF are less than 5, there is no multicolinearity problem and none of the variables should be removed. Also the mean is close to 1.
8. Calculate standardized residuals and Cooks Distances for model fit2. Remove any potentially problematic units (outliers or units with high influence). (2 p)
mydata$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
As I can see in the histogram of standarized residuals, there is no value outside +-3 so it doesn’t exist any outlier to remove from the data.
Regarding the Cook distance which gives us information about the exitance of high impact variables, I can say that there is one unit that makes the histogram non continuous so I need to remove it. Now I have to find the units.
Now we have to find the units:
head(mydata[order(-mydata$CooksD),], 6)## # A tibble: 6 × 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.061
## 6 40 2 2400 0 1 39 NO YES 1.09 0.038Only the first unit should 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, unionmydata <- mydata %>%
filter(!ID == 38)
#I removed the unit with age = 5 and distance = 45fit2 <- lm(Price ~ Age + Distance,
data = mydata)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -604.92 -229.63 -56.49 192.97 599.35
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2456.076 73.931 33.221 < 2e-16 ***
## Age -6.464 3.159 -2.046 0.044 *
## Distance -22.955 2.786 -8.240 2.52e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 276.1 on 81 degrees of freedom
## Multiple R-squared: 0.4838, Adjusted R-squared: 0.4711
## F-statistic: 37.96 on 2 and 81 DF, p-value: 2.339e-12mydata$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
We can see again a gap in the histogram so we need to eliminate the unit after the gap again
head(mydata[order(-mydata$CooksD),], 6)## # A tibble: 6 × 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.072
## 4 37 3 2540 1 1 22 YES YES 1.47 0.056
## 5 8 26 2300 1 1 25 YES YES 1.82 0.052
## 6 8 2 2820 1 0 58 YES NO 1.70 0.039library(dplyr)
mydata <- mydata %>%
filter(!ID == 55)fit2 <- lm(Price ~ Age + Distance,
data = mydata)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -627.27 -212.96 -46.23 205.05 578.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2490.112 76.189 32.684 < 2e-16 ***
## Age -7.850 3.244 -2.420 0.0178 *
## Distance -23.945 2.826 -8.473 9.53e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.5 on 80 degrees of freedom
## Multiple R-squared: 0.4968, Adjusted R-squared: 0.4842
## F-statistic: 39.49 on 2 and 80 DF, p-value: 1.173e-12mydata$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
head(mydata[order(-mydata$CooksD),], 10)## # A tibble: 10 × 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.05
## 6 40 2 2400 0 1 39 NO YES 1.04 0.038
## 7 8 2 2820 1 0 58 YES NO 1.64 0.038
## 8 10 1 2810 0 0 57 NO NO 1.58 0.033
## 9 18 1 2800 1 0 2 YES NO 1.76 0.031
## 10 18 1 2800 1 1 61 YES YES 1.76 0.031I have to remove the unit that is after the gap because there is still with high impact
library(dplyr)
mydata <- mydata %>%
filter(!ID == 53)library(dplyr)
mydata <- mydata %>%
filter(!ID == 33)fit2 <- lm(Price ~ Age + Distance,
data = mydata)mydata$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
Finally, the histogram is continuous so we can say that there is no units of high impact in the dataset. I have to continue the analysis with 81 observations.
9. Check for potential heteroskedasticity with scatterplot between standarized residuals and standardized fitted values. Explain the findings. (1 p)
mydata$StdFitted <- scale(fit2$fitted.values)
library(car)
scatterplot(y = mydata$StdResid, x = mydata$StdFitted,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
By looking at the scatterplot I can infer that there is heteroskedasticity because the variances of the errors are different through the plot.
10. Are standardized residuals distributed normally? Show the graph and formally test it. Explain the findings. (1 p)
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(mydata$StdResid)##
## Shapiro-Wilk normality test
##
## data: mydata$StdResid
## W = 0.94297, p-value = 0.001286The histogram does not show clear evidence of the standardized residuals being normally distributed. That’s why I did the shapiro test and I conclude that I can reject the null hypothesis of this test at p = 0.001286. H0 = standardized residuals are normally distributed.
This violate one of the assumptions of the test but It is not a huge problem because the sample size is big enough.
11. Estimate the fit2 again without potentially excluded units and show the summary of the model. Explain all coefficients. (2 p)
library(readxl)
mydata <- read_xlsx("./Apartments.xlsx")
mydata$ParkingF <- factor(mydata$Parking,
levels = c("0", "1"),
labels = c("NO", "YES"))
mydata$BalconyF <- factor(mydata$Balcony,
levels = c("0", "1"),
labels = c("NO", "YES"))
mydata$ID <- seq(1, nrow(mydata))
fit2 <- lm(Price ~ Age + Distance,
data = mydata)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata)
##
## 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-11- B0: If the Apartment is new and it is in the city center, it means that there is no Age and no distance to the center, the average price is 2460.101. (p < 0.01)
- B1: When we increase in one year the age of the apartment, the price will decree on average by 7.934 assuming distance remains constant. (p = 0.02)
- B2: When we increase in one kilometer the distance to the city center, the price will decree on average by 20.667 assuming age remains constant. (p < 0.01)
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 = mydata)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: there is no significant difference between fit2 and fit3. The additional terms does not improve the model.
H1: there is significant difference between fit2 and fit3. The additional terms in fit3 improves the model.
Looking at the results, I can conclude that I reject H0 at p = 0.02. fit3 provides a significantly better fit to the data than 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 = mydata)
##
## 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 values of Age, Distance and Balcony, the price of the apartment with parking will be on average 196.168 more expensive than the apartment without parking. (p = 0.003)
B4: Given the values of Age, Distance and Parking, the price of the apartment with balcony will be on average 1.935 more expensive than the apartment without Balcony. Note: it is not necessary to explain this coefficient because there is no significant evidence that Balcony can explain Price.
F-statistic hypothesis:
H0: \(\rho^2 = 0\) .I can not explain anything with the model
H1: \(\rho^2 \neq 0\) . I found linear relationship using this model to fit the data.
In this case, I reject H0 at p < 0.001. So, the dependent variable can be explained by linear relationship of Age, Distance, Parking and Balcony.
15. Save fitted values and calculate the residual for apartment ID2. (1 p)
mydata$Fitted <- fitted.values(fit3)
residualID2 <- mydata[mydata$ID == 2, "Price"] - mydata[mydata$ID == 2, "Fitted"]
residualID2## Price
## 1 442.5889