我可以使用什么测试来比较两个或多个回归模型的斜率?

机器算法验证 r 数据可视化 多元分析 假设检验
2022-01-16 19:21:51

我想测试两个变量对一个预测变量的响应差异。这是一个最小的可重现示例。

library(nlme) 
## gls is used in the application; lm would suffice for this example
m.set <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "setosa")
m.vir <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "virginica")
m.ver <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "versicolor")

我可以看到斜率系数不同:

m.set$coefficients
(Intercept) Petal.Width 
  4.7771775   0.9301727
m.vir$coefficients
(Intercept) Petal.Width 
  5.2694172   0.6508306 
m.ver$coefficients
(Intercept) Petal.Width 
   4.044640    1.426365 

我有三个问题:

  1. 如何测试斜率之间的差异?
  2. 如何测试残差之间的差异?
  3. 呈现这些比较的简单、有效的方法是什么?

一个相关问题,比较两个回归模型中的变量系数的方法,建议使用虚拟变量重新运行模型以区分斜率,是否有允许使用独立数据集的选项?

4个回答

要使用 R 代码回答这些问题,请使用以下内容:
1. 如何测试斜率之间的差异?
答案:检查来自 Petal.Width by Species 交互作用的 ANOVA p 值,然后使用 lsmeans::lstrends 比较斜率,如下所示。

library(lsmeans)
m.interaction <- lm(Sepal.Length ~ Petal.Width*Species, data = iris)
anova(m.interaction)
 Analysis of Variance Table

 Response: Sepal.Length
                      Df Sum Sq Mean Sq  F value Pr(>F)    
 Petal.Width           1 68.353  68.353 298.0784 <2e-16 ***
 Species               2  0.035   0.017   0.0754 0.9274    
 Petal.Width:Species   2  0.759   0.380   1.6552 0.1947    
 Residuals           144 33.021   0.229                    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Obtain slopes
m.interaction$coefficients
m.lst <- lstrends(m.interaction, "Species", var="Petal.Width")
 Species    Petal.Width.trend        SE  df   lower.CL upper.CL
 setosa             0.9301727 0.6491360 144 -0.3528933 2.213239
 versicolor         1.4263647 0.3459350 144  0.7425981 2.110131
 virginica          0.6508306 0.2490791 144  0.1585071 1.143154

# Compare slopes
pairs(m.lst)
 contrast                 estimate        SE  df t.ratio p.value
 setosa - versicolor    -0.4961919 0.7355601 144  -0.675  0.7786
 setosa - virginica      0.2793421 0.6952826 144   0.402  0.9149
 versicolor - virginica  0.7755341 0.4262762 144   1.819  0.1669

2. 如何测试残差之间的差异?
如果我理解这个问题,您可以将 Pearson 相关性与 Fisher 变换(也称为“Fisher r-to-z”)进行比较,如下所示。

library(psych)
library(data.table)
iris <- as.data.table(iris)
# Calculate Pearson's R
m.correlations <- iris[, cor(Sepal.Length, Petal.Width), by = Species]
m.correlations
# Compare R values with Fisher's R to Z
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("setosa", "versicolor"), .N])
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="virginica", V1], 
         n = iris[Species %in% c("setosa", "virginica"), .N])
paired.r(m.correlations[Species=="virginica", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("virginica", "versicolor"), .N])

3. 什么是呈现这些比较的简单、有效的方式?
“我们使用线性回归来比较每个物种的萼片长度与花瓣宽度的关系。我们没有发现I. Setosa (B = 0.9)、I. Versicolor (B = 1.4), 也不是I. Virginica (B = 0.6); F (2, 144) = 1.6, p = 0.19。Fisher 的 r-to-z 比较表明I. Setosa (r = 0.28) 的 Pearson 相关性是显着低于 (p = 0.02) 比I. Versicolor (r = 0.55)。类似地,I. Virginica (r = 0.28) 的相关性明显弱于 (p = 0.02) 观察到的I. Versicolor。”

最后,始终可视化您的结果!

plotly_interaction <- function(data, x, y, category, colors = col2rgb(viridis(nlevels(as.factor(data[[category]])))), ...) {
  # Create Plotly scatter plot of x vs y, with separate lines for each level of the categorical variable. 
  # In other words, create an interaction scatter plot.
  # The "colors" must be supplied in a RGB triplet, as produced by col2rgb().

  require(plotly)
  require(viridis)
  require(broom)

  groups <- unique(data[[category]])

  p <- plot_ly(...)

  for (i in 1:length(groups)) {
    groupData = data[which(data[[category]]==groups[[i]]), ]
    p <- add_lines(p, data = groupData,
                   y = fitted(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                   x = groupData[[x]],
                   line = list(color = paste('rgb', '(', paste(colors[, i], collapse = ", "), ')')),
                   name = groups[[i]],
                   showlegend = FALSE)
    p <- add_ribbons(p, data = augment(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                     y = groupData[[y]],
                     x = groupData[[x]],
                     ymin = ~.fitted - 1.96 * .se.fit,
                     ymax = ~.fitted + 1.96 * .se.fit,
                     line = list(color = paste('rgba','(', paste(colors[, i], collapse = ", "), ', 0.05)')), 
                     fillcolor = paste('rgba', '(', paste(colors[, i], collapse = ", "), ', 0.1)'),
                     showlegend = FALSE)
    p <- add_markers(p, data = groupData, 
                     x = groupData[[x]], 
                     y = groupData[[y]],
                     symbol = groupData[[category]],
                     marker = list(color=paste('rgb','(', paste(colors[, i], collapse = ", "))))
  }
  p <- layout(p, xaxis = list(title = x), yaxis = list(title = y))
  return(p)
}

plotly_interaction(iris, "Sepal.Length", "Petal.Width", "Species")

虹膜图

如何测试斜率之间的差异?

包括一个物种的假人,让它与交互,看看这个假人是否重要。为萼片长度,为踏板宽度,为三个物种的虚拟变量。比较模型PiLiPiS1,S2,S3

E(Li)=β0+β1Pi

使用允许对每个物种的影响不同的模型:Pi

E(Li)=α0+α1S2+α2S3+α4Pi+α5PiS2+α6PiS3

GLS 估计器是 MLE,第一个模型是第二个模型的子模型,因此您可以在此处使用似然比检验。可以使用该logLik函数提取可能性,并且测试的自由度将为,因为您已删除参数以到达子模型。44

什么是呈现比较的简单有效的方法?

我认为最吸引人的方法是将每个物种的回归线都绘制在同一轴上,可能带有基于标准误差的误差线。这将使物种之间的差异(或无差异)及其与的关系非常明显。Pi

编辑:我注意到正文中添加了另一个问题。所以,我要添加一个答案:

如何测试残差之间的差异?

为此,您需要对数据集进行分层并拟合单独的模型,因为我建议的基于交互的模型将限制残差方差在每个组中都相同。如果你拟合单独的模型,这个约束就会消失。在这种情况下,您仍然可以使用似然比检验(较大模型的似然现在是通过将三个独立模型的似然相加来计算的)。“空”模型取决于您要与之比较的内容

  • 如果您只想测试方差,同时保留主要影响,那么“空”模型应该是我上面写的具有交互作用的模型。测试的自由度为2

  • 如果您想与系数一起测试方差,那么空模型应该是我上面写的第一个模型。测试的自由度为6

我同意前面的建议。您应该为每个数据集拟合一个带有虚拟变量的多元回归模型。这将允许您测试截距是否不同。如果您还想知道斜率是否不同,那么您还需要包括虚拟变量和相关变量之间的交互作用。数据是独立的这一事实没有问题。请注意,如果它们既是独立的并且(例如)不同的物种,那么您将无法判断您发现的差异是由于不同的物种还是不同的数据集,因为它们完全混淆了。但是,没有测试/ get-out-of-jail-free 卡可以让您解决该问题,而无需收集新样本并再次运行您的研究。

此链接提供了对该问题的简单回答:

https://stat.ethz.ch/pipermail/r-sig-teaching/2011q4/000387.html