I have found my data set in R with “data(Wong)” function, then I have used “write.csv(Wong,”Wong.csv”” function to save Wong data set into current working directory.
mydata <- read.table("./Wong.csv", header=TRUE, sep=",")
head(mydata)
## id days duration sex age piq viq
## 1 3358 30 4 0 20.67077 87 89
## 2 3535 16 17 0 55.28816 95 77
## 3 3547 40 1 0 55.91513 95 116
## 4 3592 13 10 0 61.66461 59 73
## 5 3728 19 6 0 30.12731 67 73
## 6 3790 13 3 0 57.06229 76 69
Explanation of variables:
id: Patient ID number
days: Number of days post coma at which IQs were measured
duration: Duration of the coma in days
sex: 0-male, 1-female
age: Age in years at the time of injury
piq: Performance in points - (i.e., mathematical) IQ
viq: Performance in points - verbal IQ
# Changing sex variable into categorical factor
mydata$sex <- factor(mydata$sex,
levels = c(0, 1),
labels = c("Male", "Female"))
head(mydata)
## id days duration sex age piq viq
## 1 3358 30 4 Male 20.67077 87 89
## 2 3535 16 17 Male 55.28816 95 77
## 3 3547 40 1 Male 55.91513 95 116
## 4 3592 13 10 Male 61.66461 59 73
## 5 3728 19 6 Male 30.12731 67 73
## 6 3790 13 3 Male 57.06229 76 69
# Checking if any n/a in the dataset
anyNA(mydata)
## [1] FALSE
Explanation: There are no n/a. Now we can perform further analysis. *any() checks if at least one of those is TRUE.
library(psych)
summary(mydata[,-1])
## days duration sex age
## Min. : 13.0 Min. : 0.0 Male :260 Min. : 6.513
## 1st Qu.: 59.0 1st Qu.: 1.0 Female: 71 1st Qu.:21.737
## Median : 150.0 Median : 7.0 Median :26.877
## Mean : 481.9 Mean : 14.3 Mean :31.853
## 3rd Qu.: 416.0 3rd Qu.: 16.0 3rd Qu.:40.923
## Max. :11628.0 Max. :255.0 Max. :80.033
## piq viq
## Min. : 50.00 Min. : 64.00
## 1st Qu.: 77.00 1st Qu.: 85.00
## Median : 87.00 Median : 94.00
## Mean : 87.56 Mean : 94.96
## 3rd Qu.: 97.00 3rd Qu.:105.00
## Max. :133.00 Max. :132.00
Explanation:
Mean for “duration”: On average patients were in coma for 14.3 days.
Median for “days”: Median is 150 days, it means that 50% of the patients were assessed on 150th day post-coma or before, while the other 50% were assessed after 150 days.
Maximum for “PIQ”: The maximum PIQ score in the sample was 133 points.
hist(mydata$piq,
main = "Distribution of Performance IQ (PIQ)",
xlab = "Performance IQ Score",
ylab = "Number of Patients",
col = "lightblue",
border = "white")
Explanation: As we can see from the histogram, the distribution of PIQ scores is slightly asymmetrical to the right or positively skewed.
library(ggplot2)
ggplot(mydata, aes(x = sex, y = piq, fill = sex)) +
geom_boxplot() +
labs(title = "Boxplot of PIQ by Gender",
x = "Gender",
y = "Performance IQ") +
theme_minimal()
Explanation of boxplot: The boxplots suggest that the PIQ distribution is roughly symmetric for both genders, both having high value outliers. We can say that both boxplots are fairly simillar, only what I can see is that there is less variability in points for men. I will further use additional describtive statistics to analyze graph better.
Using describtive statistics for better interpretation of boxplot: Median:
median(mydata$piq[mydata$sex == "Female"])
## [1] 87
median(mydata$piq[mydata$sex == "Male"])
## [1] 87
Explanation: Both male and female patients have a median PIQ of 87, meaning that 50% of the patients in each group scored 87 points and below while the other 50% have scored above 87 points.
Mean:
round(mean(mydata$piq[mydata$sex == "Female"]),2)
## [1] 89.18
round(mean(mydata$piq[mydata$sex == "Male"]),2)
## [1] 87.11
Explanation: On average females have scored 2.07 points higher than males.
Standard Deviation:
round(sd(mydata$piq[mydata$sex == "Female"]),2)
## [1] 18
round(sd(mydata$piq[mydata$sex == "Male"]),2)
## [1] 14.26
Explanation: On the other hand there is difference in variability between Males and Females. Male PIQ scores are less spread out, indicating their scores are more clustered around the mean.
Conclusion: On average, females achieved higher PIQ scores compared to males; however, their scores exhibit greater variability.
library(readxl)
# Importing dataset
mydata <- read_xlsx("~/Desktop/IMB/Bootcamp - R/IMB Bootcamp/R data/R - Task 2/Business School.xlsx")
mydata <- as.data.frame(mydata)
head(mydata)
## Student ID Undergrad Degree Undergrad Grade MBA Grade
## 1 1 Business 68.4 90.2
## 2 2 Computer Science 70.2 68.7
## 3 3 Finance 76.4 83.3
## 4 4 Business 82.6 88.7
## 5 5 Finance 76.9 75.4
## 6 6 Computer Science 83.3 82.1
## Work Experience Employability (Before) Employability (After) Status
## 1 No 252 276 Placed
## 2 Yes 101 119 Placed
## 3 No 401 462 Placed
## 4 No 287 342 Placed
## 5 No 275 347 Placed
## 6 No 254 313 Placed
## Annual Salary
## 1 111000
## 2 107000
## 3 109000
## 4 148000
## 5 255500
## 6 103500
library(ggplot2)
ggplot(mydata, aes(x = `Undergrad Degree`, fill = `Undergrad Degree`)) +
geom_bar() +
labs(title = "Distribution of Undergraduate Degrees",
x = "Undergraduate Degree",
y = "Number of Students") +
theme_minimal()
Answer: Graph is showing that Business Undergrad degree is the most common in the current generation of MBA students.
library(ggplot2)
ggplot(mydata, aes(x = `Annual Salary`)) +
geom_histogram(binwidth = 10000, fill = "skyblue", color = "black") +
labs(title = "Distribution of Annual Salary",
x = "Annual Salary",
y = "Number of Students") +
theme_minimal()
library(psych)
describe(mydata$`Annual Salary`)
## vars n mean sd median trimmed mad min max
## X1 1 100 109058 41501.49 103500 104600.2 25945.5 20000 340000
## range skew kurtosis se
## X1 320000 2.22 9.41 4150.15
Answer: The mean annual salary is about 109,058, while the median is slightly lower at 103,500. This gap between mean and median as well as graphical distribution tell us that we have a positively skewed distribution — a few students earn very high salaries (up to 340,000), pulling the mean upwards.
# t-test
t.test(mydata$`MBA Grade`, mu = 74, alternative = "two.sided")
##
## One Sample t-test
##
## data: mydata$`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.04055
Answer: P-value is 0.00915, we can reject null hypothesis at 1% significance level and mean of this year class is 76.04. This means that this year average MBA grade is statistically significantly higher than previous year average grade of 74.
# Cohen's d - standardized measure of effect size
mean_val <- mean(mydata$`MBA Grade`, )
sd_val <- sd(mydata$`MBA Grade`, )
mu0 <- 74
cohen_d <- (mean_val - mu0) / sd_val
round(cohen_d,2)
## [1] 0.27
Answer: Cohen’s* d is a standardized measure of effect size that expresses the mean difference in units of standard deviation. Cohen’s d = 0.27, indicates a small effect.
Conclusion: Results suggests that although the grade improvement is statistically significant, the practical difference from the hypothesized mean of 74 is relatively modest.
*Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed., pp. 19–25). Hillsdale, NJ: Lawrence Erlbaum Associates.
library(readxl)
mydata <- read_xlsx("./Apartments.xlsx")
mydata <- as.data.frame(mydata)
mydata$ID <- 1:nrow(mydata) #An ID column was added to index the rows for easier referencing
head(mydata)
## Age Distance Price Parking Balcony ID
## 1 7 28 1640 0 1 1
## 2 18 1 2800 1 0 2
## 3 7 28 1660 0 0 3
## 4 28 29 1850 0 1 4
## 5 18 18 1640 1 1 5
## 6 28 12 1770 0 1 6
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
ID: Apartment ID number
mydata$Parking <- factor(mydata$Parking,
levels = c(0, 1),
labels = c("No", "Yes"))
mydata$Balcony <- factor(mydata$Balcony,
levels = c(0, 1),
labels = c("No", "Yes"))
head(mydata)
## Age Distance Price Parking Balcony ID
## 1 7 28 1640 No Yes 1
## 2 18 1 2800 Yes No 2
## 3 7 28 1660 No No 3
## 4 28 29 1850 No Yes 4
## 5 18 18 1640 Yes Yes 5
## 6 28 12 1770 No Yes 6
t.test(mydata$Price, mu = 1900)
##
## One Sample t-test
##
## data: mydata$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.941
Answer: Since p < 0.01, we reject the null hypothesis at the 1% significance level. The test shows evidence that the mean apartment price per m2 is not equal to 1900€.
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.03401
Answers:
Regression coefficient: If the age of an apartment goes up by 1 year, the apartment price per m2 goes down by 8.98 €, on average.This negative coefficient is statistically significant at 5% level (p = 0.034).
Coefficient of determination: R2 = 0.053 means that only 5.3% of the variability in apartment prices is explained by the linear effect of age. The rest of variability is explained by other factors not included in this simple model (such as location, size, parking, balcony, etc.).
round(-sqrt(0.05302),2)
## [1] -0.23
library(car, carData)
scatterplotMatrix(mydata[ , c(-4, -5, -6)],
smooth = FALSE)
Answer: Based on scatterplot matrix there is no problem with multicolinearity, because the scatterplot between distance and age shows points widely scattered with almost no slope.
fit2 <- lm(Price ~ Age + Distance, data = mydata)
vif(fit2)
## Age Distance
## 1.001845 1.001845
mean(vif(fit2))
## [1] 1.001845
Answer: Since both predictors have VIF close to 1, there is no multicolinearity problem in this model. Multicolinearity would be a problem if VIF greater than 5.
mydata$StdResid <- round(rstandard(fit2), 2)
mydata$CooksD <- round(cooks.distance(fit2), 2)
head(mydata[order(mydata$StdResid),], 10)
## Age Distance Price Parking Balcony ID StdResid CooksD
## 53 7 2 1760 No Yes 53 -2.15 0.07
## 13 12 14 1650 No Yes 13 -1.50 0.01
## 72 12 14 1650 No No 72 -1.50 0.01
## 20 13 8 1800 No No 20 -1.38 0.01
## 35 14 16 1660 No Yes 35 -1.26 0.01
## 36 24 5 1830 Yes No 36 -1.19 0.01
## 54 30 17 1560 No No 54 -1.10 0.01
## 5 18 18 1640 Yes Yes 5 -1.07 0.01
## 28 18 19 1620 Yes No 28 -1.07 0.01
## 64 18 18 1640 Yes Yes 64 -1.07 0.01
head(mydata[order(-mydata$StdResid),], 10)
## Age Distance Price Parking Balcony ID StdResid CooksD
## 38 5 45 2180 Yes Yes 38 2.58 0.32
## 33 2 11 2790 Yes No 33 2.05 0.07
## 2 18 1 2800 Yes No 2 1.78 0.03
## 61 18 1 2800 Yes Yes 61 1.78 0.03
## 58 8 2 2820 Yes No 58 1.66 0.04
## 57 10 1 2810 No No 57 1.60 0.03
## 22 37 3 2540 Yes Yes 22 1.58 0.06
## 25 8 26 2300 Yes Yes 25 1.57 0.03
## 11 12 5 2700 Yes No 11 1.55 0.02
## 27 16 1 2750 Yes No 27 1.55 0.02
head(mydata[order(-mydata$CooksD),], 10)
## Age Distance Price Parking Balcony ID StdResid CooksD
## 38 5 45 2180 Yes Yes 38 2.58 0.32
## 55 43 37 1740 No No 55 1.44 0.10
## 33 2 11 2790 Yes No 33 2.05 0.07
## 53 7 2 1760 No Yes 53 -2.15 0.07
## 22 37 3 2540 Yes Yes 22 1.58 0.06
## 39 40 2 2400 No Yes 39 1.09 0.04
## 58 8 2 2820 Yes No 58 1.66 0.04
## 2 18 1 2800 Yes No 2 1.78 0.03
## 25 8 26 2300 Yes Yes 25 1.57 0.03
## 31 45 21 1910 No Yes 31 0.89 0.03
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")
Explanation: Since no standardized residuals exceed ±3, there are no clear outliers, but observation 38 has a Cook’s Distance of 0.32 indicating high influence on the model, so I have decided to remove it.
Removing ID 38:
library(dplyr)
mydata <- mydata %>%
filter(!ID == "38")
library(ggplot2)
fit2 <- lm(Price ~ Age + Distance, data = mydata)
mydata$StdFittedValues <- scale(fit2$fitted.values)
scatterplot(y = mydata$StdResid, x = mydata$StdFittedValues,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
Explanation: The scatterplot of standardized residuals against standardized fitted values shows no clear pattern. The residuals appear randomly distributed around zero with roughly constant variance. This suggests that the assumption of homoskedasticity holds for our regression model.
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.9488, p-value = 0.00219
Answer: The histogram shows that standardized residuals are roughly normally distributed but not perfectly normal. The Shapiro–Wilk test (p = 0.00219 < 0.05) indicates that we reject the null hypothesis of normality, meaning the residuals are not normally distributed.
Conclusion: Although the assumption of normality is formally violated, the sample size is sufficiently large for the Central Limit Theorem to apply. Therefore, based on the theory we covered, this violation is not a major concern and the regression results remain reliable.
fit2 <- 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-12
Explanation:
Intercept: Intercept of 2456.1 is the predicted price of an apartment when both Age and Distance are equal to 0. It has no explanatory meaning in this case.
Age: Holding Distance constant, each additional year of apartment age decreases the price by 6.46 €/m2, on average. The effect is statistically significant at the 5% level (p=0.044)
Distance: Holding Age constant, each additional kilometer away from the city center decreases the price by about 22.96 €/m2, on average. This effect is statistically significant at 0.1% level, since p=2.52e-12.
Coefficient of determination: R2 = 0.484 means that 48.4% of the variability in apartment prices is explained by the linear effects of age and distance. Compared to the model with only age (R2 = 0.053), the inclusion of distance improves explanatory power, indicating that distance is an important predictor of apartment prices.
fit3 <- lm(Price ~ Age + Distance + Parking + Balcony, data = mydata)
anova(fit2,fit3)
## Analysis of Variance Table
##
## Model 1: Price ~ Age + Distance
## Model 2: Price ~ Age + Distance + Parking + Balcony
## 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 ' ' 1
Answer: The p-value (0.03051) is below 0.05 significance level, which means that fit3 explains significantly more variance in apartment prices than fit2. Adding Parking and Balcony improves the model and model fit3 should be preferred over fit2.
summary(fit3)
##
## Call:
## lm(formula = Price ~ Age + Distance + Parking + Balcony, data = mydata)
##
## 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 ***
## ParkingYes 167.531 62.864 2.665 0.00933 **
## BalconyYes -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-12
Explanation of regression coefficient for both categorical variables:
Parking: Apartments with parking are priced, on average, 167.53 €/m2 higher compared to apartments without parking, holding other variables constant. This effect is statistically significant at 1% level (p=0.00933).
Balcony: The coefficient for Balcony is not statistically significant (p = 0.80), which means we cannot conclude that having a balcony leads to any difference in apartment prices compared to apartments without a balcony.
F-Statistics part:
Null hypothesis (H0): All slope coefficients are equal to 0 (None of the predictors explain any variation in apartment prices)
Alternative hypothesis (H1): At least one slope coefficient is not equal to 0
Since the F-statistic is 22.04 with p < 0.001, we reject H0 at 0.1% significance level and conclude that at least one predictor explain a significant amount of variation in apartment prices.
mydata$fitted_fit3 <- fitted(fit3)
mydata$resid_fit3 <- resid(fit3)
mydata[mydata$ID == 2, c("Price", "fitted_fit3", "resid_fit3")]
## Price fitted_fit3 resid_fit3
## 2 2800 2372.197 427.8029
Explanation: For apartment ID2, the actual price is 2800€, while the model fit3 predicts a price of 2372.20 €. The residual is therefore +427.80, meaning that the apartment is priced 427.80 € higher than what the model predicts. This positive residual indicates that the model underestimated the apartment’s price.