Lecture 31 Stepwise models
31.1 Automated Variable Selection Procedures
Automated variable selection procedures work along the following lines:
Choose a model to start with (often the model with no covariates, or the model with all covariates included).
Test to see if there is an advantage in adding or removing covariates.
Repeat adding/removing variables until there is no advantage in changing the model.
31.1.1 Forward Model Selection
In this technique:
All one-X variable models are computed, and the one which has the smallest p-value is chosen. (provided any of the simple linear regression models are significant at the 0.05% significance level, or whatever significance level is being used. )
Then all two-variable models which include the first chosen X-variable are computed. The variable whose regression coefficient is most significant (smallest p-value) is chosen (provided the p-value < 0.05)
Then all three X-variable models which include the two variables previously chosen are computed … We continue to add variables to the model until there are no more significant ones to add.
Then we stop.
Note the criterion for “significance” can be chosen by the user: e.g. One may choose to only include variables that have P-value < the usual significance level \(\alpha\)=0.05.
In other software (e.g. Minitab) the default entry level for for forward selection is P-value < \(\alpha\)= 0.25. This is because we don’t want to ignore any potentially useful variables. It is better to first find the possible predictors, and then weed them out, than to not know about them.
For computing the P-value for entering a variable, we have two things we can look at. On the one hand we can look to see if the t- ratio for that variable’s regression coefficient is significant. (For example, when df=60, the t-ratio will be significant at \(\alpha=0.05\) if \(|t| > 2\)).
Equivalently, we since we are comparing two nested models, we can use the anova() function (e.g. anova(\(model_1,model_2\)), and check if the F statistic is significant with P-value< 0.05.
The two approaches are equivalent since, when we add just one variable, \(F = t^2\) where \(F \sim F_{1,df}\).
31.1.2 Backwards Variable Selection
This starts with model containing all possible explanatory variables (and possibly interactions).
For each variable in turn, we investigate the effect of removing the variable (or interaction, if any) from current model.
Remove the least informative variable / interaction, unless this term is nonetheless supplying significant information about the response.
Go to step 2. Stop only if all variables/interactions in the current model are important.
31.1.3 Stepwise (direction =“both”)
The phrase Stepwise variable selection is usually taken to mean that we start with a forward selection, but after each variable is entered into the model we check whether any other variables have fallen below the significance criterion \(\alpha\), and if so we remove it from the model.
One approach is to use F tests (or equivalently t tests).
31.2 Example: Hours of sunshine
We consider the Climate data, but this time we will use the variable Sun as the response, and try to explain it in terms of the other variables (main effects only). We will use a backwards elimination approach. Initially we will do this by hand for illustration purposes.
31.2.0.1 Backwards Elimination
## climate = read.csv("Climate.csv", header = TRUE, row.names = 1)
climate.lm1 = lm(Sun ~ ., data = climate)
summary(climate.lm1)
Call:
lm(formula = Sun ~ ., data = climate)
Residuals:
Min 1Q Median 3Q Max
-252.41 -91.95 -2.43 94.71 254.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.703e+03 3.809e+03 -1.497 0.1459
Lat -5.104e+00 2.938e+01 -0.174 0.8634
Long 3.959e+01 1.947e+01 2.033 0.0520 .
MnJanTemp 7.200e+01 4.581e+01 1.572 0.1276
MnJlyTemp -1.887e+01 3.152e+01 -0.599 0.5544
Rain -8.373e-02 4.716e-02 -1.775 0.0871 .
Height 2.609e-01 2.833e-01 0.921 0.3652
Sea 2.064e+02 9.785e+01 2.109 0.0443 *
NorthIsland -1.856e+02 1.355e+02 -1.370 0.1820
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 158.7 on 27 degrees of freedom
Multiple R-squared: 0.5843, Adjusted R-squared: 0.4612
F-statistic: 4.745 on 8 and 27 DF, p-value: 0.001021We would drop the variable with highest P-value, in this case Lat. NB only drop one variable at a time. We can use the update() function to save typing.
Call:
lm(formula = Sun ~ Long + MnJanTemp + MnJlyTemp + Rain + Height +
Sea + NorthIsland, data = climate)
Residuals:
Min 1Q Median 3Q Max
-254.409 -89.690 -3.238 90.960 252.261
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.048e+03 3.193e+03 -1.894 0.0686 .
Long 3.966e+01 1.913e+01 2.073 0.0475 *
MnJanTemp 7.767e+01 3.157e+01 2.460 0.0203 *
MnJlyTemp -1.671e+01 2.847e+01 -0.587 0.5619
Rain -8.055e-02 4.271e-02 -1.886 0.0697 .
Height 2.894e-01 2.266e-01 1.277 0.2120
Sea 2.078e+02 9.581e+01 2.169 0.0387 *
NorthIsland -1.777e+02 1.254e+02 -1.417 0.1674
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 155.9 on 28 degrees of freedom
Multiple R-squared: 0.5839, Adjusted R-squared: 0.4799
F-statistic: 5.613 on 7 and 28 DF, p-value: 0.000405The highest P-value is for MnJlyTemp, so we update to drop this one as well. Note instead of outputting the whole summary we can use the function drop1() to inform us of our choices in terms of dropping further variables.
Single term deletions
Model:
Sun ~ Long + MnJanTemp + Rain + Height + Sea + NorthIsland
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 689322 368.96
Long 1 98319 787641 371.76 4.1363 0.05121 .
MnJanTemp 1 140484 829807 373.64 5.9102 0.02147 *
Rain 1 142964 832286 373.74 6.0145 0.02044 *
Height 1 76044 765366 370.73 3.1992 0.08413 .
Sea 1 125566 814888 372.98 5.2826 0.02894 *
NorthIsland 1 107455 796777 372.17 4.5207 0.04212 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The F test indicates that Height has p-value =0.08, so we can drop that one as well.
Single term deletions
Model:
Sun ~ Long + MnJanTemp + Rain + Sea + NorthIsland
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 765366 370.73
Long 1 129226 894592 374.34 5.0653 0.03190 *
MnJanTemp 1 91939 857305 372.81 3.6037 0.06731 .
Rain 1 102527 867893 373.25 4.0188 0.05409 .
Sea 1 54835 820201 371.22 2.1494 0.15303
NorthIsland 1 106633 871999 373.42 4.1797 0.04977 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now we are told that Sea has p-value= 0.15, so we drop this as well.
Single term deletions
Model:
Sun ~ Long + MnJanTemp + Rain + NorthIsland
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 820201 371.22
Long 1 124840 945040 374.32 4.7184 0.037609 *
MnJanTemp 1 209006 1029207 377.39 7.8995 0.008493 **
Rain 1 74019 894219 372.33 2.7976 0.104471
NorthIsland 1 165899 986099 375.85 6.2702 0.017751 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally we drop Rain, which has p-value= 0.10. Then we stop.
Single term deletions
Model:
Sun ~ Long + MnJanTemp + NorthIsland
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 894219 372.33
Long 1 159187 1053407 376.22 5.6966 0.023079 *
MnJanTemp 1 297291 1191511 380.66 10.6387 0.002633 **
NorthIsland 1 222029 1116248 378.31 7.9454 0.008202 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = Sun ~ Long + MnJanTemp + NorthIsland, data = climate)
Residuals:
Min 1Q Median 3Q Max
-349.06 -90.02 -15.92 98.84 392.75
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7504.65 3290.81 -2.280 0.02938 *
Long 47.33 19.83 2.387 0.02308 *
MnJanTemp 85.06 26.08 3.262 0.00263 **
NorthIsland -296.09 105.04 -2.819 0.00820 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 167.2 on 32 degrees of freedom
Multiple R-squared: 0.4535, Adjusted R-squared: 0.4023
F-statistic: 8.853 on 3 and 32 DF, p-value: 0.0002034We see that annual hours of Sunshine are greater in places with high Longitude (i.e. eastern side of NZ), high mean January temperature and, perhaps surprisingly that hours are apparently lower in the North Island (but note NorthIsland and Longitude are correlated so this “apparently” may not be correct. )
31.2.0.2 Forward selection
We try using forward selection instead. Do we finish with the same model?
climate.null = lm(Sun ~ 1, data = climate)
add1(climate.null, scope = ~Lat + Long + MnJanTemp + MnJlyTemp + Height + Rain +
Sea + NorthIsland, data = climate, test = "F")Single term additions
Model:
Sun ~ 1
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 1636408 388.08
Lat 1 220053 1416355 384.88 5.2824 0.027817 *
Long 1 344588 1291820 381.57 9.0694 0.004877 **
MnJanTemp 1 506057 1130350 376.76 15.2218 0.000429 ***
MnJlyTemp 1 182734 1453674 385.82 4.2740 0.046380 *
Height 1 90941 1545467 388.02 2.0007 0.166321
Rain 1 316258 1320150 382.35 8.1451 0.007302 **
Sea 1 148252 1488156 386.66 3.3871 0.074450 .
NorthIsland 1 95788 1540620 387.91 2.1140 0.155131
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1climate.s1 = update(climate.null, . ~ . + MnJanTemp)
add1(climate.s1, scope = ~Lat + Long + MnJanTemp + MnJlyTemp + Height + Rain + Sea +
NorthIsland, data = climate, test = "F")Single term additions
Model:
Sun ~ MnJanTemp
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 1130350 376.76
Lat 1 43712 1086639 377.34 1.3275 0.25753
Long 1 14102 1116248 378.31 0.4169 0.52295
MnJlyTemp 1 22898 1107452 378.03 0.6823 0.41472
Height 1 40 1130310 378.76 0.0012 0.97294
Rain 1 132170 998181 374.29 4.3696 0.04437 *
Sea 1 50284 1080066 377.12 1.5364 0.22390
NorthIsland 1 76944 1053407 376.22 2.4104 0.13007
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1climate.s2 = update(climate.s1, . ~ . + Rain)
add1(climate.s2, scope = ~Lat + Long + MnJanTemp + MnJlyTemp + Height + Rain + Sea +
NorthIsland, data = climate, test = "F")Single term additions
Model:
Sun ~ MnJanTemp + Rain
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 998181 374.29
Lat 1 330 997850 376.27 0.0106 0.91869
Long 1 12081 986099 375.85 0.3921 0.53566
MnJlyTemp 1 16 998165 376.29 0.0005 0.98234
Height 1 1193 996988 376.24 0.0383 0.84612
Sea 1 87195 910985 373.00 3.0629 0.08968 .
NorthIsland 1 53140 945040 374.32 1.7994 0.18923
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = Sun ~ MnJanTemp + Rain, data = climate)
Residuals:
Min 1Q Median 3Q Max
-353.07 -84.98 -30.65 86.74 383.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1032.83976 343.96323 3.003 0.00507 **
MnJanTemp 62.40295 19.12695 3.263 0.00257 **
Rain -0.08159 0.03903 -2.090 0.04437 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 173.9 on 33 degrees of freedom
Multiple R-squared: 0.39, Adjusted R-squared: 0.353
F-statistic: 10.55 on 2 and 33 DF, p-value: 0.0002869You should find that the model chosen by forward selection is entirely different to the model by backwards elimination.
Usually the models are the same, but in this example they are not the same.
The reason for the difference is that sometimes two variables are both significant when together in a regression, but are not significant enough to enter one-at-a-time in forwards selection. That is why backwards elimination is usually the better approach.
31.3 Model Selection using AIC
An alternative approach to selecting a good model is to define a measure of the quality of a model, then choose the model of the highest quality.
The so called information criteria can be used as measures of model quality.
We will now examine information criteria in general, the Akaike information criterion (AIC) in particular, and consider their use in model selection.
31.4 Information Criteria
A good regression model (i.e. with a good choice of predictors) will:
- fit the data well;
- be simple (have as few predictors as possible).
Hence quality of a model can be measured by a quantity of the form \[\mbox{Measure of badness of fit} + k \times \mbox{Number of predictors} + \mbox{constant}\]
where k can be adjusted to control the relative importance of the two aspects of a model’s quality just given.
Note that any additive constant is of no real importance, since in comparing the quality of different models (our real aim) the constant term cancels out.
31.5 Akaike Information Criterion
Several information criteria in statistics are of the form of our simple equation for a “measure of badness of fit” given above.
One such measure is the Akaike Information Criterion (AIC). For a linear model with unknown error variance this is defined by \[AIC = n \log (RSS/n) + 2p + \mbox{constant}.\]
We can choose between models by selecting the one with smaller AIC.
N.B. This is not the only criterion in common use, but is perhaps the most common. Different software present different criteria.
31.6 Back to the Climate data
We have found two different models, Climate.lm6 and Climate.s2. These not only have different numbers of variables, but are not even nested! However we can compare them using the adjusted R-square, or the function AIC(). Both criteria suggest that the model lm6 is to be preferred.
[1] 0.4023174[1] 0.3530486[1] 476.4903[1] 478.4497An alternative to AIC() is the extractAIC() command. The output from extractAIC() command is number of regression
parameters (p+1) and the model’s AIC.
[1] 4.0000 372.3268[1] 3.0000 374.286131.7 AIC in Stepwise Variable Selection
We have seen that sequential variable selection can be conducted using a sequence of F tests.
An alternative is to compare models using AIC.
Specifically, at each step of the algorithm we change to the model with smallest AIC (from adding or dropping a single variable), and stop only when we cannot reduce AIC by single term additions or deletions.
This procedure is implemented using the step() command in R.
climate.null = lm(Sun ~ 1, data = climate)
climate.step <- step(climate.null, scope = ~Lat + Long + MnJanTemp + MnJlyTemp +
Height + Rain + Sea + NorthIsland, data = climate, direction = "both")Start: AIC=388.08
Sun ~ 1
Df Sum of Sq RSS AIC
+ MnJanTemp 1 506057 1130350 376.76
+ Long 1 344588 1291820 381.57
+ Rain 1 316258 1320150 382.35
+ Lat 1 220053 1416355 384.88
+ MnJlyTemp 1 182734 1453674 385.82
+ Sea 1 148252 1488156 386.66
+ NorthIsland 1 95788 1540620 387.91
+ Height 1 90941 1545467 388.02
<none> 1636408 388.08
Step: AIC=376.76
Sun ~ MnJanTemp
Df Sum of Sq RSS AIC
+ Rain 1 132170 998181 374.29
+ NorthIsland 1 76944 1053407 376.22
<none> 1130350 376.76
+ Sea 1 50284 1080066 377.12
+ Lat 1 43712 1086639 377.34
+ MnJlyTemp 1 22898 1107452 378.03
+ Long 1 14102 1116248 378.31
+ Height 1 40 1130310 378.76
- MnJanTemp 1 506057 1636408 388.08
Step: AIC=374.29
Sun ~ MnJanTemp + Rain
Df Sum of Sq RSS AIC
+ Sea 1 87195 910985 373.00
<none> 998181 374.29
+ NorthIsland 1 53140 945040 374.32
+ Long 1 12081 986099 375.85
+ Height 1 1193 996988 376.24
+ Lat 1 330 997850 376.27
+ MnJlyTemp 1 16 998165 376.29
- Rain 1 132170 1130350 376.76
- MnJanTemp 1 321969 1320150 382.35
Step: AIC=373
Sun ~ MnJanTemp + Rain + Sea
Df Sum of Sq RSS AIC
+ Height 1 94788 816197 371.04
+ MnJlyTemp 1 52902 858083 372.84
<none> 910985 373.00
+ Long 1 38986 871999 373.42
- Sea 1 87195 998181 374.29
+ NorthIsland 1 16393 894592 374.34
+ Lat 1 1945 909040 374.92
- Rain 1 169081 1080066 377.12
- MnJanTemp 1 214016 1125001 378.59
Step: AIC=371.04
Sun ~ MnJanTemp + Rain + Sea + Height
Df Sum of Sq RSS AIC
<none> 816197 371.04
+ NorthIsland 1 28556 787641 371.76
+ MnJlyTemp 1 22773 793424 372.02
+ Long 1 19420 796777 372.17
+ Lat 1 10788 805409 372.56
- Height 1 94788 910985 373.00
- Sea 1 180791 996988 376.24
- Rain 1 226572 1042769 377.86
- MnJanTemp 1 284136 1100333 379.79
Call:
lm(formula = Sun ~ MnJanTemp + Rain + Sea + Height, data = climate)
Residuals:
Min 1Q Median 3Q Max
-268.06 -104.31 1.62 90.67 344.42
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 872.1199 356.7010 2.445 0.02037 *
MnJanTemp 63.9085 19.4541 3.285 0.00253 **
Rain -0.1123 0.0383 -2.934 0.00625 **
Sea 203.2272 77.5551 2.620 0.01348 *
Height 0.3860 0.2034 1.897 0.06712 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 162.3 on 31 degrees of freedom
Multiple R-squared: 0.5012, Adjusted R-squared: 0.4369
F-statistic: 7.788 on 4 and 31 DF, p-value: 0.0001822Note AIC is more liberal than the F test, and has allowed a non-significant variable in.
The syntax below shows how to do a backwards elimination. This also results in including a non-significant variable. However this has higher adjusted R-square.
climate.full = lm(Sun ~ Lat + Long + MnJanTemp + MnJlyTemp + Height + Rain + Sea +
NorthIsland, data = climate)
climate.step2 <- step(climate.full, scope = ~1, data = climate, direction = "backward")Start: AIC=372.48
Sun ~ Lat + Long + MnJanTemp + MnJlyTemp + Height + Rain + Sea +
NorthIsland
Df Sum of Sq RSS AIC
- Lat 1 760 680941 370.52
- MnJlyTemp 1 9027 689208 370.95
- Height 1 21370 701551 371.59
<none> 680181 372.48
- NorthIsland 1 47272 727453 372.90
- MnJanTemp 1 62239 742420 373.63
- Rain 1 79398 759579 374.45
- Long 1 104127 784308 375.61
- Sea 1 112082 792263 375.97
Step: AIC=370.52
Sun ~ Long + MnJanTemp + MnJlyTemp + Height + Rain + Sea + NorthIsland
Df Sum of Sq RSS AIC
- MnJlyTemp 1 8381 689322 368.96
<none> 680941 370.52
- Height 1 39668 720609 370.56
- NorthIsland 1 48846 729788 371.01
- Rain 1 86487 767428 372.82
- Long 1 104518 785460 373.66
- Sea 1 114410 795351 374.11
- MnJanTemp 1 147165 828106 375.56
Step: AIC=368.96
Sun ~ Long + MnJanTemp + Height + Rain + Sea + NorthIsland
Df Sum of Sq RSS AIC
<none> 689322 368.96
- Height 1 76044 765366 370.73
- Long 1 98319 787641 371.76
- NorthIsland 1 107455 796777 372.17
- Sea 1 125566 814888 372.98
- MnJanTemp 1 140484 829807 373.64
- Rain 1 142964 832286 373.74
Call:
lm(formula = Sun ~ Long + MnJanTemp + Height + Rain + Sea + NorthIsland,
data = climate)
Residuals:
Min 1Q Median 3Q Max
-245.527 -94.201 -9.587 98.988 253.293
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.743e+03 3.114e+03 -1.844 0.0754 .
Long 3.808e+01 1.872e+01 2.034 0.0512 .
MnJanTemp 7.138e+01 2.936e+01 2.431 0.0215 *
Height 3.527e-01 1.972e-01 1.789 0.0841 .
Rain -9.205e-02 3.753e-02 -2.452 0.0204 *
Sea 1.740e+02 7.571e+01 2.298 0.0289 *
NorthIsland -2.188e+02 1.029e+02 -2.126 0.0421 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 154.2 on 29 degrees of freedom
Multiple R-squared: 0.5788, Adjusted R-squared: 0.4916
F-statistic: 6.641 on 6 and 29 DF, p-value: 0.0001691Question: Which of these models is preferred?
31.7.1 Epilogue (can be omitted)
If I wanted to explore further, I would investigate possible interactions.
It is also worth thinking about whether the variables make sense. Is there a logical order for the variables? Here Lat, Long, Height, Sea and NorthIsland are all geographic variables, which may determine the weather.
On the other hand Sun, MnJanTemp, MnJlyTemp and Rain are in a sense all interconnected response weather variables, even though we are using Sun as the sole response for the regression.
So it makes logical sense to first let the geographic variables explain as much of the variation in Sun hours as possible.
climate.null = lm(Sun ~ 1, data = climate)
climate.geog <- step(climate.null, scope = ~Lat * Long * Height * Sea * NorthIsland,
data = climate, direction = "both")Start: AIC=388.08
Sun ~ 1
Df Sum of Sq RSS AIC
+ Long 1 344588 1291820 381.57
+ Lat 1 220053 1416355 384.88
+ Sea 1 148252 1488156 386.66
+ NorthIsland 1 95788 1540620 387.91
+ Height 1 90941 1545467 388.02
<none> 1636408 388.08
Step: AIC=381.57
Sun ~ Long
Df Sum of Sq RSS AIC
+ Sea 1 157638 1134182 378.88
+ NorthIsland 1 100309 1191511 380.66
<none> 1291820 381.57
+ Height 1 59088 1232732 381.88
+ Lat 1 1618 1290202 383.52
- Long 1 344588 1636408 388.08
Step: AIC=378.88
Sun ~ Long + Sea
Df Sum of Sq RSS AIC
+ Long:Sea 1 241297 892885 372.27
+ NorthIsland 1 78624 1055558 378.30
<none> 1134182 378.88
+ Lat 1 3989 1130193 380.76
+ Height 1 782 1133400 380.86
- Sea 1 157638 1291820 381.57
- Long 1 353974 1488156 386.66
Step: AIC=372.27
Sun ~ Long + Sea + Long:Sea
Df Sum of Sq RSS AIC
+ NorthIsland 1 83808 809077 370.72
<none> 892885 372.27
+ Lat 1 8100 884785 373.94
+ Height 1 5468 887417 374.05
- Long:Sea 1 241297 1134182 378.88
Step: AIC=370.72
Sun ~ Long + Sea + NorthIsland + Long:Sea
Df Sum of Sq RSS AIC
<none> 809077 370.72
+ Long:NorthIsland 1 25954 783124 371.55
+ Lat 1 13602 795476 372.11
+ Height 1 10363 798715 372.26
- NorthIsland 1 83808 892885 372.27
+ Sea:NorthIsland 1 12 809065 372.72
- Long:Sea 1 246480 1055558 378.30
Call:
lm(formula = Sun ~ Long + Sea + NorthIsland + Long:Sea, data = climate)
Residuals:
Min 1Q Median 3Q Max
-360.74 -98.25 -24.93 91.57 315.67
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3895.08 3559.07 -1.094 0.28221
Long 34.00 20.76 1.638 0.11158
Sea -10765.30 3544.72 -3.037 0.00481 **
NorthIsland -173.62 96.89 -1.792 0.08291 .
Long:Sea 62.85 20.45 3.073 0.00439 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 161.6 on 31 degrees of freedom
Multiple R-squared: 0.5056, Adjusted R-squared: 0.4418
F-statistic: 7.925 on 4 and 31 DF, p-value: 0.0001603
In the plot, green locations are for the South Island and blue for the North Island, triangles represent inland places and circles are places by the sea. The plot indicates, that, for coastal towns/cities, locations on the eastern side of each island (higher longitude) tend to have more sun than locations on the western side. But for inland places the longitude makes little difference.
Next we can try entering the other weather variables.
climate.geog2 <- step(climate.geog, scope = ~Lat * Long * Height * Sea * NorthIsland *
MnJanTemp * MnJlyTemp * Rain, data = climate, direction = "both")Start: AIC=370.72
Sun ~ Long + Sea + NorthIsland + Long:Sea
Df Sum of Sq RSS AIC
+ Rain 1 124430 684647 366.71
+ MnJanTemp 1 72184 736893 369.36
<none> 809077 370.72
+ Long:NorthIsland 1 25954 783124 371.55
+ MnJlyTemp 1 18553 790524 371.89
+ Lat 1 13602 795476 372.11
+ Height 1 10363 798715 372.26
- NorthIsland 1 83808 892885 372.27
+ Sea:NorthIsland 1 12 809065 372.72
- Long:Sea 1 246480 1055558 378.30
Step: AIC=366.71
Sun ~ Long + Sea + NorthIsland + Rain + Long:Sea
Df Sum of Sq RSS AIC
+ Sea:Rain 1 86953 597694 363.82
+ Lat 1 58802 625845 365.48
<none> 684647 366.71
+ MnJanTemp 1 31200 653447 367.03
+ Long:Rain 1 30603 654044 367.07
- NorthIsland 1 56519 741166 367.57
+ Long:NorthIsland 1 20092 664555 367.64
+ MnJlyTemp 1 5540 679107 368.42
+ Sea:NorthIsland 1 2819 681828 368.56
+ Height 1 1879 682768 368.61
+ NorthIsland:Rain 1 362 684284 368.69
- Rain 1 124430 809077 370.72
- Long:Sea 1 172658 857305 372.81
Step: AIC=363.82
Sun ~ Long + Sea + NorthIsland + Rain + Long:Sea + Sea:Rain
Df Sum of Sq RSS AIC
+ MnJlyTemp 1 51554 546140 362.58
+ Lat 1 42961 554734 363.14
+ Long:Rain 1 41679 556015 363.22
+ Long:NorthIsland 1 33248 564447 363.76
<none> 597694 363.82
+ Height 1 27761 569934 364.11
+ MnJanTemp 1 13844 583850 364.98
+ Sea:NorthIsland 1 2753 594941 365.66
+ NorthIsland:Rain 1 174 597521 365.81
- NorthIsland 1 85351 683046 366.63
- Sea:Rain 1 86953 684647 366.71
- Long:Sea 1 226771 824465 373.40
Step: AIC=362.58
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Long:Sea +
Sea:Rain
Df Sum of Sq RSS AIC
+ Lat 1 89204 456937 358.16
+ MnJanTemp 1 53611 492529 360.86
- NorthIsland 1 9068 555209 361.17
+ Long:NorthIsland 1 42456 503684 361.66
+ Long:Rain 1 30258 515882 362.52
<none> 546140 362.58
+ NorthIsland:MnJlyTemp 1 24010 522131 362.96
- MnJlyTemp 1 51554 597694 363.82
+ Sea:MnJlyTemp 1 7980 538160 364.05
+ Long:MnJlyTemp 1 4503 541637 364.28
+ Height 1 3647 542494 364.33
+ NorthIsland:Rain 1 559 545581 364.54
+ MnJlyTemp:Rain 1 449 545692 364.55
+ Sea:NorthIsland 1 437 545703 364.55
- Sea:Rain 1 132966 679107 368.42
- Long:Sea 1 257769 803909 374.49
Step: AIC=358.16
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain
Df Sum of Sq RSS AIC
+ Lat:Rain 1 121655 335282 349.01
+ Long:Rain 1 99082 357855 351.36
+ Lat:NorthIsland 1 48651 408286 356.10
+ MnJlyTemp:Rain 1 31897 425040 357.55
+ NorthIsland:Rain 1 30746 426191 357.65
<none> 456937 358.16
- NorthIsland 1 41886 498823 359.31
+ NorthIsland:MnJlyTemp 1 8794 448143 359.46
+ Long:MnJlyTemp 1 4328 452609 359.81
+ MnJanTemp 1 3726 453210 359.86
+ Height 1 3423 453514 359.89
+ Sea:NorthIsland 1 3149 453788 359.91
+ Lat:Sea 1 1593 455344 360.03
+ Lat:MnJlyTemp 1 1264 455673 360.06
+ Lat:Long 1 887 456050 360.09
+ Sea:MnJlyTemp 1 276 456660 360.13
+ Long:NorthIsland 1 78 456859 360.15
- Lat 1 89204 546140 362.58
- MnJlyTemp 1 97797 554734 363.14
- Sea:Rain 1 140172 597109 365.79
- Long:Sea 1 225668 682605 370.61
Step: AIC=349.01
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + Rain:Lat
Df Sum of Sq RSS AIC
+ Lat:NorthIsland 1 19311 315971 348.88
<none> 335282 349.01
+ NorthIsland:Rain 1 16737 318545 349.17
+ NorthIsland:MnJlyTemp 1 6120 329161 350.35
+ MnJlyTemp:Rain 1 5558 329724 350.41
+ Long:MnJlyTemp 1 5108 330173 350.46
+ Sea:NorthIsland 1 4915 330366 350.48
+ Lat:Sea 1 3721 331561 350.61
+ Long:Rain 1 1416 333865 350.86
+ MnJanTemp 1 1347 333934 350.87
+ Height 1 1177 334105 350.88
+ Long:NorthIsland 1 324 334957 350.98
+ Sea:MnJlyTemp 1 12 335270 351.01
+ Lat:MnJlyTemp 1 4 335278 351.01
+ Lat:Long 1 3 335279 351.01
- NorthIsland 1 40442 375723 351.11
- MnJlyTemp 1 91599 426880 355.71
- Rain:Lat 1 121655 456937 358.16
- Sea:Rain 1 170192 505474 361.79
- Long:Sea 1 207963 543244 364.38
Step: AIC=348.88
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + Rain:Lat + NorthIsland:Lat
Df Sum of Sq RSS AIC
+ Lat:Long 1 28919 287052 347.42
+ Long:MnJlyTemp 1 25040 290931 347.90
+ NorthIsland:MnJlyTemp 1 17143 298828 348.87
<none> 315971 348.88
+ NorthIsland:Rain 1 16377 299593 348.96
- NorthIsland:Lat 1 19311 335282 349.01
+ Long:NorthIsland 1 14384 301586 349.20
+ Lat:MnJlyTemp 1 7081 308889 350.06
+ Long:Rain 1 6474 309497 350.13
+ Lat:Sea 1 5732 310239 350.22
+ MnJlyTemp:Rain 1 4737 311233 350.33
+ Sea:NorthIsland 1 2318 313652 350.61
+ MnJanTemp 1 2256 313714 350.62
+ Height 1 570 315401 350.81
+ Sea:MnJlyTemp 1 537 315433 350.81
- MnJlyTemp 1 62420 378391 353.37
- Rain:Lat 1 92315 408286 356.10
- Sea:Rain 1 98373 414344 356.63
- Long:Sea 1 201579 517550 364.64
Step: AIC=347.42
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + Rain:Lat + NorthIsland:Lat + Long:Lat
Df Sum of Sq RSS AIC
+ Long:Rain 1 30879 256173 345.32
+ NorthIsland:MnJlyTemp 1 25372 261680 346.09
<none> 287052 347.42
+ Lat:MnJlyTemp 1 13974 273078 347.62
+ NorthIsland:Rain 1 11041 276010 348.01
+ Long:MnJlyTemp 1 10232 276819 348.11
- Long:Lat 1 28919 315971 348.88
+ Sea:MnJlyTemp 1 3840 283211 348.94
+ MnJlyTemp:Rain 1 3553 283498 348.97
+ MnJanTemp 1 3287 283765 349.01
+ Lat:Sea 1 450 286601 349.36
+ Long:NorthIsland 1 251 286800 349.39
+ Sea:NorthIsland 1 10 287042 349.42
+ Height 1 4 287048 349.42
- NorthIsland:Lat 1 48227 335279 351.01
- Sea:Rain 1 59365 346417 352.19
- MnJlyTemp 1 65830 352882 352.85
- Rain:Lat 1 72159 359211 353.49
- Long:Sea 1 154025 441077 360.88
Step: AIC=345.32
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + Rain:Lat + NorthIsland:Lat + Long:Lat + Long:Rain
Df Sum of Sq RSS AIC
- Rain:Lat 1 1044 257217 343.47
+ NorthIsland:MnJlyTemp 1 24780 231393 343.66
<none> 256173 345.32
+ NorthIsland:Rain 1 13012 243161 345.45
+ Long:MnJlyTemp 1 8636 247537 346.09
+ Lat:Long:Rain 1 8232 247941 346.15
+ Lat:MnJlyTemp 1 7144 249029 346.31
+ Sea:MnJlyTemp 1 6434 249739 346.41
+ Lat:Sea 1 5367 250806 346.56
+ MnJlyTemp:Rain 1 4246 251927 346.72
- Sea:Rain 1 25770 281943 346.77
+ Height 1 2056 254117 347.03
+ Long:NorthIsland 1 1883 254290 347.06
+ Sea:NorthIsland 1 804 255369 347.21
+ Long:Sea:Rain 1 244 255929 347.29
+ MnJanTemp 1 19 256154 347.32
- Long:Rain 1 30879 287052 347.42
- Long:Lat 1 53324 309497 350.13
- MnJlyTemp 1 62295 318468 351.16
- NorthIsland:Lat 1 77692 333865 352.86
- Long:Sea 1 150162 406335 359.93
Step: AIC=343.47
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + NorthIsland:Lat + Long:Lat + Long:Rain
Df Sum of Sq RSS AIC
+ NorthIsland:MnJlyTemp 1 24039 233178 341.94
<none> 257217 343.47
+ NorthIsland:Rain 1 13686 243531 343.50
+ Long:MnJlyTemp 1 8743 248475 344.22
+ Lat:MnJlyTemp 1 7800 249418 344.36
+ Sea:MnJlyTemp 1 5783 251434 344.65
+ MnJlyTemp:Rain 1 5240 251977 344.73
+ Lat:Sea 1 4219 252998 344.87
+ Long:NorthIsland 1 1809 255409 345.22
+ Height 1 1429 255789 345.27
+ Lat:Rain 1 1044 256173 345.32
+ Sea:NorthIsland 1 430 256787 345.41
+ Long:Sea:Rain 1 145 257072 345.45
+ MnJanTemp 1 22 257195 345.47
- Sea:Rain 1 37755 294972 346.40
- Long:Lat 1 59878 317096 349.00
- MnJlyTemp 1 63286 320503 349.39
- NorthIsland:Lat 1 100638 357855 353.36
- Long:Rain 1 101994 359211 353.49
- Long:Sea 1 150964 408182 358.09
Step: AIC=341.94
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + NorthIsland:Lat + Long:Lat + Long:Rain + NorthIsland:MnJlyTemp
Df Sum of Sq RSS AIC
+ NorthIsland:Rain 1 12877 220301 341.89
<none> 233178 341.94
+ Lat:Sea 1 7745 225432 342.72
+ Long:NorthIsland 1 6997 226181 342.84
+ Sea:NorthIsland 1 6825 226353 342.87
+ Lat:MnJlyTemp 1 5342 227835 343.10
+ MnJlyTemp:Rain 1 4096 229082 343.30
+ Long:MnJlyTemp 1 4077 229101 343.30
+ Long:Sea:Rain 1 3978 229200 343.32
- NorthIsland:MnJlyTemp 1 24039 257217 343.47
+ Lat:Rain 1 1785 231393 343.66
- Long:Sea 1 27178 260356 343.91
+ Sea:MnJlyTemp 1 192 232986 343.91
+ MnJanTemp 1 90 233088 343.92
+ Height 1 30 233147 343.93
- Sea:Rain 1 49345 282523 346.85
- Long:Lat 1 68319 301497 349.19
- Long:Rain 1 86253 319431 351.27
- NorthIsland:Lat 1 121606 354784 355.05
Step: AIC=341.89
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + NorthIsland:Lat + Long:Lat + Long:Rain + NorthIsland:MnJlyTemp +
NorthIsland:Rain
Df Sum of Sq RSS AIC
<none> 220301 341.89
- NorthIsland:Rain 1 12877 233178 341.94
+ Sea:NorthIsland 1 4887 215413 343.08
+ Lat:Sea 1 4881 215420 343.09
- NorthIsland:MnJlyTemp 1 23230 243531 343.50
+ Lat:MnJlyTemp 1 2305 217995 343.51
+ Long:MnJlyTemp 1 1710 218591 343.61
+ Height 1 513 219788 343.81
+ Long:NorthIsland 1 448 219853 343.82
+ MnJlyTemp:Rain 1 410 219891 343.83
+ Long:Sea:Rain 1 409 219892 343.83
+ MnJanTemp 1 325 219976 343.84
+ Sea:MnJlyTemp 1 82 220218 343.88
+ Lat:Rain 1 48 220253 343.88
- Long:Sea 1 30383 250683 344.54
- Sea:Rain 1 50609 270910 347.34
- Long:Lat 1 74420 294721 350.37
- Long:Rain 1 86103 306403 351.77
- NorthIsland:Lat 1 132580 352880 356.85
Call:
lm(formula = Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp +
Lat + Long:Sea + Sea:Rain + NorthIsland:Lat + Long:Lat +
Long:Rain + NorthIsland:MnJlyTemp + NorthIsland:Rain, data = climate)
Residuals:
Min 1Q Median 3Q Max
-126.31 -70.88 -10.95 67.42 131.41
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.166e+05 4.193e+04 -2.781 0.01089 *
Long 7.410e+02 2.511e+02 2.950 0.00740 **
Sea -5.978e+03 3.484e+03 -1.716 0.10025
NorthIsland -9.466e+03 2.522e+03 -3.754 0.00110 **
Rain 1.136e+01 3.945e+00 2.880 0.00870 **
MnJlyTemp -6.906e+01 2.497e+01 -2.766 0.01128 *
Lat 2.689e+03 1.025e+03 2.623 0.01552 *
Long:Sea 3.489e+01 2.003e+01 1.742 0.09550 .
Sea:Rain 1.340e-01 5.961e-02 2.248 0.03492 *
NorthIsland:Lat 2.128e+02 5.848e+01 3.639 0.00145 **
Long:Lat -1.674e+01 6.140e+00 -2.726 0.01233 *
Long:Rain -6.801e-02 2.319e-02 -2.932 0.00771 **
NorthIsland:MnJlyTemp 4.444e+01 2.917e+01 1.523 0.14198
NorthIsland:Rain 1.949e-01 1.719e-01 1.134 0.26900
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 100.1 on 22 degrees of freedom
Multiple R-squared: 0.8654, Adjusted R-squared: 0.7858
F-statistic: 10.88 on 13 and 22 DF, p-value: 9.523e-07The AIC has let in some interactions that are non-significant. We will weed them out based on F tests.
Single term deletions
Model:
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + NorthIsland:Lat + Long:Lat + Long:Rain + NorthIsland:MnJlyTemp
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 233178 341.94
Long:Sea 1 27178 260356 343.91 2.6808 0.115176
Sea:Rain 1 49345 282523 346.85 4.8672 0.037631 *
NorthIsland:Lat 1 121606 354784 355.05 11.9949 0.002108 **
Long:Lat 1 68319 301497 349.19 6.7388 0.016157 *
Long:Rain 1 86253 319431 351.27 8.5078 0.007762 **
NorthIsland:MnJlyTemp 1 24039 257217 343.47 2.3712 0.137239
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1climate.geog4 = update(climate.geog3, . ~ . - NorthIsland:MnJlyTemp)
drop1(climate.geog4, test = "F")Single term deletions
Model:
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
Sea:Rain + NorthIsland:Lat + Long:Lat + Long:Rain
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 257217 343.47
MnJlyTemp 1 63286 320503 349.39 5.9050 0.0229471 *
Long:Sea 1 150964 408182 358.09 14.0859 0.0009808 ***
Sea:Rain 1 37755 294972 346.40 3.5228 0.0727414 .
NorthIsland:Lat 1 100638 357855 353.36 9.3902 0.0053234 **
Long:Lat 1 59878 317096 349.00 5.5870 0.0265311 *
Long:Rain 1 101994 359211 353.49 9.5167 0.0050676 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Single term deletions
Model:
Sun ~ Long + Sea + NorthIsland + Rain + MnJlyTemp + Lat + Long:Sea +
NorthIsland:Lat + Long:Lat + Long:Rain
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 294972 346.40
MnJlyTemp 1 33639 328611 348.29 2.8510 0.1037531
Long:Sea 1 117556 412527 356.48 9.9633 0.0041334 **
NorthIsland:Lat 1 210085 505057 363.76 17.8055 0.0002814 ***
Long:Lat 1 97397 392369 354.67 8.2547 0.0081728 **
Long:Rain 1 97456 392428 354.68 8.2598 0.0081559 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Single term deletions
Model:
Sun ~ Long + Sea + NorthIsland + Rain + Lat + Long:Sea + NorthIsland:Lat +
Long:Lat + Long:Rain
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 328611 348.29
Long:Sea 1 122977 451588 357.73 9.7301 0.0043970 **
NorthIsland:Lat 1 197823 526434 363.25 15.6520 0.0005239 ***
Long:Lat 1 78526 407136 354.00 6.2130 0.0193795 *
Long:Rain 1 97817 426428 355.67 7.7394 0.0099233 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = Sun ~ Long + Sea + NorthIsland + Rain + Lat + Long:Sea +
NorthIsland:Lat + Long:Lat + Long:Rain, data = climate)
Residuals:
Min 1Q Median 3Q Max
-165.59 -78.14 -14.90 85.53 199.25
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.013e+05 4.356e+04 -2.325 0.028140 *
Long 6.511e+02 2.600e+02 2.504 0.018887 *
Sea -7.948e+03 2.579e+03 -3.082 0.004824 **
NorthIsland -9.092e+03 2.247e+03 -4.047 0.000414 ***
Rain 8.846e+00 3.239e+00 2.732 0.011174 *
Lat 2.541e+03 1.069e+03 2.376 0.025140 *
Long:Sea 4.642e+01 1.488e+01 3.119 0.004397 **
NorthIsland:Lat 2.156e+02 5.451e+01 3.956 0.000524 ***
Long:Lat -1.592e+01 6.386e+00 -2.493 0.019379 *
Long:Rain -5.295e-02 1.903e-02 -2.782 0.009923 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 112.4 on 26 degrees of freedom
Multiple R-squared: 0.7992, Adjusted R-squared: 0.7297
F-statistic: 11.5 on 9 and 26 DF, p-value: 4.964e-07The above is the final model. It has a high \(R^2\) but is probably a bit too complicated for such a small dataset. A further step would be to weed further variables out using \(R^2_{pred}\).