GLM模型显著交互项分析：lsmeans Tukey检验Contrast Error问题咨询

Let's break down how to fix that contrast error you're hitting with lsmeans after fitting your Poisson GLM with significant interaction terms. First, let's recap your model for context:

model <- glm(formula = ParticleCount ~ ParticlePresent + AlgaePresent + ParticleTypeSize + ParticlePresent:ParticleTypeSize + AlgaePresent:ParticleTypeSize, 
             family = poisson(link = "log"), 
             data = PCB)

And your significant interaction results:

Df Deviance AIC LRT Pr(>Chi)
666.94 1013.8
ParticlePresent:ParticleTypeSize 6 680.59 1015.4 13.649 0.033818 *
AlgaePresent:ParticleTypeSize 6 687.26 1022.1 20.320 0.002428 **

Most Likely Cause: Empty Factor Combinations (Zero Observations) #

The #1 reason for contrast errors with lsmeans is missing data in some factor level combinations (i.e., "empty cells" in your data). Tukey tests require valid mean estimates for every combination you're trying to contrast, and lsmeans can't calculate those if there are no observations for a group.

Step 1: Check for Empty Cells

First, verify which combinations have no data:

# Check ParticlePresent × ParticleTypeSize combinations
table(PCB$ParticlePresent, PCB$ParticleTypeSize)

# Check AlgaePresent × ParticleTypeSize combinations
table(PCB$AlgaePresent, PCB$ParticleTypeSize)

Look for cells with a value of 0—those are the problematic groups causing the error.

Step 2: Fix Empty Cells in lsmeans

You have two straightforward ways to handle this:

Option 1: Exclude Empty Combinations Directly

Use the exclude parameter to tell lsmeans to ignore groups with no data. Replace the example combinations below with the actual empty groups from your table output:

library(lsmeans)

# For ParticlePresent:ParticleTypeSize interaction
pp_pts_lsm <- lsmeans(model, ~ ParticlePresent:ParticleTypeSize, 
                      exclude = c("No:Small", "Yes:ExtraLarge")) # Replace with your empty groups
pairs(pp_pts_lsm, adjust = "tukey")

# For AlgaePresent:ParticleTypeSize interaction
ap_pts_lsm <- lsmeans(model, ~ AlgaePresent:ParticleTypeSize, 
                      exclude = c("Absent:Medium")) # Replace with your empty groups
pairs(ap_pts_lsm, adjust = "tukey")

Option 2: Specify Only Valid Combinations

If you have multiple empty groups, you can dynamically generate a list of valid (non-empty) combinations and pass them to lsmeans with the at parameter:

# Get valid ParticlePresent × ParticleTypeSize combinations
valid_pp_pts <- expand.grid(
  ParticlePresent = unique(PCB$ParticlePresent),
  ParticleTypeSize = unique(PCB$ParticleTypeSize)
)
# Filter out combinations with zero observations
valid_pp_pts <- valid_pp_pts[apply(valid_pp_pts, 1, function(x) {
  nrow(subset(PCB, ParticlePresent == x[1] & ParticleTypeSize == x[2])) > 0
}), ]

# Calculate lsmeans only for valid groups
pp_pts_lsm <- lsmeans(model, ~ ParticlePresent:ParticleTypeSize, at = valid_pp_pts)
pairs(pp_pts_lsm, adjust = "tukey")

# Repeat for AlgaePresent:ParticleTypeSize
valid_ap_pts <- expand.grid(
  AlgaePresent = unique(PCB$AlgaePresent),
  ParticleTypeSize = unique(PCB$ParticleTypeSize)
)
valid_ap_pts <- valid_ap_pts[apply(valid_ap_pts, 1, function(x) {
  nrow(subset(PCB, AlgaePresent == x[1] & ParticleTypeSize == x[2])) > 0
}), ]

ap_pts_lsm <- lsmeans(model, ~ AlgaePresent:ParticleTypeSize, at = valid_ap_pts)
pairs(ap_pts_lsm, adjust = "tukey")

Other Potential Fixes #

If empty cells aren't the issue, check these next:

1. Ensure Your Variables Are Factors

lsmeans relies on factor levels to calculate means. If any of your predictor variables are stored as character or numeric vectors, convert them to factors first:

# Convert to factors if needed
PCB$ParticlePresent <- as.factor(PCB$ParticlePresent)
PCB$AlgaePresent <- as.factor(PCB$AlgaePresent)
PCB$ParticleTypeSize <- as.factor(PCB$ParticleTypeSize)

# Refit the model and re-run lsmeans
model <- glm(formula = ParticleCount ~ ParticlePresent + AlgaePresent + ParticleTypeSize + ParticlePresent:ParticleTypeSize + AlgaePresent:ParticleTypeSize, 
             family = poisson(link = "log"), 
             data = PCB)

2. Adjust the Scale of Your Contrasts

By default, lsmeans calculates contrasts on the response scale (i.e., back-transformed from the log link). If you want contrasts on the log link scale instead, add type = "link":

pp_pts_lsm <- lsmeans(model, ~ ParticlePresent:ParticleTypeSize, type = "link")
pairs(pp_pts_lsm, adjust = "tukey")

内容的提问来源于stack exchange，提问作者tiptothetop