下午好，我對下面的功能有疑問。任務是在 R 中開發一個函式，用于計算線性回歸的 beta 結果的異方差穩健置信區間。
正如我試圖這樣做的那樣，我的函式不回傳任何輸出。控制臺在嘗試從中獲得一些結果后根本不做任何事情。我真的很爭論為什么特別是如果我通過代碼的最后兩行手動計算它，它作業得很好。即使您沒有必要的 data.frames，也許您可??以查看我的代碼并告訴我它有什么問題或提出解決我的問題的替代方法:)
為了清楚起見：系數的原始眾多值（每個使用所有 200 個資料點）是 c(463.2121, 139.5762)，stdHC 是 lm 模型給出的 c(74.705054, 5.548689)，對于我使用的 HC-robust 標準誤差包三明治。
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)
uj5u.com熱心網友回復：
我相信以下內容可為您提供所需的輸出。
# 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
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/406818.html
標籤：