Lecture 31 (lab 8)—March 7, 2006

What was covered?

• R topics
  • the * and : notation for designating interaction terms in regression models
  • obtaining the variance-covariance matrix of regression parameters
  • defining dummy regressors to test contrasts of interest
• regression with categorical variables
• analysis of covariance
• testing linear combinations of parameters
• constructing orthogonal contrasts

R functions and commands demonstrated

• anova is used for carrying out partial F-tests, likelihood ratio tests, and the construction of ANOVA tables
• class identifies the class from which an object inherits its methods
• contrasts displays the dummy coding scheme defined for a factor. It can also be used to specify new dummy coding schemes.
• contr.treatment can be used to redefine the contrasts for a factor variable
• factor is used to declare a categorical variable as a factor. R in turn creates a set of dummy regressors for use in regression models
• ls.diag is a function that returns a number of quantities useful in diagnosing problems with least squares models
• pt is the cumulative distribution function for the t-distribution
• vcov returns the scaled variance-covariance matrix of the parameter estimates of a regression model
• %*% is the operator for matrix multiplication

R function options

• corr= (argument to summary) when set to FALSE will suppress the display of parameter correlations
• labels= (argument to factor) can be used to specify labels for factor levels
• base= (argument to contr.treatment) can be used to reassign the baseline level for treatment contrast definition

The bird-bat data set

• We examine the bird-bat data set of Speakman and Racey (1991) that we've been discussing in lecture. Once again the basic problem is to compare the energy usage of different vertebrates in flight while controlling for systematic biases that exist among the different groups. The method we employ is analysis of covariance and we use it to illustrate the different ways that categorical variables can be entered into a regression model.

> #obtain data
> bats<-read.table( 'http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab8/birdbat.txt', header=TRUE)
> bats
                     species  mass type energy  n
1           Pteropus gouldii 779.0    1  43.70  1
2     Pteropus poliocephalus 628.0    1  34.80  4
3    Hypsignathus monstrosus 258.0    1  23.30  1
4             Eidolon helvum 315.0    1  22.40  1
5        Meliphaga virescens  24.3    2   2.46  1
6    Melipsittacus undulatus  35.0    2   3.93  1
7           Sturnus vulgaris  72.8    2   9.15  6
8            Falco spaverius 120.0    2  13.80  3
9          Falco tinnunculus 213.0    2  14.60  6
10         Corvus ossifragus 275.0    2  22.80  2
11           Larus atricilla 370.0    2  26.20  2
12             Columba livia 384.0    2  25.90  4
13             Columba livia 442.0    2  29.50  6
14             Columba livia 412.0    2  43.70  6
15             Columba livia 330.0    2  34.00  6
16        Corvus crytoleucos 480.0    2  27.80  2
17     Phyllostomas hastatus  93.0    3   8.83  1
18          Plecotus auritus   8.0    3   1.35  3
19 Pipistrellus pipistrellus   6.7    3   1.12 17
20          Plecotus auritus   7.7    3   1.02  3

> attach(bats)

• I attach the data set so that we can access the variables directly without having to reference the data set by name each time.
• The variable n contained in the data frame is the sample size on which the reported mass and energy values are based. When n > 1 the reported values are means.
  • The use of averages here, while not invalid, severely limits the statistical power of this study. The aggregation was done to avoid dealing with the hierarchical structure of these data, i.e., the fact that the individual bats do not form a random sample but instead form what might best be described as a cluster sample. Near the end of the course we'll learn how such hierarchical data can be analyzed directly without aggregation.
• The variable type refers to the type of animal—non-echolocating bat, bird, or echo-locating bat. Observe that it is numeric. The numbers just refer to the categories and should be treated as labels. R doesn't know this and so when the data set was read in with the read.table function, R treated type as an ordinary numeric variable.

> class(type)
[1] "integer"

• As we've seen in lecture, when a categorical variable is used in regression it needs to be represented as a set of dummy regressors. If the levels of the categories of type had been entered as character values in the raw data set, R would have done this conversion for us automatically when the data were read in. Since the values are numbers, we'll have to do this ourselves using the factor function.

> type.f<-factor(type)
> class(type.f)
[1] "factor"
> type
[1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3
> type.f
[1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3
Levels: 1 2 3

• Observe that after declaring type to be a factor, R recognizes that the values represent levels of a categorical variable rather than quantitative values. To see the dummy coding that was used, use the contrasts function.

> contrasts(type.f)
  2 3
1 0 0
2 1 0
3 0 1

• By default R has selected the first category numerically to be the baseline value. We can obtain more informative labeling of the categories with the labels option of the factor function.

> type.f<-factor(bats$type,labels=c('nonecho bats','birds','echo bats'))
> levels(type.f)
[1] "nonecho bats" "birds" "echo bats"
> contrasts(type.f)
             birds echo bats
nonecho bats     0         0
birds            1         0
echo bats        0         1

• The baseline level of a factor can be redefined. For a treatment contrast (which is what we've been using) this can be done with the contr.treatment function. This function is applied to the levels of our factor variable and the base= argument is used to redefine which level is the baseline. In the code below I change the baseline group to be the second group, birds.

> contrasts(type.f)<-contr.treatment(levels(type.f),base=2)
> contrasts(type.f)
             nonecho bats echo bats
nonecho bats            1         0
birds                   0         0
echo bats               0         1

• We can reset the contrasts back to its earlier setting by assigning the value NULL to the contrasts.

> contrasts(type.f)<-NULL
> contrasts(type.f)
             birds echo bats
nonecho bats     0         0
birds            1         0
echo bats        0         1

Examining the data

• At issue is whether echolocation is an energy-intensive activity. Since bats typically echolocate while flying, we need to tease apart energy use due to flying from energy use due to echolocating. One way to address this is to determine if non-echolocating flyers have a lower energy expenditure than echolocators do. Before doing any formal statistical analysis it is always a good idea to examine the relationships graphically. Observe that because of the aggregation used by the authors, we don't have a lot of data with which to work.

> table(type.f)
type.f
nonecho bats birds echo bats
           4    12         4

• When there is so little data we should plot it directly rather than summarize it using means, error bars, and the like. A compromise display that I like is to use box plots to summarize the distributional patterns but also superimpose the raw data, jittered to minimize overlap, on top of the box plots. The box plot in Fig. 1 compares the energy use during flight for our three groups of flyers.

> boxplot(energy~type.f, ylab='Rate of energy use')
> points(jitter(type), energy, pch=16, col=2, cex=.7)

• The jittering is done randomly so it may take a couple of tries to get an optimal display. Observe that the jitter function was applied to the x-variable so that the jittering occurs in the x-direction (rather than the y-direction) so as to not distort the relationship of interest.
• We can formally test for a difference among the three groups by fitting a regression model and then applying the anova function to the model that is produced.

> anova(lm(energy~type.f))
Analysis of Variance Table

Response: bats$energy
          Df  Sum Sq Mean Sq F value   Pr(>F)
type.f     2 1644.87  822.44  6.7275 0.007042 **
Residuals 17 2078.24  122.25
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The significant result for the variable type.f means that not all of the groups means are the same. (Note: The box plots suggest that the variability within the three groups is not the same. The F-test reported in the ANOVA table is highly sensitive to such heteroscedasticity so the significant result should be accepted with caution.) To see where the differences lie we can display the coefficient estimates with the summary function. (I use the option corr=FALSE to suppress some of the output.)

> summary(lm(energy~type.f),corr=FALSE)

Call:
lm(formula = bats$energy ~ type.f)

Residuals:
      Min       1Q   Median      3Q      Max
-18.69333 -7.45250 -0.04167 5.97417 22.54667

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)      31.050      5.528   5.617 3.08e-05 ***
type.fbirds      -9.897      6.384  -1.550 0.13948
type.fecho bats -27.970      7.818  -3.578 0.00232 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.06 on 17 degrees of freedom
Multiple R-Squared: 0.4418, Adjusted R-squared: 0.3761
F-statistic: 6.728 on 2 and 17 DF, p-value: 0.007042

• The model we've fit is

• The estimates of β₁ and β₂ correspond to the entries labeled type.fbirds and type.fecho bats in the output. From the coding scheme for the dummy regressors we saw displayed in the output of the contrast function, it follows that the intercepts for the three groups are given by the following.

  non-echo bats: β₀
  birds:               β₀ + β₁
  echo bats:        β₀ + β₂

  • Therefore a test H₀: β₁ = 0 is a test of whether birds are different from non-echo bats. From the output we see that this test is not significant. Hence we find no evidence that the energy usages of birds and non-echo bats are significantly different.
  • A test H₀: β₂ = 0 is a test of whether echo bats are different from non-echo bats. From the output we see that this test is significant. Hence we have evidence that the energy usages of echo bats and non-echo bats are different.



Fig. 2 Body mass by animal type

• Since these data do not arise from an experiment in which animals were randomly allocated to groups, we need to be careful in drawing conclusions from any statistical tests we carry out. We first need to convince ourselves that some other factor is not responsible for the differences we see. For instance, if the groups differ systematically with respect to another variable that is known to affect energy usage, then we don't really know whether it's this variable or group membership that is responsible for the observed differences. Such a variable is called a confounding variable (also an effect modifier).
• One such candidate for a confounder is body mass. Fig. 2 compares the distribution of body mass among the three groups of animals.

> boxplot(mass~type.f,ylab='Body mass')
> points(jitter(type),mass, pch=16, col=2, cex=.7)

• In Fig. 2 we see that body mass shows the same relationship with animal type as does energy usage. The animal group with the greatest energy usage (non-echo bats) also has the biggest representatives in the sample. The group with the lowest energy usage (echo bats) also has the smallest individuals. We can test formally for significant differences in body mass between the groups.

> anova(lm(mass~type.f))
Analysis of Variance Table

Response: mass
          Df Sum Sq Mean Sq F value Pr(>F)
type.f     2 434599  217300  7.5001 0.004624 **
Residuals 17 492542   28973
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• So the differences are statistically significant and mass is clearly a confounding variable. It has the potential of obscuring the relationship between energy use in flight and echolocation.

The mass-energy relationship

• To explore the relationship between body mass and energy use further, I create a scatter plot of the two variables. To obtain different symbol types and colors for the three animal types I employ a loop.

> plot(mass, energy, type='n', xlab='mass', ylab='energy')
> for(i in 1:3)
points(mass[type==i], energy[type==i], col=i, pch=14+i)
> legend(500, 10, as.character(unique(type.f)), pch=15:17, col=1:3, cex=rep(.8,3), bty='n')

• Each time through the loop the value of i changes. To take advantage of this fact I make the choice of color, print symbol, and what's plotted depend upon i too. Since the variable type is coded numerically 1, 2, and 3, it can be used to subset the data set differently each time through the loop.

• The plot exhibits a curvilinear trend with a moderate degree of heteroscedasticity (with the variability in energy use increasing with mass). The traditional "solution" for both of these maladies is to transform the response variable using a log or square root transform.
• Energy use is usually related to body mass through a specific type of nonlinear equation called an allometric equation. This is an equation of the form

  

  Here b is called the allometric parameter. If b = 1 then the equation is called isometric. For our data y is energy usage and x is body mass.

  • Note: the rest of the civilized world would probably call this equation a power equation, but the term allometric equation has become ensconced in the biological vernacular particularly when referring to the relative change in size of different portions of the same organism.



Fig. 4  The log(mass) versus log(energy) relationship

• Historically allometric equations were fit by first linearizing them using a log-log transformation. This yields the following equation.

As was discussed in Lecture 23 when comparing the Arrhenius and log-Arrhenius models for the species-area relation, a power equation and its log-transformed equivalent are very different models. They yield different estimates of the parameters, make different assumptions about the error distributions of the residuals, and on the original scale estimate a different type of mean (arithmetic versus geometric mean). Without doing a little more work it's not obvious to me which is the better equation to use with these data. But because Speakman and Racey (1991) used the log-transformed version, so will we. The graph of the log-log transformed data is shown in Fig. 4.

> plot(log(bats$mass), log(bats$energy), type='n', xlab='mass', ylab='energy')
> for(i in 1:3)
points(log(bats$mass[bats$type==i]), log(bats$energy[bats$type==i]), col=i, pch=14+i)
> legend(2, 3.5, as.character(unique(type.f)), pch=15:17, col=1:3, cex=rep(.8,3), bty='n')

• The transformed data are an improvement over the untransformed data. The trend looks more linear and the heteroscedasticity has mostly disappeared. The plot also reveals one of the typical problems of observational and quasi-experimental studies. Since randomization was not employed (thus invoking the gods of statistics), systematic biases occur in the characteristics of our observational units.
• In this study body mass was not randomly assigned to the "treatment" groups—animal types. As a result echolocating bats came small and non-echolocating bats came big. Of course mass is not something that can be randomly assigned to an animal (although perhaps an effort could have been made to find a greater size diversity among the members of the study). In fact non-echolocating bats comprise the megachiroptera a taxonomic group that includes the flying foxes that are among the world's largest bats.
• The bottom line is that body size is undermining the goal of this study because it's hopelessly confounded with the variable of interest (animal category type). But there is a solution to this problem. Even though we didn't control for the bias through randomization and good experimental design we can still try to control for it after the fact by applying appropriate statistical methods.
• Regression is the standard technique for doing this. What we'll do is essentially correct for the bias by subtracting off (statistically) the effect due to mass and only then use the animal's category to account for whatever unexplained variability remains. We do this by fitting a single regression model that includes both body mass and animal category as predictors. Because body mass appears in the model largely to correct for a bias in the sample, we generally refer to it as a covariate rather than a predictor.
• In order for this bias-reduction method to work, we need to be able to make the same correction for all animals. Thus we need to be able to use the same slope in the regression of energy use on body size for all three animal groups.
  • If separate regression lines for the three animal groups are parallel, then we can use the vertical distance between the lines to measure the effect of the animal group category on energy use. Distances between the lines reflect differences in energy use between the groups controlling for size. Because the lines are parallel, it doesn't matter at what body mass we choose to measure this distance. Hence the notion of a "category type" effect on energy use is well-defined.
  • If on the other hand the lines are not parallel, then it would matter what body mass we use for measuring the distance between the lines. In other words, there would be no single "category type" effect. The effect would vary depending on the size of the animal.
• So our first step is to determine if a single slope suffices for all three animal categories. We begin by fitting a separate slope model and see if it can be simplified to a common slope model. This is Model 3 (non-parallel lines) from lecture 29. Let y = energy, z = mass, and x₁ and x₂ be the the two dummy variables that code for animal type, then Model 3 is the following.

• There are two different ways we can fit this model in R.
  • Method 1: Specify all the terms explicitly.

  > model3a<-lm(log(energy)~log(mass) + type.f + type.f:log(mass), data=bats)

What the various terms correspond to in the model equation is indicated below.

• Method 2: Use R's shortcut notation. Because it is standard to fit the main effects, log(mass) and type.f, any time the interaction term, type.f:log(mass), is specified, R offers the shortcut notation type.f*log(mass) for fitting all of these terms simultaneously.

> model3b<-lm(log(energy)~log(mass)*type.f, data=bats)

• Next I use the fitted model to display the regression lines for the three animal groups.

> contrasts(type.f)
             birds echo bats
nonecho bats     0         0
birds            1         0
echo bats        0         1
> coef3<-coef(model3b)
> coef3
    (Intercept)               log(mass)                 type.fbirds
     -0.2024476               0.5897821                  -1.3783902
type.fecho bats   log(mass):type.fbirds   log(mass):type.fecho bats
     -1.2680677               0.2455883                   0.2148750

From the output we see that the fitted line can be written as follows.

Based on the coding shown in the contrasts output, we can construct the lines as follows.

Fig. 5  Model 3—non-parallel lines model

> #nonecho bats x1=x2=0
> abline(c(coef3[1], coef3[2]), lty=2, col=1)
> #birds: x1=1, x2=0
> abline(c(coef3[1]+coef3[3], coef3[2]+coef3[5]), lty=2, col=2)
> #echo bats: x1=0, x2=1
> abline(c(coef3[1]+coef3[4], coef3[2]+coef3[6]), lty=2, col=3)

• There are two troublesome features revealed in the plot.
  • Observe how tightly clustered are the data for three of the echolocating bats and how the fourth observation (Phyllostomas hastatus) is far removed from this cluster. It's pretty clear from the fitted line for the echolocating bats that Phyllostomas hastatus has high leverage and is highly influential. The line has been forced to go through this point! If this point were missing or altered, we'd have a very different picture. Thus we're faced with the unsettling fact that the conclusions of this study will probably hinge on the value of a single observation.
  • The minimal overlap between the three groups in the scatter plot is also worrisome. To compare the distances between separate regression lines (assuming the parallelism hypothesis holds) we will essentially have to assume that parallelism holds over regions for which we have no data. It's a cardinal rule in regression analysis to not extrapolate beyond the range of the data. But to compare the regression lines of the two groups of bats, for example, we would have to do exactly that.

Testing the non-parallel lines model

• Based on the dummy coding scheme used for animal type, Model 3, the non-parallel lines model, is equivalent to the following three separate equations for each animal type.

• Thus to test whether the non-parallel lines model (Model 3) can be reduced to the parallel lines model (Model 2), we need to test . This is a simultaneous test of whether both interaction terms can be dropped from the model. If we fail to reject this hypothesis, then we can assume that all three animal groups have a common slope, β₁.
• The most obvious way of carrying out this test is to fit the common slope model and then compare it to the non-parallel lines model with the anova function of R.

> model2<-lm(log(energy)~log(mass)+type.f,data=bats)
> anova(model2,model3b)
Analysis of Variance Table

Model 1: log(energy) ~ log(mass) + type.f
Model 2: log(energy) ~ log(mass) * type.f
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     16 0.55332
2     14 0.50487  2   0.04845 0.6718 0.5265

• The test carried out by anova here is called a partial F-test. A partial-F test is asymptotically equivalent to the likelihood ratio test we've discussed before. Generally speaking, when a scale parameter is estimated (as it is here with two normal models), an F-test based on an F-distribution is generally preferred over the likelihood ratio test based on a chi-squared distribution. Since the reported p-value (p = 0.5265) is greater than .05, we would fail to reject the null hypothesis of a common slope. Thus it is safe for us to reduce Model 3 to Model 2.
• There is an even more direct test of Model 3. If we apply the anova function to Model 3 alone, we obtain the following.

> anova(model3b)
Analysis of Variance Table

Response: log(energy)
                 Df  Sum Sq Mean Sq  F value Pr(>F)
log(mass)         1 29.3919 29.3919 815.0382 8.265e-14 ***
type.f            2  0.0296  0.0148   0.4100 0.6713
log(mass):type.f  2  0.0484  0.0242   0.6718 0.5265
Residuals        14  0.5049  0.0361
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Observe that the line of the ANOVA table labeled log(mass):type.f reports an F-statistic and p-value identical to what we obtained when comparing Model 3 with Model 2 using a partial-F test. When the anova function is applied to a regression model in which a categorical variable is a predictor, the tests that are printed test the categorical variable as a construct. Each line of the ANOVA table reports a partial-F test in which two models are compared—the current model and a model in which the term identified in that line is omitted. Since it is the interaction term that causes the slopes in Model 3 to be different among the three groups, a test of the interaction term as a construct is a test of Model 3 versus Model 2.

• Note that the partial F-test obtained from ANOVA is not equivalent to the individual t-tests (Wald tests) obtained with the summary function.

> summary(model3b,corr=F)

Call:
lm(formula = log(energy) ~ log(mass) * type.f, data = bats)

Residuals:
      Min        1Q    Median       3Q      Max
-0.251524 -0.126431 -0.009536 0.081244 0.328402

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)
(Intercept)               -0.2024     1.2613  -0.161   0.8748
log(mass)                  0.5898     0.2061   2.861   0.0126 *
type.fbirds               -1.3784     1.2952  -1.064   0.3053
type.fecho bats           -1.2681     1.2854  -0.987   0.3406
log(mass):type.fbirds      0.2456     0.2134   1.151   0.2691
log(mass):type.fecho bats  0.2149     0.2236   0.961   0.3529
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1899 on 14 degrees of freedom
Multiple R-Squared: 0.9832, Adjusted R-squared: 0.9771
F-statistic: 163.4 on 5 and 14 DF, p-value: 6.696e-12

The column labeled t value tests the individual coefficients of the regression model. Thus the last two lines are individual tests of the interaction terms. Each is a test of one dummy variable assuming that all the rest of the variables (including the other dummy variable!) are nonzero. These two tests, individually or taken together, are clearly not equivalent to a test of Model 3 versus Model 2 which requires a test in which both coefficients are set simultaneously equal to zero.

• A graph of the parallel lines model is shown in Fig. 6.

> plot(log(bats$mass), log(bats$energy), type='n', xlab='mass', ylab='energy')
> for(i in 1:3)
points(log(bats$mass[bats$type==i]), log(bats$energy[bats$type==i]),col=i, pch=14+i)
> legend(2,3.5, as.character(unique(type.f)), pch=15:17, col=1:3, cex=rep(.8,3), bty='n')
> coef(model2)->coef2
> coef2
(Intercept)  log(mass) type.fbirds type.fecho bats
-1.57636019 0.81495749  0.10226192      0.07866368
> abline(c(coef2[1],coef2[2]), lty=2, col=1)
> abline(c(coef2[1]+coef2[3], coef2[2]), lty=2,col=2)
> abline(c(coef2[1]+coef2[4], coef2[2]), lty=2,col=3)

Testing the parallel lines model

• Because we can assume the lines are parallel, we can adjust for the differences in masses between the groups. The vertical distances between the parallel lines measure the difference in energy expenditure for members of the three groups adjusting for their differences in body mass. To test whether these distances differ significantly from zero requires yet another partial F test. This time we compare the no interaction model to a model in which there is a common intercept. If the single intercept model turns out to be adequate then we can safely conclude that after correcting for differences in body mass the groups do not significantly differ in energy usage during flight.
• Based on the dummy coding scheme used for animal type, Model 2, the parallel lines model, is equivalent to the following three separate equations for each animal type.

• Thus to test whether the parallel lines model (Model 2) can be reduced to the coincident lines model (Model 1), we need to test . This is a simultaneous test of whether both regressors associated with the animal type categorical variable can be dropped from the model. If we fail to reject this hypothesis, then we can safely assume that all three animal groups have a common intercept, β₀.
• As before we can test this hypothesis by comparing Model 1 to Model 2 with the anova function. Alternatively we can apply the anova function to Model 2 alone.

> #compare Model 1 with Model 2
> model1<-lm(log(energy)~log(mass),data=bats)
> anova(model1,model2)
Analysis of Variance Table

Model 1: log(energy) ~ log(mass)
Model 2: log(energy) ~ log(mass) + type.f
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     18 0.58289
2     16 0.55332  2   0.02957 0.4276 0.6593
> #alternatively apply anova function to Model 2 alone
> anova(model2)
Analysis of Variance Table

Response: log(energy)
          Df  Sum Sq Mean Sq  F value    Pr(>F)
log(mass)  1 29.3919 29.3919 849.9108 2.691e-15 ***
type.f     2  0.0296  0.0148   0.4276 0.6593
Residuals 16  0.5533  0.0346
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• Since the p-value of this test (0.6593) is not significant, we conclude that once we control for mass there are no significant energy use differences in flying between the groups.
• Notice what a different conclusion this is from our initial inspection of the box plots! We've used regression to make the groups comparable, and having done so we found that there are no significant differences between them. So we're done. There's nothing more to do. It's time to send the paper off to Nature to be published.

Testing linear combinations of regression parameters

• Based on our comparison of Model 2 with Model 1, we concluded that there is no difference in energy use between the animal types after having controlled for differences in mass. If we examine the tests of the individual coefficient estimates of model 2, we see that they support this conclusion.

> summary(model2,corr=FALSE)

Call:
lm(formula = log(energy) ~ log(mass) + type.f, data = bats)

Residuals:
     Min       1Q   Median      3Q     Max
-0.23224 -0.12199 -0.03637 0.12574 0.34457

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    -1.57636    0.28724  -5.488 4.96e-05 ***
log(mass)       0.81496    0.04454  18.297 3.76e-12 ***
type.fbirds     0.10226    0.11418   0.896 0.384
type.fecho bats 0.07866    0.20268   0.388 0.703
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.186 on 16 degrees of freedom
Multiple R-Squared: 0.9815, Adjusted R-squared: 0.9781
F-statistic: 283.6 on 3 and 16 DF, p-value: 4.464e-14

• The line labeled type.fbirds is a test of . If we fail to reject this hypothesis we can then conclude that non-echo bats (the baseline group) and birds have the same intercept β₀. Thus the test is a test of whether the energy use of birds is different from that of non-echo bats, controlling for mass. Based on the p-value, 0.384, we conclude that birds and non-echo bats of the same mass do not differ in energy usage.
• The line labeled type.fecho bats is a test of . If we fail to reject this hypothesis we can then conclude that non-echo bats (the baseline group) and echo bats have the same intercept β₀. Thus the test is a test of whether the energy use of echo bats is different from that of non-echo bats, controlling for mass. Based on the p-value, 0.703, we conclude that echo bats and non-echo bats of the same mass do not differ in energy usage.

• What if we are interested in whether birds differ from echo bats? This test does not appear in the output. It is not necessary to refit the model with a different coding scheme to test this hypothesis. It can be tested based on output from this model. From the Model 2 equations, this is a test of

• In lecture 30 we learned that this test can be carried out with a Wald-type test (actually a t-test).

where .

• The coefficient estimates needed for the numerator can be extracted from the model object with the coef function. The variance and covariance terms needed for the denominator are part of the variance-covariance matrix of the parameter estimates. This matrix can be obtained in three different ways.
  • Using the vcov function

  > vcov(model2)
  (Intercept)       log(mass)  type.fbirds   type.fecho         bats
  (Intercept)      0.08250476 -0.012105044 -0.019207033 -0.050561463
  log(mass)       -0.01210504  0.001983939  0.001730953  0.006869742
  type.fbirds     -0.01920703  0.001730953  0.013037676  0.014639321
  type.fecho bats -0.05056146  0.006869742  0.014639321  0.041078883

• From the scaled covariance list element obtained with the ls.diag function.

> ls.diag(model2)$cov.scaled
(Intercept)       log(mass)  type.fbirds   type.fecho         bats
(Intercept)      0.08250476 -0.012105044 -0.019207033 -0.050561463
log(mass)       -0.01210504  0.001983939  0.001730953  0.006869742
type.fbirds     -0.01920703  0.001730953  0.013037676  0.014639321
type.fecho bats -0.05056146  0.006869742  0.014639321  0.041078883

• Scaling the unscaled covariance matrix that is returned by the summary function.

> names(summary(model2))
 [1] "call"          "terms" "residuals" "coefficients"  "aliased"
 [6] "sigma"         "df"    "r.squared" "adj.r.squared" "fstatistic"
[11] "cov.unscaled"
> summary(model2)$cov.unscaled * summary(model2)$sigma^2
(Intercept)       log(mass)  type.fbirds   type.fecho         bats
(Intercept)      0.08250476 -0.012105044 -0.019207033 -0.050561463
log(mass)       -0.01210504  0.001983939  0.001730953  0.006869742
type.fbirds     -0.01920703  0.001730953  0.013037676  0.014639321
type.fecho bats -0.05056146  0.006869742  0.014639321  0.041078883

• Variances lie on the diagonal of this matrix and covariances are off the diagonal. To obtain the terms we need we can just extract them from the matrix.

> varcov.mat<-vcov(model2)
> var.b2minusb3<-varcov.mat[3,3]+varcov.mat[4,4]-2*varcov.mat[3,4]
> var.b2minusb3
[1] 0.02483792

• An easier approach is to construct the appropriate quadratic form. It's easier because we don't actually need to know the formula for the variance. For this all we need are the coefficients of the linear combination.

Using this we obtain the variance we need as follows. (The operator %*% is matrix multiplication.)

> c(0,0,1,-1) %*% varcov.mat %*% c(0,0,1,-1)
[,1]
[1,] 0.02483792

• Thus the t-statistic for our test is calculated as follows.

> coef(model2)
 (Intercept)  log(mass) type.fbirds type.fecho bats
-1.57636019 0.81495749   0.10226192      0.07866368
> t.stat<-(coef(model2)[3]-coef(model2)[4]) / sqrt(c(0,0,1,-1) %*% varcov.mat %*% c(0,0,1,-1) )
> t.stat
[,1]
[1,] 0.1497345

• The t-distribution requires an additional parameter called degrees of freedom that is contained in the model object as $df.residual. The two-tailed p-value for our test statistic is calculated as follows.

> model2$df.residual
[1] 16
> 2*(1-pt(t.stat,model2$df.residual))
[,1]
[1,] 0.8828453

• Clearly the result is not significant. Fig. 7 shows the reference t-distribution, the critical region, and the observed value of the test statistic. Since the observed value of the test statistic does not lie in the critical region, we fail to reject the null hypothesis. We have no evidence to suggest that birds and echo bats differ in energy usage, controlling for mass.

> #code for Fig. 7
> plot(c(-4, 4), c(0, .41), type='n', axes=FALSE, xlab='', ylab='')
> lines(seq(-4, 4, .1), dt(seq(-4, 4, .1), model2$df.residual))
> axis(1)
> abline(h=0, lty=1, col='grey80')
> points(t.stat, 0, pch=4, col=2, cex=2)
> segments(qt(.975, model2$df.residual), dt(qt(.975, model2$df.residual),
model2$df.residual), qt(.975, model2$df.residual), 0)
> segments(-qt(.975, model2$df.residual), dt(qt(.975, model2$df.residual),
model2$df.residual), -qt(.975, model2$df.residual), 0)
> text(3.1, .06, 'critical region', cex=.8)
> text(-3.1, .06, 'critical region', cex=.8)
> text(t.stat, .025, 'observed value', col=2, cex=.8)
> text(1.3, .3, substitute(t[x], list(x=model2$df.residual)), cex=1.5)

Testing user-specified contrasts

• As discussed in lecture 30, it is to your advantage to use coding schemes for dummy regressors that actually test hypotheses you care about. In lecture 30 we replaced the default contrasts selected by R with two other contrasts: a contrast of birds with the average of the bats, and a contrast of the two kinds of bats. We code these into R as follows.

> type.f2<-factor(type)
> contrasts(type.f2)<-cbind(c(1,-2,1),c(1,0,-1))
> contrasts(type.f2)
  [,1] [,2]
1    1    1
2   -2    0
3    1   -1

• Technically for the bat data set these are contrasts but not orthogonal contrasts. That's because the number of observations in the three groups differs. True orthogonal contrasts would require that we weight the contrast coefficients by the sample size. The tests one obtains by doing this are not very interpretable, so I don't try to make the contrasts orthogonal.
• I fit the common slope model, Model 2.

> model2a<-lm(log(energy)~log(mass)+type.f2,data=bats)

• The choice of coding scheme for the dummy variables has no effect on the basic results. If we compare the ANOVA output of this model with that from the original model 2 in which the default coding scheme was used, the results are identical.
  

> anova(model2a)
Analysis of Variance Table

Response: log(energy)
           Df  Sum Sq Mean Sq  F value Pr(>F)
log(mass)   1 29.3919 29.3919 849.9108 2.691e-15 ***
type.f2     2  0.0296  0.0148   0.4276 0.6593
Residuals  16  0.5533  0.0346
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(model2)
Analysis of Variance Table

Response: log(energy)
           Df  Sum Sq Mean Sq  F value    Pr(>F)
log(mass)   1 29.3919 29.3919 849.9108 2.691e-15 ***
type.f      2  0.0296  0.0148   0.4276 0.6593
Residuals  16  0.5533  0.0346
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• But if we use the summary function on the model, the individual tests we get are tests of the contrasts that were specified.

> summary(model2a,corr=FALSE)

Call:
lm(formula = log(energy) ~ log(mass) + type.f2, data = bats)

Residuals:
     Min       1Q   Median      3Q     Max
-0.23224 -0.12199 -0.03637 0.12574 0.34457

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.51605    0.21274  -7.126 2.40e-06 ***
log(mass)    0.81496    0.04454  18.297 3.76e-12 ***
type.f21    -0.02098    0.03103  -0.676 0.509
type.f22    -0.03933    0.10134  -0.388 0.703
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.186 on 16 degrees of freedom
Multiple R-Squared: 0.9815, Adjusted R-squared: 0.9781
F-statistic: 283.6 on 3 and 16 DF, p-value: 4.464e-14

• Observe that neither contrast is statistically significant.

Cited Reference

• Speakman, J. R. and P. A. Racey. 1991. No cost for echolocation for bats in flight. Nature 350: 421–423.

Course Home Page