Week 3 tutorial for Experimental Design

Group 1 (Motor Trend Car Road Tests)

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)}

Describe and summarize your assigned data set.

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.

Use histograms to observe distributions, look for outliers and zeroes

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).

Graph your data and explore the relationship between car weight (wt) and miles per gallon (mpg)?

#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.

Conclusion: linear regression model was statistically significant, F(1,30) = 91.38, p < 0.001, indicating that vehicle weight significantly contributed to the fuel efficiency

Do cars with more cylinders (cyl) consume more fuel?

# 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()

Run Kruskal-wallis Test

# 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"
)
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"
)
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.