1.1 Coefficient of determination

R-squared \(\left( R^2 \right)\), or the coefficient of determination, is a key statistical measure used to evaluate the goodness of fit for a linear model. It indicates the proportion of the variance in the dependent variable that can be explained by the independent variables in the model. To calculate \(R^2\), we first need to evaluate specific sums of squares based on the linear model:

\[ \text{Total Sum of Squares} \quad TSS = \sum_{i = 1}^n(Y_i - \bar Y)^2, \\ \text{Sum of Squares due to Regression} \quad SSR = \sum_{i = 1}^n(\hat Y_i - \bar{\hat Y})^2 = \sum_{i = 1}^n(\hat Y_i - \bar{ Y})^2, \\ \text{Sum of Squared Errors} \quad SSE = \sum_{i = 1}^n(Y_i - \hat Y_i)^2, \]

where \(n\) represents the size of the dataset, \(Y_i\) and \(\hat Y_i\) denote the \(i\)-th outcome and the fitted value, respectively. The means of all outcomes and all fitted values are represented as:

\[ \bar{Y} = \sum_{i = 1}^n Y_i, \qquad \bar{\hat Y} = \sum_{i = 1}^n \hat Y_i. \]

It can be shown that the equality \(\bar{Y} = \bar{\hat Y}\) holds.

In model 1 on \(\texttt{fev}\) dataset we can calculate these sums of squares as shown in the code below.

y <- fev$FEV # outcome
e <- residuals(model1) # residuals
y.hat <- fitted(model1) # fitted values

tss <- sum((y-mean(y))^2) # total sum of squares
ssr <- sum((y.hat-mean(y.hat))^2) # sum of squares due to regression
ssr <- sum((y.hat-mean(y))^2) # sum of squares due to regression
sse <- sum(e^2) # sum of squared errors

The formula shows that \(R^2\) takes values within the interval \([0, 1]\). A higher \(R^2\) value suggests a better fit of the model to the data. For example, an \(R^2\) of 0.80 means that 80% of the variability in the dependent variable can be explained by the model.

In model 1 on \(\texttt{fev}\) dataset we can calculate \(R^2\) as shown in the code below.

R2 <- 1-sse/tss # R^2
R2 <- ssr/tss # R^2

\(R^2\) is also available in the summary of the model on the prelast line denoted as Multiple R-squared.

summary(model1)

## 
## Call:
## lm(formula = FEV ~ Height + Sex, data = fev)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6763 -0.2505  0.0001  0.2347  2.0722 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.390263   0.180082 -29.932  < 2e-16 ***
## Height       0.051272   0.001167  43.933  < 2e-16 ***
## SexMale      0.125123   0.033801   3.702 0.000232 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4265 on 651 degrees of freedom
## Multiple R-squared:  0.7587, Adjusted R-squared:  0.758 
## F-statistic:  1024 on 2 and 651 DF,  p-value: < 2.2e-16

\(R^2\) can be also extracted from the summary using the following code.

summary(model1)$r.squared

## [1] 0.7587368

A disadventage of \(R^2\) is that it does not account for the number of predictors in the model. A model with more predictors will certainly have a higher \(R^2\), even if those predictors do not contribute meaningfully to the model. Because of that, we don’t use \(R^2\) to compare nested models.

This variant adjusts for the number of predictors in the model, providing a more accurate measure when comparing models with different numbers of predictors. The formula shows that models with more predictors are penalized for its complexity. \(R^2_{adj}\) decreases if the additional predictors do not significantly improve the model. This helps to prevent overfitting.

In model 1 on \(\texttt{fev}\) dataset we can calculate \(R^2_{adj}\) as shown in the code below.

R2adj <- 1-sse/(nrow(fev)-summary(model1)$df[1])/tss*(nrow(fev)-1) # R^2_{adj}

\(R^2_{adj}\) is also available in the summary of the model on the prelast line denoted as Adjusted R-squared.

summary(model1)

## 
## Call:
## lm(formula = FEV ~ Height + Sex, data = fev)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6763 -0.2505  0.0001  0.2347  2.0722 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.390263   0.180082 -29.932  < 2e-16 ***
## Height       0.051272   0.001167  43.933  < 2e-16 ***
## SexMale      0.125123   0.033801   3.702 0.000232 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4265 on 651 degrees of freedom
## Multiple R-squared:  0.7587, Adjusted R-squared:  0.758 
## F-statistic:  1024 on 2 and 651 DF,  p-value: < 2.2e-16

\(R^2_{adj}\) can be also extracted from the summary using the following code.

summary(model1)$adj.r.squared

## [1] 0.7579956

The following table summarizes \(R^2\) and \(R^2_{adj}\) for Model 1 and Model 2 on \(\texttt{fev}\) data. We can see that \(R^2\) for Model 2 is higher, as it will always be for nested models with more predictors. However, \(R^2_{adj}\) is lower for Model 2, which indicates that the predictors added to Model 2 do not significantly improve the goodness of fit. We can conclude that Model 1 has better performance than Model 2.

1.2 Information criteria

We might also use information criteria to compare quality of model fit. Information criteria help assess how well a model explains the data, balancing goodness of fit with model complexity. They are used to attribute scores to different regression models.

When comparing multiple models, those with lower \(AIC\) or \(BIC\) scores are generally considered better fits to the data. These criteria help prevent overfitting, as they reward models that are simpler. While both \(AIC\) and \(BIC\) are used for model selection, \(BIC\) applies stricter penalties for adding more parameters. This makes \(BIC\) more suitable when you want to choose simpler models.

\(AIC\) and \(BIC\) can be computed as showed in the following code.

AIC_mod1 <- AIC(model1)
BIC_mod1 <- BIC(model1)

The following table summarizes \(AIC\) and \(BIC\) for Model 1 and Model 2 on \(\texttt{fev}\) data. We can see that both criteria are in favor of the Model 1.

3.1 Tasks

As an environmental scientist or sustainability consultant, you are tasked with analyzing this dataset and addressing the following inquiries:

For this analysis, it is imperative to fit linear models incorporating various predictors, compare their goodness of fit, and determine which model is the most appropriate for modelling the relationship between the \(CO_2\) emissions and other characteristics of the vehicle.

3.2 Data preprocessing

First, load the data set, rename the variables and change appropriate ones to type factor. Also delete the Cunsumption_mpg variable (you are working with L units). Check, whether there are any missing values. How can we handle missing data?

data <- read.csv(...)

colnames(data) <- c('Make', 'Model', 'Class', 'Engine_size',
                    'Cylinders', 'Transmission', 'Fuel_type',
                    'Consumption_city', 'Consumption_hwy',
                    'Consumption_comb', 'Cunsumption_mpg',
                    'CO2')
data <- data |> 
  select(!(Cunsumption_mpg)) |>
  mutate(
    ... = factor(...),
    ...
  )

Consider the variables Make, Model and Class. Are these variables suitable for modelling using linear regression?

One potential approach is to cluster the car makes into similar groups, artificial intelligence can be helpful here.

HINT: Our prompt for ChatGPT:

Hi, I have a list of car makes, can you cluster them so 
the clusters contain the most similar car makes 
(for example cluster named luxury car contains aston martin,
bentley, rolls-royce). Cluster all the listed car makes below
and do not add other make which is not listed below. 
Clusters will be disjoint and their union will equal to the whole
list below. Write the car makes in one cluster in capital (SUBARU),
wrap up each car make in quotation marks ("SUBARU") 
and separate them by a comma (,).
Here is the list of car makes: 
ACURA, ALFA ROMEO, ASTON MARTIN , AUDI, BENTLEY, BMW, BUICK,
CADILLAC, CHEVROLET, CHRYSLER, DODGE, FIAT, FORD, GMC, HONDA, 
HYUNDAI, INFINITI, JAGUAR, JEEP, KIA, LAMBORGHINI, LAND ROVER,
LEXUS, LINCOLN, MASERATI, MAZDA, MERCEDES-BENZ, MINI, MITSUBISHI,
NISSAN, PORSCHE, RAM, ROLLS-ROYCE, SCION, SMART, SRT, SUBARU, 
TOYOTA, VOLKSWAGEN, VOLVO, GENESIS, BUGATTI.

According to the answer, cluster the car makes into similar groups, defining new variable Type. It is possible that certain makes do not fit within any of the defined clusters. If appropriate, either remove these makes from the analysis or assign them to the cluster that is most similar, based on your judgment.

Delete unuseful variables. Functions filter(), mutate(), select() and case_when() can be helpful.

data <- data |> filter(Make != ...) |>
  mutate(Type = case_when(
    Make %in% c("...", ...) ~ "...",
    ...
  )
  ) |>
  select(!c(...))

Investigate the Fuel_type variable. Use table() function. If needed, delete appropriate rows.

data <- data |>
  filter(...) |>
  mutate(... = factor(....))

3.3 Explanatory data analysis

• Identification and handling outliers.

• Visualization of significant relationships among the variables.

• Examining the correlations between the variables.

The visualization provided below may serve as a source of inspiration. It is also recommended to employ two-dimensional exploratory data analysis (especially including outcome variable) to further enhance your insights.

The correlation structure of the data can be visualized using corrplot.mixed() function from corrplot package.

Hi, how can I make ggplot in R interactive?

3.4 Model comparison

model <- lm(CO2 ~ Engine_size + Cylinders + 
              Fuel_type + Consumption_comb + Type, 
            data = data)

Fit a range of models to the \(CO_2\) emissions data and evaluate their performance based on the criteria outlined in the theoretical framework of this seminar. Can you identify a model that outperforms the one provided by your fellow researcher? Keep in mind that adding more covariates does not necessarily lead to improved model performance.

Using the models you have explored throughout your analysis, address the primary questions of your main task.