F-оценка и стандартизованная бета для матрицы ковариации, скорректированной на гетероскедастичность (hccm) в R

Question

У меня есть несколько моделей регрессии, которые провалили тесты 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 

Любые советы о том, как сообщить о них так они как можно более согласовываются со стандартами summary() и Anova(), выходными из R и car, а функция std_beta() из пакета sjmisc?

Answer 1

Если у кого-то еще есть этот вопрос, вот мое решение. Это не особенно элегантно, но работает.

Я просто использовал функцию для 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. Надеюсь, это немного сэкономит.