HW 3
Stephanie Yang
1. Import the dataset Apartments.xlsx
library(readxl)
mydata <- read_xlsx("./Apartments.xlsx")
#View(mydata)Description:
- 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
2. 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"))
#View(mydata)3. Test the hypothesis H0: Mu_Price = 1900 eur. What can you conclude? (1 p)
#Check for Normality from a graph
library(ggplot2)
ggplot(mydata, aes(x = Price)) +
geom_histogram(binwidth = 100, colour = "black") +
ylab("Frequency") +
xlab("Price")
The graph doesn’t appear to be normal. We must conduct a Shapiro-Wilk test for normality.
#Shapiro-Wilk test for normality
shapiro.test(mydata$Price)##
## Shapiro-Wilk normality test
##
## data: mydata$Price
## W = 0.94017, p-value = 0.0006513Hypothesis for Shapiro-Wilk Test:
- H0: variable is normally distributed
- H1: variable is not normally distributed
Based off the results of the Shapiro-Wilk test, we reject the null hypothesis meaning the variable is not normally distributed at p < 0.001. This means we must conduct a Wilcoxon Signed Rank Test
median(mydata$Price)## [1] 1950wilcox.test(mydata$Price,
mu = 1900,
correct = FALSE)##
## Wilcoxon signed rank test
##
## data: mydata$Price
## V = 2328, p-value = 0.02828
## alternative hypothesis: true location is not equal to 1900library(effectsize)
effectsize(wilcox.test(mydata$Price,
mu = 1900,
correct = FALSE))## r (rank biserial) | 95% CI
## --------------------------------
## 0.27 | [0.04, 0.48]
##
## - Deviation from a difference of 1900.interpret_rank_biserial(0.27, rules = "funder2019")## [1] "medium"
## (Rules: funder2019)Interpretation: Based on the sample data, we found that the median price per m2 of apartments in the Vienna area is 1950 and is different from the hypothesis of 1900 per m2 (𝑝 < 0.03, 𝑟 = 0.27 – medium effect).
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_coef1 <- sqrt(summary(fit1)$r.squared)
cor_coef1## [1] 0.230255Regression Coefficient
If age of an apartment is increased by 1 year, on average, the price per m2 of the apartment will decrease by 8.98 (p < 0.04).
Coefficient of Correlation
The Pearson Correlation Coefficient is 0.23 which means there is a positive, weak linear relationship between age of an apartment and the price per m2.
Coefficient of Determination
The Coefficient of Determination is .05302 which means that 5.302% of variability of price per m2 can be explained by linear effect of age of the apartment.
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)## Loading required package: carDatascatterplotMatrix(mydata[ , -c(4,5,6,7)], smooth=FALSE)
Based off the scatterplot matrix, there does not seem to be a problem with multicolinearity. The scatterplots are not perfectly correlated and they do not form perfect lines which means there is not multicolinearity.
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.001845The vif for both age and distance are under 5. The average vif is also about 1 meaning that there is no colinearity between price, age, and distance.
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")
shapiro.test(mydata$StdResid)##
## Shapiro-Wilk normality test
##
## data: mydata$StdResid
## W = 0.95303, p-value = 0.003645hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
head(mydata[order(mydata$StdResid),], 3) #Three units with lowest value of stand. residuals## # A tibble: 3 × 9
## Age Distance Price Parking Balcony ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 7 2 1760 0 1 No Yes -2.15 0.066
## 2 12 14 1650 0 1 No Yes -1.50 0.013
## 3 12 14 1650 0 0 No No -1.50 0.013head(mydata[order(-mydata$StdResid),], 3) #Three units with highest value of stand. residuals## # A tibble: 3 × 9
## Age Distance Price Parking Balcony ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 5 45 2180 1 1 Yes Yes 2.58 0.32
## 2 2 11 2790 1 0 Yes No 2.05 0.069
## 3 18 1 2800 1 0 Yes No 1.78 0.03head(mydata[order(-mydata$CooksD),], 6) #Six units with highest value of Cooks distance## # A tibble: 6 × 9
## Age Distance Price Parking Balcony ParkingF BalconyF StdResid CooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 5 45 2180 1 1 Yes Yes 2.58 0.32
## 2 43 37 1740 0 0 No No 1.44 0.104
## 3 2 11 2790 1 0 Yes No 2.05 0.069
## 4 7 2 1760 0 1 No Yes -2.15 0.066
## 5 37 3 2540 1 1 Yes Yes 1.58 0.061
## 6 40 2 2400 0 1 No Yes 1.09 0.038which.max(mydata$CooksD)## 38
## 38cleanedData <- mydata[-38, ]It looks like there are no outliers since there are no points with a standardized residule below -3 or above 3. The results of the Shapiro-Wilk test show that the standardized residuals are not normally distributed, (reject H0 at p < 0.004) but we have a larger sample of 84 units so we can move on. However, there was a gap in the histogram of Hooks distances so I found which point that was using the which.max function and removed that row from the data.
#Finding the new fit2
fit2_1 <- lm(Price ~ Age + Distance,
data = cleanedData)
summary(fit2_1)##
## Call:
## lm(formula = Price ~ Age + Distance, data = cleanedData)
##
## 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-12#Calculating adjusted std. residuals and hooks distances
cleanedData$StdResid2 <- round(rstandard(fit2_1), 3)
cleanedData$CooksD2 <- round(cooks.distance(fit2_1), 3)
hist(cleanedData$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(cleanedData$StdResid2)##
## Shapiro-Wilk normality test
##
## data: cleanedData$StdResid2
## W = 0.95649, p-value = 0.006355hist(cleanedData$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
9. Check for potential heteroskedasticity with scatterplot between standarized residuals and standrdized fitted values. Explain the findings. (1 p)
cleanedData$StdFitted <- scale(fit2_1$fitted.values)
library(car)
scatterplot(y = cleanedData$StdResid2, x = cleanedData$StdFitted,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riversols_test_breusch_pagan(fit2_1)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## ---------------------------------
## Response : Price
## Variables: fitted values of Price
##
## Test Summary
## -----------------------------
## DF = 1
## Chi2 = 2.927455
## Prob > Chi2 = 0.08708469Looking at the scatterplot, it looks like there might be slight heteroskedasticity because as the standardized fitted values increase, the standardized residuals start to become more spread out but only by a little bit. There doesn’t seem to be a lot of heteroskedasticity.
10. Are standardized residuals ditributed normally? Show the graph and formally test it. Explain the findings. (1 p)
hist(cleanedData$StdResid2,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(cleanedData$StdResid2)##
## Shapiro-Wilk normality test
##
## data: cleanedData$StdResid2
## W = 0.95649, p-value = 0.006355Interpretation
Looking at the scatterplot, it seems like the standardized residuals are not normally distributed. To formally test for normality, I conducted a Shapiro-Wilk test with the following hypothesis:
- H0: Distribution is normally distributed
- H1: Distribution is not normally distributed
The results of the Shapiro-Wilk test show that we reject the null hypothesis at p < .007 meaning the residuals are not normally distributed. This is a violation of the assumptions of linear regression; however, there is a large sample size (n=84) so it is ok.
11. Estimate the fit2 again without potentially excluded units and show the summary of the model. Explain all coefficients. (2 p)
summary(fit2_1)##
## Call:
## lm(formula = Price ~ Age + Distance, data = cleanedData)
##
## 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-12Explanation of coefficients:
- Age: Holding all else constant, if age of an apartment is increased by 1 year, on average, the price per m2 of the apartment will decrease by 6.46 (p < 0.05).
- Distance: Holding all else constant, if the distance of an apartment is increased by 1 km, on average, the price per m2 of the apartment will decrease by 22.96 (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.
contrasts(cleanedData$ParkingF)## Yes
## No 0
## Yes 1contrasts(cleanedData$BalconyF)## Yes
## No 0
## Yes 1fit3 <- lm(Price ~ Age + Distance + ParkingF + BalconyF,
data = cleanedData)
summary(fit3)##
## Call:
## lm(formula = Price ~ Age + Distance + ParkingF + BalconyF, data = cleanedData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -473.21 -192.37 -28.89 204.17 558.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2329.724 93.066 25.033 < 2e-16 ***
## Age -5.821 3.074 -1.894 0.06190 .
## Distance -20.279 2.886 -7.026 6.66e-10 ***
## ParkingFYes 167.531 62.864 2.665 0.00933 **
## BalconyFYes -15.207 59.201 -0.257 0.79795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.5 on 79 degrees of freedom
## Multiple R-squared: 0.5275, Adjusted R-squared: 0.5035
## F-statistic: 22.04 on 4 and 79 DF, p-value: 3.018e-1213. With function anova check if model fit3 fits data better than model fit2. (1 p)
anova(fit2_1, 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 81 6176767
## 2 79 5654480 2 522287 3.6485 0.03051 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Results
The anova function is testing the following hypothesizes:
- H0: Δrho^2 = 0
- H1: Δrho^2 > 0
We reject the null hypothesis at p < 0.04 meaning that the second model (fit3) is better.
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 = cleanedData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -473.21 -192.37 -28.89 204.17 558.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2329.724 93.066 25.033 < 2e-16 ***
## Age -5.821 3.074 -1.894 0.06190 .
## Distance -20.279 2.886 -7.026 6.66e-10 ***
## ParkingFYes 167.531 62.864 2.665 0.00933 **
## BalconyFYes -15.207 59.201 -0.257 0.79795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.5 on 79 degrees of freedom
## Multiple R-squared: 0.5275, Adjusted R-squared: 0.5035
## F-statistic: 22.04 on 4 and 79 DF, p-value: 3.018e-12Interpretation of regression coefficients
- Parking: Given all else is constant, the price per m2 of an apartment with parking, on average, is 167.53 Euros higher than apartments without parking (p < 0.01).
- Balcony: Given all else is constant, the price per m2 of an apartment with a balcony, on average, is 15.21 Euros lower than apartments without a balcony (p < 0.8).
F-statistics hypothesis
- H0 : 𝛽1 = 𝛽2 = ⋯ = 𝛽5 = 0
- H1 : at least one 𝛽𝑗 is different from 0
15. Save fitted values and calculate the residual for apartment ID2. (1 p)
cleanedData$StdFitted2 <- scale(fit3$fitted.values)
scatterplot(y = cleanedData$StdResid2, x = cleanedData$StdFitted2,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
cleanedData[2,13]## # A tibble: 1 × 1
## StdFitted2[,1]
## <dbl>
## 1 1.29fitted_val <- 2370.91
actual_val <- 2329.724 - 5.821 * 18 - 20.279 * 1 + 167.531 * 1 - 15.207 * 0
residual_val <- actual_val - fitted_val
residual_val## [1] 1.288