R 中计算线性回归的异方差稳健置信区间的函数

Question

下午好，我对下面的功能有疑问。任务是在 R 中开发一个函数，计算线性回归贝塔结果的异方差稳健置信区间。

正如我尝试这样做的那样，我的函数没有 return 任何输出。在尝试从中获取一些结果后，控制台根本不做任何事情。我真的争辩为什么特别是如果我通过我的代码的最后两行手动计算它一切正常。即使您没有必要的 data.frames，也许您可以看看我的代码并告诉我它有什么问题或提出解决我的问题的替代方法:)

为清楚起见：系数的原始数值（每个使用所有 200 个数据点）为 c(463.2121, 139.5762)，stdHC 为 c(74.705054, 5.548689)，如 lm 模型和 HC-稳健的标准错误我使用包三明治。

my_CI <- function (mod, level = 0.95)
{
  `%>%` <- magrittr::`%>%`
  standard_deviation <- stderrorHC
  Margin_Error <- abs(qnorm((1-0.95)/2))*standard_deviation 
  df_out <- data.frame(stderrorHC, mod,Margin_Error=Margin_Error,
                       'CI lower limit'=(mod - Margin_Error),
                       'CI Upper limit'=(mod + Margin_Error)) %>%
    return(df_out)
}

my_CI(mod, level = 0.95) #retrieving does not return any results for me

Definitions:
women <- read.table("women.txt")
men <- read.table("men.txt")
converged <- merge(women, men, all = TRUE)
level <- c(0.95, 0.975)
modell <- lm(formula = loan ~ education, data = converged)
mod <- modell$coefficients
vcov <- vcovHC(modell, type = "HC1")
stderrorHC <- sqrt(diag(vcov))

mod - abs(qnorm((1-level[1])/2))*stderrorHC 
mod + abs(qnorm((1-level[1])/2))*stderrorHC

补充：这是原始数据集中的一些数据。我只包含了十个数据点，因此在这种情况下我们需要在 t 分布上构建置信区间。

dataMenEductaion <- c(12, 17, 16, 11, 20, 20 , 11, 19, 15, 16)
dataMenLoan <- c(2404.72, 3075.313, 2769.543, 2009.295, 3105.121, 4269.216
                   2213.730, 4025.136, 2605.191, 2760.186)
dataWomenEducation <- c(12, 14, 16, 19 , 12, 19, 20, 17, 16, 10)
dataWomenLoan <- c(1920.667, 2278.255, 2296.804, 2977.048, 1915.740, 3557.991, 
                   3336.683, 2923.040, 2628.351, 1918.218)

Answer 1

我相信以下内容可以为您提供所需的输出。

# install.packages('sandwich')
library(sandwich) # contains vcovHC()

# data
df <- data.frame(education = c(12, 17, 16, 11, 20, 20, 11, 19, 15, 16,
                              12, 14, 16, 19 , 12, 19, 20, 17, 16, 10),
                loan = c(2404.72, 3075.313, 2769.543, 2009.295, 3105.121, 4269.216,
                         2213.730, 4025.136, 2605.191, 2760.186,
                         1920.667, 2278.255, 2296.804, 2977.048, 1915.740, 3557.991, 
                         3336.683, 2923.040, 2628.351, 1918.218))
df$sex <- factor(gl(2, nrow(df)/2, labels = c('males', 'females')))

# linear model
fit <- lm(loan ~ education + sex, data = df)
coefs <- fit$coefficients
vcov <- vcovHC(fit, type = "HC1")
stderrorHC <- sqrt(diag(vcov))

# function to compute robust SEs
my_CIs <- function (coefs, level = 0.95) {
  standard_deviation <- stderrorHC
  Margin_Error <- abs( qnorm( (1-level)/ 2) ) * standard_deviation 
  df_out <- data.frame(stderrorHC, coefs, Margin_Error = Margin_Error,
                       'CI lower limit' = (coefs - Margin_Error),
                       'CI Upper limit' = (coefs + Margin_Error))
  return(df_out)
}

输出

> my_CIs(coefs = coefs)
stderrorHC     coefs Margin_Error CI.lower.limit CI.Upper.limit
(Intercept)  295.86900  160.3716    579.89259      -419.5210      740.26416
education     23.64313  176.0111     46.33968       129.6714      222.35073
sexfemales   132.07169 -313.2632    258.85576      -572.1189      -54.40743

R 中计算线性回归的异方差稳健置信区间的函数

Function in R that computes heteroskedasticity-robust confidence intervals for a linear regression

r

function

linear-regression

confidence-interval

dataframe