In Westchester County, north of New York City, Consolidated Edison bills residential customers for electricity on a monthly basis. The company wants to predict residential usage, in order to plan purchases of fuel and budget revenue flow. The data includes information on usage (in kilowatt-hours per day) and average monthly temperature for 55 consecutive months for an all-electric home. Data on consumption of electricity and the temperature in Westchester County, NY.
attach(elusage)
head(elusage)
## Month Year Usage Temperature
## 1 8 1989 24.828 73
## 2 9 1989 24.688 67
## 3 10 1989 19.310 57
## 4 11 1989 59.706 43
## 5 12 1989 99.667 26
## 6 1 1990 49.333 41
splot(Usage, Temperature)
The scatterplot does not look all that bad but running the regression and checking the residual vs.fits plot shows that a linear model is not good for this dataset:
slr(Usage, Temperature, RFplot = TRUE)
## The least squares regression equation is:
## Usage = 116.716 - 1.365 Temperature
## R^2 = 78.05%
In our current discussion we are concerned with the shape of the model only. The other assumptions (normal residuals and equal variance) we can worry about only after we have a good model. So for the time being I will only print the residual vs. fits plot.
so let’s try the transformations:
slr(log(Usage), Temperature, RFplot = TRUE)
## The least squares regression equation is:
## log(Usage) = 5.335 - 0.032 Temperature
## R^2 = 80.58%
does not look bad
slr(Usage, log(Temperature), RFplot = TRUE)
## The least squares regression equation is:
## Usage = 320.358 - 70.314 log(Temperature)
## R^2 = 82.91%
does not look bad
slr(log(Usage), log(Temperature), RFplot = TRUE)
## The least squares regression equation is:
## log(Usage) = 9.92 - 1.599 log(Temperature)
## R^2 = 81.12%
does not look bad
What do we do now? All three of the transformations have a good residual vs fits plot, and so are good models. But which one is best?
We can choose using the
Say we are have two “good” models. How can we decide which is better? One measure of “quality of fit” is \(R^2\), the Coefficient of Determination. It is defined as
\[ R^2 = \text{cor}(\text{Observed Values, Predicted Values})^2 100\% \]
This definition immediately leads us to the following properties \(R^2\):
\(R^2=0\): no relationship between observed and predicted values.
\(R^2=100\): perfect relationship between observed and predicted values.
\(R^2\) of model A is greater than \(R^2\) of model B, then model A is better than model B.
\(R^2\) is part of the output of the slr routine.
If the relationship is linear there is a close connection between Pearson’s correlation coefficient and \(R^2\):
\[ R^2 = (r)^2100\% \]
Example Draft data:
detach(elusage)
attach(draft)
cor(Draft.Number, Day.of.Year)
## [1] -0.2260414
slr(Draft.Number, Day.of.Year, return.result = TRUE)
## Intercept Day.of.Year (Intercept) x
## 225.009 -0.226 5.110 0.000 0.000
\[ r^2 * 100 \% =(- 0.226)^2*100 \% = 0.051*100\% = 5.1 \% = R^2 \]
The big advantage of \(R^2\) is that it also works if the relationship is not linear, whereas Pearsons’s correlation coefficient does not.
Back to the electricity data. We find
exponential model \(R^2=80.6\%\)
y vs log x model \(R^2=82.9\%\)
power model \(R^2=81.1\%\)
so the y vs. log x model is best.
But how about a polynomial model? Unfortunately we cannot use \(R^2\) to choose between two polynomial models, say the quadratic and the cubic model, because the model with the higher power will never have a smaller \(R^2\)! For example a cubic model can never give a worse fit than the quadratic model.
Example: Elusage:
Model | \(R^2\) |
---|---|
Linear | \(78.05\%\) |
Quadratic | \(84.67\%\) |
Cubic | \(84.72\%\) |
… | … |
power 10 | \(\ge 84.72\%\) |
The reason for this is simple: Say we find the best quadratic model, which is \[ \text{Usage} = 196.7 - 4.640\text{ Temperature} + 0.03073\text{ Temperature}^2 \]
Now we add the cubic term Temperature3 as a predictor. One (of many) cubic models is
\[ \text{Usage} = 196.7 - 4.640\text{ Temperature} + 0.03073\text{ Temperature}^2+ 0.0\text{ Temperature}^3 \]
this is of course the same as the quadratic model above, so it has \(R^2=84.67\%\). Only the least squares cubic model is the best cubic model, so it’s \(R^2\) cannot be smaller (and because of statistical fluctuation usually will be even a bit higher, even if the cubic term is not necessary).
and this one is perfect, it has \(R^2=100\%\).
Actually, it is always possible to find a polynomial model which fits the data set perfectly, that is it has \(R^2=100\%\)!
But: we want our models to fit the relationship, not the random fluctuations in the dataset.
A model should be parsimoneous, that is as simple as possible. This is in agreement with one of the fundamental principles of science:
Ockham’s razor, named after William of Ockham
Ockham’s razor is the principle that “entities must not be multiplied beyond necessity” (entia non sunt multiplicanda praeter necessitatem). The popular interpretation of this principle is that the simplest explanation is usually the correct one.
For our problem this means: Use the polynomial model of lowest degree that can’t be improved statisticaly significantly by adding another power.
Let’s consider again the quadratic and the cubic models:
slr(Usage, Temperature, polydeg = 2, ndigit = 5)
## The least squares regression equation is:
## Usage = 196.71532 - 4.64046 Temperature +0.03073Temperature^2
## R^2 = 84.69%
slr(Usage, Temperature, polydeg = 3, ndigit = 5)
## The least squares regression equation is:
## Usage = 212.99485 - 5.68932 Temperature +0.05198Temperature^2 -0.00014Temperature^3
## R^2 = 84.72%
So the cubic model is better than the quadratic one (in terms of R2), but is it statistically significantly better? Here is a test:
quad.model <- slr(Usage, Temperature, polydeg = 2, return.model = TRUE)
cube.model <- slr(Usage, Temperature, polydeg = 3, return.model = TRUE)
nested.models.test(quad.model, cube.model)
## H0: both models are equally good.
## p value= 0.7516
and so the answer is no, the cubic is not better. So we should use the quadratic one!
Warning as always a GOOD MODEL is one with a good residual vs. fits plot. It can happen that both the quadratic and the cubic are bad models and this test fails to reject the null because they are equally bad!
Note: if we have two models, one of which is a special case of the other, we say we have nested models.
Example: quadratic and cubic
Example: y vs x and y vs x, z
In all of those cases the model with more predictors will NEVER have a smaller \(R^2\), so using \(R^2\) would always lead to the model with all the terms, which may not be best.
In choosing the best model (from our short list) proceed as follows: Model is “good” = no pattern in the Residual vs. Fits plot
If the linear model is not good, proceed as follows
check the transformation models and see which of these (if any) are good.
find the best polynomial model using method described above.
Choose as the best of the good models in 2) and 3) the one which has the highest \(R^2\).
Back to Electricity usage. We have found:
best transformation model is y vs log of x with \(R^2=82.9\%\)
best polynomial model is the quadratic with \(R^2=84.7\%\)
so best overall is quadratic.
An important comment
Having a high \(R^2\) is desirable, but neither necessary for a model to be good, nor an indication that a model is good:
Example 1970’s draft data:
slr(Draft.Number, Day.of.Year, RFplot = TRUE)
## The least squares regression equation is:
## Draft.Number = 225.009 - 0.226 Day.of.Year
## R^2 = 5.11%
the linear model is a good one, even though it has a very low \(R^2=5.1\%\)
Example fabric wear data:
attach(fabricwear)
slr(Wear, Speed, RFplot = TRUE )
## The least squares regression equation is:
## Wear = -6.947 + 0.274 Speed
## R^2 = 88.58%
the linear model is bad, even though it has a fairly high \(R^2=88.6\%\)
For this data set we have seen previously that none of the log transformations work. How about the polynomial models?
quad.model <- slr(Wear, Speed, polydeg = 2, return.model = TRUE)
cube.model <- slr(Wear, Speed, polydeg = 3, return.model = TRUE)
nested.models.test(quad.model, cube.model)
## H0: both models are equally good.
## p value= 0.3765
and so the cubic model is not better than the quadratic model. So we should use the quadratic one!