Get frequencies (absolute and relative) of levels of a categorical variable from incidence binary data by combination of columns factors

Question

I would like to have the frequencies of each levels of a categorical variable (row vector) denoting ecological type (3 levels: H,F,T) of a set of 93 herbaceous plants for the observed species present (=1) conditioning by sites (3 levels: A,B,C), habitats (3 levels: 1,2,3,4) and years (3 levels: 1,2,3).

I know the procedure is passed by tapply(), but the messy thing come from the logic operator for linking levels of the categorical variable (H,F,T) for the present species (=1) accross all of the species conditioning by combination of columns factors.

This could be summarized by a 12 x 3 contingency table indicating the numbers of each ecological types (3) of species per sites (3) and habitats (4).

Ex of my data (each habitat contain 20 lines): for each species (Sp1 to Sp93) 0 for absent and 1 for present. Vector "type" contain ecological type for each species.

Site,Habitat,Year,Sp1,Sp2,Sp3,Sp4,Sp5,Sp6,...,Sp93

type= c(H,H,F,T,F,T,H,....T) # vector of length 93

Thank you in advance.

---

I hope this would help describe my data objects better.

data = read.csv(file = "Veg_06.csv", header = TRUE)

data = data[1:240, -c(1,4:7)]

Ilot # Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ... each level has 4 sublevels (from "Site") with 20 lines each, adding up to 80 lines by levels.

Site # Factor w/ 4 levels "Am","Av","CP","CS": 2 2 2 2 2 2 2 2 2 2 ...

Sp # int [1:240] 0 0 0 0 0 0 0 0 0 0 ... either "0" or "1" for absence or presence of species.

veg # Factor w/ 3 levels "H","F","T": 3 3 2 2 3 1 2 1 2 1 ... categorical factor indicating type of species.

Answer 1

First off, I would recommend http://vita.had.co.nz/papers/tidy-data.pdf, Hadley Wickham's paper on Tidy Data, for some ideas on how to organize the data to be better suited to analysis. In essence, we think of each row as a single observation.

It sounds like fundamentally, your data is a collection of year, site, habitat, quadrant(? maybe line, not sure from the description), species with the observation point being that species was observed in that site, habitat, quadrant, and year. For simplicity, a row is present if the species is present.

In addition, there's the concept of type, which is associated with each species.

Analyzing and contingency table

Putting aside the question of how to get your data into this form, let's assume that we have the data in the form described above.

> raw <- expand.grid(species=1:93, quadrant=1:20, habitat=1:4, site=1:3, year=1:3)
> head(raw)
  species quadrant habitat site year
1       1        1       1    1    1
2       2        1       1    1    1
3       3        1       1    1    1
4       4        1       1    1    1
5       5        1       1    1    1
6       6        1       1    1    1

And let's take a small sample and a large sample

> set.seed(100); d.small <- raw[sample(nrow(raw),20), ]
> set.seed(100); d.large <- raw[sample(nrow(raw),1000), ]

We can use the ftable function to get this into a state that we want, the 12x4 contingency table, as

> ftable(habitat ~ year + site, data=d.small)
          habitat 1 2 3 4
year site
1    1            0 0 1 0
     2            0 0 1 1
     3            0 1 1 1
2    1            2 1 1 0
     2            1 1 0 2
     3            0 0 1 0
3    1            2 0 0 1
     2            0 1 0 1
     3            0 0 0 0

This will count the same species twice if it occurs in two different quadrants of the site/habitat mixture. We can discard the habitat and unique-ify to get the count across all of them

> ftable(habitat ~ year + site , data=unique(d.small[c('species', 'habitat','year','site')]))

Transforming (tidying the source data)

To transform the data as it stands into a form like this is tricky in vanilla R. With the tidyr package it gets easier (reshape does very similar things as well)

> onerow <- data.frame(year=1, site=1, habitat=2, quadrant=3, sp1=0, sp2=1,sp3=0,sp4=0,sp5=1)
> onerow
  year site habitat quadrant sp1 sp2 sp3 sp4 sp5
1    1    1       2        3   0   1   0   0   1

Here I'm making assumptions about what your data look like that seem reasonable

> subset(gather(onerow, species, present, -(year:quadrant)), present==1)
  year site habitat quadrant species present
2    1    1       2        3     sp2       1
5    1    1       2        3     sp5       1
> subset(gather(onerow, species, present, -(year:quadrant)), present==1, select=-present)
  year site habitat quadrant species
2    1    1       2        3     sp2
5    1    1       2        3     sp5

And now you can proceed with the analysis above.

Merging in the species type data

Looking at your description a little closer, I think you also want to merge in a parallel vector of species type information.

> set.seed(100); sp.type <- data.frame(species=1:93, type=factor(sample(1:4, 93, replace=T)))
> merge(d.small, sp.type)
   species quadrant habitat site year type
1        6       16       4    2    3    2
2       27        9       2    2    2    4
3       27        8       4    2    1    4
4       32       18       1    2    2    4
5       33       18       1    1    2    2
6       45       14       4    2    2    3
7       49        6       2    3    1    1
8       54        3       3    2    1    2
9       55        2       1    1    3    3
10      56        2       4    3    1    2
11      56        1       3    1    1    2
12      57        7       2    1    2    1
13      62       18       4    2    2    3
14      70       19       1    1    2    3
15      77        2       3    3    1    4
16      80        7       3    1    2    1
17      81       17       1    1    3    2
18      82        5       2    2    3    3
19      86        9       4    1    3    3
20      87       10       3    3    2    3

And now you can use the subset, unique, and ftable approach above to get the data you need.

Answer 2

Assuming you had a dataframe with (among other things) the columns named: "sites", "habitats", "years":

dfrm <- data.frame( sites = sample( LETTERS[1:3], 20, replace=TRUE),
            habitats= sample( factor(1:4), 20, replace=TRUE),
            years = sample( factor(paste("Y",1:4, sep="_")), 20, replace=TRUE) )

Then this will give you an additional factor-mode column that encodes the various levels of each row.

dfrm$three.way.inter <- with(dfrm, interaction(sites, habitats, years))

If you want non-populated levels then do nothing else. If you want possible levels that have no instances, then use drop=TRUE. Then you can analyze these within individual levels of the three classification variables.