Mivšek HW final
2025-09-27
================================================================
Matej Mivšek
Student ID: 19607754
Statistics take home exam
================================================================
Task 1
data(mtcars) # load the dataset
mydata <- mtcars # copy it into your variable
head(mydata) # first 6 rows## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 11. - Explain the data set
Data set has 32 observations and 11 variables. Data set was acquired from R’s own “mtcars” data set.
Population: All automobiles that were available on the US market in 1973–74 (or more broadly, all cars that could have been studied from that era).
Unit: An individual automobile model (e.g., Mazda RX4, Ford Pantera L, etc.).
Sample: The 32 specific automobile models included in the dataset. This is a non-random sample, chosen by the magazine for testing and it is not necessarily representative of the entire population of cars.
Explanation of the variables:
# Create table
mtcars_info <- data.frame(
Variable = c("mpg","cyl","disp","hp","drat","wt","qsec","vs","am","gear","carb"),
Unit = c("Miles per US gallon","Number of cylinders","Displacement (cu.in.)",
"Gross horsepower","Rear axle ratio","Weight (1000 lbs)","1/4 mile time",
"Engine type","Transmission","Number of forward gears","Number of carburetors"),
Numeric_Categorical = c("Ratio","Ordinal/Discrete","Ratio","Ratio","Ratio","Ratio","Ratio",
"Binary","Binary","Ordinal","Ordinal"))
# View table
mtcars_info## Variable Unit Numeric_Categorical
## 1 mpg Miles per US gallon Ratio
## 2 cyl Number of cylinders Ordinal/Discrete
## 3 disp Displacement (cu.in.) Ratio
## 4 hp Gross horsepower Ratio
## 5 drat Rear axle ratio Ratio
## 6 wt Weight (1000 lbs) Ratio
## 7 qsec 1/4 mile time Ratio
## 8 vs Engine type Binary
## 9 am Transmission Binary
## 10 gear Number of forward gears Ordinal
## 11 carb Number of carburetors OrdinalThe variables are explained in the table format above for better readability (I used function data.frame to create a table,I had to enter each value of the table seperatly, making sure that the order of the data is right)
2. - Data Manipulation
mydata$TransmissionFactor <- factor(mydata$am,
levels = c(0, 1),
labels = c("Automatic", "Manual"))
head(mydata, 3)## mpg cyl disp hp drat wt qsec vs am gear carb TransmissionFactor
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 Manual
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 Manual
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 ManualHere I created a new factor variable, so that it is easier to tell, what transmission each car has (automatic or manual). I converted numeric into a categorical variable (factor)
mydata$Efficiency <- mydata$mpg / mydata$hp
mydata$Car <- rownames(mtcars)
mydata[1:3, c("Car", "Efficiency")]## Car Efficiency
## Mazda RX4 Mazda RX4 0.1909091
## Mazda RX4 Wag Mazda RX4 Wag 0.1909091
## Datsun 710 Datsun 710 0.2451613Here I created a completely new variable (ratio) called “Efficiency ratio”.
It shows us miles per gallon per unit of horsepower. The idea is to see, how fuel-efficient the car is relative to engine power.
I also created a new variable “Car”, so I could display just car name with the newly created Efficiency ration in the table
3. - Descriptive statistics
# Map old variable names to readable names
readable_names <- c("Miles per Gallon", "Horsepower", "Weight (1000 lbs)")
original_vars <- c("mpg", "hp", "wt")
# Create descriptive statistics table
descriptive_stats <- data.frame(
Variable = readable_names,
Mean = round(sapply(mydata[, original_vars], mean), 3),
Median = round(sapply(mydata[, original_vars], median), 3),
SD = round(sapply(mydata[, original_vars], sd), 3),
Q1 = round(sapply(mydata[, original_vars], function(x) quantile(x, 0.25)), 3),
Q3 = round(sapply(mydata[, original_vars], function(x) quantile(x, 0.75)), 3),
Min = round(sapply(mydata[, original_vars], min), 3),
Max = round(sapply(mydata[, original_vars], max), 3))
descriptive_stats## Variable Mean Median SD Q1 Q3 Min Max
## mpg Miles per Gallon 20.091 19.200 6.027 15.425 22.80 10.400 33.900
## hp Horsepower 146.688 123.000 68.563 96.500 180.00 52.000 335.000
## wt Weight (1000 lbs) 3.217 3.325 0.978 2.581 3.61 1.513 5.424First I created “original_vars” and “readable_names” to make the table more readable.
This code creates a descriptive statistics table for selected variables in mydata. It uses sapply() to apply functions like mean(), median(), sd(), min(), max(), and quantile() to each variable, and round(…, 3) to round results to three decimals. Finally, data.frame() organizes all the statistics into a neat table with readable variable names.
Analysis of descriptive statistics
I choose to analyse three different variables. Here are the explanations for selected values:
Mean (20.091 for Miles per Gallon) – The average fuel efficiency across all cars in the sample is 20.091 mpg
Median (123.000 for Horsepower) – The median horsepower is 123. This means that half of the cars have engine power below (or equal) to123 hp and half have above 123 hp.
Standard Deviation (65.563 for Horsepower) – The standard deviation of 65.563 shows that horsepower values are quite spread out around the mean. A higher SD indicates greater variability in engine power across the sample of cars.
Q1 (2.581) and Q3 (3.610 for Weight) – The 1st quartile (2.581) means that 25% of cars weigh less (or equal) than 2.581 (1000 lbs), and the 3rd quartile (3.610) means that 75% of cars weigh less than 3.610 (1000 lbs).
Minimum (10.400) and Maximum (33.900 for Miles per Gallon) – The lightest and heaviest fuel efficiency values in the sample are 10.4 and 33.9 mpg
4. - Graph
# Histogram for Miles per Gallon
hist(mydata$mpg,
main = "Distribution of Miles per Gallon",
xlab = "Miles per Gallon",
col = "skyblue",
breaks = 10)
# Histogram for Horsepower
hist(mydata$hp,
main = "Distribution of Horsepower",
xlab = "Horsepower",
col = "lightgreen",
breaks = 10)
# Histogram for Weight
hist(mydata$wt,
main = "Distribution of Car Weight",
xlab = "Weight (1000 lbs)",
col = "lightcoral",
breaks = 10)
Miles per Gallon (mpg): The distribution of mpg is right-skewed, with most cars achieving between 15 and 22 mpg.
Horsepower (hp): Horsepower is left-skewed, with the majority of cars concentrated between 75 and 175 hp.
Weight (wt): Car weight is unimodal and slightly right-skewed, with most vehicles weighing between 2.5 and 3.5 thousand pounds.
# Scatterplot: Horsepower vs. Weight
plot(mydata$wt, mydata$hp,
main = "Horsepower vs. Weight",
xlab = "Weight (1000 lbs)",
ylab = "Horsepower",
pch = 19, col = "darkblue")
# Fit linear regression
model_hp <- lm(hp ~ wt, data = mydata)
# Add regression line
abline(model_hp, col = "red", lwd = 2)
# Scatterplot: Miles per Gallon vs. Weight with regression line
plot(mydata$wt, mydata$mpg,
main = "MPG vs. Weight",
xlab = "Weight (1000 lbs)",
ylab = "Miles per Gallon",
pch = 19, col = "darkred")
# Fit linear regression
model <- lm(mpg ~ wt, data = mydata)
# Add regression line
abline(model, col = "blue", lwd = 2)
Horsepower vs. Weight: Heavier cars show a moderate positive correlation with horsepower, since larger engines are needed to power heavier vehicles.
MPG vs. Weight: There is a strong negative correlation, as heavier cars consume more fuel and thus achieve lower miles per gallon.
================================================================
Task 2
library(readxl)
mydata2 <- read_xlsx("./Business School.xlsx")
mydata2 <- as.data.frame(mydata2)
head(mydata2)## Student ID Undergrad Degree Undergrad Grade MBA Grade Work Experience Employability (Before)
## 1 1 Business 68.4 90.2 No 252
## 2 2 Computer Science 70.2 68.7 Yes 101
## 3 3 Finance 76.4 83.3 No 401
## 4 4 Business 82.6 88.7 No 287
## 5 5 Finance 76.9 75.4 No 275
## 6 6 Computer Science 83.3 82.1 No 254
## Employability (After) Status Annual Salary
## 1 276 Placed 111000
## 2 119 Placed 107000
## 3 462 Placed 109000
## 4 342 Placed 148000
## 5 347 Placed 255500
## 6 313 Placed 103500Context
Population: The population of interest is all MBA students at this business school
Unit: The unit of observation is an individual student.
Variables:
variables <- data.frame(
Variable = c("Student ID", "Undergrad Degree", "Undergrad Grade", "MBA Grade",
"Work Experience", "Employability (Before)", "Employability (After)",
"Status", "Annual Salary"),
Description = c("Identifier of each student",
"Undergraduate field of study (Business, Finance, etc.)",
"Undergraduate grade (%)",
"MBA grade (%)",
"Whether the student had prior work experience (Yes/No)",
"Employability score before MBA",
"Employability score after MBA",
"Placement outcome (Placed / Not placed)",
"Salary after MBA (currency units)"),
Scale = c("Nominal",
"Nominal",
"Ratio (continuous)",
"Ratio (continuous)",
"Nominal (binary)",
"Ratio (continuous)",
"Ratio (continuous)",
"Nominal (binary)",
"Ratio (continuous)")
)
variables## Variable Description Scale
## 1 Student ID Identifier of each student Nominal
## 2 Undergrad Degree Undergraduate field of study (Business, Finance, etc.) Nominal
## 3 Undergrad Grade Undergraduate grade (%) Ratio (continuous)
## 4 MBA Grade MBA grade (%) Ratio (continuous)
## 5 Work Experience Whether the student had prior work experience (Yes/No) Nominal (binary)
## 6 Employability (Before) Employability score before MBA Ratio (continuous)
## 7 Employability (After) Employability score after MBA Ratio (continuous)
## 8 Status Placement outcome (Placed / Not placed) Nominal (binary)
## 9 Annual Salary Salary after MBA (currency units) Ratio (continuous)As in the Task 1, I presented the variables in a table format, so it is more easily readable.
1. - Graph the distribution of undergrad degrees
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphaggplot(mydata2, aes(x = `Undergrad Degree`)) +
geom_bar(fill = "steelblue") +
labs(title = "Distribution of Undergraduate Degrees",
x = "Undergraduate Degree",
y = "Count")
Which degree is most common?
Business (n=35 students) is the most common undergraduate degree.
2. - Show descriptive statistics
# Descriptive statistics
summary(mydata2$`Annual Salary`)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 20000 87125 103500 109058 124000 340000# Histogram of Annual Salary
library(ggplot2)
ggplot(mydata2, aes(x = `Annual Salary`)) +
geom_histogram(binwidth = 20000, fill = "steelblue", color = "white", alpha = 0.8) +
labs(title = "Distribution of Annual Salary",
x = "Annual Salary",
y = "Count") +
theme_minimal(base_size = 13)
The average annual salary of the students is about 109,000.
The median annual salary is slightly lower, at 103,500, indicating a right-skewed distribution.
Most students earn between 87,000 (Q1) and 124,000 (Q3).
The lowest reported salary is 20,000, while the highest reaches 340,000.
The large difference between the mean and maximum suggests the presence of a few outliers with very high salaries.
3. - Testing hypothesis
t.test(mydata2$`MBA Grade`, mu = 74)##
## One Sample t-test
##
## data: mydata2$`MBA Grade`
## t = 2.6587, df = 99, p-value = 0.00915
## alternative hypothesis: true mean is not equal to 74
## 95 percent confidence interval:
## 74.51764 77.56346
## sample estimates:
## mean of x
## 76.04055H₀: μ(MBA Grade) = 74
H₁: μ(MBA Grade) ≠ 74
We can reject the H0 at p = 0.009, concluding that the mean MBA grade (M = 76.04) is significantly higher than 74. The effect size (Cohen’s d ≈ 0.27) is small, indicating the difference is statistically significant but practically modest.
================================================================
Task 3
================================
Import the dataset Apartments.xlsx
================================
library(readxl)
mydata3 <- read_xlsx("./ApartmentsV2.xlsx")
mydata3 <- as.data.frame(mydata3)
head(mydata3)## Age Distance Price Parking Balcony
## 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
================================
Change categorical variables into factors.
================================
# Convert categorical variables into factors
mydata3$ParkingFactor <- factor(mydata3$Parking,
levels = c(0, 1),
labels = c("No", "Yes"))
mydata3$BalconyFactor <- factor(mydata3$Balcony,
levels = c(0, 1),
labels = c("No", "Yes"))
# Check structure
str(mydata3)## 'data.frame': 85 obs. of 7 variables:
## $ Age : num 7 18 7 28 18 28 14 18 22 25 ...
## $ Distance : num 28 1 28 29 18 12 20 6 7 2 ...
## $ Price : num 1640 2800 1660 1850 1640 1770 1850 1970 2270 2570 ...
## $ Parking : num 0 1 0 0 1 0 0 1 1 1 ...
## $ Balcony : num 1 0 0 1 1 1 1 1 0 0 ...
## $ ParkingFactor: Factor w/ 2 levels "No","Yes": 1 2 1 1 2 1 1 2 2 2 ...
## $ BalconyFactor: Factor w/ 2 levels "No","Yes": 2 1 1 2 2 2 2 2 1 1 ...Here I changed the categorical variables into numerical, so we can use them in the regression model.
================================
Test the hypothesis H0: Mu_Price = 1900 eur. What can you conclude?
================================
t.test(mydata3$Price, mu = 1900)##
## 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.941h0: μ(Price) = 1900 € h1: μ(Price) ≠ 1900 €
Since p < 0.05, we reject H₀ at p=0,00431. The average apartment price is significantly different from 1900 EUR. The 95% confidence interval suggests that the true mean price is between about 1937 and 2100 EUR, therefore higher than 1900 EUR.
================================
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.
================================
# Fit simple linear regression model: Price ~ Age
fit1 <- lm(Price ~ Age, data = mydata3)
# View regression summary
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.03401# Compute correlation between Price and Age
cor(mydata3$Price, mydata3$Age)## [1] -0.230255Regression coefficient: The slope is –8.98, which means that if age increases by 1 year, the apartment price decreases on average by about 8,975 EUR. The coefficient is statistically significant at p = 0,034.
Coefficient of correlation (r = –0.23): There is a weak negative linear relationship between age and price.
Coefficient of determination (R² = 0.053): Only about 5.3% of the variation in price is explained by the age of the apartment. The explanatory power of the model is low.
================================
Show the scaterplot matrix between Price, Age and Distance. Based on the matrix determine if there is potential problem with multicolinearity.
================================
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:psych':
##
## logit# Clean up column names (removes spaces if any)
names(mydata3) <- trimws(names(mydata3))
# Now plot scatterplot matrix
scatterplotMatrix(
mydata3[, c("Price", "Age", "Distance")],
main = "Scatterplot Matrix: Price, Age and Distance",
smooth = FALSE)
The scatteprlot matrix shows us a (strong) linear negative relationship between Price and Age. The relationship between Price and Age is less linear and not as strong, but it is still negative,
There is also no clear relation between Age and Distance. Thus, there is no potential multicollinearity problem. To further support this statement, I will run another analysis!.
# Load the required package
library(Hmisc)##
## Attaching package: 'Hmisc'## The following object is masked from 'package:psych':
##
## describe## The following objects are masked from 'package:base':
##
## format.pval, units# Run rcorr on numeric matrix
rcorr(as.matrix(mydata3[ , c(-4,-5,-6,-7)]))## Age Distance Price
## Age 1.00 0.04 -0.23
## Distance 0.04 1.00 -0.63
## Price -0.23 -0.63 1.00
##
## n= 85
##
##
## P
## Age Distance Price
## Age 0.6966 0.0340
## Distance 0.6966 0.0000
## Price 0.0340 0.0000The correlation analysis shows no multicollinearity among predictors, as all correlations are well below |0.9|. Distance has a semi-strong negative correlation with Price (r = -0.63). Age shows a weak negative correlation (r = -0.23). Overall, Distance is the most important predictors, and multicollinearity is not a concern.
================================
Estimate the multiple regression function: Price = f(Age, Distance). Save it in object named fit2.
================================
# Multiple regression: Price ~ Age + Distance
fit2 <- lm(Price ~ Age + Distance, data = mydata3)
# Show summary
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata3)
##
## 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-11Explanation of the output:
The intercept (2460 EUR) is the expected price when age and distance are both zero.
If age increases by 1 year, the price decreases on average by 7,9 EUR, holding distance constant (ceteris paribus).
If distance increases by 1 km, the price decreases on average by 20,7 EUR, holding age constant (ceteris paribus)
Both age (p = 0,016) and distance (p < 0,001) are statistically significant predictors.
The model explains about 44% of the variation in prices (R² = 0,44).
================================
Chech the multicolinearity with VIF statistics. Explain the findings.
================================
library(car)
# Variance Inflation Factors for fit2
vif(fit2)## Age Distance
## 1.001845 1.001845# Calculating the average VIF
mean(vif(fit2))## [1] 1.001845The VIF values for both Age and Distance are about 1,00, which is the minimum possible value. This means that the explanatory variables are completely independent of each other, and there is no multicollinearity problem in the regression model. Therefore also average VIF is close to one, which is a desirable outcome.
================================
Calculate standardized residuals and Cooks Distances for model fit2. Remove any potentially problematic units (outliers or units with high influence).
================================
mydata3$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata3$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
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.95303, p-value = 0.003645Upon visual inspection, we can see that residuals are slightly asymmetrically distributed to the right. Since all the values are between -3 and +3, we have no problem with outliers. This would be the problem for datasets with less then 30 units.
To also formally test the distribution, we need to use the Shapiro-Wilks normality test. The p-value is less then 0,05, therefore we can reject the null hypothesis that the errors are normally distributed. Therefore, it means that the errors in the population are most likely not normally distributed, as is required by the assumption.
hist(mydata3$CooksD <- round(cooks.distance(fit2), 3),
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
To study Cook’s distance, we first display the histogram of Cook’s distances. We observe that there are some units with the large gap inbetween them. They are most probably the units with high impact, so we will remove them.
head(mydata3[order(-mydata3$CooksD),], 6) #Six units with highest value of Cooks distance## Age Distance Price Parking Balcony ParkingFactor BalconyFactor StdResid CooksD
## 38 5 45 2180 1 1 Yes Yes 2.577 0.320
## 55 43 37 1740 0 0 No No 1.445 0.104
## 33 2 11 2790 1 0 Yes No 2.051 0.069
## 53 7 2 1760 0 1 No Yes -2.152 0.066
## 22 37 3 2540 1 1 Yes Yes 1.576 0.061
## 39 40 2 2400 0 1 No Yes 1.091 0.038To identify this units, I showed the top six units with the highest Cooks distance. The unit with the highest Cook’s distance is the unit 38 (the only one with price 2180!). In addition to that, unit 55 (the only one with unique price of 1740), also has a slighlty higer value, so I decided to remove this one as well.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:Hmisc':
##
## src, summarize## The following object is masked from 'package:car':
##
## recode## 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, unionmydata3 <- mydata3 %>%
filter(!Price %in% c(2180, 1740))After identifying the units, I removed them using the inverted filter function (!). Again, because this is the only unit with price=2180, I could filter by that variable. I did the same for unit with price 1740. Again, I could only do this, because I checked that these were the only units with that price, so I used the price as a unique identifying value.
Now I need to recheck Cook’s distance an reevaluate the regression model.
# Reevaluatin the regression model
fit2 <- lm(Price ~ Age + Distance, data = mydata3)
# Rechecking residuals and Cook's distance
mydata3$StdResid <- round(rstandard(fit2), 3) #Standardized residuals
mydata3$CooksD <- round(cooks.distance(fit2), 3) #Cooks distances
head(mydata3[order(mydata3$StdResid),], 3) #Three units with lowest value of stand. residuals## Age Distance Price Parking Balcony ParkingFactor BalconyFactor StdResid CooksD
## 52 7 2 1760 0 1 No Yes -2.346 0.084
## 13 12 14 1650 0 1 No Yes -1.515 0.014
## 70 12 14 1650 0 0 No No -1.515 0.014head(mydata3[order(-mydata3$CooksD),], 6) #Six units with highest value of Cooks distance## Age Distance Price Parking Balcony ParkingFactor BalconyFactor StdResid CooksD
## 33 2 11 2790 1 0 Yes No 2.173 0.084
## 52 7 2 1760 0 1 No Yes -2.346 0.084
## 22 37 3 2540 1 1 Yes Yes 1.566 0.065
## 25 8 26 2300 1 1 Yes Yes 1.853 0.054
## 31 45 21 1910 0 1 No Yes 1.073 0.050
## 38 40 2 2400 0 1 No Yes 1.044 0.038hist(mydata3$CooksD <- round(cooks.distance(fit2), 3),
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
After running the Cook’s distance histogram again, there were now more units coming out as units with large impact. I decided to remove two more units (units with Cook’s Distance of 0,084 - after second calcualtion)
library(dplyr)
mydata3 <- mydata3 %>%
filter(!Price %in% c(2790, 1760))After removing the two more units, I will check Cook’s distances again.
fit2 <- lm(Price ~ Age + Distance, data = mydata3) # refit after any filtering
cd <- round(cooks.distance(fit2), 3)
mydata3$CooksD <- cd # lengths now match
hist(mydata3$CooksD, xlab="Cook's distance", ylab="Frequency",
main="Histogram of Cook's distances")
Now there is finally no more breaks in the histogram, so I can proceede with the exercise.
================================
Check for potential heteroskedasticity with scatterplot between standarized residuals and standrdized fitted values. Explain the findings.
================================
# Reevaluatin the regression model
fit2 <- lm(Price ~ Age + Distance, data = mydata3)
mydata3$StdFitted <- scale(fit2$fitted.values)
library(car)
scatterplot(y = mydata3$StdResid, x = mydata3$StdFitted,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
# Load the olsrr package
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## rivers# Perform Breusch-Pagan test for heteroskedasticity
ols_test_breusch_pagan(fit2)##
## 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 = 1.64762
## Prob > Chi2 = 0.1992832Because I removed another unit with strong impact, I first refitted the model.
H0: The error variance is constant (homoskedasticity). H1: The error variance is not constant (heteroskedasticity).
Result: The Breusch–Pagan test yields χ² = 1.65, p = 0.199 > 0.05, so we cannot to reject H₀ and conclude there is no evidence of heteroskedasticity. I will thus assume homoskedasticity.
Graphical check: The residual–fitted plot shows random scatter without a clear funnel shape, supporting the conclusion that variance appears roughly constant.
================================
Are standardized residuals ditributed normally? Show the graph and formally test it. Explain the findings.
================================
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")
shapiro.test(mydata3$StdResid)##
## Shapiro-Wilk normality test
##
## data: mydata3$StdResid
## W = 0.94297, p-value = 0.001286H0: The standardized residuals are normally distributed. H1: The standardized residuals are not normally distributed.
The histogram of standardized residuals looks roughly symmetric but not perfectly bell-shaped. The formal Shapiro–Wilk test confirms this: W=0.943, p=0.0013. Since the p-value is below 0.0, we reject the null hypothesis of normality. This means the standardized residuals are not normally distributed**, even though the histogram might visually suggest approximate symmetry.
================================
Estimate the fit2 again without potentially excluded units and show the summary of the model. Explain all coefficients.
================================
fit2 <- lm(Price ~ Age + Distance, data = mydata3)
summary(fit2)##
## Call:
## lm(formula = Price ~ Age + Distance, data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -408.72 -217.93 -38.86 220.60 506.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2492.715 75.544 32.997 < 2e-16 ***
## Age -7.496 3.169 -2.365 0.0205 *
## Distance -24.574 2.700 -9.100 6.9e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 259.4 on 78 degrees of freedom
## Multiple R-squared: 0.5324, Adjusted R-squared: 0.5204
## F-statistic: 44.4 on 2 and 78 DF, p-value: 1.335e-13Explanation of the coeficients:
For every additional year of Age, the apartment price decreases on average by 7.5 EUR, holding Distance constant (ceteris paribus).
For every additional unit of Distance, the apartment price decreases on average by 24.6 EUR, holding Age constant (ceteris paribus).
The intercept (2492.7 EUR) is the expected price when Age and Distance are both zero.
Both predictors are statistically significant (p<0.05), meaning Age and Distance have a real effect on Price.
The model explains about 53% of the variation in prices (R2=0.53), which is a moderate fit.
The residual standard error (259.4) indicates the typical deviation of observed prices from the predicted values.
================================
Estimate the linear regression function Price = f(Age, Distance, Parking and Balcony). Be careful to correctly include categorical variables. Save the object named fit3.
================================
# Fit multiple linear regression model with factors
fit3 <- lm(Price ~ Age + Distance + ParkingFactor + BalconyFactor, data = mydata3)
# Show model summary
summary(fit3)##
## Call:
## lm(formula = Price ~ Age + Distance + ParkingFactor + BalconyFactor,
## data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -403.16 -204.08 -41.24 239.52 528.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2372.693 93.344 25.419 < 2e-16 ***
## Age -6.906 3.118 -2.215 0.0298 *
## Distance -22.254 2.839 -7.837 2.25e-11 ***
## ParkingFactorYes 137.479 60.854 2.259 0.0267 *
## BalconyFactorYes 17.234 57.099 0.302 0.7636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 254.3 on 76 degrees of freedom
## Multiple R-squared: 0.5621, Adjusted R-squared: 0.539
## F-statistic: 24.38 on 4 and 76 DF, p-value: 5.29e-13I was careful to include the Factor variables and not the original categorical variabvles.
================================
With function anova check if model fit3 fits data better than model fit2.
================================
# Compare fit2 and fit3 using ANOVA
anova(fit2, fit3)## Analysis of Variance Table
##
## Model 1: Price ~ Age + Distance
## Model 2: Price ~ Age + Distance + ParkingFactor + BalconyFactor
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 78 5248925
## 2 76 4915944 2 332981 2.5739 0.08287 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1H0: Model fit3 does not fit the data better than fit2.
H1: Model fit3 fits the data better than fit2.
The null hypothesis states that adding Parking and Balcony does not improve the model fit compared to using only Age and Distance. The alternative hypothesis states that the extended model (fit3) explains more variation than the simpler model (fit2). The ANOVA result shows F=2.57,p=0.083. Since the p-value is above 0.05, we do not reject the null hypothesis, meaning fit3 does not significantly improve the fit over fit2, although there is weak evidence at the 10% level.
================================
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 + ParkingFactor + BalconyFactor,
## data = mydata3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -403.16 -204.08 -41.24 239.52 528.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2372.693 93.344 25.419 < 2e-16 ***
## Age -6.906 3.118 -2.215 0.0298 *
## Distance -22.254 2.839 -7.837 2.25e-11 ***
## ParkingFactorYes 137.479 60.854 2.259 0.0267 *
## BalconyFactorYes 17.234 57.099 0.302 0.7636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 254.3 on 76 degrees of freedom
## Multiple R-squared: 0.5621, Adjusted R-squared: 0.539
## F-statistic: 24.38 on 4 and 76 DF, p-value: 5.29e-13Explanation of the output:
Apartments with Parking are priced about 137 EUR higher on average than those without, significant at the 5% level (p = 0.027), ceteris paribus (other variables unchanged)
Having a Balcony adds about 17 EUR on average, but this effect is not statistically significant (p = 0.764), ceteris paribus (other variables unchanged)
The model explains about 56% of the variation in apartment prices (Adjusted R² = 0.54).
F-test hypotheses:
H0: All regression coefficients (except the intercept) are equal to 0.
H1: At least one predictor has a non-zero effect on Price.
The highly significant F-statistic (F(4,76) = 24.38, p < 0.001) leads us to reject H₀, confirming that the model explains a substantial portion of the variation in Price.
================================
Save fitted values and claculate the residual for apartment ID2.
================================
# Save fitted values and residuals in the dataset
mydata3$fitted_fit3 <- fitted(fit3)
mydata3$residual_fit3 <- residuals(fit3)
# Show fitted value and residual for Apartment ID 2
mydata3[2, c("Price", "fitted_fit3", "residual_fit3")]## Price fitted_fit3 residual_fit3
## 2 2800 2363.611 436.3894The observed price of apartment ID 2 is 2800 EUR.
The fitted value from model fit3 is 2363.6 EUR.
The residual, calculated as observed minus fitted, is 436.4 EUR.
This means the model underestimates the price of apartment 2 by about 436 EUR.