R:在线性组合中添加常量，glht()

所以我试图复制我在 Hill、Griffiths 和 Lim 的《计量经济学原理》中看到的 stata 函数。我要复制的函数在 stata 中是这样的；

lincom _cons + b_1 * [arbitrary value] - c

这是零假设 H0 : B0 + B1*X = C

我可以在没有常数的情况下测试假设，但我想在测试参数的线性组合时添加常数。我浏览了 glht() 的包文档，但它只有一个示例，其中他们取出了常量。我复制了这个例子，保持常数，但我不确定当你有一个矩阵 K 和一个常数时如何测试线性组合。作为引用，这里是他们的例子；

### multiple linear model, swiss data
lmod <- lm(Fertility ~ ., data = swiss)

### test of H_0: all regression coefficients are zero
### (ignore intercept)

### define coefficients of linear function directly
K <- diag(length(coef(lmod)))[-1,]
rownames(K) <- names(coef(lmod))[-1]
K

### set up general linear hypothesis
glht(lmod, linfct = K)

我不太擅长创建假装数据集，但这是我的尝试。

library(multcomp)
test.data = data.frame(test.y = seq(200,20000,1000),
                   test.x = seq(10,1000,10))
test.data$test.y = sort(test.data$test.y + rnorm(100, mean = 10000, sd = 100)) -
  rnorm(100, mean = 5733, sd = 77)
test.lm = lm(test.y ~ test.x, data = test.data)

# to view the name of the coefficients
coef(test.lm)

# this produces an error. How can I add this intercept?
glht(test.lm,
 linfct = c("(Intercept) + test.x = 20"))

根据文档，似乎有两种方法可以解决这个问题。我可以使用函数 diag() 构造一个矩阵，然后我可以在 linfct = 参数中使用它，或者我可以使用字符串。这种方法的问题是，我不太知道如何使用 diag() 方法，同时还包括常量(等式的右侧)；在字符串方法的情况下，我不确定如何添加拦截。

任何和所有的帮助将不胜感激。

这是我正在处理的数据。这最初是在一个 .dta 文件中，所以我为可怕的格式道歉。根据我上面提到的那本书，这是 food.dta 文件。

structure(list(food_exp = structure(c(115.22, 135.98, 119.34,
114.96, 187.05, 243.92, 267.43, 238.71, 295.94, 317.78, 216,
240.35, 386.57, 261.53, 249.34, 309.87, 345.89, 165.54, 196.98,
395.26, 406.34, 171.92, 303.23, 377.04, 194.35, 213.48, 293.87,
259.61, 323.71, 275.02, 109.71, 359.19, 201.51, 460.36, 447.76,
482.55, 438.29, 587.66, 257.95, 375.73), label = "household food expenditure per week", format.stata = "%10.0g"),
income = structure(c(3.69, 4.39, 4.75, 6.03, 12.47, 12.98,
14.2, 14.76, 15.32, 16.39, 17.35, 17.77, 17.93, 18.43, 18.55,
18.8, 18.81, 19.04, 19.22, 19.93, 20.13, 20.33, 20.37, 20.43,
21.45, 22.52, 22.55, 22.86, 24.2, 24.39, 24.42, 25.2, 25.5,
26.61, 26.7, 27.14, 27.16, 28.62, 29.4, 33.4), label = "weekly  household income", format.stata = "%10.0g")), .Names = c("food_exp","income"), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -40L))

最佳答案

让我们从您的书中加载数据，然后检查我们的结果，以确保我们得到相同的结果。我以这种方式向您提供答案部分是因为它帮助我准确了解您想要什么，部分是为了让您相信这里的等效性。

我的部分困惑在于您的 lincom 示例的语法。你的语法可能是正确的，我不知道，但根据它的外观，我认为你在做一些不同的事情，因此引用你的书真的很有帮助。

首先，让我们加载数据并运行第 115 页上的线性模型:

install.packages("devtools")  # if not already installed
library(devtools)
install_git("https://github.com/ccolonescu/PoEdata")

library(PoEdata)   # loads the package in memory
library(multcomp)  # for hypo testing
data(food)         # loads the data set of interest

# EDA
summary(food)

# Model
mod <- lm(food_exp ~ income, data = food)
summary(mod) # Note: same results as PoE 4th ed. Pg 115 (other than rounding)

到现在为止还挺好。在 Pg。第 4 版的第 115 条显示了相同的回归模型，除了一些细微的舍入差异。

接下来，本书以家庭收入为 20(即 2,000 美元/周)为条件，计算每周食品支出的点估计:

# Point estimate
predict(mod, newdata = data.frame(income=20))

同样，我们得到了完全相同的结果。顺便说一句，您还可以在 Wiley 的 Using Stata for Principles of Econometrics 4th ed. 一书的免费样本中在 Stata 中看到相同的结果。

最后，我们准备进行假设检验。如前所述，我想确保我可以完全复制 Stata 所拥有的东西。您好心提供了代码，但我对您的语法有点困惑。

幸运的是，我们已经很幸运了。虽然第 4 版 Stata 指南的预览仅通过第 2 章，但荷兰的一所大学的 Economics and Business faculty 能够免费获得旧版的部分 DRM，因此我们可以引用:



最后看到我们可以像这样在 R 中复制它:

# Hyothesis Test
summary(glht(mod,  linfct = c("income = 15")))

不要被不同的输出格式所迷惑。它在 R 代码中显示的 estimate 只是回归模型中 income 的系数(“b2”)。它不会根据假设检验而改变，而在 Stata 输出中，它们执行“b2 - 15”(在 R 中是 mod$coefficients[2]-15 )。

改变的是 t ( t value ) 和 p ( Pr(>|t|) ) 值。请注意，这些来自 R 的测试统计数据与来自 Stata 的数据相匹配。

还有一个 H0 为收入 = 7.5 的例子，让我们看看 R 和 Stata 中的 t 值为 1.29，p 值为 0.203:

summary(glht(mod,  linfct = c("income = 7.5")))

您还可以使用 confint() 获得置信区间。

最后，我认为您正在查看您书中的第 3.6.4 节(第 117 页)，其中一位高管想要测试假设 income 为 20($2000/wk)，food_exp > 250:



我们可以计算 R 中的 t 值:

t = sum((mod$coefficients[1] + 20*mod$coefficients[2]-250)/sqrt(vcov(mod)[1] + 20^2 * vcov(mod)[4] + 2 * 20 *  vcov(mod)[2]))
t
[1] 2.652613

其中公式与本书前 2 页相同。

我们甚至可以把它变成一个自定义函数(适用于简单的线性回归，意味着只有 1 个自变量):

hypo_tester <- function(expenditure, income_per_week_hundreds, mod){
  t = sum((mod$coefficients[1] +
             income_per_week_hundreds*mod$coefficients[2]-expenditure)/sqrt(vcov(mod)[1] +
             income_per_week_hundreds^2 * vcov(mod)[4] + 2 * income_per_week_hundreds *  vcov(mod)[2]))
  return(t)
}

hypo_tester(250, 20, mod)
[1] 2.652613
hypo_tester(200, 20, mod)
[1] 6.179193
hypo_tester(300, 20, mod)
[1] -0.8739668

关于R:在线性组合中添加常量，glht()，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/46018122/