R - 多元多重回归中系数等价性测试

Question

我正试图进行多元多重回归分析。幸运的是，我找到了一个非常好的页面，展示了如何在Stata中进行这样的分析：

http://www.ats.ucla.edu/stat/stata/dae/mvreg.htm

问题是我正在使用R，虽然我已经学会了如何在R中进行基本的多元多重回归模型，但我仍然不确定如何查看每个因变量的系数是否不同（如链接所示）。有人知道如何在R中计算此分析吗？看起来检查同一自变量是否对每个因变量产生不同的影响是一项非常有用的工具，我很想能够做到！

更新：以下是我已经完成的关于自己数据的可复现示例：

# data
data(mtcars)

# fitting a multivariate multiple regression where mpg and cyl are predicted by wt and hp
car.mod <- lm(cbind(mpg,cyl) ~ wt + hp,data=mtcars)

# see if there is a multivariate effect of wt and hp
summary(manova(car.mod),test="W")

# Get coefficients for each dependent variable
summary(car.mod)

在这个例子中，我想知道如何测试“mpg”和“cyl”上的“wt”的等价性。显然，在Stata中可以使用“test”命令实现。

- user2917781

---

查看R标签信息页面以获取大量R特定资源，其中包括一些涵盖此主题的资源。不过你现在有什么？你需要添加一个可重现的示例。 - alistaire

好的，完成了！ - user2917781

Answer 1

据我所知，目前没有相应的包可以实现这个功能，因此我会采用手动解决方案。手动解决方案是：

z = (b1 - b2) / (b1^2 + b2^2)^(1/2).

以下是（粗糙的）代码。提取系数和标准误差可能有更优雅的解决方案。

# extract coefficients and SEs
coef_mpg <- coef(summary(car.mod))[[1]]
coef_cyl <- coef(summary(car.mod))[[2]]

# test equality of wt coefficients
z_wt <- (coef_mpg[2, 1] - coef_cyl[2, 1]) / (coef_mpg[2, 2]^2 + coef_cyl[2, 2]^2)^(1/2)
p_wt <- 2*pnorm(-abs(z_wt)) 
p_wt

但是我更希望有一个基于 bootstrap 的解决方案，它做出更少的假设。

require(boot)
b_b <- function(data=mtcars, indices) {
  d <- data[indices, ]
  model_mpg <- lm(mpg ~ wt + hp, data=d)
  model_cyl <- lm(cyl ~ wt + hp, data=d)
  return(coef(model_mpg)["wt"] - coef(model_cyl)["wt"])
}
results <- boot(data=mtcars, statistic=b_b, R=1000)
results

- Richard Herron

Answer 2

假设您有两个结果和多个预测变量，以下是一个示例。该公式由Clifford Clogg等人（1995）提出，并被Ray Paternoster等人（1998）引用。

下面是一些数据和R语言脚本，以说明@Richard Herron的过程。

df1 = data.frame(
  estimate = c(15.2418519, 2.2215987, 0.3889724, 0.5289710),
  `std.error` = c(1.0958919, 0.2487793, 0.1973446, 0.1639074),
  row.names = c('(Intercept)', 'psychoticism', 'extraversion', 'neuroticism')
  ); df1

df2 = data.frame(
  estimate = c(17.2373874, 0.8350460, -0.3714803, 1.0382513),
  `std.error` = c(1.0987151, 0.2494201, 0.1978530, 0.1643297),
  row.names = c('(Intercept)', 'psychoticism', 'extraversion', 'neuroticism')
); df2

所以在上面的数据示例中，我们假设您拥有一个多元多重回归模型。df1和df2涉及两个结果变量及其相关预测变量（这些变量相同）。

下面的脚本迭代相关数据，并在数据框中输出系数的比较。我已经对脚本进行了注释，因此很容易跟随并了解正在发生的事情。

compare_coefs <- function(.data1, .data2){
  
  # import map_dbl() and pluck() functions from the purrr library
  import::here(map_dbl, pluck, .from = purrr)

  # extract the coef and std error from df1
  b1 = map_dbl(.data1[-1, 1], pluck)
  se1 = map_dbl(.data1[-1, 2], pluck)

  # extract the coef and std error from df2
  b2 = map_dbl(.data2[-1, 1], pluck)
  se2 = map_dbl(.data2[-1, 2], pluck)
  
  # calculate the difference per Clogg et al (1995) formula 
  # as cited by Paternoster et al. (1998)
  b = b1 - b2
  s1 = se1^2
  s2 = se2^2
  sc = s1 + s2
  v = b / sqrt(sc)
  
  data.frame(diff=b, zdiff=v, `p-value`=format(2*pnorm(-abs(v)), scientific=FALSE))
  
  }

# call the function
compare_coefs(df1, df2)

- GSA