Interpreting Linear Prediction Models

Data exploration and preparation

The prediction task is the following: Can we predict the ozone level given solar radiation, wind speed, and temperature? To see whether the assumptions of the linear model are appropriate for the data at hand, we will compute the correlation between the variables:

# scatterplot matrix 
plot(ozone)

# pairwise variable correlations
cors <- cor(ozone)
print(cors)

##              Ozone    Solar.R       Wind       Temp
## Ozone    1.0000000  0.3483417 -0.6124966  0.6985414
## Solar.R  0.3483417  1.0000000 -0.1271835  0.2940876
## Wind    -0.6124966 -0.1271835  1.0000000 -0.4971897
## Temp     0.6985414  0.2940876 -0.4971897  1.0000000

# which variables are highly correlated, exclude self-correlation
cor.idx <- which(abs(cors) >= 0.5 & cors != 1, arr.ind = TRUE)
cor.names <- paste0(colnames(ozone)[cor.idx[,1]], "+", 
            colnames(ozone)[cor.idx[,2]], ": ", round(cors[cor.idx], 2))
print(cor.names)

## [1] "Wind+Ozone: -0.61" "Temp+Ozone: 0.7"   "Ozone+Wind: -0.61"
## [4] "Ozone+Temp: 0.7"

Since ozone is involved in two linear interactions, namely:

• Ozone has a positive correlation with temperature
• Ozone has a negative correlation with wind

This suggests that it should be possible to form a linear model predicting the ozone level using the remaining features.

Splitting into training and test sets

We will take 70% of samples for training and 30% for testing:

set.seed(123)
N.train <- ceiling(0.7 * nrow(ozone))
N.test <- nrow(ozone) - N.train
trainset <- sample(seq_len(nrow(ozone)), N.train)
testset <- setdiff(seq_len(nrow(ozone)), trainset)

The residuals

The first piece of information we obtain is on the residuals. The term residual comes from the residual sum of squares (RSS), which is defined as

\[\begin{align*} r_i &= y_i - f(x_i) \\ \text{RSS} &= \sum_{i=1}^n r_i^2 \end{align*}\]

where the residual \(r_i\) is defined as the difference between observed and predicted values, \(f(x_i)\), from the observed value, \(y_i\).

The residual median value suggest that the model generally predicts slightly higher values of ozone than observed. The large maximum value, however, indicates that some outlier predictions are also much too low. Looking at a numbers can be a bit abstract, so let us plot the predictions of the model against the observed values:

res <- residuals(model)
# ensure that x- and y-axis have the same range
pred <- model$fitted.values
obs <- ozone[trainset, "Ozone"]
# determine maximal range
val.range <- range(pred, obs)
plot(obs, pred, 
    xlim = val.range, ylim = val.range,  
    xlab = "Observed ozone concentration", 
    ylab = "Predicted ozone concentration",
    main = "Residuals of the linear model for the training data")
# show ideal prediction as a diagonal
abline(0,1, col = "red")
# add residuals to the plot
segments(obs, pred, obs, pred + res)

The coefficients

Now that we understand the residuals, let’s take a look at the coefficients. We can use the coefficients function to obtain the fitted coefficients of the model:

coefficients(model)

##  (Intercept)      Solar.R         Temp         Wind 
## -74.18928801   0.04506998   1.84999289  -3.34510459

Note that the intercept of the model has quite a low value. This is the value that the model will predict in the case that all independent values are zero. The low coefficient for Solar.R indicates that solar radiation does not play an important role for predicting ozone levels, which does not come as a surprise as it did not show a large correlation with the ozone level in our exploratory analysis. The coefficient for Temp indicates that ozone levels are high if the temperature is high (since ozone will form faster). The coefficient for Wind tells us that ozone levels will be low for fast winds (since ozone will be blown away).

The other values associated with the coefficients give information about the statistical certainty of the estimates.

model.summary$coefficients

##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -74.18928801 28.65267764 -2.589262 1.157483e-02
## Solar.R       0.04506998  0.02843762  1.584871 1.172601e-01
## Temp          1.84999289  0.31695917  5.836692 1.317712e-07
## Wind         -3.34510459  0.80164181 -4.172817 8.103673e-05

• Std. Error is the standard error of the coefficient estimate
• t value indicates the value of the coefficient in terms of the standard error
• Pr(>|t|) is the p-value for the t-test, which indicates the significance of the test statistic

The standard error

The standard error of the coefficients is defined as the standard deviation of the feature variance:

\[\hat{\sigma} = \sqrt{\text{diag}\left(\sigma_{\text{err}}^2 (X^T X)^{-1}\right)}\]

where \({X^T X}^{-1}\) is the inverse of the design matrix \(X\), which includes the intercept, and where \(\sigma_{\text{err}}^2\) is the variance of the error.

In R, the standard errors of the model estimates can be computed via:

# compute variance-covariance matrix of model estimates
var.matrix <- vcov(model)
print(var.matrix)

##              (Intercept)       Solar.R         Temp          Wind
## (Intercept) 820.97593617  0.0787079320 -8.596270586 -16.239437823
## Solar.R       0.07870793  0.0008086984 -0.002707639  -0.002198751
## Temp         -8.59627059 -0.0027076392  0.100463112   0.131551223
## Wind        -16.23943782 -0.0021987513  0.131551223   0.642629584

# use variances from variance-covariance matrix
std.err <- sqrt(diag(vcov(model)))
print(std.err)

## (Intercept)     Solar.R        Temp        Wind 
## 28.65267764  0.02843762  0.31695917  0.80164181

Now, you are probably wondering where the values from vcov are coming from. It is defined as the variance-covariance matrix of the design matrix scaled by the variance of the error:

scaled.var.matrix <- model.summary$cov.unscaled * model.summary$sigma^2
print(scaled.var.matrix)

##              (Intercept)       Solar.R         Temp          Wind
## (Intercept) 820.97593617  0.0787079320 -8.596270586 -16.239437823
## Solar.R       0.07870793  0.0008086984 -0.002707639  -0.002198751
## Temp         -8.59627059 -0.0027076392  0.100463112   0.131551223
## Wind        -16.23943782 -0.0021987513  0.131551223   0.642629584

The variance that is used to scale to unscaled variance-covariance matrix is the estimated variance of the error, which is defined by

\[ \sigma_{\text{err}}^2 = \frac{1}{n-p} \text{RSS}\]

The cov.unscaled parameter is simply the variance-covariance matrix of all features \(\beta_i\) (including the intercept) of the design matrix:

# include intercept as a feature via 'model.matrix'
X <- model.matrix(model) # design matrix
# solve X^T %*% X = I to find the inverse of X^T * X
unscaled.var.matrix <- solve(crossprod(X), diag(4))
print(paste("Is this the same?", isTRUE(all.equal(unscaled.var.matrix, model.summary$cov.unscaled, check.attributes = FALSE))))

## [1] "Is this the same? TRUE"

# divide model coefficients by standard errors
t.vals <- coef(model) / std.err
print(t.vals)

## (Intercept)     Solar.R        Temp        Wind 
##   -2.589262    1.584871    5.836692   -4.172817

model.df <- df.residual(model)

Further statistics

The following additional statistics are provided by the summary function: the multiple R-squared, the adjusted R-squared, and the F-statistic. A nice summary of these statistics is available at StackExchange.

rss <- sum(res^2)
mean.rss <- rss / model.df
res.std.error <- sqrt(mean.rss)
print(res.std.error)

## [1] 22.24615

r.squared <- cor(pred, obs)^2
print(r.squared)

## [1] 0.6219815

n <- length(trainset) # number of samples
p <- length(coefficients(model)) - 1 # number of features w/o intercept
r.squared.adj <- 1 - (1 - r.squared) * ((n - 1) / (n - p - 1))
print(r.squared.adj)

## [1] 0.6066564

# create intercept-only model
null.model <- lm(formula = Ozone ~ 1, ozone, trainset)
print(null.model)

## 
## Call:
## lm(formula = Ozone ~ 1, data = ozone, subset = trainset)
## 
## Coefficients:
## (Intercept)  
##       44.13

# plot the incercept-only prediction function
plot(ozone$Ozone[trainset], null.model$fitted.values)

rss <- function(model) {
    return(sum(model$residuals^2))
}
f.statistic <- function(model1, model2) {
    n <- length(model1$fitted.values)
    n2 <- length(model2$fitted.values)
    if (n != n2) {
        stop("Both models should have the same number of training samples!")
    }
    p1 <- length(coefficients(model1))
    p2 <- length(coefficients(model2))
    enum <- (rss(model1) - rss(model2)) / (p2 - p1)
    denom <- rss(model2) / (n - p2)
    out <- enum/denom
    return(out)
}
# compare the intercept-only model and the model with 3 features
f.statistic(null.model, model)

## [1] 40.58587