Members: MOHAMAD ZAHID BIN BAHAROM, ELOISE JINGYI YU, LIM HUI EN
# load the packages needed
library(ggplot2)
library(dplyr)
# Load the ‘mtcars’ data set
{data(mtcars)}
The mtcars dataset contains 32 observations and 11
variables.
| Variable | Description |
|---|---|
mpg |
Miles per US gallon |
cyl |
Number of cylinders |
disp |
Displacement (cu. in.) |
hp |
Gross horsepower |
drat |
Rear axle ratio |
wt |
Weight (1000 lbs) |
qsec |
1/4 mile time |
vs |
Engine (0 = V-shaped, 1 = straight) |
am |
Transmission (0 = automatic, 1 = manual) |
gear |
Number of forward gears |
carb |
Number of carburetors |
summary(mtcars)
## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000
We are looking at mpg, cyl and wt variables. Based on the summary, there is no NAs and zeros in the 3 variables we are interested in.
Now let’s use histograms to look for normality, outliers/data entry errors, zeroes.
hist(mtcars$mpg, breaks = 25)
hist (mtcars$cyl, breaks = 20)
hist (mtcars$wt, breaks = 25)
-
mpg looks a little right skewed aka not normal. -
cyl seemed to have 3 groups. Maybe we can treat this as a
single categorical variable with > 2 groups. - wt appear
to have normal distribution.
We should also check for any violations of the 4 assumptions of linear model (linearity, homoscedasticity, independence and normality).
Linearity is seldom meaningful to check when x is a categorical variable. We can assume independence is not violated as each entry is unique.
We’ll use diagnostic plot for homoscedasticity check.
data_mod <- lm (mpg ~ wt + cyl, data = mtcars)
plot(data_mod)
The first plot (Residuals vs Fitted values) suggests a mild homoscedasticity (slight unequal variance).
#Q2 - linear models, report the stats.
#RQ = is a heavier car associated with lower fuel efficiency?
#Fit the linear model
model_1 <- lm(mpg ~ wt, data = mtcars) #check the linear model between mpg and wt
#Visualise the relationship between car weight (wt) and miles per gallon (mpg)
plot(mtcars$wt, mtcars$mpg,
xlab = "Weight (1,000 lbs)",
ylab = "Fuel efficiency (mpg)",
main = "Fuel Efficiency and Car Weight",
pch = 19,
col= "steelblue")
abline(model_1, col = "red", lwd = 2)
#There is a negative linear relationship between car weight and MPG, indicating that heavier cars generally have lower fuel efficiency.
summary(model_1) #summary of the model
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
#extract the important values from the model
coef(model_1) #check the intercept and slope of the model
## (Intercept) wt
## 37.285126 -5.344472
confint(model_1) #95% confidence interval (CI) between these variables
## 2.5 % 97.5 %
## (Intercept) 33.450500 41.119753
## wt -6.486308 -4.202635
summary(model_1)$r.squared #R2 values
## [1] 0.7528328
#Y = a + bX, regression equation
#y = mpg
#x = wt
#intercept (a) = 37.29
#slope (b) = 5.34
#fitted equation : mpg = 37.29 - 5.34(wt)
#display the coefficients result of the model
coefficient_results <- data.frame(
Estimate = summary(model_1)$coefficients[, "Estimate"],
SE = summary(model_1)$coefficients[, "Std. Error"],
Lower_95_CI = confint(model_1)[, 1],
Upper_95_CI = confint(model_1)[, 2],
t_value = summary(model_1)$coefficients[, "t value"],
p_value = summary(model_1)$coefficients[, "Pr(>|t|)"])
print(coefficient_results, digits = 4) #reduce the decimal digits
## Estimate SE Lower_95_CI Upper_95_CI t_value p_value
## (Intercept) 37.285 1.8776 33.450 41.120 19.858 8.242e-19
## wt -5.344 0.5591 -6.486 -4.203 -9.559 1.294e-10
statistically significant difference between mpg and wt (p < 0.001) 95% confidence interval (CI) = a plausible value range for the true population effect based on the sample size, it has lower and upper limit 95% CI suggests that additional 1000 pounds of car weight is associated with a true average decrease of 4.20 to 6.49 mpg (estimated decrease 5.34)
#display the overall results of the model
overall_results <- data.frame(
N = nobs(model_1), #sample size
R_squared = summary(model_1)$r.squared, #R2
Adjusted_R_squared = summary(model_1)$adj.r.squared, #adjusted R2
Residual_SE = summary(model_1)$sigma, #residual standard error
Degrees_of_freedom = summary(model_1)$df[2], #degree of freedom
F_statistic = unname(summary(model_1)$fstatistic["value"]), #F-statistic
Model_p_value = pf(
summary(model_1)$fstatistic["value"], #p-value indicating significant differences
summary(model_1)$fstatistic["numdf"],
summary(model_1)$fstatistic["dendf"],
lower.tail = FALSE))
print(overall_results, digits = 4) #reduce the decimal digits
## N R_squared Adjusted_R_squared Residual_SE Degrees_of_freedom
## value 32 0.7528 0.7446 3.046 30
## F_statistic Model_p_value
## value 91.38 1.294e-10
R2 = 0.7528 (75.3% variation was observed in mpg among these cars) residual standard error suggests that observed values typically vary roughly 3.05 mpg around the fitted line.
# we approach this relationship by using two different figures; (1) box plot and (2) violin plot
# boxplot using ggplot
ggplot(mtcars,aes(x = as.factor(cyl), y = mpg)) + #Fix: treat cylinders as a category
geom_boxplot(width = 0.5,outlier.shape = 1) + #to plot box plot, with 0.5 width and shape 1
labs(title = "Do cars with more cylinders consume more fuel?", #the title of the figure
subtitle = "More cylinders contribute to low fuel consumptions", #conclusion of the figure
x = "Cylinders",
y = "Miles per Gallon (mpg)") +
theme_classic(base_size = 14)
# use factor(__) for Cly, as one 3 cly values 4,6 and 8, treat as Categorical groups
# colour based on Cly value
# geom_violin() - trim = FALSE so do not cut of at min and max value
# overlay boxplot and points onto violin plot
# Create the violin plot with an overlaid boxplot
ggplot(mtcars, aes(x = factor(cyl), y = mpg, fill = factor(cyl))) +
geom_violin(trim = FALSE, alpha = 0.3) +
geom_boxplot(width = 0.15, fill = "white", color = "black") +
geom_jitter(width = 0.05, alpha = 0.3, color = "black") +
labs(
title = "Cylinders vs Miles per Gallon",
x = "Number of Cylinders",
y = "Miles Per Gallon (MPG)",
fill = "Cylinders"
) +
theme_minimal()
# use bulit in kruskal.test() function
# Kruskal-Wallis test
kruskal_result <- kruskal.test(mpg ~ factor(cyl), data = mtcars)
kruskal_result
##
## Kruskal-Wallis rank sum test
##
## data: mpg by factor(cyl)
## Kruskal-Wallis chi-squared = 25.746, df = 2, p-value = 2.566e-06
# run anova to see
anova_result <- aov(mpg ~ factor(cyl), data = mtcars)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(cyl) 2 824.8 412.4 39.7 4.98e-09 ***
## Residuals 29 301.3 10.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#presenting results in table for kruskal.test
kruskal_table <- data.frame(
Chi_squared = kruskal_result$statistic,
df = kruskal_result$parameter,
p_value = kruskal_result$p.value)
knitr::kable(
kruskal_table,
digits = 6,
caption = "Kruskal-Wallis test results for MPG among cylinder groups"
)
| Chi_squared | df | p_value | |
|---|---|---|---|
| Kruskal-Wallis chi-squared | 25.74616 | 2 | 3e-06 |
#presenting results in table for anova_result
anova_table <- as.data.frame(summary(anova_result)[[1]])
knitr::kable(
anova_table,
digits = 3,
caption = "One-way ANOVA results for MPG among cylinder groups"
)
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| factor(cyl) | 2 | 824.785 | 412.392 | 39.698 | 0 |
| Residuals | 29 | 301.263 | 10.388 | NA | NA |
The Kruskal–Wallis test was more suitable because the variance
differed between the cylinder groups
and the data were not fully normally distributed.
Both tests showed a significant difference in MPG among the cylinder groups (Kruskal–Wallis: χ²(2) = 25.746, p < 0.001, ANOVA: F(2, 29) = 39.7, p < 0.001). Therefore, the null hypothesis that MPG does not differ among the 4, 6, and 8 cylinder groups was rejected.
Conclusion #### Overall, cars with more cylinders generally had lower MPG, indicating that they were less fuel-efficient and consumed more fuel.