17.6 – ANCOVA – analysis of covariance

Introduction

Analysis of covariance (ANCOVA) is intended to help with analysis of designs with categorical treatment variables on some response (dependent) variable, but a known confounding variable is also present. Thus, the researcher is also likely to know of additional ratio scale variables that covary with the response variable and, moreover, must be included in the experimental design in some way.

Take for example the well-known relationship between body size and whole-animal metabolic rate as measured by rates of carbon dioxide production or rates of oxygen consumption for aerobic organisms. We may be interested in how addition or blocking of stress hormones affects resting metabolism; we may be interested in comparing men and women for activity metabolism, and so on. We’d like to know if the regressions were the same (eg, metabolic rate covaried with body mass in the same way — that is, the slope of the relationship was the same).

This situation arises frequently in biology. For example, we might want to know if male and female birds have different mean field metabolic rates, in which case we might be tempted to use a one-way ANOVA or t-test (since one factor with two levels). However, if males and females also differ for body size, then any differences we might see in metabolic rate could be due to differences in metabolic rate are confounded by differences in the covariable body size. We already discussed one approach to correction: calculate a ratio. Thus, a logical approach to correcting or normalizing for the covariation would be to divide body mass (units of kilograms) into metabolic rate (eg, volume of oxygen, O2, consumed), and make comparisons, say, among different species, on mass-specific trait (\frac{ml \ O_{2}}{hours\cdot mass}). However, because the regression between mass and metabolic rate is allometric, i.e., not equal to one, the ratio does not, in fact normalize for body mass.  We made this point in Chapter 6.2, and remarked that analysis of covariance ANCOVA was a solution.

ANCOVA allows you to test for mean differences in traits like metabolic rate between two or more groups, but only after first accounting for covariation due to another variable (eg, body size). However, ANCOVA makes the assumption that relationship between the covariable and the response variable is the same in the two groups. This is the same as saying that the regression slopes are the same. We discussed how to use t-test to test hypothesis of equal slopes between regression models in Chapter 17.5, but a more elegant way is to include this in your model.

Example

We return to our sample of 13 tadpoles (Rana pipiens), hatched in the laboratory. I’ve repeated the data set in this page, scroll or click here.

Our linear model was \dot{V} O_{2}\sim \text{Body.mass} and a scatterplot of the data set is shown in Figure 1 (a repeat of Figure 3 from Chapter 17.5, but now points identified to developmental group).

Oxygen consumption of tadpoles by body weight

Figure 1. Scatterplot of oxygen consumption by R. pipiens tadpoles vs body mass (g) by developmental group (Gosner stages I or II).

The project looked at whether metabolism as measured by oxygen consumption was consistent across two developmental stages. Metamorphosis in frogs and other amphibians represents profound reorganization of the organism as the tadpole moves from water to air. Thus, we would predict some cost as evidenced by change in metabolism associated with later stages of development. Figure 2 shows a box plot of tadpole oxygen consumption by Gosner (1960) developmental stage (Figure 2 is a repeat of Figure 4 Chapter 17.5).

boxplot oxygen consumption of tadploes by developmental stage

Figure 2. Copy of Figure 4, Chapter 17.5; boxplot of oxygen consumption of R. pipiens tadpoles by Gosner developmental stages.

Looking at Figure 2 we see a trend consistent with our prediction; developmental stage may be associated with increased metabolism. However, older tadpoles also tend to be larger, and the plot in Figure 2 does not account for that. Thus, body mass is a confounding variable in this example. There are several options for analysis here (eg, ANCOVA), but one way to view this is to compare the slopes for the two developmental stages. While this test does not compare the means, it does ask a related question: is there evidence of change in rate of oxygen consumption relative to body size between the two developmental stages? The assumption that the slopes are equal is a necessary step for conducting the ANCOVA.

So, divide the data set (Table 1) into two groups by developmental stage (12 tadpoles could be assigned to one of two developmental stages; one was at a lower Gosner stage than the others and so is dropped from the subset.

Table 1. \text{Body.mass} and \dot{V} O_{2} of 12 R. pipiens tadpole larvae by Gosner developmental stage.

Gosner stage I

\text{Body.mass} \dot{V} O_{2}
1.76 109.41
1.88 329.06
1.95 82.35
2.13 198
2.26 607.7

Gosner stage II

Body mass \dot{V} O_{2}
2.28 362.71
2.35 556.6
2.62 612.93
2.77 514.02
2.97 961.01
3.14 892.41
3.79 976.97

The slopes and standard errors we obtained in Chapter 17.5 were

Gosner Stage I Gosner stage II
slope 750.0 399.9
standard error of slope 444.6 111.2

Make a plot (Figure 3).

scatterplot oxygen consumption of anuran tadpoles by body mass

Figure 3. Scatterplot with best-fit regression lines of \dot{V} O_{2} by \text{Body.mass} for Gosner State I (closed circle, solid line) and Gosner Stage II (open circle, dashed line) R. pipiens tadpoles.

R code for plot in Figure 3.

#Used Rcmdr scatterplot(), then modified code
scatterplot(VO2~Body.mass | Gosner, regLine=FALSE, smooth=FALSE, 
boxplots=FALSE, xlab="Body mass (g)", ylab="Oxygen consumption (ml/h)", 
main="", cex=1.4, cex.axis=1.5, cex.lab=1.5, pch=c(19,19), by.groups=TRUE,
col=c("blue","red"), grid=FALSE, 
legend=list(coords="bottomright"), data=Tadpoles)
#Get regression equations for groups, subset by Gosner
abline(lm(VO2~Body.mass, data=Stage01), lty=1, lwd=2, col="blue")
abline(lm(VO2~Body.mass, data=Stage02), lty=1, lwd=2, col="red")

Returning to the important question, are the two slopes statistically indistinguishable (Ho: bI = bII), where bI is the slope for the Gosner Stage I subset and bII is the slope for the Gosner Stage II subset? We look at the plot, and since the lines cross, we tend to see a difference. Of course, we need to consider that our perception of slope differences may simply be chance, especially because the sample size is small. Proceed to test.

R code

The ANCOVA is a new ANOVA model where the factor variables are adjusted or corrected for the effects of the continuous variable.

R code for ANCOVA example, crossed or interaction model.

tadpole.1 <- lm(VO2 ~ Body.mass*Gosner, data=example.Tadpole)
summary(tadpole.1)
Anova(tadpole.1, type="II")

Output

summary(tadpole.1)

Call:
lm(formula = VO2 ~ Body.mass * Gosner, data = example.Tadpole)

Residuals:
    Min      1Q   Median      3Q      Max 
-167.80 -117.93    13.81   94.66   214.65

Coefficients:
                       Estimate   Std. Error   t value    Pr(>|t|) 
(Intercept)             -1231.6       783.0     -1.573    0.1544 
Body.mass                 750.0       390.7      1.919    0.0912 .
Gosner[T.II]              790.4       859.2      0.920    0.3845 
Body.mass:Gosner[T.II]   -350.1       409.5     -0.855    0.4174 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 155.8 on 8 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.821, Adjusted R-squared: 0.7539 
F-statistic: 12.23 on 3 and 8 DF, p-value: 0.002336

The coefficients for the first factor (GII) and then the differences in the coefficient for the second factor. You can just add the second coefficient to the first so they’re on the same scale.

Anova Table (Type II tests)

Response: VO2
                  Sum Sq   Df    F value     Pr(>F) 
Body.mass         330046    1    13.6030   0.006146 **
Gosner              5630    1     0.2321   0.642908 
Body.mass:Gosner   17736    1     0.7310   0.417423 
Residuals 194102 8 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Suggests interaction is not significant, i.e., the slopes are not different.

We can then proceed to check to see if the intercepts are different, now that we’ve confirmed no significant difference in slope.

R code for ANCOVA as additive model

tadpole.2 <- lm(VO2 ~ Body.mass + Gosner, data=example.Tadpole)
summary(tadpole.2)
Anova(tadpole.2, type="II")

Output

> summary(tadpole.2)

Call:
lm(formula = VO2 ~ Body.mass + Gosner, data = example.Tadpole)

Residuals:
    Min       1Q    Median      3Q      Max 
-163.12  -125.53    -20.27   83.71   228.56

Coefficients:
                   Estimate   Std. Error     t value    Pr(>|t|) 
(Intercept)         -595.37       239.87      -2.482     0.03487 * 
Body.mass            431.20       115.15       3.745     0.00459 **
Gosner[T.II]          64.96       132.83       0.489     0.63648 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 153.4 on 9 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8047, Adjusted R-squared: 0.7613 
F-statistic: 18.54 on 2 and 9 DF, p-value: 0.0006432

Anova(tadpole.2, type="II")
Anova Table (Type II tests)

Response: VO2
              Sum Sq    Df    F value      Pr(>F) 
Body.mass     330046     1    14.0221   0.004593 **
Gosner          5630     1     0.2392   0.636482 
Residuals     211839     9 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note 1. There is no test of interaction in the added model. This model would be appropriate IF the slopes are equal.

Instead of the interaction, or as added, try a nested model, with body mass nested within stage.

tadpole.3 <- lm(VO2 ~ Body.mass/Gosner, data=example.Tadpole)

summary(tadpole.3)

Call:
lm(formula = VO2 ~ Body.mass/Gosner, data = example.Tadpole)

Residuals:
    Min       1Q    Median      3Q      Max 
-168.66  -131.14    -20.28   90.33   225.36

Coefficients:
                 Estimate      Std. Error    t value Pr(>|t|) 
(Intercept)           -575.10      319.51     -1.800   0.1054 
Body.mass              423.65      162.50      2.607   0.0284 *
Body.mass:Gosner[T.II]  21.95       63.73      0.344   0.7384 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 154.4 on 9 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8021, Adjusted R-squared: 0.7581 
F-statistic: 18.24 on 2 and 9 DF, p-value: 0.0006823

 

Gets the true coefficient (nested lm() version).

The two test different hypotheses.

lm(VO2 ~ Body.mass * Gosner) tests whether or not the regression has a nonzero slope.

lm(VO2 ~ Body.mass */Gosner) test whether or not the slopes and intercepts from different factors are statistically significant.

Questions

  1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
  2. An OLS approach was used for the analysis of tadpole oxygen consumption body mass. Consider the RMA approach — would that be a more appropriate regression model? Explain why or why not.
  3. Consider an experiment in which you plan to administer a treatment that has a carry-over effect. For example, Compare and contrast “crossed” and “nested” designs.
  4. True or False. The nested design option for the ANCOVA assumes the slopes for the two groups of tadpoles for the regression line of \dot{V} O_{2} by \text{Body.mass} are equal. Explain your choice.
  5. Metabolic rates like oxygen consumption over time are well-known examples of allometric relationships. That is, many biological variables (eg, \dot{V} O_{2} is related as aMb, where M is body mass, slope b is scaling exponent), and best evaluated on log-log scale. Repeat the analysis above on log10-transformed \dot{V} O_{2} and \text{Body.mass} for
    • crossed design (eg, tadpole.1 model)
    • added design (eg, tadpole.2 model)
    • nested design (eg, tadpole.3 model)
  6. Create the plot and add the fitted lines from crossed design to the plot.

Note 2. About log-transform of a variable. The most straight-forward tact is to create two new variables. For example,

lgVO2 <- log10(VO2)

Another option is to transform the variables within the call to lm() function. For example, try

lm(log10(VO2) ~ log10(Body.mass ), data=example.Tadpole)

Hint: don’t forget to attach your data set to avoid having to call the variable as, for example, example.Tadpole$VO2

Quiz Chapter 17.6

ANCOVA - analysis of covariance

Data sets

Oxygen consumption, \dot{V}O_{2}, of Anuran tadpoles, dataset=example.Tadpole

GosnerBody massVO2
NA1.46170.91
I1.76109.41
I1.88329.06
I1.9582.35
I2.13198
II2.28362.71
I2.26607.7
II2.35556.6
II2.62612.93
II2.77514.02
II2.97961.01
II3.14892.41

Gosner refers to Gosner (1960), who developed a criteria for judging metamorphosis staging.


Chapter 17 contents

17.7 – Regression model fit

Introduction

In Chapter 17.5 and 17.6 we introduced the example of tadpoles body size and oxygen consumption. We ran a simple linear regression, with the following output from R

RegModel.1 <- lm(VO2~Body.mass, data=example.Tadpole)

summary(RegModel.1)

Call:
lm(formula = VO2 ~ Body.mass, data = example.Tadpole)

Residuals:
    Min      1Q    Median       3Q       Max 
-202.26 -126.35     30.20    94.01    222.55

Coefficients:
                Estimate     Std. Error    t value    Pr(>|t|) 
(Intercept)      -583.05         163.97     -3.556     0.00451 ** 
Body.mass         444.95          65.89      6.753   0.0000314 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 145.3 on 11 degrees of freedom
Multiple R-squared: 0.8057, Adjusted R-squared: 0.788 
F-statistic: 45.61 on 1 and 11 DF, p-value: 0.00003144

You should be able to pick out the estimates of slope and intercept from the table (intercept was -583 and slope was 445). Additionally, as part of your interpretation of the model, you should be able to report how much variation in VO2 was explained by tadpole body mass (coefficient of determination, R^2, was 0.81, which means about 81% of variation in oxygen consumption by tadpoles is explained by knowing the body mass of the tadpole.

What’s left to do? We need to evaluate how well our model fits the data, i.e., we evaluate regression model fit. By fit we mean, how well does our model agree to the raw data? A poor fit and the model predictions are far from the raw data; a good fit, the model accurately predicts the raw data. In other words, if we draw a line through the raw data, a good fit line crosses through most of the points — a poor fit the line rarely passes through the cloud of points.

How to judge fit? This we can do by evaluating the error components relative to the portion of the model that explains the data. Additionally, we can perform a number of diagnostics of the model relative to the assumptions we made to perform linear regression. These diagnostics form the subject of Chapter 17.8. Here, we ask how well does the model,

    \begin{align*} \dot{V}O_{2} = b_{0}+b_{1}\left ( Body \ mass \right ) \end{align*}

fit the data?

Model fit statistics

The second part of fitting a model is to report how well the model fits the data. The next sections apply to this aspect of model fitting. The first area to focus on is the magnitude of the residuals: the greater the spread of residuals, the less well a fitted line explains the data.

In addition to the output from lm() function, which focuses on the coefficients, we typically generate the ANOVA table also.

Anova(RegModel.1, type="II")
Anova Table (Type II tests)

Response: VO2
                Sum Sq   Df    F value       Pr(>F) 
Body.mass       962870    1     45.605   0.00003144 ***
Residuals       232245   11 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Standard error of regression

S, the Residual Standard Error (aka Standard error of regression), is an overall measure to indicate the accuracy of the fitted line: it tells us how good the regression is in predicting the dependence of response variable on the independent variable. A large value for S indicates a poor fit. One equation for S is given by

    \begin{align*} S = \sqrt{\frac{SS_{residual}}{n-2}} \end{align*}

In the above example, S = 145.3 (underlined, bold in regression output above). We can see how if SS_{residual} is large, S will be large indicating poor fit of the linear model to the data. However, by itself S is not of much value as a diagnostic as it is difficult to know what to make of 145.3, for example. Is this a large value for S? Is it small? We don’t have any context to judge S, so additional diagnostics have been developed.

Coefficient of determination

R^2, the coefficient of determination, is also used to describe model fit. R^2 can take on values between 0 and 1 (0% to 100%). A good model fit has a high R^2 value. In our example above, R^2 = 0.8057 or 80.57%. One equation for R^2 is given by

    \begin{align*} R^2 = \frac{SS_{regression}}{SS_{total}} \end{align*}

A value of R^2 close to 1 means that the regression “explains” nearly all of the variation in the response variable, and would indicate the model is a good fit to the data. As we mentioned in Chapter 17.2, the coefficient of determination, R^2, is the squared value of r, the product moment correlation.

Adjusted R-squared

Before moving on we need to remark on the difference between R^2 and adjusted R^2. For Simple Linear Regression there is but one predictor variable, X; for multiple regression there can be many additional predictor variables. Without some correction,  R^2 will increase with each additional predictor variables. This doesn’t mean the model is more useful, however, and in particular, one cannot compare  R^2 between models with different numbers of predictors. Therefore, an adjustment is used so that the coefficient of determination remains a useful way to assess how reliable a model is and to permit comparisons of models. Thus, we have the Adjusted \bar{R}^2, which is calculated as

    \begin{align*} \bar{R}^2 = 1 - \frac{SS_{residual}}{SS_{total}} \cdot \frac{DF_{total}}{DF_{residual}} \end{align*}

In our example above, Adjusted R 2 = 0.3806 or 38.06%.

Which should you report? Adjusted R^2 ,because it is independent of the number of parameters in the model.

Both \bar{R}^2 and S are useful for regression diagnostics, a topic which we will discuss next (Chapter 17.8).

Questions

  1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
  2. As you may know, the US census is used to reapportion congress and also to determine the number of electoral college votes. In honor of the election for US President that’s just days away, in the next series of questions in this Chapter and subsequent sections of Chapter 17 and 18, I’ll ask you to conduct a regression analysis on the electoral college. For starters, make the regression of Electoral votes on 2010 census population. (Ignore for now the other columns, just focus on POP_2019 and Electoral.) Report the
    • regression coefficients (slope, intercept)
    • percent of the variation in electoral college votes explained by the regression ( R^2 ).
  3. Make a scatterplot and add the regression line to the plot

Quiz Chapter 17.7

Regression model fit

Data set

USA population year 2010 and year 2019 with Electoral College counts

StateRegionDivisionPOP_2010POP_2019Electoral
AlabamaSouthEast South Central477973649031859
AlaskaWestPacific7102317315453
ArizonaWestMountain6392017727871711
ArkansasSouthWest South Central291591830178046
CaliforniaWestPacific372539563951222355
ColoradoWestMountain502919657587369
ConnecticutNortheastNew England357409735652877
DelawareSouthSouth Atlantic8979349828953
District of ColumbiaSouthSouth Atlantic6017237057493
FloridaSouthSouth Atlantic188013102147773729
GeorgiaSouthSouth Atlantic96876531061742316
HawaiiWestPacific136030114158724
IdahoWestMountain156758217870654
IllinoisMidwestEast North Central128306321267182120
IndianaMidwestEast North Central6483802673221911
IowaMidwestWest North Central304635531550706
KansasMidwestWest North Central285311829133146
KentuckySouthEast South Central433936744676738
LouisianaSouthWest South Central453337246487948
MaineNortheastNew England132836113442124
MarylandSouthSouth Atlantic5773552604568010
MassachusettsNortheastNew England6547629689250311
MichiganMidwestEast North Central9883640988363516
MinnesotaMidwestWest North Central5303925563963210
MississippiSouthEast South Central296729729761496
MissouriMidwestWest North Central5988927613742810
MontanaWestMountain98941510687783
NebraskaMidwestWest North Central182634119344085
NevadaWestMountain270055130801566
New HampshireNortheastNew England131647013597114
New JerseyNortheastMid-Atlantic8791894888219014
New MexicoWestMountain205917920968295
New YorkNortheastMid-Atlantic193781021945356129
North CarolinaSouthSouth Atlantic95354831048808415
North DakotaMidwestWest North Central6725917620623
OhioMidwestEast North Central115365041168910018
OklahomaSouthWest South Central375135139569717
OregonWestPacific383107442177377
PennsylvaniaNortheastMid-Atlantic127023791280198920
Rhode IslandNortheastNew-England105256710593614
South CarolinaSouthSouth-Atlantic462536451487149
South DakotaMidwestWest-North-Central8141808846593
TennesseeSouthEast-South-Central6346105682917411
TexasSouthWest-South-Central251455612899588138
UtahWestMountain276388532059586
VermontNortheastNew-England6257416239893
VirginiaSouthSouth-Atlantic8001024853551913
WashingtonWestPacific6724540761489312
West VirginiaSouthSouth-Atlantic185299417921475
WisconsinMidwestEast-North-Central5686986582243410
WyomingWestMountain5636265787593

Chapter 17 contents

17.3 – Estimation of linear regression coefficients

Introduction

In discussing correlations we made the distinction between inference, testing the statistical significance of the estimate, and the process of getting the estimate of the parameter itself. Estimating parameters is possible for any data set; whether or not the particular model is a good and useful model is another matter. Statisticians speak about the fit of a model… that a model with one or more independent predictor variables explains a substantial amount of the variation in the dependent variable, that it describes the relationship between the predictors and the dependent variable without bias.

Note 1. We introduced the concept of statistical fit in our discussion of Chi-square, goodness of fit, Chapter 9.1.

A number of tools have been developed to assess model fit. For starters, I’ll list just two ways you can approach whether or not a linear model fits your data ore requires some intervention on your part.

Assess fit of a linear model

Recall our R output from the regression of Number of Matings on Body Mass from the bird data set. We used the linear model function.

LinearModel.1 <- lm(Matings ~ Body.Mass, data=bird_matings)
summary(LinearModel.1)
Call: lm(formula = Matings ~ Body.Mass, data = bird_matings)
Residuals:       Min       1Q    Median       3Q      Max
            -2.29237 -1.34322 .-0.03178..1.33792 .2.70763
Coefficients:
             Estimate  Std. Error  t value  Pr(>|t|)
(Intercept)   -8.4746      4.6641   -1.817    0.1026 
Body.Mass      0.3623      0.1355    2.673    0.0255 * 
--- 
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.776 on 9 degrees of freedom 
Multiple R-squared: 0.4425, Adjusted R-squared: 0.3806 
F-statistic: 7.144 on 1 and 9 DF, p-value: 0.02551

Request R to print the ANOVA table.

Anova(RegModel.1, type="II")
Anova Table (Type II tests)
Response: Matings 
            Sum Sq   Df   F value    Pr(>F) 
Body.Mass   22.528    1    7.1438   0.02551 * 
Residuals   28.381    9

With a little rounding we have the following statistical model

    \begin{align*} Y_{i} = -8.5 + 0.36 \cdot X_{i} \end{align*}

and in English, Number of matings equals Body.mass multiplied by 0.36 then subtract 8.5; the intercept was -8.5, the slope was 0.36.

Note 2. The results of the regression analyses were stored in the object called “LinearModel.1“. This is a nice feature of Rcmdr — it automatically provides an object name for you. Note that with each successive run of the linear model function via Rcmdr that it will change the object name by adding numbers successively. For example, after LinearModel.1 the next run of lm() in Rcmdr will automatically be called “LinearModel.2” and so on. In your own work you may specify the names of the objects directly or allow Rcmdr to do it for you, but do keep track of the object names!

From the R output we see that the estimate of the slope was +0.36, statistically different from zero (p = 0.025). The intercept was -8.5, but not statistically significant (p = 0.103), which means the intercept may be zero.

As a general rule, if you make an estimate of a parameter or coefficient, then you should provide a confidence interval of the form.

estimate + critical value X standard error of the estimate

Note 3. Reminder: Approximate 95% CI can be obtained by + twice the standard error for the coefficient.

For confidence interval of regression slope we have a couple of options in R

Option 1.

confint(LinearModel.1, 'Body.Mass', level=0.95)

               2.5 %    97.5 %
Body.Mass 0.05565899 0.6689173

Option 2.

Goal: Extract the coefficients from the output of the linear model and calculate the approximate SE with nine degrees of freedom. This is the big advantage of saving output from functions as objects. Typically, much more information is about the results are available, and, additionally, can be retrieved for additional use. Extracting coefficients from the objects is the best option, but does come with a learning curve.  Let’s get started.

First, what information is available in the linear model output beyond the default information? To find out, use the names() function

names(LinearModel.1)

R output

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"

Another way is to use the summary() function call.

summary(LinearModel.1)$coefficients
              Estimate Std. Error   t value   Pr(>|t|)
(Intercept) -8.4745763  4.6640856 -1.816986 0.10259256
Body.Mass    0.3622881  0.1355472  2.672781 0.02550595

How can we get just the standard error for the slope? Note that the estimates are reported in a 2X4 matrix like so

1,1 1,2 1,3 1,4
2,1 2,2 2,3 2,4

Therefore, to get the standard error for the slope we identify that it is stored in cell 2,2 of the matrix and we write

summary(LinearModel.1)$coefficients[2,2]

which returns

[1] 0.1355472

Let’s use this information to calculate confidence intervals

slp=summary(LinearModel.1)$coefficients[2,1]
slpErrs=summary(LinearModel.1)$coefficients[2,2] 
slp + c(-1,1)*slpErrs*qt(0.975, 9)

where qt() is the quantile function for the t distribution and “9” is the degrees of freedom from the regression. Results follow

coef + c(-1,1)*errs*qt(0.975, 9)
[1] 0.05565899 0.66891728

And for the intercept

int=summary(LinearModel.1)$coefficients[1,1]
intErrs=summary(LinearModel.1)$coefficients[1,2] 
int + c(-1,1)*intErrs*qt(0.975, 9)

Results

int + c(-1,1)*intErrs*qt(0.975, 9)
[1] -19.025471   2.076318

In conclusion, part of fitting a model includes reporting the estimates of the coefficients (model parameters). And, in general, when estimation is performed, reporting of suitable confidence intervals are expected.

Extract additional statistics from R’s linear model function

The summary() function is used to report the general results from ANOVA and linear model function output in R software, but additional functions can be used to extract the rest of the output, e.g., coefficient of determination. To complete our example of extracting information from the summary() function, we next turn to summary.lm() function to see what is available.

At the R prompt type and submit

summary.lm(LinearModel.1)

returns the following R output

Call:
lm(formula = Matings ~ Body.Mass, data = bird_matings)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.29237 -1.34322 -0.03178  1.33792  2.70763 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -8.4746     4.6641  -1.817   0.1026  
Body.Mass     0.3623     0.1355   2.673   0.0255 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.776 on 9 degrees of freedom
Multiple R-squared:  0.4425,    Adjusted R-squared:  0.3806 
F-statistic: 7.144 on 1 and 9 DF,  p-value: 0.02551

Looks exactly like the output from summary(). Let’s look at what is available in the summary.lm() function

names(summary.lm(LinearModel.1))
 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

We see some information we got from summary(), e.g., “coefficients”. If we interrogate the name coefficients like so

summary.lm(LinearModel.1)$coefficients

we get

summary.lm(LinearModel.1)$coefficients
              Estimate Std. Error   t value   Pr(>|t|)
(Intercept) -8.4745763  4.6640856 -1.816986 0.10259256
Body.Mass    0.3622881  0.1355472  2.672781 0.02550595

which, again, is a 2X4 matrix (see above)

So to get the standard error for the slope we identify that it is stored in cell 2,2 of the matrix and call it LinearModel.1$coefficients[2,2].

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.

2. What does a small standard error of the slope for a SLR model tell us?

3. For slope, intercept, for just about any statistical estimate, why are we obligated to provide a confidence interval?

Quiz Chapter 17.3

Estimation of linear regression coefficients


Chapter 17 contents

17 – Linear Regression

Introduction

Regression is a toolkit for developing models of cause and effect between one ratio scale data type dependent response variables, and one (simple linear regression) or more or more (multiple linear regression) ratio scale data type independent predictor variables. By convention the dependent variable(s) is denoted by Y, the independent variable(s) represented by X1, X2, \cdots Xn for n independent variables. Like ANOVA, linear regression is simply a special case of the general linear model, first introduced in Chapter 12.7.

Components of a statistical model

Regression statistical methods return model estimates of the intercept and slope coefficients, plus statistics of regression fit (eg, R2, aka “R-squared,” the coefficient of determination).

Chapter 17.1 – 17.9 cover the simple linear model

    \begin{align*} Y_{i} = \alpha + \beta X_{i} + \epsilon_{i} \end{align*}

Chapter 18.1 – 18.5 cover the multiple regression linear model

    \begin{align*} Y_{i} = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \cdots + \beta_{n}X_{n} + \epsilon_{i} \end{align*}

where α or β0 represent the Y-intercept and β or β1, β2, … βn represent the regression slopes.

Regression and correlation test linear hypotheses

We state that the relationship between two variables is linear (the alternative hypothesis) or it is not (the null hypothesis). The difference? Correlation is a test of linear association (are variables correlated, we ask?), imply possible causation, but are not sufficient evidence for causation: we do not imply that one variable causes another to vary, even if the correlation between the two variables is large and positive, for example. Correlations are used in statistics on data sets not collected from explicit experimental designs incorporated to test specific hypotheses of cause and effect.

Linear regression, however, is to cause and effect as correlation is to association. With regression and ANOVA, we are indeed making a case for a particular understanding of the cause of variation in a response variable: modeling cause and effect is the goal. Regression, ANOVA, and other general linear models are designed to permit the statistician to control for the effects of confounding variables provided the causal variables themselves are uncorrelated.

When to use correlation and when to apply linear regression to data set? For two ratio scale variables, one can always apply either approach. Correlation is always appropriate if both variables are measured in the course of the study whereas a regression modeling approach would be implied where one variable was the result of manipulation by the researcher.

Assumptions of linear regression

The key assumption in linear regression is that a straight line indeed is the best fit of the relationship between dependent and independent variables. The additional assumptions of parametric tests (Chapter 13) also hold. In Chapter 18 we conclude with an extension of regression from one to many predictor variables and the special and important topic of correlated predictor variables or multicollinearity.

Build a statistical model, make predictions

In our exploration of linear regression we begin with simple linear regression, also called ordinary least squares regression, starting with one predictor variable. Practical aspects of model diagnostics are presented. Regression may be used to describe or to provide a predictive statistical framework. In Chapter 18 we conclude with an extension of regression from one to many predictor variables. We conclude with a discussion of model selection. Throughout, use of Rcmdr and R have multiple ways to analyze linear regression models are presented; we will continue to emphasize the general linear model approach, but note that use of linear model in Rcmdr provides a number of default features that are conveniently available.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.

Suggested homework

Homework 11: Simple linear regression.

Quiz

In addition to the quiz on this page, a total of 81 questions among the several subchapters, a mix of true or false and multiple choice question format.

Quiz Chapter 17

Linear Regression

References

Linear regression is a huge topic; references I include are among my favorite on the subject, but are only a small and incomplete sampling. For simplicity, I merged references for Chapter 17 and Chapter 18 into one page at References and suggested readings (Ch17 & 18)


Chapter 17 contains