coefplot.gllvm按因子而非水平绘图及参考水平问题咨询

问题背景 #

使用GLLVM包分析季节、水团对浮游生物分布的影响时，调用coefplot.gllvm()默认输出每个因子水平（如season_factorSpring）的独立子图，希望合并为季节、水团两个分组图；同时发现每个因子缺失一个水平，该水平为R因子处理时默认的参考水平，需解决显示问题。

解决方案 #

1. 按因子分组合并绘图 #

coefplot.gllvm()默认按单个参数（因子水平）拆分绘图，通过group参数可将同属一个因子的水平归为一组，生成合并后的子图：

# 查看模型系数的名称，确认因子水平前缀
names(coef(fit_env9)$params)

# 构建分组向量：将season_factor开头的参数归为"季节"组，water_factor开头的归为"水团"组
group_vec <- ifelse(grepl("season_factor", names(coef(fit_env9)$params)), 
                    "季节", "水团")

# 调用绘图函数，指定group参数，同时调整布局为1行2列
coefplot.gllvm(fit_env9, cex.ylab = 0.7, mar = c(4,9,2,1), 
               mfrow = c(1,2), group = group_vec)

2. 显示因子的参考水平 #

R中因子默认以第一个水平为参考水平，模型仅输出其他水平相对于参考水平的差值系数，因此参考水平不会显示。要显示所有水平，可通过设置总和对比（sum-to-zero） 修改因子的对比方式：

# 重新设置因子并指定总和对比
ALLData$water_factor <- factor(ALLData$water_type, levels = c("NW", "STW", "SAW"))
contrasts(ALLData$water_factor) <- contr.sum(3)  # 3个水平的总和对比

ALLData$season_factor <- factor(ALLData$season, levels = c("Spring", "Summer", "Autumn"))
contrasts(ALLData$season_factor) <- contr.sum(3)

# 重新拟合模型
fit_env9 <- gllvm(y, data = ALLData, family = "negative.binomial", num.lv = 1,
                  formula = ~ season_factor + water_factor, seed = 1234)

# 再次绘图，此时所有因子水平都会显示，同时按分组合并
coefplot.gllvm(fit_env9, cex.ylab = 0.7, mar = c(4,9,2,1), mfrow = c(1,2),
               group = c(rep("季节", 3), rep("水团", 3)))

额外提示 #

无需将协变量数据data_ev转换为矩阵，直接传入原始数据框给gllvm()的data参数，模型会自动处理因子类型，代码更简洁可靠。

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