3

In my RDA triplot I would like to display 'sites', 'species' and their constraints which in my case are Field and Trt. The problem is that not all levels of the constraints are displayed in the plot. There are two levels of each factor.

My RDA code is:

Dummy.rda <- rda(species.rda ~ Field + Trt,RDA.env, scale=TRUE)

summary(Dummy.rda, scaling=3)  #Here I see only one level of each reported in:Biplot scores for constraining variables. However all levels appear in: Centroids for factor constraints

anova.cca(Dummy.rda, step=100, by='margin') # degrees of freedom are correct for both factors (df=1)

plot(Dummy.rda, scaling = 3) #This displays all levels of Field and Trt but only one of each has an arrow

plot(Dummy.rda, display = "species", xlim = xlims, ylim = ylims, 
       scaling = 3)
text(Dummy.rda, scaling = 3, display = "bp")  # I want to customize the RDA plot, but this 'text' only displays 1 level of each of Field and Trt.
Didzis Elferts
  • 95,661
  • 14
  • 264
  • 201
Steve Crittenden
  • 31
  • 1
  • 3
  • The question is a bout the text function ? If it is , Can you show what do you want to plot or a sample of your data to reproduce it? – agstudy Dec 27 '12 at 12:12
  • 1
    try ```display = 'cn' ``` in your last text call. This will plot the factor centroids instead of arrows. – EDi Dec 27 '12 at 15:27
  • The question is why all levels of constraints are not displayed. The code above displays 'sites' and 'species' properly, I would like to use text() to add arrows to represent the environmental constraints. In my case Trt has 2 levels and Field has 2 levels, but with the script above only one level is displayed on the triplot. – Steve Crittenden Dec 28 '12 at 15:18
  • 1
    @SteveCrittenden Sorry to come to this late, but what you are seeing is the standard way factor variables are managed in a regression. one of the dummy variables created from the factor has to be dropped as it is linearly dependent upon the other dummy variables. Hence we have to leave it out of the analysis as an explicit variable - its effect is included. See my Answer for a more detailed discussion/explanation. – Gavin Simpson May 10 '13 at 16:11

1 Answers1

4

The missing levels are because you are trying to view the factor variables as if they were continuous ones - strictly they should not be displayed as biplot arrows I guess. Anyway, just as in regression with dummy variables, one of the levels of the factor cannot be included because it is linearly dependent on the dummy variables for the remaining levels in the model matrix. Consider this example:

require("vegan")
data(dune)
data(dune.env)

mod <- rda(dune ~ Management, data = dune.env)

> model.matrix(mod)
   ManagementHF ManagementNM ManagementSF
2             0            0            0
13            0            0            1
4             0            0            1
16            0            0            1
6             1            0            0
1             0            0            1
8             1            0            0
5             1            0            0
....<truncated>

What you see in the output from model.matrix() are the variables that went into the ordination. Notice that there are three variables in the model matrix but there are four levels in the Management factor:

> with(dune.env, levels(Management))
[1] "BF" "HF" "NM" "SF"

The convention in R is to use the first level of a factor as the reference level. In a regression this would be included in the intercept, but we don't have one of those in RDA. Notice how in the first row of the model.matrix() output, all values are 0; this indicates that that row was in the BF management group. But as only three variables went into the model proper, we can only represent them by three biplot arrows - that is just the way the maths works.

What we can do is plot the group centroids and that is what is shown in the summary() output you refer to and which can be extracted using scores():

> scores(mod, display = "cn")
                   RDA1       RDA2
ManagementBF -1.2321245  1.9945903
ManagementHF -1.1847246  0.7128810
ManagementNM  2.1149031  0.4258579
ManagementSF -0.5115703 -2.0172205
attr(,"const")
[1] 6.322924

Hence to add the centroids to the existing plot do:

text(mod, scaling = 3, display = "cn")

No matter what you do, you can't add a biplot arrow for the reference group.

I hope that explains the behaviour you are seeing?

Gavin Simpson
  • 170,508
  • 25
  • 396
  • 453