Factor levels in R: reference groups and empty sites

dplyr

data wrangling

data cleaning

regression

ecology tutorial

Factor levels in R decide which group a model compares against and which sites vanish from a count. Set them from the design, keep empty sites, catch stray NAs.

Author

Tidy Ecology

Published

2026-09-27

A fire experiment: quadrats that were left alone, quadrats that were burnt, quadrats that were mown, and the dry biomass cut from each at the end of the season. You fit lm(biomass ~ treatment), open the coefficient table, and there is a row for treatmentcontrol but no row for burnt. Nothing failed. R has made a choice about your treatment column, and the table is an answer to a question you did not ask.

A factor is R’s way of storing a grouping column: a set of labels plus a list of the allowed values, the levels, in a fixed order. That list decides three things that show up in ecological results. Its first entry is the group every model coefficient is measured against. Its members are the groups that a count or a summary is allowed to report, including groups with nothing in them. And any value that is not on the list stops being a value at all. Each of the three main sections below runs one of these on a small simulated data set, shows the number that comes out wrong without any warning, and then fixes it. Two further sections add the cases next to them: an empty cell in a crossed design, where R does warn but quietly, and a grouping column stored as numbers, which never becomes a factor at all.

The short answer. Turn every grouping column into a factor yourself, with the levels written out from your study design rather than taken from the data. Put the group you want to compare against first. Include the sites that produced no records. Then check that no value dropped out on the way.

This post assumes you have met group_by() and summarise(), which Summarising ecological data by group introduces. The data are simulated, and every seed is in the code.

library(dplyr)
library(ggplot2)
options(scipen = 6)   # print small p-values as decimals

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          axis.text        = element_text(colour = te_body),
          strip.text       = element_text(colour = te_ink, face = "bold"))
}

The first level is the reference group

Eight quadrats per treatment, with true mean biomass fixed in the code at 120 g for control, 90 g for burnt and 95 g for mown: mowing is built to remove almost as much as fire does. The treatment column is plain text, which is what read.csv() and data.frame() give you: since R 4.0.0 neither turns text into factors unless you ask.

set.seed(1012)
treat_means <- c(control = 120, burnt = 90, mown = 95)    # true means, g per quadrat
n_per <- 8

plots <- data.frame(treatment = rep(names(treat_means), each = n_per))
plots$biomass <- round(rnorm(nrow(plots), mean = treat_means[plots$treatment], sd = 15), 1)

stopifnot(is.character(plots$treatment),
          is.character(read.csv(text = "treatment\nburnt")$treatment))
head(plots, 3)

  treatment biomass
1   control   107.9
2   control   117.5
3   control   111.5

lm() needs groups, so it converts the text column to a factor on the way in. factor() sorts the distinct values to make the levels, the first sorted level becomes the reference, and “burnt” comes before “control” in the alphabet.

stopifnot(identical(levels(factor(plots$treatment)), c("burnt", "control", "mown")))

fit_default <- lm(biomass ~ treatment, data = plots)
cf_def <- coef(summary(fit_default))
round(cf_def, 4)

                 Estimate Std. Error t value Pr(>|t|)
(Intercept)       88.9500     5.1737 17.1927   0.0000
treatmentcontrol  37.8000     7.3167  5.1662   0.0000
treatmentmown     -4.4875     7.3167 -0.6133   0.5463

Read the table the way the model means it. The intercept, 89.0 g, is the mean of the burnt quadrats. treatmentcontrol is control minus burnt, +37.8 g. And treatmentmown is mown minus burnt: -4.5 g with a p-value of 0.55. Read as “the effect of mowing”, that row says mowing did very little. It says nothing of the sort. It compares two disturbance treatments with each other, and the comparison the experiment was built for, mown against untouched, is not in the table at all.

The fix is one line. Make the column a factor and name the reference, either with relevel() or by writing out all the levels in the order you want (useful later for plots and tables):

plots$treatment <- relevel(factor(plots$treatment), ref = "control")
# same result: factor(plots$treatment, levels = c("control", "burnt", "mown"))
levels(plots$treatment)

[1] "control" "burnt"   "mown"

fit_control <- lm(biomass ~ treatment, data = plots)
round(coef(summary(fit_control)), 4)

               Estimate Std. Error t value Pr(>|t|)
(Intercept)    126.7500     5.1737 24.4988        0
treatmentburnt -37.8000     7.3167 -5.1662        0
treatmentmown  -42.2875     7.3167 -5.7796        0

Now the intercept is the control mean, 126.8 g, and each row is a treatment against it: burnt -37.8 g, mown -42.3 g with a p-value below 0.0001. The same word, mown, sat on a row with a small estimate and a large p-value a moment ago, and now sits on a row with a large negative estimate. relevel() only works on a factor; on a text column it stops with an error, which is why the factor() call comes first.

Nothing about the model itself has changed. The coefficients are differences between the same three group means, taken from a different starting group. (How a factor becomes an intercept plus differences is in t-tests and ANOVA as linear models in R.)

all.equal(fitted(fit_default), fitted(fit_control))

[1] TRUE

max_fit_diff <- max(abs(fitted(fit_default) - fitted(fit_control)))
round(c(logLik_default = as.numeric(logLik(fit_default)),
        logLik_control = as.numeric(logLik(fit_control)),
        r2_default = summary(fit_default)$r.squared,
        r2_control = summary(fit_control)$r.squared), 4)

logLik_default logLik_control     r2_default     r2_control 
      -96.8517       -96.8517         0.6575         0.6575

grp_mean <- tapply(plots$biomass, plots$treatment, mean)
round(grp_mean, 2)

control   burnt    mown 
 126.75   88.95   84.46

The fitted values of the two models agree to within floating-point rounding. The log-likelihood is -96.8517 both times and R-squared is 0.6575 both times. The reference level changes what each coefficient means; it does not change the fit, the predictions or any test of the treatment as a whole (the next section shows one exception for the predictions: an empty cell in a crossed design). That is why the mistake is easy to miss: every diagnostic you run on the model looks the same. Contrasts and post-hoc comparisons in R goes on to other codings (sum-to-zero, planned contrasts) and to comparing every pair.

Two panels of quadrat biomass. Each shows the eight quadrats per treatment as pale green points, a dark green diamond at each treatment mean, a dashed horizontal line at the reference mean and vertical rust arrows labelled with the coefficient. Left panel, reference burnt near 89 g: arrows of plus 37.8 to control and minus 4.5 to mown. Right panel, reference control near 127 g: arrows of minus 37.8 to burnt and minus 42.3 to mown. — Figure 1: The same biomass data, the same three group means and the same fitted model, read against two reference levels. The dashed line is the intercept and each arrow is one coefficient, drawn from the reference mean to the group mean. With burnt as the reference, the mown arrow is short; with control as the reference it is long and points down.

The same levels set the order of groups in plots and tables: alphabetical unless you say otherwise. To sort groups by a statistic instead, for bars in order of their mean, reorder() builds the level order from the data:

levels(reorder(plots$treatment, plots$biomass, FUN = mean))

[1] "mown"    "burnt"   "control"

Unless a cell of the design is empty

The reference level did not change any prediction above because every treatment had quadrats in it, so every prediction was the mean of quadrats that were cut. With two grouping columns crossed, that can fail. A second survey cuts biomass in three habitats in three seasons, four quadrats for each habitat and season, except that the wetland was under water in spring and its quadrats were never cut. The true means are fixed in the code, with a wetland that is poor in spring and peaks in autumn.

set.seed(1015)
cell_mean <- rbind(forest    = c(spring = 60, summer =  80, autumn =  70),   # true means, g
                   grassland = c(spring = 90, summer = 140, autumn = 100),
                   wetland   = c(spring = 70, summer = 130, autumn = 150))
survey <- expand.grid(quadrat = 1:4, season = colnames(cell_mean),
                      habitat = rownames(cell_mean), stringsAsFactors = FALSE)
survey <- subset(survey, !(habitat == "wetland" & season == "spring"))   # flooded, never cut
survey$biomass <- round(rnorm(nrow(survey), cell_mean[cbind(survey$habitat, survey$season)],
                              sd = 15), 1)
table(survey$habitat, survey$season)

           
            autumn spring summer
  forest         4      4      4
  grassland      4      4      4
  wetland        4      0      4

The model biomass ~ habitat * season asks for a separate mean in every habitat and season, the interaction. Fit it twice: once with the text columns as they are, so forest and autumn are the references because they come first in the alphabet, and once with grassland and summer set as the references by hand. Then ask both fits for the cell that was never cut.

fit_alpha <- lm(biomass ~ habitat * season, data = survey)

survey_hand <- survey
survey_hand$habitat <- factor(survey$habitat, levels = c("grassland", "forest", "wetland"))
survey_hand$season  <- factor(survey$season, levels = c("summer", "spring", "autumn"))
fit_hand <- lm(biomass ~ habitat * season, data = survey_hand)

wet_spring <- data.frame(habitat = "wetland", season = "spring")
pred_alpha <- suppressWarnings(predict(fit_alpha, wet_spring))   # the warning is shown below
pred_hand  <- suppressWarnings(predict(fit_hand, wet_spring))
round(c(forest_autumn_first = as.numeric(pred_alpha),
        grassland_summer_first = as.numeric(pred_hand)), 2)

   forest_autumn_first grassland_summer_first 
                132.85                  64.15

Both fits give the same fitted value for every quadrat that was cut, and they still predict 132.85 g and 64.15 g for wetland in spring, a cell built at 70 g. Neither number comes from a wetland quadrat in spring, because there are none. Each is put together from three cells that were cut: the wetland mean in the reference season, plus the reference habitat in spring, minus the reference habitat in the reference season. For the alphabetical fit that is wetland in autumn plus forest in spring minus forest in autumn:

cell_avg <- function(h, s) mean(survey$biomass[survey$habitat == h & survey$season == s])
three_cells <- cell_avg("wetland", "autumn") + cell_avg("forest", "spring") - cell_avg("forest", "autumn")
round(c(three_cells = three_cells, alphabetical_fit = as.numeric(pred_alpha)), 2)

     three_cells alphabetical_fit 
          132.85           132.85

It is the prediction of the alphabetical fit. For the fit set by hand the three cells are wetland in summer plus grassland in spring minus grassland in summer. This is arithmetic, not something the simulation found, and it means the gap between the two predictions is set by how differently the habitats change between seasons, the very thing the interaction was meant to estimate.

R does flag the problem, in three places that are easy to scroll past:

coef(fit_alpha)[is.na(coef(fit_alpha))]

habitatwetland:seasonspring 
                         NA

grep("singularities", capture.output(summary(fit_alpha)), value = TRUE)

[1] "Coefficients: (1 not defined because of singularities)"

tryCatch(predict(fit_alpha, wet_spring), warning = conditionMessage)

[1] "prediction from rank-deficient fit; attr(*, \"non-estim\") has doubtful cases"

An NA coefficient means R dropped a column it could not estimate. Fitting a senescence model to individual records meets one of the same kind, where age is tied to female and season, and the slope R does report is what one arbitrary constraint happens to produce. Since R 4.3.0, predict() on an lm() fit takes rankdeficient = "NA" and returns NA for such a row instead of a number, which is the honest answer here: wetland in spring was not measured.

If a number is needed anyway, fit biomass ~ habitat + season without the interaction. It predicts 104.9 g whatever the reference levels, against the 70 g the cell was built at, because it assumes the wetland differs from the other habitats by the same amount in every season. That is an assumption to state in the methods, and in this simulation it is false by construction. Multilevel post-stratification of records fills the empty cells of a survey frame with a main-effects model and with a multilevel model, and measures what each costs.

Sites with no records drop out of a count

A season of camera trapping for pine marten at twenty sites, stored the way detection data usually is: one row per detection. A site where the camera ran all season and saw nothing has no row. Each site’s expected number of detections is drawn at random, so some sites end up empty.

set.seed(1013)
site_ids  <- sprintf("S%02d", 1:20)                 # the survey design
site_rate <- exp(rnorm(length(site_ids), log(1.5), 1))
n_det     <- rpois(length(site_ids), site_rate)

records <- data.frame(site = rep(site_ids, times = n_det))
c(detections = nrow(records), sites_surveyed = length(site_ids),
  sites_with_none = sum(n_det == 0))

     detections  sites_surveyed sites_with_none 
             50              20               8

Count detections per site and summarise:

per_site_text <- count(records, site)
c(sites_in_table = nrow(per_site_text),
  mean_detections = mean(per_site_text$n),
  share_detected  = mean(per_site_text$n > 0))

 sites_in_table mean_detections  share_detected 
      12.000000        4.166667        1.000000

The table has 12 sites, not 20. The 8 sites with no detections are not in records, so no function working from records alone can know they exist. The mean comes out at 4.17 detections per site, and the share of sites where the marten was detected comes out at 1.00: every site in the table has at least one detection, because that is how a site gets into the table.

Turning the column into a factor does not help if the levels come from the same records, since factor() can only list the values it sees: nlevels(factor(records$site)) is 12. The levels have to come from the design, the list of sites you actually surveyed. Even then, count() and group_by() drop empty levels by default, so the argument .drop = FALSE is the second half of the fix:

records$site <- factor(records$site, levels = site_ids)

per_site_default <- count(records, site)
per_site_all     <- count(records, site, .drop = FALSE)

stopifnot(nrow(per_site_default) == sum(n_det > 0),
          nrow(per_site_all) == length(site_ids),
          nrow(count(data.frame(site = as.character(records$site)), site,
                     .drop = FALSE)) == sum(n_det > 0))
c(default = nrow(per_site_default), drop_false = nrow(per_site_all))

   default drop_false 
        12         20

The last line of the stopifnot() is the other half of the rule: .drop = FALSE does nothing for a text column, because a text column has no list of levels to keep. You need both, a factor with the design’s levels and .drop = FALSE. Base R’s table() keeps empty levels without being asked, and group_by() takes the same argument:

sum(table(records$site) == 0)

[1] 8

per_site <- records |>
  group_by(site, .drop = FALSE) |>
  summarise(detections = n(), .groups = "drop")

rbind(zeros_dropped = c(sites = nrow(per_site_text), mean = mean(per_site_text$n),
                        detected = mean(per_site_text$n > 0)),
      zeros_kept    = c(sites = nrow(per_site), mean = mean(per_site$detections),
                        detected = mean(per_site$detections > 0)))

              sites     mean detected
zeros_dropped    12 4.166667      1.0
zeros_kept       20 2.500000      0.6

With the empty sites back, the mean is 2.50 detections per site and the share of sites with a detection is 0.60. Dropping the zeros inflated the mean by 67 per cent and turned a species detected at 12 of 20 sites into one detected at every site in the table. No simulation is needed for the first number: the same 50 detections are divided by 12 sites instead of 20, so the mean is always inflated by the number of surveyed sites over the number with a detection, minus one (here 20/12 - 1).

Bar chart of detections for sites S01 to S20. Twelve sites have green bars, the tallest at 14 detections for S05. Eight sites, S02, S07, S09, S12 and S17 to S20, have no bar and a red cross on the zero line. A dashed red horizontal line sits at 4.17 and a solid dark green line at 2.50. — Figure 2: Detections per site for the twenty simulated camera sites. Red crosses mark the sites with no detections, which have no rows in the detection table. The dashed line is the mean over the sites that count() returns by default; the solid line is the mean over all twenty sites.

The same zeros can come back through a join instead, when the site list lives in its own table; Joining ecological tables without losing zeros takes that route. Either way the site list has to exist somewhere, and it can be lost on the way through a file: What a CSV loses shows an unsampled site disappearing from the level set when a factor is written to CSV and read back.

Levels set by hand turn mismatches into NA

Writing the levels out is the fix in both sections above, and it has a trap of its own. Here is a habitat column for 150 vegetation plots, typed by several people over a season. Most entries are clean, and a few carry a capital letter, a trailing space or a two-word spelling.

set.seed(1014)
n_rec    <- 150
true_hab <- sample(c("forest", "grassland", "wetland"), n_rec, replace = TRUE,
                   prob = c(0.5, 0.3, 0.2))
typed <- true_hab
slip  <- runif(n_rec) < 0.15                         # entries typed differently
variants <- list(forest = c("Forest", "forest "), grassland = c("Grassland", "grassland "),
                 wetland = c("wet land", "Wetland"))
typed[slip] <- vapply(true_hab[slip], function(h) sample(variants[[h]], 1), character(1))
table(typed)

typed
    forest     Forest    forest   grassland  Grassland grassland    wet land 
        62          4          7         38          3          6          1 
   wetland    Wetland 
        28          1

You want the habitats in a fixed order, so you write the levels out:

hab_levels <- c("forest", "grassland", "wetland")
habitat <- factor(typed, levels = hab_levels)
c(NA_in_typed = sum(is.na(typed)), NA_in_factor = sum(is.na(habitat)))

 NA_in_typed NA_in_factor 
           0           22

factor() matches each value against the levels exactly, character for character, and every value without a match becomes NA. It does this without a warning. 22 of the 150 plots, 15 per cent of the sheet, now have no habitat, although the column had no missing values a line earlier. The R documentation for factor() states this behaviour; it is not a bug, and nothing downstream will point back to it.

What it does to a class share depends on where the mismatches fall:

round(cbind(true = share_true, after_levels = share_hand,
            plots_lost = lost, plots_left = kept), 3)

           true after_levels plots_lost plots_left
forest    0.487        0.484         11         62
grassland 0.313        0.297          9         38
wetland   0.200        0.219          2         28

table() leaves NA out unless told otherwise, so every share is now a share of the 128 plots that survived. Wetland lost 2 of its 30 plots, fewer than its fair part of the 22 losses, so its share rose from 0.200 to 0.219 although no wetland plot was added. Grassland lost 9 of its 47, more than its fair part, and its share fell from 0.313 to 0.297. In this simulation every habitat is mistyped at the same rate, so the shares move only a little. When one habitat has a spelling of its own (one recorder who always writes “wet land”) the losses pile up in that class and its share drops with them. The counts are wrong regardless: here every habitat has fewer plots than were surveyed.

The fix is a check before and a check after. Before, ask which values have no level to go to; clean them; ask again:

setdiff(unique(typed), hab_levels)

[1] "grassland " "Forest"     "Grassland"  "forest "    "Wetland"   
[6] "wet land"

cleaned <- tolower(trimws(typed))
setdiff(unique(cleaned), hab_levels)

[1] "wet land"

cleaned[cleaned == "wet land"] <- "wetland"
habitat <- factor(cleaned, levels = hab_levels)
stopifnot(sum(is.na(habitat)) == sum(is.na(typed)))
all(table(habitat) == table(factor(true_hab, levels = hab_levels)))

[1] TRUE

trimws() and tolower() handle the spaces and capitals, and the two-word spelling needs a decision, written down in the code. Cleaning species names before you count does the same job for species names with a lookup table, which scales better than one line per variant. The lasting fix is on the sheet itself: Broman and Woo (2018) make one consistent code per category part of their first rule for spreadsheet data. The stopifnot() after the factor() call is the one to keep in every script: the number of missing values must not change when a column becomes a factor. The forcats package offers fct(), which stops with an error instead of creating NA (chapter 16 of R for Data Science shows the two side by side), if you prefer a loud failure.

One more factor trap is common enough to name: as.numeric() on a factor returns the level positions, not the values. Debugging and defensive R code measures what that does to a regression on elevation.

A grouping column stored as numbers

Field numbers, block numbers and site codes are often typed as 1, 2, 3. R reads such a column as numbers, and a model given numbers fits one straight line through them instead of a separate mean for each group. Say the fire experiment’s 24 quadrats sit in four fields, two quadrats of each treatment in each field, and the fields differ in productivity in an order that has nothing to do with their numbers.

set.seed(1016)
field_effect <- c(0, 40, -30, 20)          # g per quadrat, added in fields 1 to 4
fields <- data.frame(field = rep(1:4, each = 6),
                     treatment = rep(rep(names(treat_means), each = 2), times = 4))
fields$biomass <- round(treat_means[fields$treatment] + field_effect[fields$field] +
                          rnorm(nrow(fields), sd = 15), 1)
str(fields)

'data.frame':   24 obs. of  3 variables:
 $ field    : int  1 1 1 1 1 1 2 2 2 2 ...
 $ treatment: chr  "control" "control" "burnt" "burnt" ...
 $ biomass  : num  131.7 140.6 109.5 92.2 67.6 ...

str() shows the problem before any model does: field is int. Fit the treatment test with the field column as it is, and again as a factor:

anova(lm(biomass ~ field + treatment, data = fields))

Analysis of Variance Table

Response: biomass
          Df  Sum Sq Mean Sq F value Pr(>F)  
field      1     6.4    6.44  0.0079 0.9299  
treatment  2  4672.3 2336.15  2.8749 0.0799 .
Residuals 20 16252.1  812.61                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(lm(biomass ~ factor(field) + treatment, data = fields))

Analysis of Variance Table

Response: biomass
              Df  Sum Sq Mean Sq F value    Pr(>F)    
factor(field)  3 10664.8  3554.9 11.4392 0.0001977 ***
treatment      2  4672.3  2336.2  7.5174 0.0042336 ** 
Residuals     18  5593.8   310.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The Df column gives it away: field takes 1 degree of freedom, one slope, and factor(field) takes 3, one difference for each field after the first. A straight line through field effects of 0, 40, -30 and 20 g removes almost none of them (0.2 per cent of their spread), so they stay in the residuals and the treatment test weighs its effect against that extra noise. In this data set the treatment p-value is 0.080 with the numeric field and 0.0042 with the factor. One data set could be luck, so here is the same design drawn 1000 times: with the field effects as above, with the same four effects sorted so that they rise with the field number, and with the four effects put in a new random order for every data set:

p_treatment <- function(effects) {
  y <- treat_means[fields$treatment] + effects[fields$field] + rnorm(nrow(fields), sd = 15)
  c(numeric = anova(lm(y ~ field + treatment, data = fields))["treatment", "Pr(>F)"],
    factor  = anova(lm(y ~ factor(field) + treatment, data = fields))["treatment", "Pr(>F)"])
}
n_sim <- 1000
power_mixed  <- rowMeans(replicate(n_sim, p_treatment(field_effect)) < 0.05)
power_sorted <- rowMeans(replicate(n_sim, p_treatment(sort(field_effect))) < 0.05)
power_random <- rowMeans(replicate(n_sim, p_treatment(sample(field_effect))) < 0.05)
rbind(order_as_above = power_mixed, order_sorted = power_sorted, order_shuffled = power_random)

               numeric factor
order_as_above   0.160  0.951
order_sorted     0.953  0.954
order_shuffled   0.391  0.958

With the field effects as above, the treatment test found the effect (p below 0.05) in 0.951 of the data sets with factor(field) and in 0.160 with the numeric column. With the effects sorted, the straight line catches nearly all of them, and the factor and the numeric column find the effect in 0.954 and 0.953 of the data sets, the same within Monte Carlo error. The order used above is the worst of the 24 possible orders of these four effects (together with the same order reversed). With the order shuffled for every data set, the numeric column finds the effect in 0.391 of the data sets and the factor in 0.958. (Each share comes from 1000 data sets. Drawn again, a share would typically move by its Monte Carlo standard error, here at most 0.015.) So the damage comes from the order of the codes, and field numbers are labels: nothing makes them line up with productivity, so a numeric column loses power by an amount you cannot know in advance. The trap is for grouping columns in lm(), aov() and glm(). A random effect in lme4, (1 | field), is safe, because lme4 turns the grouping column into a factor itself.

What to check in your own data

Look at levels() of every grouping column before a model sees it, and make sure the first level is the group you want to compare against. If a column is still text when it goes into lm() or glm(), the reference is whatever comes first in the alphabet.

Run str() on the data too. A grouping column listed as int or num (field, block or site numbers) goes into lm() as one slope, not as groups; is.numeric() on every grouping column should be FALSE by the time a model sees it.

Before fitting an interaction of two grouping columns, run table() on the pair; a 0 means a combination the model cannot estimate, and any prediction for it depends on the reference levels.

Write the level list from the design (the sites surveyed, the months sampled, the treatments applied), not from unique() on the records. A site with no records can only appear in a summary if something other than the records knows about it.

When you count or summarise by group, compare the number of rows in the result with the number of groups in the design. If they differ, add .drop = FALSE and check the grouping column is a factor. With .drop = FALSE, an empty group gets n() equal to zero, but mean() of a measurement in that group has nothing to average and returns NaN, so decide what an empty site should show before you plot it.

Every time you call factor() with levels =, run setdiff(unique(x), levels) first and compare the count of NA values before and after. A difference is a spelling you have not dealt with.

Honest limits

Alphabetical order depends on the locale when labels mix capitals and lower case: factor(c("burnt", "Control")) can put either first on different machines, which is one more reason to write the levels out. With two grouping factors, .drop = FALSE keeps every combination of their levels, including combinations the design never had, so check the row count against the design there too.

References

Broman KW, Woo KH 2018 The American Statistician 72(1):2-10 (10.1080/00031305.2017.1375989)

R Core Team 2024 R documentation: Factors (https://stat.ethz.ch/R-manual/R-devel/library/base/html/factor.html)

R Core Team 2024 R documentation: Reorder levels of factor, relevel (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/relevel.html)

Wickham H, Francois R, Henry L, Muller K, Vaughan D 2023 dplyr documentation: Count the observations in each group (https://dplyr.tidyverse.org/reference/count.html)

Wickham H, Cetinkaya-Rundel M, Grolemund G 2023 R for Data Science, 2nd edition, chapter 16: Factors (https://r4ds.hadley.nz/factors.html)

--- title: "Factor levels in R: reference groups and empty sites" description: "Factor levels in R decide which group a model compares against and which sites vanish from a count. Set them from the design, keep empty sites, catch stray NAs." date: "2026-09-27 09:00" categories: [R, dplyr, data wrangling, data cleaning, regression, ecology tutorial] image: thumbnail.png image-alt: "Two panels of the same quadrat biomass data. With burnt as the reference level, arrows from the dashed reference line read +37.8 to control and -4.5 to mown; with control as the reference, they read -37.8 to burnt and -42.3 to mown." --- A fire experiment: quadrats that were left alone, quadrats that were burnt, quadrats that were mown, and the dry biomass cut from each at the end of the season. You fit `lm(biomass ~ treatment)`, open the coefficient table, and there is a row for `treatmentcontrol` but no row for burnt. Nothing failed. R has made a choice about your treatment column, and the table is an answer to a question you did not ask. A factor is R's way of storing a grouping column: a set of labels plus a list of the allowed values, the levels, in a fixed order. That list decides three things that show up in ecological results. Its first entry is the group every model coefficient is measured against. Its members are the groups that a count or a summary is allowed to report, including groups with nothing in them. And any value that is not on the list stops being a value at all. Each of the three main sections below runs one of these on a small simulated data set, shows the number that comes out wrong without any warning, and then fixes it. Two further sections add the cases next to them: an empty cell in a crossed design, where R does warn but quietly, and a grouping column stored as numbers, which never becomes a factor at all. **The short answer.** Turn every grouping column into a factor yourself, with the levels written out from your study design rather than taken from the data. Put the group you want to compare against first. Include the sites that produced no records. Then check that no value dropped out on the way. This post assumes you have met `group_by()` and `summarise()`, which [Summarising ecological data by group](../grouped-summaries-in-r/) introduces. The data are simulated, and every seed is in the code. ```{r setup} #| message: false library(dplyr) library(ggplot2) options(scipen = 6) # print small p-values as decimals te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), axis.text = element_text(colour = te_body), strip.text = element_text(colour = te_ink, face = "bold")) } ``` ## The first level is the reference group Eight quadrats per treatment, with true mean biomass fixed in the code at 120 g for control, 90 g for burnt and 95 g for mown: mowing is built to remove almost as much as fire does. The treatment column is plain text, which is what `read.csv()` and `data.frame()` give you: since R 4.0.0 neither turns text into factors unless you ask. ```{r plots} set.seed(1012) treat_means <- c(control = 120, burnt = 90, mown = 95) # true means, g per quadrat n_per <- 8 plots <- data.frame(treatment = rep(names(treat_means), each = n_per)) plots$biomass <- round(rnorm(nrow(plots), mean = treat_means[plots$treatment], sd = 15), 1) stopifnot(is.character(plots$treatment), is.character(read.csv(text = "treatment\nburnt")$treatment)) head(plots, 3) ``` `lm()` needs groups, so it converts the text column to a factor on the way in. `factor()` sorts the distinct values to make the levels, the first sorted level becomes the reference, and "burnt" comes before "control" in the alphabet. ```{r fit-default} stopifnot(identical(levels(factor(plots$treatment)), c("burnt", "control", "mown"))) fit_default <- lm(biomass ~ treatment, data = plots) cf_def <- coef(summary(fit_default)) round(cf_def, 4) ``` Read the table the way the model means it. The intercept, `r sprintf("%.1f", cf_def["(Intercept)", "Estimate"])` g, is the mean of the burnt quadrats. `treatmentcontrol` is control minus burnt, `r sprintf("%+.1f", cf_def["treatmentcontrol", "Estimate"])` g. And `treatmentmown` is mown minus burnt: `r sprintf("%+.1f", cf_def["treatmentmown", "Estimate"])` g with a p-value of `r sprintf("%.2f", cf_def["treatmentmown", "Pr(>|t|)"])`. Read as "the effect of mowing", that row says mowing did very little. It says nothing of the sort. It compares two disturbance treatments with each other, and the comparison the experiment was built for, mown against untouched, is not in the table at all. The fix is one line. Make the column a factor and name the reference, either with `relevel()` or by writing out all the levels in the order you want (useful later for plots and tables): ```{r fit-control} plots$treatment <- relevel(factor(plots$treatment), ref = "control") # same result: factor(plots$treatment, levels = c("control", "burnt", "mown")) levels(plots$treatment) fit_control <- lm(biomass ~ treatment, data = plots) round(coef(summary(fit_control)), 4) ``` ```{r control-numbers} #| include: false cf_ctl <- coef(summary(fit_control)) stopifnot(inherits(try(relevel(c("a", "b"), ref = "b"), silent = TRUE), "try-error"), cf_def["treatmentmown", "Pr(>|t|)"] > 0.05, cf_ctl["treatmentmown", "Pr(>|t|)"] < 0.0001) ``` Now the intercept is the control mean, `r sprintf("%.1f", cf_ctl["(Intercept)", "Estimate"])` g, and each row is a treatment against it: burnt `r sprintf("%+.1f", cf_ctl["treatmentburnt", "Estimate"])` g, mown `r sprintf("%+.1f", cf_ctl["treatmentmown", "Estimate"])` g with a p-value below 0.0001. The same word, mown, sat on a row with a small estimate and a large p-value a moment ago, and now sits on a row with a large negative estimate. `relevel()` only works on a factor; on a text column it stops with an error, which is why the `factor()` call comes first. Nothing about the model itself has changed. The coefficients are differences between the same three group means, taken from a different starting group. (How a factor becomes an intercept plus differences is in [t-tests and ANOVA as linear models in R](../linear-models-t-test-anova/).) ```{r same-fit} all.equal(fitted(fit_default), fitted(fit_control)) max_fit_diff <- max(abs(fitted(fit_default) - fitted(fit_control))) round(c(logLik_default = as.numeric(logLik(fit_default)), logLik_control = as.numeric(logLik(fit_control)), r2_default = summary(fit_default)$r.squared, r2_control = summary(fit_control)$r.squared), 4) grp_mean <- tapply(plots$biomass, plots$treatment, mean) round(grp_mean, 2) ``` The fitted values of the two models agree to within floating-point rounding. The log-likelihood is `r sprintf("%.4f", as.numeric(logLik(fit_control)))` both times and R-squared is `r sprintf("%.4f", summary(fit_control)$r.squared)` both times. The reference level changes what each coefficient means; it does not change the fit, the predictions or any test of the treatment as a whole (the next section shows one exception for the predictions: an empty cell in a crossed design). That is why the mistake is easy to miss: every diagnostic you run on the model looks the same. [Contrasts and post-hoc comparisons in R](../contrasts-and-post-hoc/) goes on to other codings (sum-to-zero, planned contrasts) and to comparing every pair. ```{r fp-guard-1} #| include: false stopifnot(max_fit_diff / max(abs(fitted(fit_control))) < 1e-10) ``` ```{r fig-reference} #| echo: false #| fig-width: 7.5 #| fig-height: 4.2 #| fig-cap: "The same biomass data, the same three group means and the same fitted model, read against two reference levels. The dashed line is the intercept and each arrow is one coefficient, drawn from the reference mean to the group mean. With burnt as the reference, the mown arrow is short; with control as the reference it is long and points down." #| fig-alt: "Two panels of quadrat biomass. Each shows the eight quadrats per treatment as pale green points, a dark green diamond at each treatment mean, a dashed horizontal line at the reference mean and vertical rust arrows labelled with the coefficient. Left panel, reference burnt near 89 g: arrows of plus 37.8 to control and minus 4.5 to mown. Right panel, reference control near 127 g: arrows of minus 37.8 to burnt and minus 42.3 to mown." panels <- c(burnt = "reference: burnt (alphabetical)", control = "reference: control (set by hand)") orders <- list(burnt = c("burnt", "control", "mown"), control = c("control", "burnt", "mown")) pts <- do.call(rbind, lapply(names(panels), function(r) data.frame(panel = panels[[r]], treatment = as.character(plots$treatment), x = match(as.character(plots$treatment), orders[[r]]), biomass = plots$biomass))) means <- do.call(rbind, lapply(names(panels), function(r) data.frame(panel = panels[[r]], treatment = orders[[r]], x = 1:3, y = as.numeric(grp_mean[orders[[r]]]), ref = as.numeric(grp_mean[[r]])))) arrows <- subset(means, x > 1) arrows$label <- sprintf("%+.1f", arrows$y - arrows$ref) ggplot(pts, aes(x, biomass)) + geom_jitter(width = 0.08, height = 0, colour = "#93a87f", size = 1.8, alpha = 0.8) + geom_hline(data = means[means$x == 1, ], aes(yintercept = ref), linetype = "dashed", colour = te_ink, linewidth = 0.5) + geom_segment(data = arrows, aes(x = x + 0.22, xend = x + 0.22, y = ref, yend = y), colour = te_rust, linewidth = 0.9, arrow = arrow(length = unit(0.18, "cm"))) + geom_text(data = arrows, aes(x = x + 0.3, y = (ref + y) / 2, label = label), hjust = 0, colour = te_rust, size = 3.6) + geom_point(data = means, aes(x, y), shape = 18, size = 4.2, colour = te_forest) + geom_text(data = means, aes(x = x, y = 52, label = treatment), colour = te_ink, size = 3.8) + facet_wrap(~ panel) + scale_x_continuous(breaks = NULL, limits = c(0.6, 3.75)) + scale_y_continuous(breaks = c(75, 100, 125, 150)) + labs(x = NULL, y = "biomass (g per quadrat)") + theme_datasheet() ``` The same levels set the order of groups in plots and tables: alphabetical unless you say otherwise. To sort groups by a statistic instead, for bars in order of their mean, `reorder()` builds the level order from the data: ```{r reorder} levels(reorder(plots$treatment, plots$biomass, FUN = mean)) ``` ## Unless a cell of the design is empty The reference level did not change any prediction above because every treatment had quadrats in it, so every prediction was the mean of quadrats that were cut. With two grouping columns crossed, that can fail. A second survey cuts biomass in three habitats in three seasons, four quadrats for each habitat and season, except that the wetland was under water in spring and its quadrats were never cut. The true means are fixed in the code, with a wetland that is poor in spring and peaks in autumn. ```{r survey} set.seed(1015) cell_mean <- rbind(forest = c(spring = 60, summer = 80, autumn = 70), # true means, g grassland = c(spring = 90, summer = 140, autumn = 100), wetland = c(spring = 70, summer = 130, autumn = 150)) survey <- expand.grid(quadrat = 1:4, season = colnames(cell_mean), habitat = rownames(cell_mean), stringsAsFactors = FALSE) survey <- subset(survey, !(habitat == "wetland" & season == "spring")) # flooded, never cut survey$biomass <- round(rnorm(nrow(survey), cell_mean[cbind(survey$habitat, survey$season)], sd = 15), 1) table(survey$habitat, survey$season) ``` The model `biomass ~ habitat * season` asks for a separate mean in every habitat and season, the interaction. Fit it twice: once with the text columns as they are, so forest and autumn are the references because they come first in the alphabet, and once with grassland and summer set as the references by hand. Then ask both fits for the cell that was never cut. ```{r empty-cell} fit_alpha <- lm(biomass ~ habitat * season, data = survey) survey_hand <- survey survey_hand$habitat <- factor(survey$habitat, levels = c("grassland", "forest", "wetland")) survey_hand$season <- factor(survey$season, levels = c("summer", "spring", "autumn")) fit_hand <- lm(biomass ~ habitat * season, data = survey_hand) wet_spring <- data.frame(habitat = "wetland", season = "spring") pred_alpha <- suppressWarnings(predict(fit_alpha, wet_spring)) # the warning is shown below pred_hand <- suppressWarnings(predict(fit_hand, wet_spring)) round(c(forest_autumn_first = as.numeric(pred_alpha), grassland_summer_first = as.numeric(pred_hand)), 2) ``` Both fits give the same fitted value for every quadrat that was cut, and they still predict `r sprintf("%.2f", pred_alpha)` g and `r sprintf("%.2f", pred_hand)` g for wetland in spring, a cell built at 70 g. Neither number comes from a wetland quadrat in spring, because there are none. Each is put together from three cells that were cut: the wetland mean in the reference season, plus the reference habitat in spring, minus the reference habitat in the reference season. For the alphabetical fit that is wetland in autumn plus forest in spring minus forest in autumn: ```{r three-cells} cell_avg <- function(h, s) mean(survey$biomass[survey$habitat == h & survey$season == s]) three_cells <- cell_avg("wetland", "autumn") + cell_avg("forest", "spring") - cell_avg("forest", "autumn") round(c(three_cells = three_cells, alphabetical_fit = as.numeric(pred_alpha)), 2) ``` ```{r empty-cell-checks} #| include: false pred_na <- predict(fit_alpha, wet_spring, rankdeficient = "NA") pred_add <- predict(lm(biomass ~ habitat + season, data = survey), wet_spring) stopifnot(isTRUE(all.equal(fitted(fit_alpha), fitted(fit_hand))), isTRUE(all.equal(as.numeric(pred_alpha), cell_avg("wetland", "autumn") + # three-cell identity cell_avg("forest", "spring") - cell_avg("forest", "autumn"))), isTRUE(all.equal(as.numeric(pred_hand), cell_avg("wetland", "summer") + cell_avg("grassland", "spring") - cell_avg("grassland", "summer"))), "rankdeficient" %in% names(formals(predict.lm)), is.na(pred_na), isTRUE(all.equal(pred_add, predict(lm(biomass ~ habitat + season, data = survey_hand), wet_spring))), sum(is.na(coef(fit_alpha))) == 1, # the three signs shown below length(grep("singularities", capture.output(summary(fit_alpha)))) == 1, inherits(tryCatch(predict(fit_alpha, wet_spring), warning = identity), "warning")) ``` It is the prediction of the alphabetical fit. For the fit set by hand the three cells are wetland in summer plus grassland in spring minus grassland in summer. This is arithmetic, not something the simulation found, and it means the gap between the two predictions is set by how differently the habitats change between seasons, the very thing the interaction was meant to estimate. R does flag the problem, in three places that are easy to scroll past: ```{r empty-cell-signs} coef(fit_alpha)[is.na(coef(fit_alpha))] grep("singularities", capture.output(summary(fit_alpha)), value = TRUE) tryCatch(predict(fit_alpha, wet_spring), warning = conditionMessage) ``` An `NA` coefficient means R dropped a column it could not estimate. [Fitting a senescence model to individual records](../fitting-a-senescence-model/) meets one of the same kind, where age is tied to female and season, and the slope R does report is what one arbitrary constraint happens to produce. Since R 4.3.0, `predict()` on an `lm()` fit takes `rankdeficient = "NA"` and returns `NA` for such a row instead of a number, which is the honest answer here: wetland in spring was not measured. If a number is needed anyway, fit `biomass ~ habitat + season` without the interaction. It predicts `r sprintf("%.1f", pred_add)` g whatever the reference levels, against the 70 g the cell was built at, because it assumes the wetland differs from the other habitats by the same amount in every season. That is an assumption to state in the methods, and in this simulation it is false by construction. [Multilevel post-stratification of records](../poststratifying-unstructured-records/) fills the empty cells of a survey frame with a main-effects model and with a multilevel model, and measures what each costs. ## Sites with no records drop out of a count A season of camera trapping for pine marten at twenty sites, stored the way detection data usually is: one row per detection. A site where the camera ran all season and saw nothing has no row. Each site's expected number of detections is drawn at random, so some sites end up empty. ```{r records} set.seed(1013) site_ids <- sprintf("S%02d", 1:20) # the survey design site_rate <- exp(rnorm(length(site_ids), log(1.5), 1)) n_det <- rpois(length(site_ids), site_rate) records <- data.frame(site = rep(site_ids, times = n_det)) c(detections = nrow(records), sites_surveyed = length(site_ids), sites_with_none = sum(n_det == 0)) ``` Count detections per site and summarise: ```{r count-text} per_site_text <- count(records, site) c(sites_in_table = nrow(per_site_text), mean_detections = mean(per_site_text$n), share_detected = mean(per_site_text$n > 0)) ``` The table has `r nrow(per_site_text)` sites, not `r length(site_ids)`. The `r sum(n_det == 0)` sites with no detections are not in `records`, so no function working from `records` alone can know they exist. The mean comes out at `r sprintf("%.2f", mean(per_site_text$n))` detections per site, and the share of sites where the marten was detected comes out at `r sprintf("%.2f", mean(per_site_text$n > 0))`: every site in the table has at least one detection, because that is how a site gets into the table. Turning the column into a factor does not help if the levels come from the same records, since `factor()` can only list the values it sees: `nlevels(factor(records$site))` is `r nlevels(factor(records$site))`. The levels have to come from the design, the list of sites you actually surveyed. Even then, `count()` and `group_by()` drop empty levels by default, so the argument `.drop = FALSE` is the second half of the fix: ```{r count-factor} records$site <- factor(records$site, levels = site_ids) per_site_default <- count(records, site) per_site_all <- count(records, site, .drop = FALSE) stopifnot(nrow(per_site_default) == sum(n_det > 0), nrow(per_site_all) == length(site_ids), nrow(count(data.frame(site = as.character(records$site)), site, .drop = FALSE)) == sum(n_det > 0)) c(default = nrow(per_site_default), drop_false = nrow(per_site_all)) ``` The last line of the `stopifnot()` is the other half of the rule: `.drop = FALSE` does nothing for a text column, because a text column has no list of levels to keep. You need both, a factor with the design's levels and `.drop = FALSE`. Base R's `table()` keeps empty levels without being asked, and `group_by()` takes the same argument: ```{r keep-zeros} sum(table(records$site) == 0) per_site <- records |> group_by(site, .drop = FALSE) |> summarise(detections = n(), .groups = "drop") rbind(zeros_dropped = c(sites = nrow(per_site_text), mean = mean(per_site_text$n), detected = mean(per_site_text$n > 0)), zeros_kept = c(sites = nrow(per_site), mean = mean(per_site$detections), detected = mean(per_site$detections > 0))) ``` ```{r limits-checks} #| include: false # claims in the closing sections: n() is 0 and mean() is NaN for an empty level; # two factors with .drop = FALSE give every combination of levels two <- data.frame(site = factor(c("A", "A"), levels = c("A", "B")), season = factor(c("spring", "autumn"), levels = c("spring", "autumn")), x = c(1, 2)) empty_summary <- summarise(group_by(two, site, .drop = FALSE), n = n(), m = mean(x)) stopifnot(empty_summary$n[2] == 0, is.nan(empty_summary$m[2]), nrow(count(two, site, season, .drop = FALSE)) == 4, # closed form of the inflation: surveyed sites / sites with a detection isTRUE(all.equal(mean(per_site_text$n) / mean(per_site$detections), length(site_ids) / sum(n_det > 0)))) ``` With the empty sites back, the mean is `r sprintf("%.2f", mean(per_site$detections))` detections per site and the share of sites with a detection is `r sprintf("%.2f", mean(per_site$detections > 0))`. Dropping the zeros inflated the mean by `r sprintf("%.0f", 100 * (mean(per_site_text$n) / mean(per_site$detections) - 1))` per cent and turned a species detected at `r sum(per_site$detections > 0)` of `r length(site_ids)` sites into one detected at every site in the table. No simulation is needed for the first number: the same `r nrow(records)` detections are divided by `r nrow(per_site_text)` sites instead of `r length(site_ids)`, so the mean is always inflated by the number of surveyed sites over the number with a detection, minus one (here `r length(site_ids)`/`r sum(n_det > 0)` - 1). ```{r fig-empty-sites} #| echo: false #| fig-width: 7.5 #| fig-height: 3.8 #| fig-cap: "Detections per site for the twenty simulated camera sites. Red crosses mark the sites with no detections, which have no rows in the detection table. The dashed line is the mean over the sites that count() returns by default; the solid line is the mean over all twenty sites." #| fig-alt: "Bar chart of detections for sites S01 to S20. Twelve sites have green bars, the tallest at 14 detections for S05. Eight sites, S02, S07, S09, S12 and S17 to S20, have no bar and a red cross on the zero line. A dashed red horizontal line sits at 4.17 and a solid dark green line at 2.50." per_site$empty <- per_site$detections == 0 mean_lines <- data.frame(y = c(mean(per_site_text$n), mean(per_site$detections)), what = c("mean, empty sites dropped", "mean, all sites")) ggplot(per_site, aes(site, detections)) + geom_col(fill = te_forest, width = 0.7) + geom_point(data = subset(per_site, empty), shape = 4, size = 3, stroke = 1.3, colour = te_rust) + geom_hline(data = mean_lines, aes(yintercept = y, linetype = what, colour = what), linewidth = 0.8) + scale_linetype_manual(values = c("mean, all sites" = "solid", "mean, empty sites dropped" = "dashed")) + scale_colour_manual(values = c("mean, all sites" = te_forest, "mean, empty sites dropped" = te_rust)) + labs(x = NULL, y = "detections", linetype = NULL, colour = NULL) + theme_datasheet() + theme(legend.position = "top", panel.grid.major.x = element_blank(), axis.text.x = element_text(size = 8.5)) ``` The same zeros can come back through a join instead, when the site list lives in its own table; [Joining ecological tables without losing zeros](../joining-ecological-tables/) takes that route. Either way the site list has to exist somewhere, and it can be lost on the way through a file: [What a CSV loses](../what-a-csv-loses/) shows an unsampled site disappearing from the level set when a factor is written to CSV and read back. ## Levels set by hand turn mismatches into NA Writing the levels out is the fix in both sections above, and it has a trap of its own. Here is a habitat column for 150 vegetation plots, typed by several people over a season. Most entries are clean, and a few carry a capital letter, a trailing space or a two-word spelling. ```{r habitats} set.seed(1014) n_rec <- 150 true_hab <- sample(c("forest", "grassland", "wetland"), n_rec, replace = TRUE, prob = c(0.5, 0.3, 0.2)) typed <- true_hab slip <- runif(n_rec) < 0.15 # entries typed differently variants <- list(forest = c("Forest", "forest "), grassland = c("Grassland", "grassland "), wetland = c("wet land", "Wetland")) typed[slip] <- vapply(true_hab[slip], function(h) sample(variants[[h]], 1), character(1)) table(typed) ``` You want the habitats in a fixed order, so you write the levels out: ```{r hand-levels} hab_levels <- c("forest", "grassland", "wetland") habitat <- factor(typed, levels = hab_levels) c(NA_in_typed = sum(is.na(typed)), NA_in_factor = sum(is.na(habitat))) ``` ```{r no-warning} #| include: false hab_warning <- NULL invisible(withCallingHandlers(factor(typed, levels = hab_levels), warning = function(w) hab_warning <<- conditionMessage(w))) share_true <- prop.table(table(true_hab)) share_hand <- prop.table(table(habitat)) lost <- table(factor(true_hab[is.na(habitat)], levels = hab_levels)) kept <- table(habitat) fair <- sum(is.na(habitat)) * share_true # losses in proportion to class size stopifnot(is.null(hab_warning), sum(table(habitat)) == sum(!is.na(habitat)), all(kept < table(true_hab)), lost[["wetland"]] < fair[["wetland"]], lost[["grassland"]] > fair[["grassland"]]) ``` `factor()` matches each value against the levels exactly, character for character, and every value without a match becomes `NA`. It does this without a warning. `r sum(is.na(habitat))` of the `r n_rec` plots, `r sprintf("%.0f", 100 * mean(is.na(habitat)))` per cent of the sheet, now have no habitat, although the column had no missing values a line earlier. The R documentation for `factor()` states this behaviour; it is not a bug, and nothing downstream will point back to it. What it does to a class share depends on where the mismatches fall: ```{r shares} round(cbind(true = share_true, after_levels = share_hand, plots_lost = lost, plots_left = kept), 3) ``` `table()` leaves `NA` out unless told otherwise, so every share is now a share of the `r sum(!is.na(habitat))` plots that survived. Wetland lost `r lost[["wetland"]]` of its `r sum(true_hab == "wetland")` plots, fewer than its fair part of the `r sum(is.na(habitat))` losses, so its share rose from `r sprintf("%.3f", share_true[["wetland"]])` to `r sprintf("%.3f", share_hand[["wetland"]])` although no wetland plot was added. Grassland lost `r lost[["grassland"]]` of its `r sum(true_hab == "grassland")`, more than its fair part, and its share fell from `r sprintf("%.3f", share_true[["grassland"]])` to `r sprintf("%.3f", share_hand[["grassland"]])`. In this simulation every habitat is mistyped at the same rate, so the shares move only a little. When one habitat has a spelling of its own (one recorder who always writes "wet land") the losses pile up in that class and its share drops with them. The counts are wrong regardless: here every habitat has fewer plots than were surveyed. The fix is a check before and a check after. Before, ask which values have no level to go to; clean them; ask again: ```{r fix-levels} setdiff(unique(typed), hab_levels) cleaned <- tolower(trimws(typed)) setdiff(unique(cleaned), hab_levels) cleaned[cleaned == "wet land"] <- "wetland" habitat <- factor(cleaned, levels = hab_levels) stopifnot(sum(is.na(habitat)) == sum(is.na(typed))) all(table(habitat) == table(factor(true_hab, levels = hab_levels))) ``` `trimws()` and `tolower()` handle the spaces and capitals, and the two-word spelling needs a decision, written down in the code. [Cleaning species names before you count](../cleaning-species-names/) does the same job for species names with a lookup table, which scales better than one line per variant. The lasting fix is on the sheet itself: Broman and Woo (2018) make one consistent code per category part of their first rule for spreadsheet data. The `stopifnot()` after the `factor()` call is the one to keep in every script: the number of missing values must not change when a column becomes a factor. The `forcats` package offers `fct()`, which stops with an error instead of creating `NA` (chapter 16 of R for Data Science shows the two side by side), if you prefer a loud failure. One more factor trap is common enough to name: `as.numeric()` on a factor returns the level positions, not the values. [Debugging and defensive R code](../debugging-and-defensive-code/) measures what that does to a regression on elevation. ## A grouping column stored as numbers Field numbers, block numbers and site codes are often typed as 1, 2, 3. R reads such a column as numbers, and a model given numbers fits one straight line through them instead of a separate mean for each group. Say the fire experiment's 24 quadrats sit in four fields, two quadrats of each treatment in each field, and the fields differ in productivity in an order that has nothing to do with their numbers. ```{r fields} set.seed(1016) field_effect <- c(0, 40, -30, 20) # g per quadrat, added in fields 1 to 4 fields <- data.frame(field = rep(1:4, each = 6), treatment = rep(rep(names(treat_means), each = 2), times = 4)) fields$biomass <- round(treat_means[fields$treatment] + field_effect[fields$field] + rnorm(nrow(fields), sd = 15), 1) str(fields) ``` `str()` shows the problem before any model does: `field` is `int`. Fit the treatment test with the field column as it is, and again as a factor: ```{r fields-fit} anova(lm(biomass ~ field + treatment, data = fields)) anova(lm(biomass ~ factor(field) + treatment, data = fields)) ``` ```{r fields-numbers} #| include: false p_num_one <- anova(lm(biomass ~ field + treatment, data = fields))["treatment", "Pr(>F)"] p_fac_one <- anova(lm(biomass ~ factor(field) + treatment, data = fields))["treatment", "Pr(>F)"] lin_share <- summary(lm(field_effect ~ I(1:4)))$r.squared # share of the spread a line removes perm <- expand.grid(rep(list(1:4), 4)); perm <- perm[apply(perm, 1, anyDuplicated) == 0, ] share_all <- apply(perm, 1, function(o) summary(lm(field_effect[o] ~ I(1:4)))$r.squared) stopifnot(is.integer(fields$field), p_fac_one < 0.05, p_num_one > 0.05, anova(lm(biomass ~ field + treatment, data = fields))["field", "Df"] == 1, anova(lm(biomass ~ factor(field) + treatment, data = fields))["factor(field)", "Df"] == 3, # a line through the field codes removes almost none of the effects in this order, # and nearly all of the same effects sorted; this order is the worst of 24 (with its reverse) lin_share < 0.05, summary(lm(sort(field_effect) ~ I(1:4)))$r.squared > 0.95, nrow(perm) == 24, sum(share_all <= lin_share + 1e-12) == 2) ``` The `Df` column gives it away: `field` takes 1 degree of freedom, one slope, and `factor(field)` takes 3, one difference for each field after the first. A straight line through field effects of 0, 40, -30 and 20 g removes almost none of them (`r sprintf("%.1f", 100 * lin_share)` per cent of their spread), so they stay in the residuals and the treatment test weighs its effect against that extra noise. In this data set the treatment p-value is `r sprintf("%.3f", p_num_one)` with the numeric field and `r sprintf("%.4f", p_fac_one)` with the factor. One data set could be luck, so here is the same design drawn 1000 times: with the field effects as above, with the same four effects sorted so that they rise with the field number, and with the four effects put in a new random order for every data set: ```{r fields-power} p_treatment <- function(effects) { y <- treat_means[fields$treatment] + effects[fields$field] + rnorm(nrow(fields), sd = 15) c(numeric = anova(lm(y ~ field + treatment, data = fields))["treatment", "Pr(>F)"], factor = anova(lm(y ~ factor(field) + treatment, data = fields))["treatment", "Pr(>F)"]) } n_sim <- 1000 power_mixed <- rowMeans(replicate(n_sim, p_treatment(field_effect)) < 0.05) power_sorted <- rowMeans(replicate(n_sim, p_treatment(sort(field_effect))) < 0.05) power_random <- rowMeans(replicate(n_sim, p_treatment(sample(field_effect))) < 0.05) rbind(order_as_above = power_mixed, order_sorted = power_sorted, order_shuffled = power_random) ``` ```{r fields-power-se} #| include: false all_power <- c(power_mixed, power_sorted, power_random) mc_se <- sqrt(all_power * (1 - all_power) / n_sim) stopifnot(power_mixed[["numeric"]] < power_mixed[["factor"]] - 10 * max(mc_se), abs(power_sorted[["numeric"]] - power_sorted[["factor"]]) < 2 * max(mc_se), power_random[["numeric"]] < power_random[["factor"]] - 10 * max(mc_se), power_random[["numeric"]] > power_mixed[["numeric"]] + 5 * max(mc_se)) ``` With the field effects as above, the treatment test found the effect (p below 0.05) in `r sprintf("%.3f", power_mixed[["factor"]])` of the data sets with `factor(field)` and in `r sprintf("%.3f", power_mixed[["numeric"]])` with the numeric column. With the effects sorted, the straight line catches nearly all of them, and the factor and the numeric column find the effect in `r sprintf("%.3f", power_sorted[["factor"]])` and `r sprintf("%.3f", power_sorted[["numeric"]])` of the data sets, the same within Monte Carlo error. The order used above is the worst of the 24 possible orders of these four effects (together with the same order reversed). With the order shuffled for every data set, the numeric column finds the effect in `r sprintf("%.3f", power_random[["numeric"]])` of the data sets and the factor in `r sprintf("%.3f", power_random[["factor"]])`. (Each share comes from `r n_sim` data sets. Drawn again, a share would typically move by its Monte Carlo standard error, here at most `r sprintf("%.3f", max(mc_se))`.) So the damage comes from the order of the codes, and field numbers are labels: nothing makes them line up with productivity, so a numeric column loses power by an amount you cannot know in advance. The trap is for grouping columns in `lm()`, `aov()` and `glm()`. A random effect in lme4, `(1 | field)`, is safe, because lme4 turns the grouping column into a factor itself. ## What to check in your own data Look at `levels()` of every grouping column before a model sees it, and make sure the first level is the group you want to compare against. If a column is still text when it goes into `lm()` or `glm()`, the reference is whatever comes first in the alphabet. Run `str()` on the data too. A grouping column listed as `int` or `num` (field, block or site numbers) goes into `lm()` as one slope, not as groups; `is.numeric()` on every grouping column should be `FALSE` by the time a model sees it. Before fitting an interaction of two grouping columns, run `table()` on the pair; a 0 means a combination the model cannot estimate, and any prediction for it depends on the reference levels. Write the level list from the design (the sites surveyed, the months sampled, the treatments applied), not from `unique()` on the records. A site with no records can only appear in a summary if something other than the records knows about it. When you count or summarise by group, compare the number of rows in the result with the number of groups in the design. If they differ, add `.drop = FALSE` and check the grouping column is a factor. With `.drop = FALSE`, an empty group gets `n()` equal to zero, but `mean()` of a measurement in that group has nothing to average and returns `NaN`, so decide what an empty site should show before you plot it. Every time you call `factor()` with `levels =`, run `setdiff(unique(x), levels)` first and compare the count of `NA` values before and after. A difference is a spelling you have not dealt with. ## Honest limits Alphabetical order depends on the locale when labels mix capitals and lower case: `factor(c("burnt", "Control"))` can put either first on different machines, which is one more reason to write the levels out. With two grouping factors, `.drop = FALSE` keeps every combination of their levels, including combinations the design never had, so check the row count against the design there too. ## References Broman KW, Woo KH 2018 The American Statistician 72(1):2-10 (10.1080/00031305.2017.1375989) R Core Team 2024 R documentation: Factors (https://stat.ethz.ch/R-manual/R-devel/library/base/html/factor.html) R Core Team 2024 R documentation: Reorder levels of factor, relevel (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/relevel.html) Wickham H, Francois R, Henry L, Muller K, Vaughan D 2023 dplyr documentation: Count the observations in each group (https://dplyr.tidyverse.org/reference/count.html) Wickham H, Cetinkaya-Rundel M, Grolemund G 2023 R for Data Science, 2nd edition, chapter 16: Factors (https://r4ds.hadley.nz/factors.html) ## Related tutorials - [Joining ecological tables without losing zeros](../joining-ecological-tables/) - [Contrasts and post-hoc tests](../contrasts-and-post-hoc/) - [What a CSV loses: archiving an R analysis](../what-a-csv-loses/) - [Cleaning species names before you count](../cleaning-species-names/)