# Assessing a binary decission based on continuos and multi-level categorical variables

I have been asked to generate a tool to assess if a particular new set of measurements fit within a list of already accepted ones. The problem is that there are different categorical variables with multiple levels and with some levels having only a few or one set of measurements

My idea would be to generate a linear model and then look at the residuals of the new set of measurements and then assess if there are recurrent discrepancies with the model (i.e. all the measurements have higher/lower values than expected) or individual measurements that are very different from the expected value.

The data would look something like this:

Eventually, what I would like to do is to check if a new petrol consumption measurements set is similar/different than what is expected for that model in the measurement's conditions (temperature, factory, oil brand, number of wheels).

A new entry (red) would look like this compared with existing data:

I would like to test if the new entry comes from the same distribution or (or is an outlier) from my reference dataset. The new entry can be classified as faulty/failure if it is obvious that it comes from a different 'population' - e.g. in this case the new entry should be marked as a failure because the fuel consumption at 50mph is way too high.

I would like to do those diagnostics on the whole dataset (no separation by model or manufacturer) and then do the tests with the 'subsets' - e.g. by Model, Manufacturer, "Type_of_Motorcycle, Factory_of_production, etc to be able to do diagnostics what could have been the cause of my new entry failing - e.g. the new Centre has not assembled it properly, or because the new model appearing on the market is not good enough.

Note, that the trial entry has been simulated as being of an existing model with an existing compatible oil. As the technology develops, these manufacturers will release new models which may be compatible with absolutely new oil types.

Also note, that the new models assessment of fuel consumption at a particular speed may change - e.g. 2-wheeled Kawasaki Ninja99 may be tested for fuel consumption at speed 20mph instead of 10mph and it will still be required to test if the new motorcycle behaves in the expected way (as the mechanics and principle of action hasn't changed).

So far, I tried to model this using the following:

lm( Fuel_consumption ~ Speed + Speed:Type_of_Motorbike + Manufacturer + Model + Compatible_oil_brand )


The code to generate the synthetic data set. Sorry for the long code but that's the best way I know to explain the data:

## Variables

# You have two types of motorcycles being manufactured
Type_of_Motorbike <- c( "2_wheels", "3_wheels")

# From 3 different manufacturers:
Manufacturer <- c("BMW", "Honda", "Kawasaki")

# And in total among all the 3 manufacturers you have 22 unique motorbike models
Model_BMW <- c( paste( "Series", as.character( 1:8), sep = "_" ))
Model_Kawa <- c( paste( "Ninja", as.character( 1:7), sep = "_" ))
Model_Honda <- c( paste( "Honda", as.character( 1:7), sep = "_" ))

# You have 12 types (one included later) of oil you can put into your motorbike, some bikes accept only one type of
# oil while others can be put several different types
Compatible_oil_brand <- c("oil_A", "oil_B", "oil_C","oil_D",
"oil_E", "oil_F", "oil_G","oil_H",
"oil_J", "oil_K", "oil_L")

# You have 19 different factories which produce motorbikes - some of them can assembel several
# models of motorbikes from different manufacturers
Producing_factory <- c("Vilnius", "Pamplona", "London", "Glasgow", "Kaunas", "Moletai", "Visaginas",
"Prague", "Moscow", "Warsaw", "Leeds", "Klaipeda","Valencia","Huesca",
"Sofia", "Tokyo", "Paris", "Madrid", "Zaragoza", "Voronezh","Tudela",
"SanSebastian", "Vienna", "Ostrava", "Budapest", "Toledo", "Alicante", "Berlin",
"Bristol", "Lima", "Minsk", "Bedford", "Tallinn", "Elizondo", "Osaka", "Sochi",
"BuenosAires", "Bejing", "HongKong", "Riga", "Seville", "Porto", "Nagano", "Santander")

# You do a test how good a bike is by measuring its fuel consumption per mile at different speeds:
Test_speed_mph <- c(10, 30, 50, 70 ,90)

# You will also have a temperature of engine variable which reflects the temperature at which engine was at before
# you started the fuel consumption test.

#########
set.seed(123)
motorcycles <- data.frame(Type_of_Motorbike = sample(Type_of_Motorbike, 96, replace = TRUE),
Manufacturer = c( rep( "BMW", 52 ), rep( "Kawasaki", 32 ), rep( "Honda", 12 )),
Model = c( sample(Model_BMW, 52, replace = TRUE), sample(Model_Kawa, 32, replace = TRUE), sample(Model_Honda, 12, replace = TRUE) ),
Compatible_oil_brand= c( sample(Compatible_oil_brand, 95, replace = TRUE), "oil_M" ),
Producing_factory = sample(Producing_factory, 96, replace = TRUE),
Speed_10mph = rnorm(96, mean= 15, sd=7),
Speed_30mph = rnorm(96, mean= 40, sd=9),
Speed_50mph = rnorm(96, mean= 68, sd=20),
Speed_70mph = rnorm(96, mean= 98, sd=34),
Speed_90mph = rnorm(96, mean= 145, sd=45),
Temperature_of_engine = rnorm(96, mean=25, sd=10)
)


The data would look like this:

library(tidyverse)
library(dplyr)

######################################## Visualising the data:
# install.packages("data.table")
library(data.table)
motorcycles_melted <- melt(data = motorcycles,
id.vars = c( "Manufacturer",
"Model",
"Type_of_Motorbike",
"Producing_factory",
"Compatible_oil_brand",
"Temperature_of_engine"),
variable.name = "Speed",
value.name = "Fuel_consumption" )

# converting categorical variable speed into numeric:
table(motorcycles_melted$$Speed) motorcycles_melted$$Speed <- as.character(motorcycles_melted$$Speed) table(motorcycles_melted$$Speed)

motorcycles_melted$$Speed[motorcycles_melted$$Speed=="Speed_10mph"] <- 10
motorcycles_melted$$Speed[motorcycles_melted$$Speed=="Speed_30mph"] <- 30
motorcycles_melted$$Speed[motorcycles_melted$$Speed=="Speed_50mph"] <- 50
motorcycles_melted$$Speed[motorcycles_melted$$Speed=="Speed_70mph"] <- 70
motorcycles_melted$$Speed[motorcycles_melted$$Speed=="Speed_90mph"] <- 90
motorcycles_melted$$Speed <- as.numeric(as.character(motorcycles_melted$$Speed))

ggplot(data = motorcycles_melted, aes(x=Speed, y=Fuel_consumption, color=Type_of_Motorbike)) +
geom_point() +theme_classic() +
xlab("Speed") + ylab("Fuel Consumption per Mile")


And I would get a new set of measurements like this:

#### I get a new model and I want to see if it is within the distribution. Note, the new entry
# comes from a new Factory centre - in Miami!
trial <- data.frame(Type_of_Motorbike = c( rep(x= "2_wheels", 5)),
Manufacturer = c( rep(x="BMW", 5)),
Model = c(rep(x="Series_1",5)),
Compatible_oil_brand= c(rep(x="oil_C",5)),
Producing_factory = c(rep(x="Miami",5)),
Speed = c(10, 30, 50, 70,90),
Fuel_consumption= c(10, 36, 200, 150, 100),
Temperature_of_engine = c(rep(x= 18, 5)))

summary.statistics3 <-  aggregate(Fuel_consumption ~ Speed, data=motorcycles_melted, mean)

# Each of the black dots is an individual measurement of fuel consumption per mile, there is no
# separation by manufacturer or model. The red dots are the simulated entry - fuel
#consumption of my newly entered motorbike entry at the corresponding speed. The green dot is the
# mean of all black dots - i.e. the mean fuel consumption at 10mph across all the models
# and manufacturers:

ggplot(data = motorcycles_melted, aes(x=Speed, y=Fuel_consumption)) +
geom_point() +theme_classic() +
xlab("Speed") + ylab("Fuel Consumption per Mile") +
geom_point(data = trial, aes(x = Speed, y = Fuel_consumption, shape = factor(Type_of_Motorbike), color="red")) +
theme(axis.title.x = element_text(size=14),axis.title.y = element_text(size=14), legend.position = "none" ) +
geom_point(data = summary.statistics3, aes(x = Speed, y = Fuel_consumption), color= "green")

# Another exploratory plot:
# Each of the black dots is an individual measurement of fuel consumption per mile for that particular
# model, in this case Series 12. The red dots are the simulated entry - fuel
#consumption of my newly entered motorbike entry of Series 12 at the corresponding speed.
#The green dot is the mean of all black dots - i.e. the mean fuel consumption at 10mph across
# the data for Series 12 in the reference dataset.
summary.statistics4 <- aggregate(Fuel_consumption ~ Model+ Speed, data=motorcycles_melted, mean)

ggplot( subset(motorcycles_melted, motorcycles_melted$$Model %in% (trial$$Model)), aes(x=Speed, y=Fuel_consumption)) + geom_point() +theme_classic() +
xlab("Speed") + ylab(paste ("Fuel Consumption per Mile",trial$$Model, sep = " ", collapse = NULL))+ geom_point(data = trial, aes(x = Speed, y = Fuel_consumption, shape = factor(Model), color="red")) + theme(axis.title.x = element_text(size=14),axis.title.y = element_text(size=14), legend.position = "none" ) + geom_point(data = subset(summary.statistics4, summary.statistics4$$Model %in% (trial\$Model)), aes(x = Speed, y = Fuel_consumption), color="green")

• Try to add the plots so we can have a look at the data (I or other will edit them into the question). BTW, at the moment the question reads as a wall of text and code (the impression being that anyone trying to answer it should invest a lot of time) maybe that's why you don't receive and comment or answer. – llrs Jul 29 at 10:37
• Thanks for the suggestion. I have added a more "summarized" version on top and left the code for the synthetic data set in the bottom. I hope it helps to make it look a bit less "text/code wall". – EduardoGC Jul 29 at 11:50