Linear model without intercept

We want to know what is the expression of gene 1159 under each condition and compare them. As discussed in a previous article, we can find the gene expression for each condition with a linear model without intercept.

model0 <- lm(ex[1159,] ~ condition + 0)
summary(model0)

The fitted model’s coefficients correspond to the mean expression level for each condition. Each 𝑝-value corresponds to the null hypothesis “real gene expression is zero”You can find several web pages with refreshers about 𝑝-values and statistical significance on the web.

. As expected, all 𝑝-values are small, since the gene is expressed in all conditions.

The stars in the last columns are used to show which 𝑝-values are below the traditional thresholds. Three stars means 𝑝-value below 1%, two stars represent “under 5%”, and so on.

One way to understand the model is the following diagram. Each arrow represents one comparison made when we fit the linear model.

In this basic case, each circle represent the average expression for each condition, and the arrows represent the estimation calculated by fitting the linear model. It is important to know that we can have only four arrows in this diagram, since we have four conditions, and that limits how many combinations we can take.

The expression levels are hard to interpret without context. Are they high or low? The mean values are different. Are they significantly different? In other words, are the underlying expression levels really different, or the difference of means corresponds only to the noise?

Linear model with intercept

To answer these questions it is more useful to fit a linear model with intercept, that will show the difference between a reference condition against all other conditions.

model1 <- lm(ex[1159,] ~ condition)
summary(model1)

The intercept shows the expression level of the reference condition. Again, we can represent this model with a diagram.

In this case the reference is automatically chosen to be “the first level of the condition factor”, which in this case is Blue. This is because, unless we say otherwise, the factor levels are sorted alphabetically.

levels(condition)

[1] "Blue"     "Darkness" "Red"      "White"   

This does not necessarily corresponds to the most interesting question⊕This is the key question.

. We need a way to choose the reference level.

Choosing the reference level

There are several ways to choose the reference level. We will discuss some of the easy methods. The easiest one is to use the function relevel() to create a new factor with the same levels in different order.

From the the sake of argument I will assume that a good reference level is "Darkness". Thus, the new factor will be called dark_vs_. We create and examine it as follows.

dark_vs_ <- relevel(condition, "Darkness")
levels(dark_vs_)

[1] "Darkness" "Blue"     "Red"      "White"   

We see that this new factor dark_vs_ contains the same values as the old factor condition, but the levels are in a different order. Now we can fit the linear model again.

model2<- lm( ex[1159,] ~ dark_vs_ )

summary(model2)

In the new results the intercept shows the expression level for "Darkness", and the other coefficients show the differential expression for "Blue", "Red" and "White" light.

Why doing this complicated way? We could just calculate the difference between the averages? Because of the 𝑝-values. This time the null hypothesis is “the expression levels are the same”. Small 𝑝-values show that we do not have enough evidence to reject the null hypothesis. In other words, there is no statistically significant difference between the expression levels. The gene 1159 does not respond to the light.

Using Contrast matrices

There is another approach to compare gene expressions, more flexible than just reorganizing factor levels. As usual, higher flexibility implies higher complexity.

To make the example more clear, we use the original condition factor, without reordering.

The function contr.treatment() can be used to create a contrasts matrix, like this.

contrast1 <- contr.treatment(levels(condition), base = 2)

contrast1

This function takes the names of factor’s levels and the number indicating the base level used for reference. The resulting contrasts matrix has one row and one column for each level of condition, except that there is no column for the base level. The row that only has zeros will be the base level.

Now when we fit the linear model we need to include a contrasts= option. This option takes a list, where each element has the name of a factor. There can be several factors in a model, each one with its own contrast. In our case we have only one factor.

model3b <- lm(ex[3342,] ~ condition, contrasts = list(condition=contrast1))

summary(model3b)

We see that we recovered the same result as model3 without re-leveling the factor. It is more complicated, but will have some advantages.

Making the contrast matrix manually

Once we know how the contrast matrix looks, we can make it manually. In R we use matrix(data, nrow, byrow) to create a new matrix from a vector data. If we do not include byrow=TRUE then the matrix will be filled by columns. Try it and see the difference.

The code to make a matrix like contrast1 is this

contrast2 <- matrix(
c(1,0,0,
  0,0,0,
  0,1,0,
  0,0,1), nrow=4, byrow=TRUE)
colnames(contrast2) <- c("B_D","R_D","W_D")
rownames(contrast2) <- levels(condition)

contrast2

Obviously, it produces the same model and the same coefficients as before.

Other combinations

We can make other contrast matrices, but there are technical constraints on what can be done. One way to easy avoid these contrast is to build a single comparison. For example, an important biological issue is to compare the gene expression under blue light versus under red light, and find if the difference is statistically significant.

To answer this question we build a one-column matrix with a code like this

contrast3 <- matrix(c( 1/2, 0,-1/2, 0), nrow=4)
colnames(contrast3) <- c("B-R")
rownames(contrast3) <- levels(condition)

which produces the following matrix

contrast3

Not the fitted model has only two coefficients

model5 <- lm(ex[3342,] ~ condition,
 contrasts=list(condition=contrast3))
summary(model5)

Now it will be useful to keep in mind the mean expression for each condition.

tapply(ex[3342,], condition, mean)

    Blue Darkness      Red    White 
 10.8700  13.7450  10.8675  14.5225 

Close inspection shows that coefficient conditionB-R is the difference between mean expression in "Blue" and "Red" conditions. More importantly, we get a 𝑝-value for the difference. We conclude that, although "Blue" and "Red" are different from "Darkness".

We also observe that the intercept is the mean of all expressions.

mean(ex[3342,])

[1] 12.50125

Unfortunately it is not easy to do this comparison for several conditions at the same time. For this article we will limit us to only one comparison per linear model.

All versus the average

We have the mean expression for each condition and the global mean of all expressions. It is easy to calculate the difference between the global mean and the mean of each condition.

tapply(ex[3342,], condition, mean) - mean(ex[3342,])

    Blue Darkness      Red    White 
-1.63125  1.24375 -1.63375  2.02125 

It is not so easy to see if this difference is significant. We can answer this question using contrast matrices created with the contr.sum() function.

contrast4 <- contr.sum(levels(condition))
colnames(contrast4) <- c("B","D","R")

The contr.sum() does not assign names to the columns, so we assign them manually. The resulting matrix is this

contrast4

We see that the last row is different from the rest, and indeed it is not included in the output. This is always the case when we do not have enough degrees of freedom. Since we are comparing to the average, we lose one degree of freedom. We can represent this model with the following diagram.

As before, we can have up to four arrows —that is, contrasts— given that we are limited by the number of conditions. To include the mean we must “sacrifice” one of the conditions to make space for the average. The “lost” row is replaced with a row of -1.

The fitted model has the following coefficients

model6 <- lm(ex[3342,] ~ condition,
 contrasts=list(condition=contrast4))
summary(model6)

We observe that the intercept is the global average, and the rest of the coefficients are the difference between the average and the average expression under each condition. Also, we have 𝑝-values to evaluate the statistical significance of each difference.

How would you make a contrast matrix to compare the global average to “Blue”, “Red” and “White”, omitting “Darkness”?