R 中异方差校正协方差矩阵 (hccm) 的 F 分数和标准化 Beta
F-score and standardized Beta for heteroscedasticity-corrected covariance matrix (hccm) in R
我有多个回归模型未通过 Breusch-Pagan 检验,因此我使用异方差校正协方差矩阵重新计算了方差,如下所示:coeftest(lm.model,vcov=hccm(lm.model))
。 coeftest()
来自 lmtest
包,而 hccm()
来自 car
包。
我想提供 F 分数和标准化贝塔值,但我不确定该怎么做,因为输出看起来像这样...
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.000261 0.038824 0.01 0.995
age 0.004410 0.041614 0.11 0.916
exercise -0.044727 0.023621 -1.89 0.059 .
tR -0.038375 0.037531 -1.02 0.307
allele1_num 0.013671 0.038017 0.36 0.719
tR:allele1_num -0.010077 0.038926 -0.26 0.796
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
有关如何报告这些的任何建议,以便它们尽可能与 R 和 car
的标准 summary()
和 Anova()
输出以及函数 std_beta()
保持一致] 来自 sjmisc
包?
万一其他人有这个问题,这是我的解决方案。它不是特别优雅,但它确实有效。
我只是将 std_beta 的函数用作模板,然后将标准错误的输入更改为从 std_beta()
函数派生的输入。
# This is taken from std_beta function from sj_misc package.
# =====================================
b <-coef(lm.model) # Same Estimate
b <-b[-1] # Same intercept
fit.data <- as.data.frame(stats::model.matrix(lm.model)) # Same model.
fit.data <- fit.data[, -1] # Keep intercept
fit.data <- as.data.frame(sapply(fit.data, function(x) if (is.factor(x))
to_value(x, keep.labels = F)
else x))
sx <- sapply(fit.data, sd, na.rm = T)
sy <- sapply(as.data.frame(lm.model$model)[1], sd, na.rm = T)
beta <- b * sx/sy
se <-coeftest(lm.model,vcov=hccm(lm.model))[,2] # ** USE HCCM covariance for SE **
se <- se[-1]
beta.se <- se * sx/sy
data.frame(beta = beta, ci.low = (beta - beta.se *
1.96), ci.hi = (beta + beta.se * 1.96))
对于 F 分数,我只是将 t 值平方。
我希望这可以节省一些时间。
我有多个回归模型未通过 Breusch-Pagan 检验,因此我使用异方差校正协方差矩阵重新计算了方差,如下所示:coeftest(lm.model,vcov=hccm(lm.model))
。 coeftest()
来自 lmtest
包,而 hccm()
来自 car
包。
我想提供 F 分数和标准化贝塔值,但我不确定该怎么做,因为输出看起来像这样...
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.000261 0.038824 0.01 0.995
age 0.004410 0.041614 0.11 0.916
exercise -0.044727 0.023621 -1.89 0.059 .
tR -0.038375 0.037531 -1.02 0.307
allele1_num 0.013671 0.038017 0.36 0.719
tR:allele1_num -0.010077 0.038926 -0.26 0.796
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
有关如何报告这些的任何建议,以便它们尽可能与 R 和 car
的标准 summary()
和 Anova()
输出以及函数 std_beta()
保持一致] 来自 sjmisc
包?
万一其他人有这个问题,这是我的解决方案。它不是特别优雅,但它确实有效。
我只是将 std_beta 的函数用作模板,然后将标准错误的输入更改为从 std_beta()
函数派生的输入。
# This is taken from std_beta function from sj_misc package.
# =====================================
b <-coef(lm.model) # Same Estimate
b <-b[-1] # Same intercept
fit.data <- as.data.frame(stats::model.matrix(lm.model)) # Same model.
fit.data <- fit.data[, -1] # Keep intercept
fit.data <- as.data.frame(sapply(fit.data, function(x) if (is.factor(x))
to_value(x, keep.labels = F)
else x))
sx <- sapply(fit.data, sd, na.rm = T)
sy <- sapply(as.data.frame(lm.model$model)[1], sd, na.rm = T)
beta <- b * sx/sy
se <-coeftest(lm.model,vcov=hccm(lm.model))[,2] # ** USE HCCM covariance for SE **
se <- se[-1]
beta.se <- se * sx/sy
data.frame(beta = beta, ci.low = (beta - beta.se *
1.96), ci.hi = (beta + beta.se * 1.96))
对于 F 分数,我只是将 t 值平方。 我希望这可以节省一些时间。