So, an expected pattern of results emerges for the non-dyslexic students. Those with higher IQs tend to have higher reading accuracy; this does not mean higher IQ’s cause more accurate reading because random assignment of IQ is not possible. However, for the dyslexic students, the relationship is not what one would might expect. It is slightly negative, showing that higher IQ’s are related to lower reading accuracy. What we conclude from this is that we should not expect higher IQ’s to show higher performance on a test like this.
Checking the assumptions is always recommended before getting focused on the inferences in the model. When fitting models with multiple groups, it is possible to see “groups” in the fitted values (x-axis in Residuals vs Fitted and Scale-Location plots) and that is not a problem – it is a feature of these models. You should look for issues in the residuals for each group but the residuals should overall still be normally distributed and have the same variability everywhere. It is a bit hard to see issues in Figure 8.27 because of the group differences, but note the line of residuals for the higher fitted values. This is an artifact of the upper threshold in the reading accuracy test used. As in the first year of college GPA, these observations were censored – their true score was outside the range of values we could observe – and so we did not really get a measure of how good these students were since a lot of their abilities were higher than the test could detect and they all binned up at the same value of getting all the questions correct. The relationship in this group might be even stronger if we could really observe differences in the highest level readers. We should treat the results for the non-dyslexic group with caution even though they are clearly scoring on average higher and have a different slope than the results for the dyslexic students. The normality and influence diagnostics do not suggest any major issues other than this.
par(mfrow=c(2,2), oma=c(0,0,2,0))
plot(dys_model,
sub.caption="Plot of diagnostics for Dyslexia Interaction model")
For these models, we have relaxed an earlier assumption that data were collected from only one group. In fact, we are doing specific research that is focused on questions about the differences between groups. However, these models still make assumptions that, within a specific group, the linearity assumption is met. They also assume that the variability in the residuals is the same for all observations. Sometimes it can be difficult to check the assumptions by looking at the overall diagnostic plots and it may be easier to go back to the original scatterplot or plot the residuals vs fitted values by group to fully assess the results. For example, Figure 8.28 shows a scatterplot of the residuals vs the quantitative explanatory variable by the groups. The variability in the residuals is a bit larger in the non-dyslexic group, possibly suggesting that variability in the reading test is higher for higher scoring individuals even though we couldn’t observe all of that variability because there were so many perfect scores in this group.
(ref:fig8-27) Plot of Residuals vs Fitted from interaction dyslexia data model with groups indicated.
scatterplot(residuals(dys_model)~fitted(dys_model)|dys,
data=dyslexic3, smooth=F)
If we feel comfortable enough with the assumptions to trust the inferences here (this might be dangerous), then we can consider what some of the model inferences provide us in this situation. For example, the test for \(H_0: {\color{red}{\boldsymbol{\beta_3}}}=0\) vs \(H_A: {\color{red}{\boldsymbol{\beta_3}}}\ne 0\) provides an interesting comparison. Under the null hypothesis, the two groups would have the same slope so it provides an opportunity to directly consider whether the relationship (via the slope) is different between the groups in their respective populations. We find \(t=-1.907\) which, if the assumptions are met, follows a \(t(40)\)-distribution under the null hypothesis. This test statistic has a corresponding p-value of 0.0638. So it provides some evidence of a difference in the slopes but it isn’t strong evidence of that. There are serious issues (like getting the wrong idea about directions of relationships) if we ignore a potentially important interaction and some statisticians would recommend retaining interactions even if the evidence is only moderate for its inclusion in the model. For the original research question of whether the relationships differ for the two groups, we only have marginal evidence to support that result. Possibly with a larger sample size or a reading test that only a few students could get 100% on, the researchers might have detected a more pronounced difference in the slopes for the two groups.
In the presence of a categorical by quantitative interaction, term-plots can be generated that plot the results for each group on the same display. This basically provides a plot of the “simplified” SLR models for each group. In Figure 8.29 we can see noticeable differences in the slopes and intercepts. Note that testing for differences in intercepts between groups is not very interesting when there are different slopes because if you change the slope, you have to change the intercept. The plot shows that there are clear differences in the means even though we don’t have a test to directly assess that in this complicated of a model111.
plot(allEffects(dys_model), ci.style="bands", multiline=T, lty=c(1,2))
It certainly appears in the plots that IQ has a different impact on the mean score in the two groups (even though the p-value only provided marginal evidence). To reinforce the potential dangers of forcing the same slope for both groups, consider the additive model for these data. Again, this just shifts one group off the other one, but both have the same slope. The following model summary and term-plots (Figure 8.30) suggest the potentially dangerous conclusion that can come from assuming a common slope when that might not be the case.
dys_modelR <- lm(score~ziq+dys, data=dyslexic3)
summary(dys_modelR)
##
## Call:
## lm(formula = score ~ ziq + dys, data = dyslexic3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26062 -0.05565 0.02932 0.07577 0.13217
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.89178 0.02312 38.580 < 2e-16
## ziq 0.02620 0.01957 1.339 0.188
## dys1 -0.26879 0.03905 -6.883 2.41e-08
##
## Residual standard error: 0.1049 on 41 degrees of freedom
## Multiple R-squared: 0.6858, Adjusted R-squared: 0.6705
## F-statistic: 44.75 on 2 and 41 DF, p-value: 4.917e-11
plot(allEffects(dys_modelR))
This model provides little to no evidence that IQ is related to reading score for all students (\(t_{41}=1.34\), p-value=0.188) but strong evidence of a difference in the y-intercepts (\(t_{41}=-6.88\), p-value \(<0.00001\)).
Since the IQ term has a large p-value, then we could drop it from the model – leaving a model that only includes the grouping variable:
dys_modelR2 <- lm(score~dys, data=dyslexic3)
summary(dys_modelR2)
##
## Call:
## lm(formula = score ~ dys, data = dyslexic3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25818 -0.04510 0.02514 0.09520 0.09694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.90480 0.02117 42.737 <2e-16
## dys1 -0.29892 0.03222 -9.278 1e-11
##
## Residual standard error: 0.1059 on 42 degrees of freedom
## Multiple R-squared: 0.6721, Adjusted R-squared: 0.6643
## F-statistic: 86.08 on 1 and 42 DF, p-value: 1e-11
plot(allEffects(dys_modelR2))
These results, including the term-plot in Figure 8.31, show that there is evidence of a difference in the mean reading scores between the two groups and maybe that is all these data really say… This is the logical outcome if we decide that the interaction is not important IN THIS DATA SET. In general, if the interaction is dropped, the interaction model can be reduced to considering an additive model with the categorical and quantitative predictor variables. Either or both of those variables could also be considered for removal, possibly starting with the variable with the larger p-value, leaving a string of ever-simpler models possible if large p-values are continually encountered112.
It is useful to note that the last model has returned us to the first model we encountered in Chapter 2 where we were just comparing the means for two groups. However, the researchers probably were not seeking to make the discovery that dyslexic students have a tougher time than non-dyslexic students on a reading test but sometimes that is all that the data support. The key part of this sequence of decisions was how much evidence you think a p-value of 0.06 contains…
For more than two categories in a categorical variable, the model contains more indicators to keep track of but uses the same ideas. We have to deal with modifying the intercept and slope coefficients for every deviation group so the task is onerous but relatively repetitive. The general model is:
\[\begin{array}{rl} y_i=\beta_0 &+ \beta_1x_i +\beta_2I_{\text{Level }2,i}+\beta_3I_{\text{Level }3,i} +\cdots+\beta_JI_{\text{Level }J,i} \\ &+\beta_{J+1}I_{\text{Level }2,i}\:x_i+\beta_{J+2}I_{\text{Level }3,i}\:x_i +\cdots+\beta_{2J-1}I_{\text{Level }J,i}\:x_i +\varepsilon_i.\ \end{array}\]
Specific to the audible tolerance/headache data that had four groups. The model with an interaction present is
\[\begin{array}{rl} \text{du2}_i = \beta_0 &+ \beta_1\cdot\text{du1}_i + \beta_2I_{T2,i} + \beta_3I_{T3,i} + \beta_4I_{\text{Control},i} \\ &+ \beta_5I_{T2,i}\cdot\text{du1}_i + \beta_6I_{T3,i}\cdot\text{du1}_i + \beta_7I_{\text{Control},i}\cdot\text{du1}_i+\varepsilon_i.\ \end{array}\]
Based on the following output, the estimated general regression model is
\[\begin{array}{rl} \widehat{\text{du2}}_i = 1.33 &+ 0.733\cdot\text{du1}_i - 0.236I_{T2,i} - 0.316I_{T3,i} - 1.091I_{\text{Control},i} \\ &+ 0.066I_{T2,i}\cdot\text{du1}_i + 0.199I_{T3,i}\cdot\text{du1}_i + 0.106I_{\text{Control},i}\cdot\text{du1}_i.\ \end{array}\]
Then we could work out the specific equation for each group with replacing their indicator variable in two places with 1s and the rest of the indicators with 0. For example, for the Control group:
\[\begin{array}{rll} \widehat{\text{du2}}_i &= 1.33 &+ 0.733\cdot\text{du1}_i - 0.236*0 - 0.316*0 - 1.091*1 \\ &&+ 0.066*0*\text{du1}_i + 0.199*0*\text{du1}_i+ 0.106*1*\text{du1}_i \\ \widehat{\text{du2}}_i&=1.33&+0.733\cdot\text{du1}_i - 1.091 + 0.106\cdot\text{du1}_i \\ \widehat{\text{du2}}_i&=0.239 &+ 0.839\cdot\text{du1}_i.\ \end{array}\]
head2 <- lm(du2~du1*treatment, data=Headache)
summary(head2)
##
## Call:
## lm(formula = du2 ~ du1 * treatment, data = Headache)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.8072 -1.0969 -0.3285 0.8192 10.6039
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.33157 0.66027 2.017 0.0467
## du1 0.73319 0.09969 7.355 8.53e-11
## treatmentT2 -0.23560 1.13414 -0.208 0.8359
## treatmentT3 -0.31613 0.95767 -0.330 0.7421
## treatmentControl -1.09084 0.95020 -1.148 0.2540
## du1:treatmentT2 0.06623 0.17473 0.379 0.7055
## du1:treatmentT3 0.19904 0.13350 1.491 0.1395
## du1:treatmentControl 0.10604 0.14326 0.740 0.4611
##
## Residual standard error: 2.148 on 90 degrees of freedom
## Multiple R-squared: 0.7573, Adjusted R-squared: 0.7384
## F-statistic: 40.12 on 7 and 90 DF, p-value: < 2.2e-16
Or we can let the term-plots (Figures 8.32 and
8.33) show us all four different simplified models. Here we
can see that all the slopes “look” to be pretty similar. When the interaction
model is fit and the results “look” like the additive model, there is a good
chance that we will be able to avoid all this complication and just use the
additive model without missing anything interesting.
There are two different
options for displaying interaction models. Version 1 (Figure
8.32) has a
different panel for each level of the categorical variable and Version 2
(Figure 8.33) puts all the lines on the same plot. In this
case, neither
version shows much of a difference and Version 2 overlaps so much that you
can’t see all the groups. In these situations, it can be useful to make the term-plots with multiline=T
and multiline=F
and select the version that captures the results best.
(ref:fig8-32) Term-plot for decibel tolerance interaction model (version 2). This plot is not printed in color because it is impossible to distinguish the four groups whether in color or black and white.
plot(allEffects(head2), x.var="du1", ci.style="bands")
plot(allEffects(head2), x.var="du1", multiline=T, ci.style="bands")
In situations with more than 2 levels, the \(t\)-tests for the interaction or changing y-intercepts are not informative for deciding if you really need different slopes or intercepts for all the groups. They only tell you if a specific group is potentially different from the baseline group and the choice of the baseline is arbitrary. To assess whether we really need to have varying slopes or intercepts with more than two groups we need to develop \(F\)-tests for the interaction part of the model.
8.12 F-tests for MLR models with quantitative and categorical variables and interactions
For models with multi-category \((J>2)\) categorical variables we need a method for deciding if all the extra complexity present in the additive or interaction models is necessary. We can appeal to model selection methods such as the adjusted R2 that focus on balancing model fit and complexity but interests often move to trying to decide if the differences are more extreme than we would expect by chance if there were no group differences in intercepts or slopes. Because of the multi-degree of freedom aspects of the use of indicator variables (\(J-1\) variables for a \(J\) level categorical variable), we have to develop tests that combine and assess information across multiple “variables” – even though these indicators all pertain to a single original categorical variable. ANOVA \(F\)-tests did exactly this sort of thing in the One and Two-Way ANOVA models and can do that for us here. There are two models that we perform tests in – the additive and the interaction models. We start with a discussion of the tests in an interaction setting since that provides us the first test to consider in most situations to assess evidence of whether the extra complexity of varying slopes is really needed. If we don’t “need” the varying slopes or if the plot really does have lines for the groups that look relatively parallel, we can fit the additive model and either assess evidence of the need for different intercepts or for the quantitative predictor – either is a reasonable next step. Basically this establishes a set of nested models (each model is a reduced version of another more complicated model higher in the tree of models) displayed in Figure 8.34. This is based on the assumption that we would proceed through the model, dropping terms if the p-values are large (“not significant” in the diagram) to arrive at a final model.
If the initial interaction test suggests the interaction is important, then no further refinement should be considered and that model should be explored (this was the same protocol suggested in the 2-WAY ANOVA situation, the other place where we considered interactions). If the interaction is not deemed important based on the test, then the model should be re-fit using both variables in an additive model. In that additive model, both variables can be assessed conditional on the other one. If both have small p-values, then that is the final model and should be explored further. If either the categorical or quantitative variable have large p-values, then they can be dropped from the model and the model re-fit with only one variable in it, usually starting with dropping the component with the largest p-value if both are not “small”. Note that if there is only a categorical variable remaining, then we would call that linear model a One-Way ANOVA (quantitative response and \(J\) group categorical explanatory) and if the only remaining variable is quantitative, then a SLR model is being fit. If that final variable has a large p-value in either model, it can be removed and all that is left to describe the responses is a mean-only model. Otherwise the single variable model is the final model. Usually we will not have to delve deeply into this tree of models, but it is good to consider the potential paths that an analysis could involve before it is started.
To perform the first test (after checking that assumptions are met, of
course),
we can apply the Anova
function from the car
package to an interaction
model113.
It will provide three tests, one for each variable by themselves, which are not
too interesting, and then the interaction test. This will result in an
\(F\)-statistic that, if the assumptions are met, will follow an
\(F(J-1, n-2J)\)-distribution
under the null hypothesis. This tests
the hypotheses:
\(\boldsymbol{H_0:}\) The slope for \(\boldsymbol{x}\) is the same for all \(\boldsymbol{J}\) groups in the population vs
\(\boldsymbol{H_A:}\) The slope for \(\boldsymbol{x}\) in at least one group differs from the others in the population.
This test is also legitimate in the case of a two-level categorical variable \((J=2)\) and then follows an \(F(1, n-4)\)-distribution under the null hypothesis. With \(J=2\), the p-value from this test matches the results for the \(t\)-test \((t_{n-4})\) for the single slope-changing coefficient in the model summary output. The noise tolerance study, introduced in Section 8.10, provides a situation for exploring the results in detail.
With the \(J=4\) level categorical variable (Treatment), the model for the second noise tolerance measurement (du2) as a function of the interaction between Treatment and initial noise tolerance (du1) is
\[\begin{array}{rl} \text{du2}_i = \beta_0 &+ \beta_1\cdot\text{du1}_i + \beta_2I_{T2,i} + \beta_3I_{T3,i} + \beta_4I_{\text{Control},i} \\ &+ \beta_5I_{T2,i}\cdot\text{du1}_i + \beta_6I_{T3,i}\cdot\text{du1}_i + \beta_7I_{\text{Control},i}\cdot\text{du1}_i+\varepsilon_i. \end{array}\]
We can re-write the previous hypotheses in one of two more specific ways:
\(H_0:\) The slope for du1 is the same for all four Treatment groups in the population OR
\(H_0: \beta_5=\beta_6=\beta_7=0\)
- This defines a null hypothesis that all the deviation coefficients for getting different slopes for the different treatments are 0 in the population.
\(H_A:\) The slope for du1 is NOT the same for all four Treatment groups in the population (at least one group has a different slope) OR
\(H_A:\) At least one of \(\beta_5,\beta_6,\beta_7\) is different from 0 in the population.
- The alternative states that at least one of the deviation coefficients for getting different slopes for the different Treatments is not 0 in the population.
In this situation, the results for the test of these hypotheses is in the row
labeled du1:treatment
in the Anova
output. The ANOVA table below shows
a test statistic of \(F=0.768\) with the numerator df of 3, coming from \(J-1\),
and the denominator df of 90, coming from \(n-2J=98-2*4=90\) and also provided
in the Residuals
row in the table, leading to an \(F(3, 90)\)-distribution
for the test statistic under the null hypothesis.
The p-value from this
distribution is 0.515, showing little to no evidence against
the null hypothesis. The conclusion is that there is insufficient evidence to
claim that the slope coefficient for du1 in explaining du2 is different for
at least one of the Treatment groups in the population.
Anova(head2)
## Anova Table (Type II tests)
##
## Response: du2
## Sum Sq Df F value Pr(>F)
## du1 1197.78 1 259.5908 <2e-16
## treatment 23.90 3 1.7265 0.1672
## du1:treatment 10.63 3 0.7679 0.5150
## Residuals 415.27 90
Without evidence to support an interaction, we should consider both the quantitative and categorical variables in an additive model. The ANOVA table for the additive model contains two interesting tests. One test is for the quantitative variable discussed previously. The other is for the categorical variable, assessing whether different y-intercepts are needed. The additive model here is
\[\text{du2}_i = \beta_0 + \beta_1\cdot\text{du1}_i + \beta_2I_{T2,i} + \beta_3I_{T3,i} + \beta_4I_{\text{Control},i} +\varepsilon_i.\ \]
The hypotheses assessed in the ANOVA test for treatment are:
\(H_0:\) The y-intercept for the model with du1 is the same for all four Treatment groups in the population OR
\(H_0: \beta_2=\beta_3=\beta_4=0\)
- This defines a null hypothesis that all the deviation coefficients for getting different y-intercepts for the different Treatments are 0 in the population.
\(H_A:\) The y-intercepts for the model with du1 is NOT the same for all four Treatment groups in the population (at least one group has a different y-intercept) OR
\(H_A:\) At least one of \(\beta_2,\beta_3,\beta_4\) is different from 0 in the population.
- The alternative states that at least one of the deviation coefficients for getting different y-intercepts for the different Treatments is not 0 in the population.
The \(F\)-test for the categorical variable in an additive model follows \(F(J-1, n-J-1)\)-distribution under the null hypothesis. For this example, the test statistic for Treatment follows an \(F(3, 93)\)-distribution under the null hypothesis. The observed test statistic has a value of 1.74, generating a p-value of 0.164. So we would fail to reject the null hypothesis and conclude that there is little to no evidence to support the conclusion of some difference in y-intercepts between the treatment groups, in a model with du1, in the population. We could interpret this in the fashion we used initially in MLR by stating this result as: there is little to no evidence of a difference in the mean du2 for the Treatment groups after controlling for du1.
head1 <- lm(du2~du1+treatment, data=Headache)
Anova(head1)
## Anova Table (Type II tests)
##
## Response: du2
## Sum Sq Df F value Pr(>F)
## du1 1197.8 1 261.5491 <2e-16
## treatment 23.9 3 1.7395 0.1643
## Residuals 425.9 93
In the same ANOVA table, there is a test for the du1 effect. This tests \(H_0: \beta_1=0\) vs \(H_A: \beta_1\ne 0\) in a model with different y-intercepts estimated for the different groups. If we remove this term from the model, all we are left with is different y-intercepts for the groups. A model just with different y-intercepts is typically called a One-Way ANOVA model. Here, there is evidence that the quantitative variable is needed in the model after controlling for the different y-intercepts for different treatments. Note that this interpretation retains the conditional wording regardless of whether the other variable had a small p-value or it did not. If you want an unconditional interpretation for a variable, then you will need to refit the model without the other variable(s) after deciding that they are not important.
8.13 AICs for model selection
There are a variety of techniques for selecting among a set of potential models or refining an initially fit MLR model. Hypothesis testing can be used (in the case where we have nested models either by adding or deleting a single term at a time) or comparisons of adjusted R2 across different potential models (which is valid for nested or non-nested model comparisons). Diagnostics should play a role in the models considered and in selecting among models that might appear to be similar on a model comparison criterion. In this section, a new model selection method is introduced that has stronger theoretical underpinnings, a slightly more interpretable scale, and, often, better performance in picking an optimal114 model than the adjusted R2. The measure is called the AIC (Akaike’s An Information Criterion115, (Akaike 1974). It is extremely popular, but sometimes misused, in some fields such as Ecology and has been applied in almost every other potential application area where statistical models can be compared. Burnham and Anderson (2002) have been responsible for popularizing the use of AIC for model selection, especially in Ecology. The AIC is an estimate of the distance (or discrepancy or divergence) between a candidate model and the true model, on a log-scale, based on a measure called the Kullback-Leibler divergence. The models that are closer (have a smaller distance) to the truth are better and we can compare how close two models are to the truth, picking the one that has a smaller distance (smaller AIC) as better. The AIC includes a component that is on the log-scale, so negative values are possible and you should not be disturbed if you are comparing large magnitude negative numbers – just pick the model with the smallest AIC score.
The AIC is optimized (smallest) for a model that contains the optimal balance of simplicity of the model with quality of fit to the observations. Scientists are driven to different degrees by what is called the principle of parsimony: that simpler explanations (models) are better if everything else is equal or even close to equal. In this case, it would mean that if two models are similarly good on AIC, then select the simpler of the two models since it is more likely to be correct in general than the more complicated model. The AIC is calculated as \(AIC=-2log(Likelihood)+2m\), where the Likelihood provides a measure of fit of the model (we let R calculate it for us) and gets smaller for better fitting models and \(m\) = (number of estimated \(\beta\text{'s}+1\)). The value \(m\) is called the model degrees of freedom for AIC calculations and relates to how many total parameters are estimated. Note that it is a different measure of degrees of freedom than used in ANOVA \(F\)-tests. The main things to understand about the formula for the AIC is that as \(m\) increases, the AIC will go up and that as the fit improves, the likelihood will increase (so -2log-likelihood will get smaller)116.
There are some facets of this discussion to keep in mind when comparing models. More complicated models always fit better (we saw this for the R2 measure, as the proportion of variation explained always goes up if more “stuff” is put into the model even if the “stuff” isn’t useful). The AIC resembles the adjusted R2 in that it incorporates the count of the number of parameters estimated. This allows the AIC to make sure that enough extra variability is explained in the responses to justify making the model more complicated (increasing \(m\)). The optimal model on AIC has to balance adding complexity and increasing quality of the fit. Since this measure provides an estimate of the distance or discrepancy to the “true model”, the model with the smallest value “wins” – it is top-ranked on the AIC. Note that the top-ranked AIC model will often not be the best fitting model since the best fitting model is always the most complicated model considered. The top AIC model is the one that is estimated to be closest to the truth, where the truth is still unknown…
To help with interpreting the scale of AICs, they are often reported in a table sorted from smallest to largest values with the AIC and the “delta AIC” or, simply, \(\Delta\text{AIC}\) reported. The
\[\Delta\text{AIC}=\text{AIC}_{\text{model}} - \text{AIC}_{\text{topModel}}\]
and so provides a value of 0 for the top-ranked AIC model and a measure of how much worse on the AIC scale the other models are. A rule of thumb is that a 2 unit difference on AICs \((\Delta\text{AIC}=2)\) is decent evidence of a difference in the models and more than 4 units \((\Delta\text{AIC}>4)\) is a really big difference. This is more based on experience than a distinct reason or theoretical result but seems to provide reasonable results in most situations. Often researchers will consider any models within 2 AIC units of the top model \((\Delta\text{AIC}<2)\) as indistinguishable on AICs and so either select the simplest model of the choices or report all the models with similar “support”, allowing the reader to explore the suite of similarly supported potential models . It is important to remember that if you search across too many models, even with the AIC to support your model comparisons, you might find a spuriously top model. Individual results that are found by exploring many tests or models have higher chances to be spurious and results found in this manner are difficult to reproduce when someone repeats a similar study117. For these reasons, there is a set of general recommendations that have been developed for using AICs:
Consider a suite of models (often pre-specified and based on prior research in the area of interest) and find the models with the top (in other words, smallest) AIC results.
- The suite of candidate models need to contain at least some good models. Selecting the best of a set of BAD models only puts you at the top of $%#%-mountain, which is not necessarily a good thing.
Report a table with the models considered, sorted from smallest to largest AICs (\(\Delta\text{AICs}\) from smaller to larger) that includes a count of number of parameters estimated118, the AICs, and \(\Delta\text{AICs}\).
- Remember to incorporate the mean-only model in the model selection results. This allows you to compare the top model to one that does not contain any predictors.
Interpret the top model or top models if a few are close on the AIC-scale to the top model.
DO NOT REPORT P-VALUES OR CALL TERMS “SIGNIFICANT” when models were selected using AICs.
- Hypothesis testing and AIC model selection are not compatible philosophies and testing in models selected by AICs invalidates the tests as they have inflated Type I error rates. The AIC results are your “evidence” – you don’t need anything else. If you wanted to report p-values, use them to select your model.
You can describe variables as “important” or “useful” and report confidence intervals to aid in interpretation of the terms in the selected model(s) but need to avoid performing hypothesis tests with the confidence intervals.
Remember that the selected model is not the “true” model – it is only the best model according to AIC among the set of models you provided.
Model assumptions need to be met to use AICs. They assume that the model is specified correctly up to possibly comparing different predictor variables. Perform diagnostic checks on your initial model and the top model.
When working with AICs, there are two options. Fit the models of interest and
then run the AIC
function on each model. This can be tedious, especially
when we have many possible models to consider. We can make
it easy to fit all the potential candidate models that are implied by a
complicated starting model by using the dredge
function from the MuMIn
package (Barton 2019).
The name (dredge) actually speaks to what fitting
all possible models really engages – what is called data dredging.
The term is meant to
refer to considering way too many models for your data set, probably finding
something good from the process, but maybe identifying something spurious since
you looked at so many models. Note that if you take a hypothesis testing
approach where you plan to remove any terms with large p-values in this same
situation, you are really considering all possible models as well because you
could have removed some or all model components.
Methods that consider all
possible models are probably best used in exploratory analyses where you do not
know if any or all terms should be important. If you have more specific
research questions, then you probably should try to focus on comparisons of models
that help you directly answer those questions, either with AIC or p-value methods.
The dredge
function provides an automated method of assessing all possible
simpler models based on an initial (full) model. It generates a table of AIC
results, \(\Delta\text{AICs}\), and also shows when various predictors are in or
out of the model for all reduced models possible from an initial model. For
quantitative predictors, the estimated slope is reported when that predictor is
in the model. For categorical variables and interactions with them, it just
puts a “+” in the table to let you know that the term is in the models. Note
that you must run the options(na.action = "na.fail")
code to get dredge
to work.
To explore the AICs and compare their results to the adjusted R2
that we used before for model selection, we can revisit the Snow Depth data set
with related results found in Section 8.4 and
Table 8.1. In that situation we were considering
a “full” model that included Elevation, Min.Temp, and Max.Temp as potential
predictor variables after removing two influential points. And we considered all
possible reduced models from that “full”119 model. Note
that the dredge
output adds one more model that adjusted R2
can’t consider – the mean-only model that contains no predictor variables. In
the following output it is the last model in the output (worst ranked on AIC).
Including the mean-only model in these results helps us “prove” that there is
support for having something in the model, but only if there is better support for
other models than this simplest possible model.
In reading dredge
output120 as it is constructed here, the models are sorted by
top to bottom AIC values (smallest AIC to largest). The column delta
is for
the \(\Delta\text{AICs}\) and shows a 0 for the first row, which is the
top-ranked AIC model. Here it is for the model with Elevation and Max.Temp but
not including Min.Temp. This was also the top ranked model from adjusted
R2, which is reproduced in the adjRsq
column. The AIC
is calculated using the previous formula based on the df
and logLik
columns. The df
is also a useful column for comparing models as it helps
you see how complex each model is. For example, the top model used up 4
model df (three \(\beta\text{'s}\) and the residual error variance) and the
most complex model that included four predictor variables used up 5 model df.
require(MuMIn)
options(na.action = "na.fail") #Must run this code once to use dredge
snotel2R<-snotel2[-c(9,22),]; m6 <- lm(Snow.Depth~Elevation+Min.Temp+Max.Temp, data=snotel2R)
dredge(m6, rank="AIC",
extra = c("R^2", adjRsq=function(x) summary(x)$adj.r.squared))
## Global model call: lm(formula = Snow.Depth ~ Elevation + Min.Temp + Max.Temp, data = snotel2R)
## ---
## Model selection table
## (Int) Elv Max.Tmp Min.Tmp R^2 adjRsq df logLik AIC delta
## 4 -167.50 0.02408 1.2530 0.8495 0.8344 4 -80.855 169.7 0.00
## 8 -213.30 0.02686 1.2430 0.9843 0.8535 0.8304 5 -80.541 171.1 1.37
## 2 -80.41 0.01791 0.8087 0.7996 3 -83.611 173.2 3.51
## 6 -130.70 0.02098 1.0660 0.8134 0.7948 4 -83.322 174.6 4.93
## 5 179.60 -5.0090 0.6283 0.6106 3 -91.249 188.5 18.79
## 7 178.60 -0.2687 -4.6240 0.6308 0.5939 4 -91.170 190.3 20.63
## 3 119.50 -2.1800 0.4131 0.3852 3 -96.500 199.0 29.29
## 1 40.21 0.0000 0.0000 2 -102.630 209.3 39.55
## weight
## 4 0.568
## 8 0.286
## 2 0.098
## 6 0.048
## 5 0.000
## 7 0.000
## 3 0.000
## 1 0.000
## Models ranked by AIC(x)
There are two models that are clearly favored over the others with \(\Delta\text{AICs}\) for the model with Elevation and Max.Temp of 0 and for the model with all three predictors of 1.37. The \(\Delta\text{AIC}\) for the third ranked model (contains just Elevation) is 3.51 suggesting clear support for the top model over this. The difference between the second and third ranked models also provides relatively strong support for the more complex model over the model with just Elevation. And the mean-only model had a \(\Delta\text{AIC}\) of nearly 40 – suggesting extremely strong support for the top model versus using no predictors. So we have pretty clear support for models that include the Elevation and Max.Temp variables (in both top models) and some support for also including the Min.Temp, but the top model did not require its inclusion.
We could add further explorations of the term-plots and confidence intervals for the slopes from the top or, here, possibly top two models. We would not spend any time with p-values since we already used the AIC to assess importance of the model components and they are invalid if we model select prior to reporting them. We can quickly compare the slopes for variables that are shared in the two models since they are both quantitative variables using the output. It is interesting that the Elevation and Max.Temp slopes change little with the inclusion of Min.Temp in moving from the top to second ranked model (0.02408 to 0.0286 and 1.253 to 1.243).
This was an observational study and so we can’t consider causal inferences here as discussed previously. Generally, the use of AICs does not preclude making causal statements but if you have randomized assignment of levels of an explanatory variable, it is more philosophically consistent to use hypothesis testing methods in that setting. If you went to the effort to impose the levels of a treatment on the subjects, it also makes sense to see if the differences created are beyond what you might expect by chance if the treatment didn’t matter.
8.14 Case study: Forced expiratory volume model selection using AICs
Researchers were interested in studying the effects of smoking by children on
their lung development by measuring the forced expiratory volume (FEV,
measured in Liters) in a representative sample of children (\(n=654\)) between
the ages of 3 and 19; this data set is available in the FEV
data set in the
coneproj
package (Meyer and Liao (2018), Liao and Meyer (2014)).
Measurements on the age (in years)
and height (in inches) as well as the sex and smoking status of the
children were made. We would expect both the age and height to have positive
relationships with FEV (lung capacity) and that smoking might decrease the lung
capacity but also that older children would be more likely to smoke. So the
height and age might be confounded with smoking status and smoking
might diminish lung development for older kids – resulting in a
potential interaction between age and smoking. The sex of the child might
also matter and should be considered or at
least controlled for since the response is a size-based measure. This creates
the potential for including up to four variables (age, height, sex, and
smoking status) and possibly the interaction between age and smoking status.
Initial explorations suggested that modeling the log-FEV would be more successful
than trying to model the responses on the original scale. Figure
8.35 shows the suggestion of different slopes for the
smokers than non-smokers and that there aren’t very many smokers under 9 years
old in the data set.
So we will start with a model that contains an age by smoking interaction
and include height and sex as additive terms. We are not sure if any of
these model components will be needed, so the simplest candidate model will be
to remove all the predictors and just have a mean-only model (FEV~1
). In
between the mean-only and most complicated model are many different options
where we can drop the interaction or drop the additive terms or drop the terms
involved in the interaction if we don’t need the interaction.
require(coneproj)
data(FEV)
FEV <- as.tibble(FEV)
FEV$sex <- factor(FEV$sex) #Make sex a factor
levels(FEV$sex) <- c("Female","Male") #Make sex labels explicit
FEV$smoke <- factor(FEV$smoke) #Make smoking status a factor
levels(FEV$smoke) <- c("Nonsmoker","Smoker") #Make smoking status labels explicit
require(car)
scatterplot(log(FEV)~age|smoke, data=FEV, smooth=F,
main="Plot of log(FEV) vs Age of children by smoking status",
legend=list(coords="topleft", columns=1, cex.leg=0.75))
To get the needed results, start with the full model – the most
complicated model you want to consider. It is good to check assumptions before
considering reducing the model as they rarely get better in simpler models and
the AIC is only appropriate to use if the model assumptions are reasonably
well-met. As suggested above, our “fullish” model for the log(FEV) values is
specified as log(FEV)~height+age*smoke+sex
.
(ref:fig8-35) Diagnostics for the log(FEV) model that includes height, sex, and an interaction between age and smoking status (the full model).
fm1 <- lm(log(FEV)~height+age*smoke+sex, data=FEV)
summary(fm1)
##
## Call:
## lm(formula = log(FEV) ~ height + age * smoke + sex, data = FEV)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62926 -0.08783 0.01136 0.09658 0.40751
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.919494 0.080571 -23.824 < 2e-16
## height 0.042066 0.001759 23.911 < 2e-16
## age 0.025368 0.003642 6.966 8.03e-12
## smokeSmoker 0.107884 0.113646 0.949 0.34282
## sexMale 0.030871 0.011764 2.624 0.00889
## age:smokeSmoker -0.011666 0.008465 -1.378 0.16863
##
## Residual standard error: 0.1454 on 648 degrees of freedom
## Multiple R-squared: 0.8112, Adjusted R-squared: 0.8097
## F-statistic: 556.8 on 5 and 648 DF, p-value: < 2.2e-16
par(mfrow=c(2,2), oma=c(0,0,2,0))
plot(fm1, sub.caption="Diagnostics for full FEV model")
The diagnostic plots suggest that there are a few outlying points (Figure 8.36) but they are not influential and there is no indication of violations of the normality and constant variance assumptions. If we select a different model(s), we would want to check its diagnostics and make sure that the results do not look noticeably worse than these do.
The AIC
function can be used to generate the AIC values for a
single or set
of candidate models. It will also provide the model degrees of freedom used for
each model. For example, suppose that the want to compare fm1
to a model
without the interaction term in the model, called fm1R
. You need to fit both
models and then apply the AIC
function to them with commas between the model
names:
fm1R <- lm(log(FEV)~height+age+smoke+sex, data=FEV)
AIC(fm1, fm1R)
## df AIC
## fm1 7 -658.5178
## fm1R 6 -658.6037
These results tells us that the fm1R
model (the one without the interaction)
is better on the AIC by 0.09 AIC units. Note that this model does not “fit” as well
as the full model, it is just the top AIC model – the AIC results suggest that
it is slightly closer to the truth than the more complicated model. But this
provides only an assessment of the difference between including or excluding
the interaction between age and smoking in a model with two other predictors.
We are probably also interested in whether the other terms are
needed in the model. The full suite of results from dredge provide model
comparisons that help us to assess the presence/absence of each model component
including the interaction.
require(MuMIn)
options(na.action = "na.fail") #Must run this code once to use dredge
dredge(fm1, rank="AIC",
extra = c("R^2", adjRsq=function(x) summary(x)$adj.r.squared))
## Global model call: lm(formula = log(FEV) ~ height + age * smoke + sex, data = FEV)
## ---
## Model selection table
## (Int) age hgh sex smk age:smk R^2 adjRsq df logLik
## 16 -1.944000 0.02339 0.04280 + + 0.81060 0.80950 6 335.302
## 32 -1.919000 0.02537 0.04207 + + + 0.81120 0.80970 7 336.259
## 8 -1.940000 0.02120 0.04299 + 0.80920 0.80830 5 332.865
## 12 -1.974000 0.02231 0.04371 + 0.80880 0.80790 5 332.163
## 28 -1.955000 0.02388 0.04315 + + 0.80920 0.80800 6 332.802
## 4 -1.971000 0.01982 0.04399 0.80710 0.80650 4 329.262
## 7 -2.265000 0.05185 + 0.79640 0.79580 4 311.594
## 3 -2.271000 0.05212 0.79560 0.79530 3 310.322
## 15 -2.267000 0.05190 + + 0.79640 0.79550 5 311.602
## 11 -2.277000 0.05222 + 0.79560 0.79500 4 310.378
## 30 -0.067780 0.09493 + + + 0.64460 0.64240 6 129.430
## 26 -0.026590 0.09596 + + 0.62360 0.62190 5 110.667
## 14 -0.015820 0.08963 + + 0.62110 0.61930 5 108.465
## 6 0.004991 0.08660 + 0.61750 0.61630 4 105.363
## 10 0.022940 0.09077 + 0.60120 0.60000 4 91.790
## 2 0.050600 0.08708 0.59580 0.59520 3 87.342
## 13 0.822000 + + 0.09535 0.09257 4 -176.092
## 9 0.888400 + 0.05975 0.05831 3 -188.712
## 5 0.857400 + 0.02878 0.02729 3 -199.310
## 1 0.915400 0.00000 0.00000 2 -208.859
## AIC delta weight
## 16 -658.6 0.00 0.414
## 32 -658.5 0.09 0.397
## 8 -655.7 2.87 0.099
## 12 -654.3 4.28 0.049
## 28 -653.6 5.00 0.034
## 4 -650.5 8.08 0.007
## 7 -615.2 43.42 0.000
## 3 -614.6 43.96 0.000
## 15 -613.2 45.40 0.000
## 11 -612.8 45.85 0.000
## 30 -246.9 411.74 0.000
## 26 -211.3 447.27 0.000
## 14 -206.9 451.67 0.000
## 6 -202.7 455.88 0.000
## 10 -175.6 483.02 0.000
## 2 -168.7 489.92 0.000
## 13 360.2 1018.79 0.000
## 9 383.4 1042.03 0.000
## 5 404.6 1063.22 0.000
## 1 421.7 1080.32 0.000
## Models ranked by AIC(x)
There is a lot of information in the output and some of the needed information
in the second set of rows, so we will try to point out some useful features to
consider. The left columns describe the models being estimated. For example, the
first row of results is for a model with an intercept (Int
), age (age
)
, height (hgh
), sex (sex
), and smoking(smk
). For sex and
smoking, there are “+
”s in the output row when they are included in that
model but no coefficient since they are categorical variables. There is no
interaction between age and smoking in the top ranked model. The top AIC
model has an \(\boldsymbol{R}^2=0.8106\), adjusted R2 of 0.8095,
model df=6 (from an intercept, four slopes, and the residual variance),
log-likelihood (logLik
)=335.302, an AIC=-658.6 and \(\Delta\text{AIC}\) of 0.00.
The next best model adds the interaction between age and smoking, resulting
in increases in the R2, adjusted R2, and
model df, but increasing the AIC by 0.09 units \((\Delta\text{AIC}=0.09)\).
This suggests that these two models are essentially equivalent on the AIC because
the difference is so small. The simpler model is a little bit better on AIC so
you could focus on it or on the slightly more complicated model – but you should
probably note that the evidence is equivocal for these two models.
The comparison to other potential models shows the strength of evidence in support of all the other model components. The intercept-only model is again the last in the list with the least support on AICs with a \(\Delta\text{AIC}\) of 1080.32, suggesting it is not worth considering in comparison with the top model. Comparing the mean-only model to our favorite model on AICs is a bit like the overall \(\boldsymbol{F}\)-test we considered in Section 8.7 because it compares a model with no predictors to a complicated model. Each model with just one predictor included is available in the table as well, with the top single predictor model based on height having a \(\Delta\text{AIC}\) of 43.96. So we certainly need to pursue something more complicated than SLR based on the AIC results. Closer to the top model is the third-ranked model that includes age, height, and sex. It has a \(\Delta\text{AIC}\) of 2.87 so we would say that these results present marginal support for the top two models over this model. It is the simplest model of the top three but not close enough to be considered in detail.
This table also provides the opportunity to compare the model selection results from the adjusted R2 compared to the AIC. The AIC favors the model without an interaction between age and smoking whereas the adjusted R2 favors the most complicated model considered here that included an age and smoking interaction. The AIC provides units that are more interpretable than adjusted R2 even though the scale for the AIC is a bit mysterious as distances from the unknown true model.
The top AIC model (and possibly the other similar models) can then be explored in more detail. You should not then focus on hypothesis testing in this model121. Confidence intervals and term-plots are useful for describing the different model components and making inferences for the estimated sizes of differences in the population. These results should not be used for deciding if terms are “significant” when the models are selected using measures like the AIC or adjusted R2. But you can discuss the estimated model components to go with how you arrived at having them in the model.
In this situation, the top model is estimated to be
\[\log(\widehat{\text{FEV}})_i = -1.94 + 0.043\cdot\text{Height}_i+ 0.0234\cdot\text{Age}_i -0.046I_{\text{Smoker},i}+0.0293I_{\text{Male},i}\]
based on the estimated coefficients provided below. Using these results and the term-plots (Figure 8.37) we see that in this model there are positive slopes for Age and Height on log-FEV, a negative coefficient for smoking (Smoker), and a positive coefficient for sex (Males). We could go further with interpretations such as for the age term: For a 1 year increase in age, we expect, on average, a 0.0234 log-liter increase in FEV, after controlling for the height, smoking status, and sex of the children. We can even interpret this on the original scale since this was a log(y) response model using the same techniques as in Section 7.6. If we exponentiate the slope coefficient of the quantitative variable, \(\exp(0.0234)=1.0237\). This provides the interpretation on the original FEV scale, for a 1 year increase in age, we expect 2.4% increase in the median FEV, after controlling for the height, smoking status, and sex of the children. The only difference from Section 7.6 when working with a log(y) model now is that we have to note that the model used to generate the slope coefficient had other components and so this estimate is after adjusting for them.
(ref:fig8-36) Term-plots for the top AIC model for log(FEV) that includes height, age, smoking status and sex in the model.
fm1R$coefficients
## (Intercept) height age smokeSmoker sexMale
## -1.94399818 0.04279579 0.02338721 -0.04606754 0.02931936
confint(fm1R)
## 2.5 % 97.5 %
## (Intercept) -2.098414941 -1.789581413
## height 0.039498923 0.046092655
## age 0.016812109 0.029962319
## smokeSmoker -0.087127344 -0.005007728
## sexMale 0.006308481 0.052330236
plot(allEffects(fm1R))
Like any statistical method, the AIC works better with larger sample sizes and
when the assumptions are met. It also will detect important variables in models
more easily when the effects of the predictor variables are strong. Along with
the AIC results, it is good to report the coefficients for your top estimated
model(s), confidence intervals for the coefficients and/or term-plots, and
R2. This provides a useful summary of the reasons for selecting
the model(s), information on the importance of the terms within the model, and
a measure of the variability explained by the model. The R2 is not used to
select the model, but after selection can be a nice summary of model quality.
For fm1R
, the \(R^2=0.8106\) suggesting that the selected model explains
81% of the variation in log-FEV values.
The AICs are a preferred modeling strategy in some areas such as Ecology. As
with this and many other
methods discussed in this book, it is sometimes as easy to find journal
articles with mistakes in using statistical methods as it is to find papers
doing it correctly. After completing this material, you have the potential to
have the knowledge and experience of two statistics classes and now are better
trained than some researchers that frequently use these methods. This set of tools
can be easily mis-applied. Try to make sure that you are thinking carefully
through your problem before jumping to the statistical results. Make a graph
first, think carefully about your study design and variables collected and what
your models of interest might be, what assumptions might be violated based on
the data collection story, and then start fitting models. Then check your
assumptions and only proceed on with any inference if those conditions are
reasonably well-met. The AIC provides an alternative method for selecting among
different potential models and they do not need to be nested (a requirement of
hypothesis testing methods used to sequentially simplify models). The automated
consideration of all possible models in the dredge
function should not be
considered in all situations but can be useful in a preliminary model exploration
study where no clear knowledge exists about useful models to consider. Reporting the summary of AIC results beyond just reporting the top model(s) selected for focused exploration provides the evidence to support that selection.
8.15 Chapter summary
This chapter explored the most complicated models we’re going to explore. MLR models can incorporate features of SLR and ANOVAs. The MLR’s used in this chapter highlight the flexibility of the linear modeling framework to move from two-sample mean models to multi-predictor models with interactions of categorical and quantitative variables. It is useful to use the pertinent names for the simpler models, but at this point we could have called everything we are doing fitting linear models. The power of the linear model involves being able to add multiple predictor variables to the model and handle categorical predictors using indicator variables. All this power comes with some responsibility in that you need to know what you are trying to fit and how to interpret the results provided. We introduced each scenario working from simple to the most complicated version of the models, trying to motivate when you would encounter them, and the specific details of the interpretations of each type of model. In Chapter 9, case studies are used to review the different methods discussed with reminders of how to identify and interpret the particular methods used.
When you have to make modeling decisions, you should remember the main priorities in modeling. First, you need to find a model that can address research questions of interest. Second, find a model that is trustworthy and has assumptions that are reasonably well met. Third, report the logic and evidence that was used to identify and support the model. All too often, researchers present only a final model with little information on how they arrived at it. You should be reporting the reasons for decisions made and the evidence supporting them. For example, if you were considering an interaction model and the interaction was dropped and an additive model is re-fit and interpreted, the evidence related to the interaction test should still be reported. Similarly, if a larger MLR is considered and some variables are removed, the evidence (reason) for those removals should be provided. Because of multicollinearity in models, you should never remove more than one quantitative predictor at a time or else you could remove two variables that are important but were “hiding” when both were included in the model.
8.16 Summary of important R code
There is very little “new” R code in this chapter since all these methods were
either used in the ANOVA or SLR chapters. The models are more complicated but
are built off of methods from previous chapters. In this code, y
is a
response variable, x1
, x2
, …, xK
are quantitative
explanatory variables, group
is a factor variable and the data are
in DATASETNAME
.
scatterplot(y~x1|group, data=DATASETNAME, smooth=F)
Requires the
car
package.Provides a scatterplot with a regression line for each group.
MODELNAME
<-
lm(y~ x1+x2+…+xK, data=DATASETNAME)- Estimates an MLR model using least squares with \(K\) quantitative predictors.
**MODELNAME
<-
lm(y~ x1*group, data=DATASETNAME)**- Estimates an interaction model between a quantitative and categorical variable, providing different slopes and intercepts for each group.
MODELNAME
<-
lm(y~ x1+group, data=DATASETNAME)- Estimates an additive model with a quantitative and categorical variable, providing different intercepts for each group.
summary(MODELNAME)
- Provides parameter estimates, overall \(F\)-test, R2, and adjusted R2.
par(mfrow=c(2, 2)); plot(MODELNAME)
- Provides four regression diagnostic plots in one plot.
confint(MODELNAME, level=0.95)
Provides 95% confidence intervals for the regression model coefficients.
Change
level
if you want other confidence levels.
plot(allEffects(MODELNAME))
Requires the
effects
package.Provides a plot of the estimated regression lines with 95% confidence interval for the mean.
vif(MODELNAME)
Requires the
car
package.Provides VIFs for an MLR model. Only use in additive models - not meaningful for terms involved in interactions.
predict(MODELNAME, se.fit=T)
- Provides fitted values for all observed \(x\text{'s}\) with SEs for the mean.
predict(MODELNAME, newdata=tibble(x1 = X1_NEW, x2 = X2_NEW, …, xK = XK_NEW, interval=“confidence”)
- Provides fitted value for specific values of the quantitative predictors with CI for the mean.
predict(MODELNAME, newdata=tibble(x1 = X1_NEW, x2 = X2_NEW, …, xK = XK_NEW, interval=“prediction”)
- Provides fitted value for specific values of the quantitative predictors with PI for a new observation.
Anova(MODELNAME)
Requires the
car
package.Use to generate ANOVA tables and \(F\)-tests useful when categorical variables are included in either the additive or interaction models.
AIC(MODELNAME_1, MODELNAME_2)
- Use to get AIC results for two candidate models called
MODELNAME_1
andMODELNAME_2
.
- Use to get AIC results for two candidate models called
options(na.action = “na.fail”)
dredge(FULL_MODELNAME, rank=“AIC”)Requires the
MuMIn
package.Provides AIC and delta AIC results for all possible simpler models given a full model called
FULL_MODELNAME
.
8.17 Practice problems
The original research goal for the treadmill data set used for practice problems in the last two chapters was to replace the costly treadmill oxygen test with a cheap to find running time measurement but there were actually quite a few variables measured while the run time was found – maybe we can replace the treadmill test result with a combined prediction built using a few variables using the MLR techniques. The following code will get us re-started in this situation.
treadmill <- read_csv("http://www.math.montana.edu/courses/s217/documents/treadmill.csv")
tm1 <- lm(TreadMillOx~RunTime, data=treadmill)
8.1. Fit the MLR that also includes the running pulse (RunPulse
), the
resting pulse (RestPulse
), body weight (BodyWeight
), and Age (Age
)
of the subjects using the following code. Report and interpret the R2 for
this model.
8.2. Compare R2 and the adjusted R2 to the results for the SLR model
that just had RunTime
in the model. What do these results suggest?
8.3. Interpret the estimated RunTime
slope coefficients from the SLR model
and this MLR model. Explain the differences in the estimates.
8.4. Find the VIFs for this model and discuss whether there is an issue with multicollinearity noted in these results.
8.5. Report the value for the overall \(F\)-test for the MLR model and interpret the result.
8.6. Drop the variable with the largest p-value in the MLR model and re-fit it. Compare the resulting R2 and adjusted R2 values to the others found previously.
8.7. Use the dredge
function as follows to consider some other potential
reduced models and report the top two models according to adjusted
R2 values. What model had the highest R2?
Also discuss and compare the model selection results provided by the delta AICs
here.
require(MuMIn)
options(na.action = "na.fail") #Must run this code once to use dredge
dredge(mlr1, rank="AIC",
extra=c("R^2", adjRsq=function(x) summary(x)$adj.r.squared))
8.8. For one of the models, interpret the Age slope coefficient. Remember that only male subjects between 38 and 57 participated in this study. Discuss how this might have impacted the results found as compared to a more general population that could have been sampled from.
8.9. The following code creates a new three-level variable grouping the ages into low, middle, and high for those observed. The scatterplot lets you explore whether the relationship between treadmill oxygen and run time might differ across the age groups.
treadmill$Ageb <- cut(treadmill$Age, breaks=c(37,44.5,50.5,58))
summary(treadmill$Ageb)
require(car)
scatterplot(TreadMillOx~RunTime|Ageb, data=treadmill, smooth=F, lwd=2)
Based on the plot, do the lines look approximately parallel or not?
8.10. Fit the MLR that contains a RunTime
by Ageb
interaction – do not
include any other variables. Compare the R2 and adjusted
R2 results to previous models.
8.11. Find and report the results for the \(F\)-test that assesses evidence relative to the need for different slope coefficients.
8.12. Write out the overall estimated model. What level was R using as baseline? Write out the simplified model for two of the age levels.
8.13. Fit the additive model with RunTime
and predict the mean treadmill
oxygen values for subjects with run times of 11 minutes in each of the three
Ageb
groups.
8.14. Find the \(F\)-test results for the binned age variable in the additive model. Report and interpret those results.
References
Akaike, Hirotugu. 1974. “A New Look at the Statistical Model Identification.” IEEE Transactions on Automatic Control 19: 716–23.
Barton, Kamil. 2019. MuMIn: Multi-Model Inference. https://CRAN.R-project.org/package=MuMIn.
Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference. NY: Springer.
Diez, David M, Christopher D Barr, and Mine Cetinkaya-Rundel. 2017. Openintro: Data Sets and Supplemental Functions from ’Openintro’ Textbooks. https://CRAN.R-project.org/package=openintro.
Fox, John. 2003. “Effect Displays in R for Generalised Linear Models.” Journal of Statistical Software 8 (15): 1–27. http://www.jstatsoft.org/v08/i15/.
Fox, John, and Michael Friendly. 2018. Heplots: Visualizing Hypothesis Tests in Multivariate Linear Models. https://CRAN.R-project.org/package=heplots.
Fox, John, Sanford Weisberg, and Brad Price. 2018b. CarData: Companion to Applied Regression Data Sets. https://CRAN.R-project.org/package=carData.
Liao, Xiyue, and Mary C. Meyer. 2014. “Coneproj: An R Package for the Primal or Dual Cone Projections with Routines for Constrained Regression.” Journal of Statistical Software 61 (12): 1–22. http://www.jstatsoft.org/v61/i12/.
Merkle, Ed, and Michael Smithson. 2018. Smdata: Data to Accompany Smithson & Merkle, 2013. https://CRAN.R-project.org/package=smdata.
Meyer, Mary C., and Xiyue Liao. 2018. Coneproj: Primal or Dual Cone Projections with Routines for Constrained Regression. https://CRAN.R-project.org/package=coneproj.
If you take advanced applied mathematics courses, you can learn more about the algorithms being used by
lm
. Everyone else only cares about the algorithms when they don’t work – which is usually due to the user’s inputs in these models not the algorithm itself.↩Sometimes the
effects
plots ignores the edge explanatory observations with the default display. Always check the original variable summaries when considering the range of observed values. The code modifications used to make these plots are a little annoying and can be avoided if you remember to check the original summary statistics or scatterplots for variable extrema.↩We used this same notation in the fitting the additive Two-Way ANOVA and this is also additive in terms of these variables. Interaction models are discussed later in the chapter.↩
Imagine showing up to a ski area expecting a 40 inch base and there only being 11 inches. I’m sure ski areas are always more accurate than this model in their reporting of amounts of snow on the ground…↩
The site name is redacted to protect the innocence of the reader. More information on this site, located in Beaverhead County, is available at http://www.wcc.nrcs.usda.gov/nwcc/site?sitenum=355&state=mt.↩
The
seq
function has syntax ofseq(from=startingpoint, to=endingpoint, length.out=#ofvalues_between_start_and_end)
and therep
function has syntax ofrep(numbertorepeat, #oftimes).
↩Also see Section 8.13 for another method of picking among different models.↩
This section was inspired by a similar section from De Veaux, Velleman, and Bock (2011).↩
There are some social science models where the model is fit with the mean subtracted from each predictor so all have mean 0 and the precision of the y-intercept is interesting. But even in these models, the test for the y-intercept being 0 is rarely of interest.↩
Either someone had a weighted GPA with bonus points, or more likely here, there was a coding error in the data set since only one observation was over 4.0 in the GPA data. Either way, we could remove it and note that our inferences for HSGPA do not extend above 4.0.↩
When there are just two predictors, the VIFs have to be the same since the proportion of information shared is the same in both directions. With more than two predictors, each variable can have a different VIF value.↩
We are actually making an educated guess about what these codes mean. Other similar data sets used 1 for males but the documentation on these data is a bit sparse. We proceed with a small potential that the conclusions regarding differences in
sex
are in the wrong direction.↩Some people also call them dummy variables to reflect that they are stand-ins for dealing with the categorical information.↩
That may not read how I intended…↩
This is true for additive uses of indicator variables. In Section 8.11, we consider interactions between quantitative and categorical variables which has the effect of changing slopes and intercepts. The simplification ideas to produce estimated equations for each group are used there but we have to account for changing slopes by group too.↩
When making the SATM term-plot, the categorical variable is held at the most frequently occurring value in the data set. If you drop
ci.style="lines"
from the effect plot options, it is best to copy the figures as Bitmaps or save them as an image or they may (for some reason) lose the shaded bands in some word processing programs.↩Models like this with a categorical variable and quantitative variable are often called ANCOVA or analysis of covariance models but really are just versions of our linear models we’ve been using throughout this material.↩
Note that we employed some specific options in the
legend
option to get the legend to fit on this scatterplot better. Usually you can avoid this but thecoords
option defined a location and thecolumns
option made it a two column legend.↩Standardizing quantitative predictor variables is popular in social sciences, often where the response variable is also standardized. In those situations, they generate what are called “standardized slopes” (https://en.wikipedia.org/wiki/Standardized_coefficient) that estimate the change in SDs in the response for a 1 SD change in the explanatory variable.↩
There is a way to test for a difference in the two lines at a particular \(x\) value but it is beyond the scope of this material and uses what are called “contrasts”.↩
This is an example of what is called “step down” testing for model refinement which is a commonly used technique for arriving at a final model to describe response variables. Note that each step in the process should be reported, not just the final model that only has variables with small p-values remaining in it.↩
We could also use the
anova
function to do this but usingAnova
throughout this material provides the answers we want in the additive model and it has no impact for the only test of interest in the interaction model since the interaction is the last component in the model.↩In most situations, it would be crazy to assume that the true model for a process has been obtained so we can never pick the “correct” model. In fact, we won’t even know if we are picking a “good” model, but just the best from a set of the candidate models. But we can study the general performance of methods using simulations where we know the true model and the AIC has some useful properties in identifying the correct model when it is in the candidate set of models. No such similar theory exists for the adjusted R2.↩
Most people now call this Akaike’s (pronounced ah-kah-ee-kay) Information Criterion, but he used the AIC nomenclature to mean An Information Criterion – he was not so vain as to name the method after himself in the original paper that proposed it.↩
More details on these components of the methods will be left for more advanced material - we will focus on an introduction to using the AIC measure here.↩
Reproducibility ideas are used in statistics first by making data and code used available to others (like all the code and data in this book) and second by trying to use methods that when others perform similar studies they will find similar results (from Physics, think of the famous cold-fusion experiments http://en.wikipedia.org/wiki/Cold_fusion.)↩
Although sometimes excluded, the count of parameters should include counting the residual variance as a parameter.↩
We put quotes on “full” or sometimes call it the “fullish” model because we could always add more to the model, like interactions or other explanatory variables. So we rarely have a completely full model but we do have our “most complicated that we are considering” model.↩
The options in
extra=...
are to get extra information displayed that you do not necessarily need. You can simply rundredge(m6,rank="AIC")
to get just the AIC results.↩Hypothesis testing so permeates the use of statistics that even after using AICs many researchers are pressured to report p-values for model components. Some of this could be confusion caused when people first learned these statistical methods because when we teach you statistics we show you how to use various methods, one after another, and forget to mention that you should not use every method we taught you in every analysis.↩


Activities for you to
enjoy!
Our Resorts have fantasic activities for you and your family to enjoy!
Pool/ Spa
Arcade/ Bowling Alley
Tours
Golf tournements
Tennis
Campfire/ BBQ pit
Pitures
here.


Airsoft brands we work with
Here at Ireland Airsoft, we care about our customers, for this reason, we have taken a long research and checking nearly a hundred different companies until we could find the best brands for Airsoft guns in the market. Along the way, however, we realized that the companies were split into two categories – budget brands and performance brands. Most people are looking for one or the other, so in our Shops, you are able to find great guns for every budget.
In order to make it easier for you while shopping for your perfect gun, we have broken down below all our brands into Customer satisfaction and Average price.
Top Brands by Customer satisfaction:
As you can see, from our researches, we found out that most people have a higher satisfaction with G&G, Tokyo Marui and UK Arms, but please take note that all other brands listed here are all of great quality, but we have broken down the top ones to make it easier for your purchase.
Along with customer satisfaction, budget is an important variable, as you might be looking for the best gun in the market, or you could be looking for the best bang for your buck, and if that's the case, the most cheap guns are: Wellfare, ICS, Cybergun, if you click below in "Create Chart you will be able to see it all our brands listed in a Graph per average price.
Top Brands by Average Price:


class TybTickets::Member
Public Instance Methods
stats()
click to toggle source
Retrieves membership statistics with the number of completed tasks
# File app/models/tyb_tickets/member.rb, line 9 scope :stats, ->(from_date, to_date) { joins(:tickets) .group(:member_id) .select("tyb_tickets_members.*, COUNT(*) as tickets_count") .where("closed_at BETWEEN ? AND ?", from_date, to_date) .order(tickets_count: :desc) }
tickets()
click to toggle source
Members tickets
# File app/models/tyb_tickets/member.rb, line 6 has_many :tickets


[RRR-44] Need to be able to choose resources from a list | |
| Status: | Closed |
| Project: | ReportRunner |
| Component/s: | Web Wizard |
| Affects Version/s: | 0.3 |
| Fix Version/s: | 0.3 |
| Type: | Improvement | Priority: | Major |
| Reporter: | Henri Yandell | Assignee: | Henri Yandell |
| Resolution: | Fixed | ||
| Environment: |
| Description |
| Rather than just having a report run on only one resource, it should be able to run against a list of options. Currently thinking that the resource attribute of a report would look something like: resource="JiraDS or ConfluenceDS". This won't scope well when people want to run a report that uses two resources though. *** It's a param! Doh :) Though it has to be answered before the param page can, so maybe a new tag: <resource-param name="SqlReportDS" options="JiraDS,ConfluenceDS,RollerDS"/> This scopes nicely and can even replace the resource attribute of the report tag: <resource-param name="SqlReportDS" value="JiraDS"/> would replace resource="JiraDS" in report tag. |
| Comment by bayard [ Mon, 6 Dec 2004 18:02:22 -0800 (PST) ] |
| Half coded, and now needs recoding. The <resource-param> should be a part of a ReportGroup, not a Report. The Resource should use reflection to be set on the report. So SqlReport should have setSqlReportDS() as a method. The report xml may need nothing resource-wise...ponder. The jsp choice should happen in the same place as currently, ie) after the report is chosen. |
| Comment by bayard [ Tue, 7 Dec 2004 16:42:01 -0800 (PST) ] |
| This is implemented. Bit of cleanup to do, but I'll make that a new issue. The use tag needs to die in the reportrunner.xml and should just be a resource-param without options. |


WITNESS LIST Finance Committee February 26, 2001 - 8:00A Markup Session on S.B. 1 ON: Riley, Leonard Executive Director (Texas Workers Compensation Commission), Austin, TX Markup Session on SB 1 ON: Dodson, Gay Executive Director (Texas State Board of Pharmacy) Hamilton, Billy Deputy Commissioner (E-Texas - Comptroller), Austin, TX Mathis, Benny Director (Structural Pest Control Board) SB 174 exemption from ad valorem taxation FOR: Dougherty, Stephen W. General Manager (International Distribution Corp. & Texas Warehouse Association), Friendswood, TX Frasier, A. K. President (Network Distribution & Texas Warehouse Association), Castle Hills, TX Perryman, Ray Economist (Texas Warehouse Association), Waco, TX AGAINST: Jeffers, James P. City Manager (City of Plainview), Plainview, TX Lavine, Dick Fiscal Analyst (Center for Public Policy Priorities), Austin, TX Registering, but not testifying: For: Binko, Joseph J. (Texas Warehouse Association), Houston, TX Providing written testimony: For: Kelley, Larry Manager (Wilson Warehouse Co. of Texas, Inc. & Texas Warehouse Association), Winnie, TX Martineez, Jr., Mercurio County Judge (Webb County), Laredo, TX Vetters, Bob (Texas Warehouse Association), San Antonio, TX Zuniga, Federico C. (Diamond Distribution), Laredo, TX SJR 6 proposing a constitutional amendmen FOR: Dougherty, Stephen W. General Manager (International Distribution Corp. & Texas Wrehouse Association), Friendswood, TX Frasier, A.K. President (Network Distribution & Texas Warehouse Association), Castle Hills, TX Perryman, Ray Economist (Texas Warehouse Association), Waco, TX AGAINST: Jeffers, James P. City Manager (City of Plainview), Plainview, TX Lavine, Dick Fiscal Analyst (Center for Public Policy Priorities), Austin, TX Registering, but not testifying: For: Binko, Joseph J. (Texas Warehouse Association), Houston, TX Providing written testimony: For: Kelley, Larry Manager (Wilson Warehouse Co. of Texas, Inc. & Texas Warehouse Association), Winnie, TX Martiez, Jr., Mercurio County Judge (Webb County Webb County), Laredo, TX Vetters, Bob (Texas Warehouse Association), San Antonio, TX Zuniga, Federico C. (Diamond Distribution), Laredo, TX


Wow, I made a cool Web page!
To render this as a Web page and see it in your browser, choose RubyFrontier > Publish Page.
Even better, choose RubyFrontier > Publish Site to see all three example pages.
This is an example of a page whose HTML is created “by hand”. We use SmartyPants, though, to get the curly quotes in that previous sentence.
Next: My Second Web Page


--- layout: default title: "Resume" description: "This page is a part of Kalpit Thakkar's personal website. It has Kalpit Thakkar's resume which is available for download. Want to know Kalpit Thakkar better? Read Kalpit Thakkar's resume and drop Kalpit Thakkar a mail" quote: "Dig a little deeper, you might find what you are looking for" ---


Tech
SHOP
Log in
Log in
Sign up
Shopping Cart
Home
About
Contact
Categories
Camera
Iphone
Labtop
Monitor
Phone
Powerbank
Printer
Pablet
Storage
Tv
Back bag
Name
Name can not be empty.
Email
Invalid email address.
Password
Invalid password, password must be 8 characters or more.
Phone number
Invalid phone number, phone number must be 11 digits.
Address
Address can not be empty.
Birthday
Birthday can not pass the current day.
Gender
Male
Female


Short repetition from last week about conservative forces. Discussion of conditions for conservative forces and the Earth-Sun gravitional force example. Reading suggestion: Taylor sections 4.3, 4.4 and 4.8.
Potential curves and discussion of the Earth-Sun example, analytical and numerical considerations. Reading suggestions: Taylor section 4.6, 4.8 and 4.9.
If we get time, we start with harmonic oscillations and Hooke's law. Reading suggestion: Taylor section 5.1.
Energy conservation is most convenient as a strategy for addressing problems where time does not appear. For example, a particle goes from position \( x_0 \) with speed \( v_0 \), to position \( x_f \); what is its new speed? However, it can also be applied to problems where time does appear, such as in solving for the trajectory \( x(t) \), or equivalently \( t(x) \).
$$
\begin{equation}
\boldsymbol{F}(\boldsymbol{r})=-\boldsymbol{\nabla} V(\boldsymbol{r}).
\tag{1}
\end{equation}
$$
We say a force is conservative if it satisfies the following conditions:
The energy \( E \) of a given system is defined as the sum of kinetic and potential energies,
$$
E=K+V(\boldsymbol{r}).
$$
We define the potential energy at a point \( \boldsymbol{r} \) as the negative work done from a starting point \( \boldsymbol{r}_0 \) to a final point \( \boldsymbol{r} \)
$$
V(\boldsymbol{r})=-W(\boldsymbol{r}_0\rightarrow\boldsymbol{r})= -\int_{\boldsymbol{r}_0}^{\boldsymbol{r}}d\boldsymbol{r}'\boldsymbol{F}(\boldsymbol{r}').
$$
If the potential depends on the path taken between these two points there is no unique potential.
We study a classical electron which moves in the \( x \)-direction along a surface. The force from the surface is
$$
\boldsymbol{F}(x)=-F_0\sin{(\frac{2\pi x}{b})}\boldsymbol{e}_1.
$$
The constant \( b \) represents the distance between atoms at the surface of the material, \( F_0 \) is a constant and \( x \) is the position of the electron.
This is indeed a conservative force since it depends only on position and its curl is zero, that is \( -\boldsymbol{\nabla}\times \boldsymbol{F}=0 \). This means that energy is conserved and the integral over the work done by the force is independent of the path taken.
Using the work-energy theorem we can find the work \( W \) done when moving an electron from a position \( x_0 \) to a final position \( x \) through the integral
$$
W=\int_{x_0}^x \boldsymbol{F}(x')dx' = -\int_{x_0}^x F_0\sin{(\frac{2\pi x'}{b})} dx',
$$
which results in
$$
W=\frac{F_0b}{2\pi}\left[\cos{(\frac{2\pi x}{b})}-\cos{(\frac{2\pi x_0}{b})}\right].
$$
Since this is related to the change in kinetic energy we have, with \( v_0 \) being the initial velocity at a time \( t_0 \),
$$
v = \pm\sqrt{\frac{2}{m}\frac{F_0b}{2\pi}\left[\cos{(\frac{2\pi x}{b})}-\cos{(\frac{2\pi x_0}{b})}\right]+v_0^2}.
$$
The potential energy, due to energy conservation is
$$
V(x)=V(x_0)+\frac{1}{2}mv_0^2-\frac{1}{2}mv^2,
$$
with \( v \) given by the velocity from above.
We can now, in order to find a more explicit expression for the potential energy at a given value \( x \), define a zero level value for the potential. The potential is defined, using the work-energy theorem, as
$$
V(x)=V(x_0)+\int_{x_0}^x (-F(x'))dx',
$$
and if you recall the definition of the indefinite integral, we can rewrite this as
$$
V(x)=\int (-F(x'))dx'+C,
$$
where \( C \) is an undefined constant. The force is defined as the gradient of the potential, and in that case the undefined constant vanishes. The constant does not affect the force we derive from the potential.
We have then
$$
V(x)=V(x_0)-\int_{x_0}^x \boldsymbol{F}(x')dx',
$$
which results in
$$
V(x)=-\frac{F_0b}{2\pi}\left[\cos{(\frac{2\pi x}{b})}-\cos{(\frac{2\pi x_0}{b})}\right]+V(x_0).
$$
We can now define
$$
-\frac{F_0b}{2\pi}\cos{(\frac{2\pi x_0}{b})}=V(x_0),
$$
which gives
$$
V(x)=-\frac{F_0b}{2\pi}\left[\cos{(\frac{2\pi x}{b})}\right].
$$
We have defined work as the energy resulting from a net force acting on an object (or sseveral objects), that is
$$
W(\boldsymbol{r}\rightarrow \boldsymbol{r}+d\boldsymbol{r})= \boldsymbol{F}(\boldsymbol{r})d\boldsymbol{r}.
$$
If we write out this for each component we have
$$
W(\boldsymbol{r}\rightarrow \boldsymbol{r}+d\boldsymbol{r})=\boldsymbol{F}(\boldsymbol{r})d\boldsymbol{r}=F_xdx+F_ydy+F_zdz.
$$
The work done from an initial position to a final one defines also the difference in potential energies
$$
W(\boldsymbol{r}\rightarrow \boldsymbol{r}+d\boldsymbol{r})=-\left[V(\boldsymbol{r}+d\boldsymbol{r})-V(\boldsymbol{r})\right].
$$
We can write out the differences in potential energies as
$$
V(\boldsymbol{r}+d\boldsymbol{r})-V(\boldsymbol{r})=V(x+dx,y+dy,z+dz)-V(x,y,z)=dV,
$$
and using the expression the differential of a multi-variable function \( f(x,y,z) \)
$$
df=\frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz,
$$
we can write the expression for the work done as
$$
W(\boldsymbol{r}\rightarrow \boldsymbol{r}+d\boldsymbol{r})=-dV=-\left[\frac{\partial V}{\partial x}dx+\frac{\partial V}{\partial y}dy+\frac{\partial V}{\partial z}dz \right].
$$
Comparing the last equation with
$$
W(\boldsymbol{r}\rightarrow \boldsymbol{r}+d\boldsymbol{r})=F_xdx+F_ydy+F_zdz,
$$
we have
$$
F_xdx+F_ydy+F_zdz=-\left[\frac{\partial V}{\partial x}dx+\frac{\partial V}{\partial y}dy+\frac{\partial V}{\partial z}dz \right],
$$
leading to
$$
F_x=-\frac{\partial V}{\partial x},
$$
and
$$
F_y=-\frac{\partial V}{\partial y},
$$
and
$$
F_z=-\frac{\partial V}{\partial z},
$$
or just
$$
\boldsymbol{F}=-\frac{\partial V}{\partial x}\boldsymbol{e}_1-\frac{\partial V}{\partial y}\boldsymbol{e}_2-\frac{\partial V}{\partial z}\boldsymbol{e}_3=-\boldsymbol{\nabla}V(\boldsymbol{r}).
$$
And this connection is the one we wanted to show.
The net energy, \( E=V+K \) where \( K \) is the kinetic energy, is then conserved,
$$
\begin{eqnarray}
\frac{d}{dt}(K+V)&=&\frac{d}{dt}\left(\frac{m}{2}(v_x^2+v_y^2+v_z^2)+V(\boldsymbol{r})\right)\\
\nonumber
&=&m\left(v_x\frac{dv_x}{dt}+v_y\frac{dv_y}{dt}+v_z\frac{dv_z}{dt}\right)
+\partial_xV\frac{dx}{dt}+\partial_yV\frac{dy}{dt}+\partial_zV\frac{dz}{dt}\\
\nonumber
&=&v_xF_x+v_yF_y+v_zF_z-F_xv_x-F_yv_y-F_zv_z=0.
\end{eqnarray}
$$
The same proof can be written more compactly with vector notation,
$$
\begin{eqnarray}
\frac{d}{dt}\left(\frac{m}{2}v^2+V(\boldsymbol{r})\right)
&=&m\boldsymbol{v}\cdot\dot{\boldsymbol{v}}+\boldsymbol{\nabla} V(\boldsymbol{r})\cdot\dot{\boldsymbol{r}}\\
\nonumber
&=&\boldsymbol{v}\cdot\boldsymbol{F}-\boldsymbol{F}\cdot\boldsymbol{v}=0.
\end{eqnarray}
$$
Inverting the expression for kinetic energy,
$$
\begin{equation}
v=\sqrt{2K/m}=\sqrt{2(E-V)/m},
\tag{2}
\end{equation}
$$
allows one to solve for the one-dimensional trajectory \( x(t) \), by finding \( t(x) \),
$$
\begin{equation}
t=\int_{x_0}^x \frac{dx'}{v(x')}=\int_{x_0}^x\frac{dx'}{\sqrt{2(E-V(x'))/m}}.
\tag{3}
\end{equation}
$$
Note this would be much more difficult in higher dimensions, because you would have to determine which points, \( x,y,z \), the particles might reach in the trajectory, whereas in one dimension you can typically tell by simply seeing whether the kinetic energy is positive at every point between the old position and the new position.
We will now venture into a study of a system which is energy conserving. The aim is to see if we (since it is not possible to solve the general equations analytically) we can develop stable numerical algorithms whose results we can trust!
We solve the equations of motion numerically. We will also compute quantities like the energy numerically.
We start with a simpler case first, the Earth-Sun system in two dimensions only. The gravitational force \( F_G \) on the earth from the sun is
$$
\boldsymbol{F}_G=-\frac{GM_{\odot}M_E}{r^3}\boldsymbol{r},
$$
where \( G \) is the gravitational constant,
$$
M_E=6\times 10^{24}\mathrm{Kg},
$$
the mass of Earth,
$$
M_{\odot}=2\times 10^{30}\mathrm{Kg},
$$
the mass of the Sun and
$$
r=1.5\times 10^{11}\mathrm{m},
$$
is the distance between Earth and the Sun. The latter defines what we call an astronomical unit AU.
From Newton's second law we have then for the \( x \) direction
$$
\frac{d^2x}{dt^2}=-\frac{F_{x}}{M_E},
$$
and
$$
\frac{d^2y}{dt^2}=-\frac{F_{y}}{M_E},
$$
for the \( y \) direction.
Here we will use that \( x=r\cos{(\theta)} \), \( y=r\sin{(\theta)} \) and
$$
r = \sqrt{x^2+y^2}.
$$
We can rewrite
$$
F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x,
$$
and
$$
F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y,
$$
for the \( y \) direction.
We can rewrite these two equations
$$
F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x,
$$
and
$$
F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y,
$$
as four first-order coupled differential equations
$$
\frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x,
$$
$$
\frac{dx}{dt}=v_x,
$$
$$
\frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y,
$$
$$
\frac{dy}{dt}=v_y.
$$
The four coupled differential equations
$$
\frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x,
$$
$$
\frac{dx}{dt}=v_x,
$$
$$
\frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y,
$$
$$
\frac{dy}{dt}=v_y,
$$
can be turned into dimensionless equations or we can introduce astronomical units with \( 1 \) AU = \( 1.5\times 10^{11} \).
Using the equations from circular motion (with \( r =1\mathrm{AU} \))
$$
\frac{M_E v^2}{r} = F = \frac{GM_{\odot}M_E}{r^2},
$$
we have
$$
GM_{\odot}=v^2r,
$$
and using that the velocity of Earth (assuming circular motion) is
\( v = 2\pi r/\mathrm{yr}=2\pi\mathrm{AU}/\mathrm{yr} \), we have
$$
GM_{\odot}= v^2r = 4\pi^2 \frac{(\mathrm{AU})^3}{\mathrm{yr}^2}.
$$
The four coupled differential equations can then be discretized using Euler's method as (with step length \( h \))
$$
v_{x,i+1}=v_{x,i}-h\frac{4\pi^2}{r_i^3}x_i,
$$
$$
x_{i+1}=x_i+hv_{x,i},
$$
$$
v_{y,i+1}=v_{y,i}-h\frac{4\pi^2}{r_i^3}y_i,
$$
$$
y_{i+1}=y_i+hv_{y,i},
$$
The code here implements Euler's method for the Earth-Sun system using a more compact way of representing the vectors. Alternatively, you could have spelled out all the variables \( v_x \), \( v_y \), \( x \) and \( y \) as one-dimensional arrays.
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
DeltaT = 0.001
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using Euler's method
for i in range(n-1):
# Set up the acceleration
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using Euler's forward method
v[i+1] = v[i] + DeltaT*a
r[i+1] = r[i] + DeltaT*v[i]
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
#ax.set_xlim(0, tfinal)
ax.set_ylabel('y[AU]')
ax.set_xlabel('x[AU]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunEuler")
plt.show()
We notice here that Euler's method doesn't give a stable orbit. It means that we cannot trust Euler's method. In a deeper way, as we will see in homework 5, Euler's method does not conserve energy. It is an example of an integrator which is not symplectic.
Here we present thus two methods, which with simple changes allow us to avoid these pitfalls. The simplest possible extension is the so-called Euler-Cromer method. The changes we need to make to our code are indeed marginal here. We need simply to replace
r[i+1] = r[i] + DeltaT*v[i]
in the above code with the velocity at the new time \( t_{i+1} \)
r[i+1] = r[i] + DeltaT*v[i+1]
By this simple caveat we get stable orbits. Below we derive the Euler-Cromer method as well as one of the most utlized algorithms for sovling the above type of problems, the so-called Velocity-Verlet method.
Let us repeat Euler's method. We have a differential equation
$$
\begin{equation}
y'(t_i)=f(t_i,y_i)
\tag{4}
\end{equation}
$$
and if we truncate at the first derivative, we have from the Taylor expansion
$$
\begin{equation}
y_{i+1}=y(t_i) + (\Delta t) f(t_i,y_i) + O(\Delta t^2), \tag{5}
\end{equation}
$$
which when complemented with \( t_{i+1}=t_i+\Delta t \) forms
the algorithm for the well-known Euler method.
Note that at every step we make an approximation error
of the order of \( O(\Delta t^2) \), however the total error is the sum over all
steps \( N=(b-a)/(\Delta t) \) for \( t\in [a,b] \), yielding thus a global error which goes like
\( NO(\Delta t^2)\approx O(\Delta t) \).
To make Euler's method more precise we can obviously decrease \( \Delta t \) (increase \( N \)), but this can lead to loss of numerical precision. Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like.
Euler's method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at \( y_1 \) using the velocity at \( v_0 \). A simple variation is to determine \( x_{n+1} \) using the velocity at \( v_{n+1} \), that is (in a slightly more generalized form)
$$
\begin{equation}
y_{n+1}=y_{n}+ v_{n+1}+O(\Delta t^2)
\tag{6}
\end{equation}
$$
and
$$
\begin{equation}
v_{n+1}=v_{n}+(\Delta t) a_{n}+O(\Delta t^2).
\tag{7}
\end{equation}
$$
The acceleration \( a_n \) is a function of \( a_n(y_n, v_n, t_n) \) and needs to be evaluated
as well. This is the Euler-Cromer method.
Exercise: go back to the above code with Euler's method and add the Euler-Cromer method.
Let us stay with \( x \) (position) and \( v \) (velocity) as the quantities we are interested in.
We have the Taylor expansion for the position given by
$$
x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_i+O((\Delta t)^3).
$$
The corresponding expansion for the velocity is
$$
v_{i+1} = v_i+(\Delta t)a_i+\frac{(\Delta t)^2}{2}v^{(2)}_i+O((\Delta t)^3).
$$
Via Newton's second law we have normally an analytical expression for the derivative of the velocity, namely
$$
a_i= \frac{d^2 x}{dt^2}\vert_{i}=\frac{d v}{dt}\vert_{i}= \frac{F(x_i,v_i,t_i)}{m}.
$$
If we add to this the corresponding expansion for the derivative of the velocity
$$
v^{(1)}_{i+1} = a_{i+1}= a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2)=a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2),
$$
and retain only terms up to the second derivative of the velocity since our error goes as \( O(h^3) \), we have
$$
(\Delta t)v^{(2)}_i\approx a_{i+1}-a_i.
$$
We can then rewrite the Taylor expansion for the velocity as
$$
v_{i+1} = v_i+\frac{(\Delta t)}{2}\left( a_{i+1}+a_{i}\right)+O((\Delta t)^3).
$$
Our final equations for the position and the velocity become then
$$
x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_{i}+O((\Delta t)^3),
$$
and
$$
v_{i+1} = v_i+\frac{(\Delta t)}{2}\left(a_{i+1}+a_{i}\right)+O((\Delta t)^3).
$$
Note well that the term \( a_{i+1} \) depends on the position at \( x_{i+1} \). This means that you need to calculate
the position at the updated time \( t_{i+1} \) before the computing the next velocity. Note also that the derivative of the velocity at the time
\( t_i \) used in the updating of the position can be reused in the calculation of the velocity update as well.
We can now easily add the Verlet method to our original code as
DeltaT = 0.01
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up forces, air resistance FD, note now that we need the norm of the vecto
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
rabs = sqrt(sum(r[i+1]*r[i+1]))
anew = -4*(pi**2)*r[i+1]/(rabs**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
ax.set_ylabel('y[AU]')
ax.set_xlabel('x[AU]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunVV")
plt.show()
You can easily generalize the calculation of the forces by defining a function which takes in as input the various variables. We leave this as a challenge to you.
The concept of line integrals plays an important role in our discussion of energy conservation, our definition of potentials and conservative forces.
Let us remind ourselves of some the basic elements (most of you may have seen this in a calculus course under the general topic of vector fields).
We define a path integration \( C \), that is we integrate from a point \( \boldsymbol{r}_1 \) to a point \( \boldsymbol{r}_2 \). Let us assume that the path \( C \) is represented by an arc length \( s \). In three dimension we have the following representation of \( C \)
$$
\boldsymbol{r}(s)=x(s)\boldsymbol{e}_1+y(s)\boldsymbol{e}_2+z(s)\boldsymbol{e}_3,
$$
then our integral of a function \( f(x,y,z) \) along the path \( C \) is defined as
$$
\int_Cf(x,y,z)ds=\int_a^bf\left(x(s),y(s),z(s)\right)ds,
$$
where the initial and final points are \( a \) and \( b \), respectively.
With the definition of a line integral, we can in tunrn set up the theorem of independence of integration path.
Let us define \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) to be functions which are defined and continuous in a domain \( D \) in space. Then a line integral like the above is said to be independent of path in \( D \), if for every pair of endpoints \( a \) and \( b \) in \( D \) the value of the integral is the same for all paths \( C \) in \( D \) starting from a point \( a \) and ending in a point \( b \). The integral depends thus only on the integration limits and not on the path.
An expression of the form
$$
fdx+gdy+hdz,
$$
where \( f \), \( g \) and \( h \) are functions defined in \( D \), is a called a first-order differential form
in three variables.
The form is said to be exact if it is the differential
$$
du= \frac{\partial u}{\partial x}dx+\frac{\partial u}{\partial y}dy+\frac{\partial u}{\partial z}dz,
$$
of a differentiable function \( u(x,y,z) \) everywhere in \( D \), that is
$$
du=fdx+gdy+hdz.
$$
It is said to be exact if and only if we can then set
$$
f=\frac{\partial u}{\partial x},
$$
and
$$
g=\frac{\partial u}{\partial y},
$$
and
$$
h=\frac{\partial u}{\partial z},
$$
everywhere in the domain \( D \).
In vector language the above means that the differential form
$$
fdx+gdy+hdz,
$$
is exact in \( D \) if and only if the vector function (it could be a force, or velocity, acceleration or other vectors we encounter in this course)
$$
\boldsymbol{F}=f\boldsymbol{e}_1+g\boldsymbol{e}_2+h\boldsymbol{e}_3,
$$
is the gradient of a function \( u(x,y,z) \)
$$
\boldsymbol{v}=\boldsymbol{\nabla}u=\frac{\partial u}{\partial x}\boldsymbol{e}_1+\frac{\partial u}{\partial y}\boldsymbol{e}_2+\frac{\partial u}{\partial z}\boldsymbol{e}_3.
$$
If this is the case, we can state the path independence theorem which states that with functions \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) that fulfill the above exactness conditions, the line integral
$$
\int_C\left(fdx+gdy+hdz\right),
$$
is independent of path in \( D \) if and only if the differential form under the integral sign is exact in \( D \).
This is the path independence theorem.
 You can find this in any vector analysis chapter in a mathematics textbook.
We note however that the path integral from a point \( p \) to a final point \( q \) is given by
$$
\int_p^q\left(fdx+gdy+hdz\right)=\int_p^q\left(\frac{\partial u}{\partial x}dx+\frac{\partial u}{\partial y}dy+\frac{\partial u}{\partial z}dz\right)=\int_p^qdu.
$$
Assume now that we have a dependence on a variable \( s \) for \( x \), \( y \) and \( z \). We have then
$$
\int_p^qdu=\int_{s_1}^{s_2}\frac{du}{ds}ds = u(x(s),y(s),z(s))\vert_{s=s_1}^{s=s_2}=u(q)-u(p).
$$
This last equation
$$
\int_p^q\left(fdx+gdy+hdz\right)=u(q)-u(p),
$$
is the analogue of the usual formula
$$
\int_a^bf(x)dx=F(x)\vert_a^b=F(b)-F(a),
$$
with \( F'(x)=f(x) \).
We remember that a the work done by a force \( \boldsymbol{F}=f\boldsymbol{e}_1+g\boldsymbol{e}_2+h\boldsymbol{e}_3 \) on a displacemnt \( d\boldsymbol{r} \) is
$$
W=\int_C\boldsymbol{F}d\boldsymbol{r}=\int_C(fdx+gdy+hdz).
$$
From the path independence theorem, we know that this has to result in the difference between the two endpoints only. This is exact if and only if the force is the force \( \boldsymbol{F} \) is the gradient of a scalar function \( u \). We call this scalar function, which depends only the positions \( x,y,z \) for the potential energy \( V(x,y,z)=V(\boldsymbol{r}) \).
We have thus
$$
\boldsymbol{F}(\boldsymbol{r})\propto \boldsymbol{\nabla}V(\boldsymbol{r}),
$$
and we define this as
$$
\boldsymbol{F}(\boldsymbol{r})= -\boldsymbol{\nabla}V(\boldsymbol{r}).
$$
Such a force is called a conservative force. The above expression can be used to demonstrate energy conservation.
Finally we can define the criterion for exactness and independence of path. This theorem states that if \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) are continuous functions with continuous first partial derivatives in the domain \( D \), then the line integral
$$
\int_C\left(fdx+gdy+hdz\right),
$$
is independent of path in \( D \) when
$$
\frac{\partial h}{\partial y}=\frac{\partial g}{\partial z},
$$
and
$$
\frac{\partial f}{\partial z}=\frac{\partial h}{\partial x},
$$
and
$$
\frac{\partial g}{\partial x}=\frac{\partial f}{\partial y}.
$$
This leads to the curl of \( \boldsymbol{F} \) being zero
$$
\boldsymbol{\nabla}\times\boldsymbol{F}=\boldsymbol{\nabla}\times\left(-\boldsymbol{\nabla}V(\boldsymbol{r})\right)=0!
$$
A conservative force \( \boldsymbol{F} \) is a defined as the partial derivative of a scalar potential which depends only on the position,
$$
\boldsymbol{F}(\boldsymbol{r})= -\boldsymbol{\nabla}V(\boldsymbol{r}).
$$
This leads to conservation of energy and a path independent line integral as long as the curl of the force is zero, that is
$$
\boldsymbol{\nabla}\times\boldsymbol{F}=\boldsymbol{\nabla}\times\left(-\boldsymbol{\nabla}V(\boldsymbol{r})\right)=0.
$$