非参数ANCOVA在R中报错求助：斜率差异分析与代码调试

研究需求与数据前提 #

• 分析两个区域（Area1、Area2）的蛋白浓度随时间（Age）的变化趋势，核心目标是判断二者回归斜率是否存在显著差异
• 已知两个区域的蛋白浓度均随时间线性上升，但Area1数据不服从正态分布，Area2不满足方差齐性，因此选择非参数ANCOVA作为分析方法

报错信息 #

调用sm.ancova函数时触发以下错误：

Error in optimise(f = get(fname), interval = log(c(start/8, start * 4)), :
invalid 'xmin' value

现有代码与数据 #

运行代码 #

nonpm <- sm.ancova(Area1, Area2, Age, model="parallel")

原始数据 #

Area1 <- c(19.3875, 44.689546, 39.43958333, 45.22222222, 57.15703125, 34.9421875, 32.6390625, 62.50298507, 60.5515625, 61.5140625, 72.221875, 48.778125) 

Area2 <- c(60.56111111, 62.536733, 75.42638889, 73.02777778, 144.340625, 42.0578125, 49.7578125, 79.37565399, 126.946875, 92.1765625, 130.521875, 125.7265625) 

Age <- c(119, 154, 174, 176, 183, 188, 190, 194, 230, 232, 245, 247) 

问题解析与解决方法 #

1. 非参数ANCOVA基础认知 #

非参数ANCOVA是参数ANCOVA的替代方案，无需依赖正态性、方差齐性等前提假设。它通过秩变换或非参数回归的方式，检验协变量（Age）调整后的组间差异，或组间回归斜率的差异，适配你的数据特征。

2. 修复sm.ancova报错 #

报错的核心原因是sm.ancova要求输入长格式数据，而你当前使用的是宽格式（两组数据分别作为独立变量传入）。需要先转换数据格式，再按函数要求的公式形式调用：

步骤1：转换为长格式数据

library(tidyverse)
# 构建长格式数据框
data_df <- data.frame(
  Age = rep(Age, 2),
  Concentration = c(Area1, Area2),
  Area = rep(c("Area1", "Area2"), each = length(Age))
)

步骤2：正确调用sm.ancova

注意：若要检验斜率差异，需设置model="general"（model="parallel"是假设斜率相同，仅检验截距差异）：

library(sm)
# 拟合非参数ANCOVA模型，检验斜率差异
nonpm <- sm.ancova(Concentration ~ Age * Area, data = data_df, model = "general")
# 查看分析结果
summary(nonpm)

3. 替代分析方案 #

如果sm包的方法仍存在问题，可尝试以下两种更稳定的非参数/稳健分析方法：

• 秩变换ANCOVA：先对因变量做秩变换，再执行参数ANCOVA，适配非正态、方差不齐的数据：

# 对蛋白浓度做秩变换
data_df$rank_conc <- rank(data_df$Concentration)
# 拟合模型并做方差分析
rank_ancova <- lm(rank_conc ~ Age * Area, data = data_df)
anova(rank_ancova)

• 稳健ANCOVA：使用robust包的稳健回归实现，抗异常值干扰，无需严格分布假设：

library(robust)
# 拟合稳健ANCOVA模型
robust_model <- lmRob(Concentration ~ Age * Area, data = data_df)
robust_ancova <- ancova(robust_model)
# 查看结果
summary(robust_ancova)

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