data <- matrix(scan(),ncol=4,byrow=TRUE)
 1 1 4 15.5
 1 2 2 33.9
 1 3 3 13.2
 1 4 1 29.1
 2 1 2 16.3
 2 2 3 26.6
 2 3 1 19.4
 2 4 4 22.8
 3 1 3 10.8
 3 2 1 31.1
 3 3 4 17.1
 3 4 2 30.3
 4 1 1 14.7
 4 2 4 34.
 4 3 2 19.7
 4 4 3 21.6

data

Gas1 <- data.frame(data)
colnames(Gas1) <- c("driver","car","gas","mileage")

Gas1$driver <- as.factor(Gas1$driver)    # to make sure that driver is a factor, not numeric
Gas1$gas <- as.factor(Gas1$gas)    # to make sure that gas is a factor, not numeric
Gas1$car <- as.factor(Gas1$car)    # to make sure that car is a factor, not numeric

rm(data)

Gas1

# LS analysis

gas1.lm <- lm(mileage ~ driver + car + gas, data = Gas1)

summary(gas1.lm)
anova(gas1.lm)

 # check model assumptions

par(mfrow=c(2,2))  
plot(gas1.lm)          
par(mfrow=c(1,1)) 


 # use the emmeans package to conduct pairwise tests on the treatment
 
library(emmeans)

gas1.lm.means.treatment <- emmeans(gas1.lm,"gas")
summary(gas1.lm.means.treatment)
pairs(gas1.lm.means.treatment)



rm(Gas1)