Task 1
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphadata(diamonds)mydata <- print(diamonds)## # A tibble: 53,940 × 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
## 6 0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48
## 7 0.24 Very Good I VVS1 62.3 57 336 3.95 3.98 2.47
## 8 0.26 Very Good H SI1 61.9 55 337 4.07 4.11 2.53
## 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 2.49
## 10 0.23 Very Good H VS1 59.4 61 338 4 4.05 2.39
## # … with 53,930 more rowsdim(mydata)## [1] 53940 10My data set includes 53940 observations with 10 variables.
colnames(mydata) <- c("Weight_in_carats", "Cut_Quality", "Color", "Clarity","Total_depth_percentage","Table", "Price", "Length", "Width", "Depth")
head(mydata)## # A tibble: 6 × 10
## Weight_in_carats Cut_Quality Color Clarity Total_depth_percentage Table Price Length Width Depth
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
## 6 0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48Description:
- ID: ID of a diamond
- Weight_in_carats: weight of the diamond 0.2-5.01 (1 carat = 200mg)
- Cut_Quality: quality of the cut Fair, Good, Very Good, Premium, Ideal
- Color: diamond color J (worst) to D (best)
- Clarity: measurement of how clear the diamond is I1 (worst), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (best)
- Total_depth_percentage:total depth percentage 43-79 (excludes the crown)
- Table: width of top of diamond (crown of the diamond) relative to widest point 43-95
- Price: price in US dollars $326-$18,823
- Length:length in mm 0-10.74
- Width:width in mm 0-58.9 (the farthest distance)
- Depth:depth in mm 0-31.8 (from the crown to the base)
mydata$ID <- seq.int(nrow(mydata))Adding new variable:
mydata$Cut_QualityF <- factor(mydata$Cut_Quality,
levels = c( "Fair", "Good", "Very Good", "Premium", "Ideal"),
labels = c(1,2,3,4,5))mydata <- mydata[c(11,1,2,12,3,4,5,6,7,8,9,10)]Table including just Ideal cut of the diamond:
JustIdealCut <- mydata[mydata$Cut_Quality == "Ideal", ]
head(JustIdealCut)## # A tibble: 6 × 12
## ID Weight_in_carats Cut_Quality Cut_QualityF Color Clarity Total_depth_percentage Table Price Length Width Depth
## <int> <dbl> <ord> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 0.23 Ideal 5 E SI2 61.5 55 326 3.95 3.98 2.43
## 2 12 0.23 Ideal 5 J VS1 62.8 56 340 3.93 3.9 2.46
## 3 14 0.31 Ideal 5 J SI2 62.2 54 344 4.35 4.37 2.71
## 4 17 0.3 Ideal 5 I SI2 62 54 348 4.31 4.34 2.68
## 5 40 0.33 Ideal 5 I SI2 61.8 55 403 4.49 4.51 2.78
## 6 41 0.33 Ideal 5 I SI2 61.2 56 403 4.49 4.5 2.75nrow(mydata[mydata$Cut_Quality == "Ideal", ])/ nrow(mydata)## [1] 0.3995365Ideal cut represents 39.9% of all observed diamonds in this data set.
library(psych)
describeBy(mydata$Price, mydata$Cut_Quality)##
## Descriptive statistics by group
## group: Fair
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1610 4358.76 3560.39 3282 3695.65 2183.13 337 18574 18237 1.78 3.07 88.73
## ------------------------------------------------------------------------------------------
## group: Good
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 4906 3928.86 3681.59 3050.5 3251.51 2853.26 327 18788 18461 1.72 3.04 52.56
## ------------------------------------------------------------------------------------------
## group: Very Good
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 12082 3981.76 3935.86 2648 3243.22 2855.49 336 18818 18482 1.6 2.24 35.81
## ------------------------------------------------------------------------------------------
## group: Premium
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 13791 4584.26 4349.2 3185 3822.23 3371.43 326 18823 18497 1.33 1.07 37.03
## ------------------------------------------------------------------------------------------
## group: Ideal
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 21551 3457.54 3808.4 1810 2656.14 1630.86 326 18806 18480 1.84 2.98 25.94Price described based on Cut Quality of diamonds.
Arithmetic mean of the price for Ideal cut is 3457.54 $ Half of the Ideal cut diamonds are valued at 3457.54$ or less, the other 50% are valued higher than that.
Max price for Ideal cut is 18806 $.
Premium cut has a standard deviation of 4349.2$. Compared to all other cuts, this one has the highest standard deviation, which means that prices are more spread out from the mean. The lowest standard deviation is in group Fair cuts.
sum(mydata$Price < 500)## [1] 17291729 diamonds cost less than 500€.
sum(mydata$Price >= 15000)## [1] 16561656 diamonds cost equal or more to 15000$
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:pastecs':
##
## first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(naniar)
mydata <- mydata %>%
replace_with_na(replace = list(Length = c(0.000),
Width = c(0.000),
Depth = c(0.000)))I noticed that some data had 0 value so I changed the number to na.
mydata <- tidyr::drop_na(mydata)I deleted all na values which helped in next step (summary) due to min being always 0.
summary(mydata[c(-1,-4)])## Weight_in_carats Cut_Quality Color Clarity Total_depth_percentage Table Price
## Min. :0.2000 Fair : 1609 D: 6774 SI1 :13063 Min. :43.00 Min. :43.00 Min. : 326
## 1st Qu.:0.4000 Good : 4902 E: 9797 VS2 :12254 1st Qu.:61.00 1st Qu.:56.00 1st Qu.: 949
## Median :0.7000 Very Good:12081 F: 9538 SI2 : 9185 Median :61.80 Median :57.00 Median : 2401
## Mean :0.7977 Premium :13780 G:11284 VS1 : 8170 Mean :61.75 Mean :57.46 Mean : 3931
## 3rd Qu.:1.0400 Ideal :21548 H: 8298 VVS2 : 5066 3rd Qu.:62.50 3rd Qu.:59.00 3rd Qu.: 5323
## Max. :5.0100 I: 5421 VVS1 : 3654 Max. :79.00 Max. :95.00 Max. :18823
## J: 2808 (Other): 2528
## Length Width Depth
## Min. : 3.730 Min. : 3.680 Min. : 1.07
## 1st Qu.: 4.710 1st Qu.: 4.720 1st Qu.: 2.91
## Median : 5.700 Median : 5.710 Median : 3.53
## Mean : 5.732 Mean : 5.735 Mean : 3.54
## 3rd Qu.: 6.540 3rd Qu.: 6.540 3rd Qu.: 4.04
## Max. :10.740 Max. :58.900 Max. :31.80
## sapply(mydata[c(2,9,10,11,12 )], FUN = median)## Weight_in_carats Price Length Width Depth
## 0.70 2401.00 5.70 5.71 3.53Half of the diamonds weigh up to and including 0.7 carats, others more.
str(mydata)## tibble [53,920 × 12] (S3: tbl_df/tbl/data.frame)
## $ ID : int [1:53920] 1 2 3 4 5 6 7 8 9 10 ...
## $ Weight_in_carats : num [1:53920] 0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
## $ Cut_Quality : Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
## $ Cut_QualityF : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 5 4 2 4 2 3 3 3 1 3 ...
## $ Color : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 2 6 7 7 6 5 2 5 ...
## $ Clarity : Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 2 3 5 4 2 6 7 3 4 5 ...
## $ Total_depth_percentage: num [1:53920] 61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
## $ Table : num [1:53920] 55 61 65 58 58 57 57 55 61 61 ...
## $ Price : int [1:53920] 326 326 327 334 335 336 336 337 337 338 ...
## $ Length : num [1:53920] 3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
## $ Width : num [1:53920] 3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
## $ Depth : num [1:53920] 2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...library(pastecs)
round(stat.desc(mydata[ ,c(2,8,9,10,11,12)]), 2)## Weight_in_carats Table Price Length Width Depth
## nbr.val 53920.00 53920.00 53920.00 53920.00 53920.00 53920.00
## nbr.null 0.00 0.00 0.00 0.00 0.00 0.00
## nbr.na 0.00 0.00 0.00 0.00 0.00 0.00
## min 0.20 43.00 326.00 3.73 3.68 1.07
## max 5.01 95.00 18823.00 10.74 58.90 31.80
## range 4.81 52.00 18497.00 7.01 55.22 30.73
## sum 43011.89 3098072.50 211959155.00 309049.32 309225.11 190879.30
## median 0.70 57.00 2401.00 5.70 5.71 3.53
## mean 0.80 57.46 3930.99 5.73 5.73 3.54
## SE.mean 0.00 0.01 17.17 0.00 0.00 0.00
## CI.mean.0.95 0.00 0.02 33.66 0.01 0.01 0.01
## var 0.22 4.99 15898405.35 1.25 1.30 0.49
## std.dev 0.47 2.23 3987.28 1.12 1.14 0.70
## coef.var 0.59 0.04 1.01 0.20 0.20 0.20hist(mydata$Price,
main = "Distribution of the variable Price",
ylab = "Frequency",
xlab = "Diamond price",
breaks = seq(0,20000,200))
This is positive asymmetric distribution, with a high concentration of observations below the U$5,000 mark
library(ggplot2)
ggplot(mydata, aes(x = Price)) +
geom_histogram(binwidth = 1000, colour="black", fill="gray") +
facet_wrap(~Cut_Quality, ncol = 5) +
ylab("Frequency")+
theme_bw()+
theme(axis.text.x = element_text(angle = 30))
This is a histogram of a diamond price distributed by cut. As we can see, among all cuts, the Ideal cut has the highest frequency for a specific price group.
ggplot(data=mydata, aes(x=Cut_Quality, y=Price)) +
geom_boxplot()+
xlab("Cut Quality") +
ylab("Price")
Another figure of relationship between price and Cut quality- Boxplot. We can see that no matter the cut quality mean of the price is similar. The biggest mean is in group Premium cut at 4584.26 and the lowest Ideal at 3457.54. Interestingly, the highest priced diamond is Premium cut valued at 18823. Another interesting thing is the ideal cut being the “cheapest” based on average price, which means all ideal cuts are either small diamonds or the color is not the clearest.
subset(mydata, Price == max(Price))## # A tibble: 1 × 12
## ID Weight_in_carats Cut_Quality Cut_QualityF Color Clarity Total_depth_percentage Table Price Length Width Depth
## <int> <dbl> <ord> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 27750 2.29 Premium 4 I VS2 60.8 60 18823 8.5 8.47 5.16The cut for the highest priced diamond is Premium.
ggplot(data=mydata, aes(x=Weight_in_carats, y=Price, color=Color)) +
geom_point()+
xlab("Weight in carats") +
ylab("Price")
The scatterplot shows Diamonds based on their color, weight in carats and price. Interestingly, diamonds with the purest color on average weight less than the ones with the lowest color score.
ggplot(mydata) +
geom_histogram(mapping = aes(x = Length), binwidth = 0.01)
Histogram shows different lenghts of diamonds. High concentration of diamonds is placed below the 9mm mark.Graph is multimodal.
mydata$RowWidth <- seq.int(nrow(mydata))head(mydata[order(-mydata$Width), ])## # A tibble: 6 × 13
## ID Weight_in_carats Cut_Quality Cut_QualityF Color Clarity Total_depth_per…¹ Table Price Length Width Depth RowWi…²
## <int> <dbl> <ord> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int>
## 1 24068 2 Premium 4 H SI2 58.9 57 12210 8.09 58.9 8.06 24059
## 2 49190 0.51 Ideal 5 E VS1 61.8 55 2075 5.15 31.8 5.12 49173
## 3 27416 5.01 Fair 1 J I1 65.5 59 18018 10.7 10.5 6.98 27402
## 4 27631 4.5 Fair 1 J I1 65.8 58 18531 10.2 10.2 6.72 27615
## 5 25999 4.01 Premium 4 I I1 61 61 15223 10.1 10.1 6.17 25988
## 6 26000 4.01 Premium 4 J I1 62.5 62 15223 10.0 9.94 6.24 25989
## # … with abbreviated variable names ¹Total_depth_percentage, ²RowWidthmydata <- mydata[c(-24059,-49173), ]ggplot(mydata) +
geom_histogram(mapping = aes(x = Width), binwidth = 0.01)
Histogram shows different widths of diamonds. High concentration of diamonds is placed below the 9mm mark. Graph is multimodal.
mydata$RowDepth <- seq.int(nrow(mydata))head(mydata[order(-mydata$Depth), ])## # A tibble: 6 × 14
## ID Weight_in_carats Cut_Quality Cut_QualityF Color Clarity Total_d…¹ Table Price Length Width Depth RowWi…² RowDe…³
## <int> <dbl> <ord> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int> <int>
## 1 48411 0.51 Very Good 3 E VS1 61.8 54.7 1970 5.12 5.15 31.8 48394 48393
## 2 27416 5.01 Fair 1 J I1 65.5 59 18018 10.7 10.5 6.98 27402 27401
## 3 27631 4.5 Fair 1 J I1 65.8 58 18531 10.2 10.2 6.72 27615 27614
## 4 27131 4.13 Fair 1 H I1 64.8 61 17329 10 9.85 6.43 27117 27116
## 5 23645 3.65 Fair 1 H I1 67.1 53 11668 9.53 9.48 6.38 23636 23636
## 6 26445 4 Very Good 3 I I1 63.3 58 15984 10.0 9.94 6.31 26432 26431
## # … with abbreviated variable names ¹Total_depth_percentage, ²RowWidth, ³RowDepthmydata <- mydata[c(-48392), ]ggplot(mydata) +
geom_histogram(mapping = aes(x = Depth), binwidth = 0.01)
Histogram shows different depths of diamonds. High concentration of diamonds is placed between 2,2 mm and 5,5 mm mark.Graph is multimodal.
TASK 2
mydata2 <- read.table("C:/Users/Eneja/Desktop/take home exam/Task 2/Body mass.csv", header=TRUE, sep=";", dec=",")summary(mydata2$Mass)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 49.70 60.23 62.80 62.88 64.50 83.20Arithmetic mean of the mass is 62.88. Half of the ninthgraders at the beginning of the 2021/2022 school year have body weight of 62.88 or less, the other 50% weight more than that.
Maximum body weight of the ninthgraders at the beginning of the 2021/2022 school year is 83.20. Minimum body weight of the ninthgraders at the beginning of the 2021/2022 school year is 49.70.
library(pastecs)
round(stat.desc(mydata2[ , -1]), 2)## nbr.val nbr.null nbr.na min max range sum median mean
## 50.00 0.00 0.00 49.70 83.20 33.50 3143.80 62.80 62.88
## SE.mean CI.mean.0.95 var std.dev coef.var
## 0.85 1.71 36.14 6.01 0.10Difference between maximum and minimum body weight of the ninthgraders at the beginning of the 2021/2022 school year is 33.50.
Body weigh has a standard deviation of 6.01.
library(psych)
describe(mydata2 [, -1])## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 50 62.88 6.01 62.8 62.56 3.34 49.7 83.2 33.5 0.85 2.11 0.85hist(mydata2$Mass,
main = "Distribution of the variable Mass",
ylab = "Frequency",
xlab = "Mass",
breaks = seq(40,90,5),
right = FALSE)
Histogram is unimodal, with a high concentration of body weight of the ninthgraders at the beginning of the 2021/2022 school year below 65.
library(ggplot2)
ggplot(NULL, aes(c(-4, 4)))+
geom_line(stat = "function", fun = dt, args = list (df = 49))+
xlab("Density") +
ylab("Sample estimates")+
labs(title= "Distribution of sample estimates")
qt(p = 0.025, df = 49, lower.tail = FALSE)## [1] 2.009575qt(p = 0.025, df = 49, lower.tail = TRUE)## [1] -2.009575t.test(mydata2$Mass,
mu = 59.5,
alternative = "two.sided")##
## One Sample t-test
##
## data: mydata2$Mass
## t = 3.9711, df = 49, p-value = 0.000234
## alternative hypothesis: true mean is not equal to 59.5
## 95 percent confidence interval:
## 61.16758 64.58442
## sample estimates:
## mean of x
## 62.876It is extremely unlikely that the average body weight of ninthgraders at the beginning of the 2021/2022 school year is the same as it was before online schooling in 2019/2019 school year. According to data the body weight of ninthgraders at the beginning of the 2021/2022 school year is larger then it was. We reject Ho at p = 0.0002
x1=3.9711^2
x2=3.9711^2 + 49
sqrt(x1/x2)## [1] 0.4934295Difference between the body weight of ninthgraders in 2019/2019 school year and at the beginning of the 2021/2022 school year is large.
TASK 3
#install.packages("car")
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logit#install.packages("ggplot2")
library(ggplot2)
#install.packages("reshape2")
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithsImport the dataset Apartments.xlsx
library(readxl)
mydata3 <- read_xlsx("C:/Users/Eneja/Desktop/take home exam/Task 3/Apartments.xlsx")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
Change categorical variables into factors
mydata3$ParkingF <- factor(mydata3$Parking,
levels = c(0,1),
labels = c("No", "Yes"))
mydata3$BalconyF <- factor(mydata3$Balcony,
levels = c(0,1),
labels = c("No", "Yes"))Test the hypothesis H0: Mu_Price = 1900 eur. What can you conclude?
t.test(mydata3$Price,
mu = 1900,
alternative = "two.sided")##
## One Sample t-test
##
## data: mydata3$Price
## t = 2.9022, df = 84, p-value = 0.004731
## alternative hypothesis: true mean is not equal to 1900
## 95 percent confidence interval:
## 1937.443 2100.440
## sample estimates:
## mean of x
## 2018.941It is extremely unlikely that the average price per m2 is equal to 1900. According to data the average price per m2 is larger. We reject Ho at p = 0.0002
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.
fit1 <- lm(Price ~ Age,
data = mydata3)
summary(fit1)##
## Call:
## lm(formula = Price ~ Age, data = mydata3)
##
## 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.03401Estimate of regression coefficient: If age is increased by 1 % point then price on average decreases by 8.975 % point.
Coefficient of determination: 5.3% of variablity of price of the appartemnt is explained by linear effect Age.
sqrt(summary(fit1)$r.squared)## [1] 0.230255Coefficient of correlation shows that relationship between price and age is weak.
Show the scateerplot matrix between Price, Age and Distance. Based on the matrix determine if there is potential problem with multicolinearity.
library(car)
scatterplotMatrix(mydata3[ ,c(3,1,2)],
smooth = FALSE)
The slope between age and distance is almost 0, so I assume we don’t have a problem with multicolinearity. Later this will be checked with VIF statistics.
Estimate the multiple regression function: Price = f(Age, Distance). Save it in object named fit2.
fit2 <- lm(Price ~ Age + Distance,
data = mydata3)Chech the multicolinearity with VIF statistics. Explain the findings.
vif(fit2)## Age Distance
## 1.001845 1.001845We don’t have strong connection between them since the values are below 5. Therefore we leave them in a model.
Calculate standardized residuals and Cooks Distances for model fit2. Remove any potentially problematic case (outlier or unit with big influence).
mydata3$StdResid <- round(rstandard(fit2), 3)
mydata3$CooksD <- round(cooks.distance(fit2), 3)hist(mydata3$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
mydata3$RowStdRes <- seq.int(nrow(mydata3))head(mydata3[order(mydata3$StdResid),], 3)## # A tibble: 3 × 10
## Age Distance Price Parking Balcony ParkingF BalconyF StdResid CooksD RowStdRes
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <dbl> <dbl> <int>
## 1 7 2 1760 0 1 No Yes -2.15 0.066 53
## 2 12 14 1650 0 1 No Yes -1.50 0.013 13
## 3 12 14 1650 0 0 No No -1.50 0.013 72mydata3 <- mydata3[c(-53), ]hist(mydata3$CooksD,
xlab = "Cooks Distance",
ylab = "Frequency",
main = "Histogram of Cooks distance")
mydata3$RowcooksD <- seq.int(nrow(mydata3))head(mydata3[order(-mydata3$CooksD),], 6)## # A tibble: 6 × 11
## Age Distance Price Parking Balcony ParkingF BalconyF StdResid CooksD RowStdRes RowcooksD
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <dbl> <dbl> <int> <int>
## 1 5 45 2180 1 1 Yes Yes 2.58 0.32 38 38
## 2 43 37 1740 0 0 No No 1.44 0.104 55 54
## 3 2 11 2790 1 0 Yes No 2.05 0.069 33 33
## 4 37 3 2540 1 1 Yes Yes 1.58 0.061 22 22
## 5 40 2 2400 0 1 No Yes 1.09 0.038 39 39
## 6 8 2 2820 1 0 Yes No 1.66 0.037 58 57mydata3 <- mydata3[c(-38,-54, -33, -22), ]fit2 <- lm(Price ~ Age + Distance,
data = mydata3)hist(mydata3$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
hist(mydata3$CooksD,
xlab = "Cooks Distance",
ylab = "Frequency",
main = "Histogram of Cooks distance")
mydata3$StdFittedValues <- scale(fit2$fitted.values)Check for potential heteroskedasticity with scatterplot between standarized residuals and standrdized fitted values. Explain the findings.
library(car)
scatterplot(y = mydata3$StdResid, x = mydata3$StdFittedValues,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
It is linear and homoskedastic.
Are standardized residuals ditributed normally? Show the graph and formally test it. Explain the findings.
hist(mydata3$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(mydata3$StdResid)##
## Shapiro-Wilk normality test
##
## data: mydata3$StdResid
## W = 0.93418, p-value = 0.0004761Ho: variable is normally distributed.
Based on p value, std.residuals are not normally distributed.
Estimate the fit2 again without potentially excluded cases and show the summary of the model. Explain all coefficients.
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -411.50 -203.69 -45.24 191.11 492.56
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2502.467 75.024 33.356 < 2e-16 ***
## Age -8.674 3.221 -2.693 0.00869 **
## Distance -24.063 2.692 -8.939 1.57e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 256.8 on 77 degrees of freedom
## Multiple R-squared: 0.5361, Adjusted R-squared: 0.524
## F-statistic: 44.49 on 2 and 77 DF, p-value: 1.437e-13The expected price for an apartment excluding the age and distance is 2502.467 m2. If average age of an apartment is increased by 1 year, apartment price decreases by 8.674 % point (p< 0.01) assuming that everything else is unchanged. If distance is increased by 1 % point then price on average decreases by 24.063 % point (p< 0.0001) assuming that everything else is unchanged.
53.6% of price of an apartment is explained by linear effect of Age and Distance
sqrt(summary(fit2)$r.squared)## [1] 0.732187Linear relationship between price of an apartment and both explanatory variables is strong.
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 = mydata3)
summary(fit3)##
## Call:
## lm(formula = Price ~ Age + Distance + ParkingF + BalconyF, data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -390.93 -198.19 -53.64 186.73 518.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2393.316 93.930 25.480 < 2e-16 ***
## Age -7.970 3.191 -2.498 0.0147 *
## Distance -21.961 2.830 -7.762 3.39e-11 ***
## ParkingFYes 128.700 60.801 2.117 0.0376 *
## BalconyFYes 6.032 57.307 0.105 0.9165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 252.7 on 75 degrees of freedom
## Multiple R-squared: 0.5623, Adjusted R-squared: 0.5389
## F-statistic: 24.08 on 4 and 75 DF, p-value: 7.764e-13With function anova check if model fit3 fits data better than model fit2.
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 77 5077362
## 2 75 4791128 2 286234 2.2403 0.1135H0 = fit2 is more appropriate H1 = fit 3 is more appropriate
according to p value, we cannot reject H0, therefore fit 2 fits better.
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?
summary(fit3)##
## Call:
## lm(formula = Price ~ Age + Distance + ParkingF + BalconyF, data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -390.93 -198.19 -53.64 186.73 518.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2393.316 93.930 25.480 < 2e-16 ***
## Age -7.970 3.191 -2.498 0.0147 *
## Distance -21.961 2.830 -7.762 3.39e-11 ***
## ParkingFYes 128.700 60.801 2.117 0.0376 *
## BalconyFYes 6.032 57.307 0.105 0.9165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 252.7 on 75 degrees of freedom
## Multiple R-squared: 0.5623, Adjusted R-squared: 0.5389
## F-statistic: 24.08 on 4 and 75 DF, p-value: 7.764e-13Given the age and distance, apartments with parking have on average price higher by 128.7. We could not confirm that the apartment with a balcony affects on price. p=0.9165
Save fitted values and claculate the residual for apartment ID2.
mydata3$FittedValues <- fitted.values(fit3)
residuals (fit3)[2]## 2
## 443.4026