R：测试多元回归中不同方程的系数是否相等（使用 linearhypothesis()）？
R: testing whether a coefficient is equal across the different equations in a multivariate regression (using linearhypothesis())?
我有一个关于如何在 R 中比较多元回归系数的问题。
我进行了一项调查，测量了三种不同的态度（量表变量）。我的目标是估计受访者的一些特征（年龄、性别、教育和意识形态）是否可以解释他们的（positve/negative）态度。
有人建议我进行多元多元回归，而不是三元多元多元回归。我的多变量模型的代码是：
MMR <- lm(cbind(Attitude_1, Attitude_2, Attitude_3) ~
Age + Gender + Education + Ideological_position,
data = survey)
summary(MMR)
接下来我要做的是估计 'Gender' 的系数在三个单独的模型中是否具有统计显着性。
我在 Stata (https://stats.idre.ucla.edu/stata/dae/multivariate-regression-analysis/), but I don't have a license, so I have to find an alternative in R. I know a similar question has been asked here before (R - Testing equivalence of coefficients in multivariate multiple regression) 中找到了如何执行此操作的非常清楚的说明，但答案是 R 中不存在可用于此目的的包（或函数）。因为这个答案是几年前提供的，所以我想知道同时是否实现了一些新的包或功能。
更准确地说，我想知道我是否可以使用 linearHypothesis() 函数 (https://www.rdocumentation.org/packages/car/versions/3.0-11/topics/linearHypothesis)？我已经知道这个函数允许我测试，例如，性别系数是否等于教育系数：
linearhypothesis(MMR, c("GenderFemale", "EducationHigh-educated")
我是否也可以使用这个函数来测试方程建模中的性别系数Attitude_1是否等于方程建模中的性别系数Attitude_2或Attitude_3？
如有任何帮助，我们将不胜感激！
由于问题中提供的模型不可重现（缺少输入）让我们改用此模型。
fm0 <- lm(cbind(cyl, mpg) ~ wt + hp, mtcars)
我们将讨论两种方法，使用我们的线性假设，即 cyl 和 mpg 组的截距相同，wt 斜率相同且 hp 斜率相同。
1) Mean/Variance
在这种方法中，我们仅将整个比较基于系数及其方差协方差矩阵。
library(car)
v <- vcov(fm0)
co <- setNames(c(coef(fm0)), rownames(v))
h1 <- c("cyl:(Intercept) = mpg:(Intercept)", "cyl:wt = mpg:wt", "cyl:hp = mpg:hp")
linearHypothesis(NULL, h1, coef. = co, vcov. = v)
给予：
Linear hypothesis test
Hypothesis:
cyl:((Intercept) - mpg:(Intercept) = 0
cyl:wt - mpg:wt = 0
cyl:hp - mpg:hp = 0
Model 1: restricted model
Model 2: structure(list(), class = "formula", .Environment = <environment>)
Note: Coefficient covariance matrix supplied.
Df Chisq Pr(>Chisq)
1
2 3 878.53 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
为了解释 linearHypothesis 在做什么，请注意，在这种情况下，假设矩阵是 L <- t(c(1, -1)) %x% diag(3) 并给定 v 然后作为大样本近似值，我们得到 L %*% co 分布为 N(0, L %*% v %*% t(L))在零假设下，因此 t(L %*% co) %*% solve(L %*% v %*% t(L)) %*% L %*% co 分布为具有 nrow(L) 自由度的卡方分布。
L <- t(c(1, -1)) %>% diag(3)
nrow(L) # degrees of freedom
SSH <- t(L %*% co) %*% solve(L %*% v %*% t(L)) %*% L %*% co # chisq
p <- pchisq(SSH, nrow(L), lower.tail = FALSE) # p value
2) 长格式模型
使用这种方法（不等同于上面显示的第一种方法）将 mtcars 从宽格式转换为长格式，mt2。我们在最后展示了如何使用 reshape 或 pivot_longer 来做到这一点，但现在我们将只明确地形成它。将 lhs 定义为 fm0 公式左侧的 32x2 矩阵，即 cbind(cyl, mpg)。请注意，它的列名称是 c("cyl", "mpg")。将 lhs 列逐列串入一个 64 长向量的 cyl 列，然后是 mpg 列，这为我们提供了新的因变量 y。我们还形成了一个分组变量 g。与 y 的长度相同，表示 y 的相应元素来自 lhs 中的哪一列。
定义了mt2我们就可以形成fm1了。在形成 fm1 时，我们将使用基于 fm0 sigma 值的权重向量 w 来反映 cyl 和 mpg 这两个组具有由向量 sigma(fm0) 给出的不同 sigma 值的事实。
我们在下面显示 fm0 和 fm1 模型具有相同的系数，然后 运行 linearHypothesis。
library(car)
lhs <- fm0$model[[1]]
g. <- colnames(lhs)[col(lhs)]
y <- c(lhs)
mt2 <- with(mtcars, data.frame(wt, hp, g., y))
w <- 1 / sigma(fm0)[g.]^2
fm1 <- lm(y ~ g./(wt + hp) + 0, mt2, weights = w)
# note coefficient names
variable.names(fm1)
## [1] "g.cyl" "g.mpg" "g.cyl:wt" "g.mpg:wt" "g.cyl:hp" "g.mpg:hp"
# check that fm0 and fm1 have same coefs
all.equal(c(t(coef(fm0))), coef(fm1), check.attributes = FALSE)
## [1] TRUE
h2 <- c("g.mpg = g.cyl", "g.mpg:wt = g.cyl:wt", "g.mpg:hp = g.cyl:hp")
linearHypothesis(fm1, h2)
给予：
Linear hypothesis test
Hypothesis:
- g.cyl + g.mpg = 0
- g.cyl:wt + g.mpg:wt = 0
- g.cyl:hp + g.mpg:hp = 0
Model 1: restricted model
Model 2: y ~ g./(wt + hp) + 0
Res.Df RSS Df Sum of Sq F Pr(>F)
1 61 1095.8
2 58 58.0 3 1037.8 345.95 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
若L为假设矩阵，与(1)中的L相同，只是列重新排序，q为其行数，n为mt2的行数，则SSH/q分布F(q, n-q-1) 所以我们有：
n <- nrow(mt2)
L <- diag(3) %x% t(c(1, -1)) # note difference from (1)
q <- nrow(L)
SSH <- t(L %*% coef(fm1)) %*% solve(L %*% vcov(fm1) %*% t(L)) %*% L %*% coef(fm1)
SSH/q # F value
pf(SSH/q, q, n-q-1, lower.tail = FALSE) # p value
方差分析
linearHypothesis 的替代方法是定义简化模型，然后使用方差分析比较两个模型。 mt2 和 w 来自上面。没有使用包。
fm2 <- lm(y ~ hp + wt, mt2, weights = w)
anova(fm2, fm1)
给予：
Analysis of Variance Table
Model 1: y ~ hp + wt
Model 2: y ~ g./(wt + hp) + 0
Res.Df RSS Df Sum of Sq F Pr(>F)
1 61 1095.8
2 58 58.0 3 1037.8 345.95 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
从宽到长交替计算
形成 mt2 的另一种方法是使用 reshape 将 mtcars 从宽形式重塑为长形式。
mt2a <- mtcars |>
reshape(dir = "long", varying = list(colnames(lhs)), v.names = "y",
timevar = "g.", times = colnames(lhs)) |>
subset(select = c("wt", "hp", "g.", "y"))
或使用 tidyverse（其中的行顺序不同，但只要在形成 fm1 和 w 时始终使用 mat2b，那应该无关紧要。
library(dplyr)
library(tidyr)
mt2b <- mtcars %>%
select(mpg, cyl, wt, hp) %>%
pivot_longer(all_of(colnames(lhs)), names_to = "g.", values_to = "y")
我有一个关于如何在 R 中比较多元回归系数的问题。
我进行了一项调查，测量了三种不同的态度（量表变量）。我的目标是估计受访者的一些特征（年龄、性别、教育和意识形态）是否可以解释他们的（positve/negative）态度。
有人建议我进行多元多元回归，而不是三元多元多元回归。我的多变量模型的代码是：
MMR <- lm(cbind(Attitude_1, Attitude_2, Attitude_3) ~
Age + Gender + Education + Ideological_position,
data = survey)
summary(MMR)
接下来我要做的是估计 'Gender' 的系数在三个单独的模型中是否具有统计显着性。
我在 Stata (https://stats.idre.ucla.edu/stata/dae/multivariate-regression-analysis/), but I don't have a license, so I have to find an alternative in R. I know a similar question has been asked here before (R - Testing equivalence of coefficients in multivariate multiple regression) 中找到了如何执行此操作的非常清楚的说明，但答案是 R 中不存在可用于此目的的包（或函数）。因为这个答案是几年前提供的，所以我想知道同时是否实现了一些新的包或功能。
更准确地说，我想知道我是否可以使用 linearHypothesis() 函数 (https://www.rdocumentation.org/packages/car/versions/3.0-11/topics/linearHypothesis)？我已经知道这个函数允许我测试，例如，性别系数是否等于教育系数：
linearhypothesis(MMR, c("GenderFemale", "EducationHigh-educated")
我是否也可以使用这个函数来测试方程建模中的性别系数Attitude_1是否等于方程建模中的性别系数Attitude_2或Attitude_3？
如有任何帮助，我们将不胜感激！
由于问题中提供的模型不可重现（缺少输入）让我们改用此模型。
fm0 <- lm(cbind(cyl, mpg) ~ wt + hp, mtcars)
我们将讨论两种方法，使用我们的线性假设，即 cyl 和 mpg 组的截距相同，wt 斜率相同且 hp 斜率相同。
1) Mean/Variance
在这种方法中，我们仅将整个比较基于系数及其方差协方差矩阵。
library(car)
v <- vcov(fm0)
co <- setNames(c(coef(fm0)), rownames(v))
h1 <- c("cyl:(Intercept) = mpg:(Intercept)", "cyl:wt = mpg:wt", "cyl:hp = mpg:hp")
linearHypothesis(NULL, h1, coef. = co, vcov. = v)
给予：
Linear hypothesis test
Hypothesis:
cyl:((Intercept) - mpg:(Intercept) = 0
cyl:wt - mpg:wt = 0
cyl:hp - mpg:hp = 0
Model 1: restricted model
Model 2: structure(list(), class = "formula", .Environment = <environment>)
Note: Coefficient covariance matrix supplied.
Df Chisq Pr(>Chisq)
1
2 3 878.53 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
为了解释 linearHypothesis 在做什么，请注意，在这种情况下，假设矩阵是 L <- t(c(1, -1)) %x% diag(3) 并给定 v 然后作为大样本近似值，我们得到 L %*% co 分布为 N(0, L %*% v %*% t(L))在零假设下，因此 t(L %*% co) %*% solve(L %*% v %*% t(L)) %*% L %*% co 分布为具有 nrow(L) 自由度的卡方分布。
L <- t(c(1, -1)) %>% diag(3)
nrow(L) # degrees of freedom
SSH <- t(L %*% co) %*% solve(L %*% v %*% t(L)) %*% L %*% co # chisq
p <- pchisq(SSH, nrow(L), lower.tail = FALSE) # p value
2) 长格式模型
使用这种方法（不等同于上面显示的第一种方法）将 mtcars 从宽格式转换为长格式，mt2。我们在最后展示了如何使用 reshape 或 pivot_longer 来做到这一点，但现在我们将只明确地形成它。将 lhs 定义为 fm0 公式左侧的 32x2 矩阵，即 cbind(cyl, mpg)。请注意，它的列名称是 c("cyl", "mpg")。将 lhs 列逐列串入一个 64 长向量的 cyl 列，然后是 mpg 列，这为我们提供了新的因变量 y。我们还形成了一个分组变量 g。与 y 的长度相同，表示 y 的相应元素来自 lhs 中的哪一列。
定义了mt2我们就可以形成fm1了。在形成 fm1 时，我们将使用基于 fm0 sigma 值的权重向量 w 来反映 cyl 和 mpg 这两个组具有由向量 sigma(fm0) 给出的不同 sigma 值的事实。
我们在下面显示 fm0 和 fm1 模型具有相同的系数，然后 运行 linearHypothesis。
library(car)
lhs <- fm0$model[[1]]
g. <- colnames(lhs)[col(lhs)]
y <- c(lhs)
mt2 <- with(mtcars, data.frame(wt, hp, g., y))
w <- 1 / sigma(fm0)[g.]^2
fm1 <- lm(y ~ g./(wt + hp) + 0, mt2, weights = w)
# note coefficient names
variable.names(fm1)
## [1] "g.cyl" "g.mpg" "g.cyl:wt" "g.mpg:wt" "g.cyl:hp" "g.mpg:hp"
# check that fm0 and fm1 have same coefs
all.equal(c(t(coef(fm0))), coef(fm1), check.attributes = FALSE)
## [1] TRUE
h2 <- c("g.mpg = g.cyl", "g.mpg:wt = g.cyl:wt", "g.mpg:hp = g.cyl:hp")
linearHypothesis(fm1, h2)
给予：
Linear hypothesis test
Hypothesis:
- g.cyl + g.mpg = 0
- g.cyl:wt + g.mpg:wt = 0
- g.cyl:hp + g.mpg:hp = 0
Model 1: restricted model
Model 2: y ~ g./(wt + hp) + 0
Res.Df RSS Df Sum of Sq F Pr(>F)
1 61 1095.8
2 58 58.0 3 1037.8 345.95 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
若L为假设矩阵，与(1)中的L相同，只是列重新排序，q为其行数，n为mt2的行数，则SSH/q分布F(q, n-q-1) 所以我们有：
n <- nrow(mt2)
L <- diag(3) %x% t(c(1, -1)) # note difference from (1)
q <- nrow(L)
SSH <- t(L %*% coef(fm1)) %*% solve(L %*% vcov(fm1) %*% t(L)) %*% L %*% coef(fm1)
SSH/q # F value
pf(SSH/q, q, n-q-1, lower.tail = FALSE) # p value
方差分析
linearHypothesis 的替代方法是定义简化模型，然后使用方差分析比较两个模型。 mt2 和 w 来自上面。没有使用包。
fm2 <- lm(y ~ hp + wt, mt2, weights = w)
anova(fm2, fm1)
给予：
Analysis of Variance Table
Model 1: y ~ hp + wt
Model 2: y ~ g./(wt + hp) + 0
Res.Df RSS Df Sum of Sq F Pr(>F)
1 61 1095.8
2 58 58.0 3 1037.8 345.95 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
从宽到长交替计算
形成 mt2 的另一种方法是使用 reshape 将 mtcars 从宽形式重塑为长形式。
mt2a <- mtcars |>
reshape(dir = "long", varying = list(colnames(lhs)), v.names = "y",
timevar = "g.", times = colnames(lhs)) |>
subset(select = c("wt", "hp", "g.", "y"))
或使用 tidyverse（其中的行顺序不同，但只要在形成 fm1 和 w 时始终使用 mat2b，那应该无关紧要。
library(dplyr)
library(tidyr)
mt2b <- mtcars %>%
select(mpg, cyl, wt, hp) %>%
pivot_longer(all_of(colnames(lhs)), names_to = "g.", values_to = "y")