r - 哪个用于稳健线性回归的函数/包可与 glmulti 一起使用(即，行为类似于 glm)？

Question

背景:使用 glmulti 进行多模型推理

glmulti是一个 R 函数/包，用于一般线性模型的自动模型选择，在给定因变量和一组预测变量的情况下构造所有可能的一般线性模型，通过经典 glm 拟合它们函数，然后允许多模型推理(例如，使用从 AICc、BIC 导出的模型权重)。 glmulti 理论上也适用于任何其他返回系数、模型的对数似然和自由参数数量(也许还​​有其他信息？)的函数，其格式与 glm 确实如此。

我的目标:具有稳健错误的多模型推理

我想使用glmulti对定量因变量的误差进行稳健建模，以防止异常值的影响。

例如，我可以假设线性模型中的误差分布为 t distribution而不是正态分布。凭借其峰度参数，t 分布可以具有重尾，因此对异常值更加稳健(与正态分布相比)。

但是，我并不致力于使用 t 分布方法。我对任何返回对数似然的方法都很满意，因此可以与 glmulti 中的多模型方法一起使用。但这意味着，不幸的是我无法使用 R 中众所周知的稳健线性模型(例如，来自 robust 的 lmRob 或来自 robustbase 的 lmrob )，因为它们确实不在对数似然框架下运行，因此无法与glmulti一起使用。

问题:我找不到适用于 glmulti 的稳健回归函数

我发现在对数似然框架下运行的唯一稳健的 R 线性回归函数是 heavyLm(来自 heavy 包)；它使用 t 分布对误差进行建模。不幸的是，heavyLm 不能与 glmulti 一起使用(至少不是开箱即用的)，因为它没有用于 loglik 的 S3 方法(可能还有其他方法)的东西)。

举例说明:

library(glmulti)
library(heavy)

使用数据集stackloss

head(stackloss)

正则高斯线性模型:

summary(glm(stack.loss ~ ., data = stackloss))

使用 glm 的默认高斯链接函数通过 glmulti 进行多模型推理

stackloss.glmulti <- glmulti(stack.loss ~ ., data = stackloss, level=1, crit=bic)
print(stackloss.glmulti)
plot(stackloss.glmulti)

具有 t 分布误差的线性模型(默认为 df=4)

summary(heavyLm(stack.loss ~ ., data = stackloss))

使用glmulti调用heavyLm作为拟合函数进行多模型推理

stackloss.heavyLm.glmulti <- glmulti(stack.loss ~ .,
 data = stackloss, level=1, crit=bic, fitfunction=heavyLm)

出现以下错误:

    Initialization...
    Error in UseMethod("logLik") :
    no applicable method for 'logLik' applied to an object of class "heavyLm".

如果我定义以下函数，

logLik.heavyLm <- function(x){x$logLik}

glmulti 可以得到对数似然，但随后出现下一个错误:

    Initialization...
    Error in .jcall(molly, "V", "supplyErrorDF", 
    as.integer(attr(logLik(fitfunc(as.formula(paste(y,  : 
    method supplyErrorDF with signature ([I)V not found

问题:哪个用于稳健线性回归的函数/包可以与 glmulti 一起使用(即，行为类似于 glm)？

可能有一种方法可以定义更多函数来让 heavyLm 与 glmulti 一起工作，但在开始这个旅程之前，我想问是否有人

• 知道一个稳健的线性回归函数，该函数 (a) 在对数似然框架下运行，并且 (b) 行为类似于 glm(因此可以与 glmulti 一起使用) -开箱即用)。
• heavyLm 已与 glmulti 合作。

非常感谢任何帮助!

Answer 1

这是使用 heavyLm 的答案。尽管这是一个相对较老的问题，但在使用 HeavyLm 时仍然会出现您提到的相同问题(即错误消息 Error in .jcall(molly, "V", "supplyErrorDF"...)。

问题是 glmulti 需要模型的自由度，作为您需要提供的属性传递给函数 logLik.heavyLm 返回值的属性;有关详细信息，请参阅函数 logLik 的文档。此外，事实证明，您还需要提供一个函数来返回用于拟合模型的数据点的数量，因为信息标准(AIC、BIC 等)也取决于该值。这是通过下面代码中的函数 nobs.heavyLm 完成的。

这是代码:

nobs.heavyLm <- function(mdl) mdl$dims[1] # the sample size (number of data points)

logLik.heavyLm <- function(mdl) {
    res <- mdl$logLik
    attr(res, "nobs") <- nobs.heavyLm(mdl) # this is not really needed for 'glmulti', but is included to adhere to the format of 'logLik'
    attr(res, "df") <- length(mdl$coefficients) + 1 + 1 # I am also considering the scale parameter that is estimated; see mdl$family
    class(res) <- "logLik"
    res
}

当与您提供的代码放在一起时，会产生以下结果:

Initialization...
TASK: Exhaustive screening of candidate set.
Fitting...
Completed.

> print(stackloss.glmulti)
glmulti.analysis
Method: h / Fitting: glm / IC used: bic
Level: 1 / Marginality: FALSE
From 8 models:
Best IC: 117.892471265874
Best model:
[1] "stack.loss ~ 1 + Air.Flow + Water.Temp"
Evidence weight: 0.709174196998897
Worst IC: 162.083142797858
2 models within 2 IC units.
1 models to reach 95% of evidence weight.

因此在 2 个 BIC 单位阈值内生产 2 个模型。

但重要的是:我不确定上面的自由度表达式是否严格正确。对于标准线性模型，自由度等于 p + 1，其中 p 是模型中参数的数量，额外参数(+ 1 >) 是“误差”方差(用于计算可能性)。在上面的函数 logLik.heavyLm 中，我不清楚是否还应该将 heavyLm 估计的“尺度参数”算作额外的自由度，并且因此 p + 1 + 1，如果似然性也是该参数的函数，就会出现这种情况。不幸的是，我无法证实这一点，因为我无法访问 heavyLm 引用的引用文献(Dempster 等人的论文，1980)。因此，我正在计算尺度参数，从而提供模型复杂性的(稍微)保守估计，从而惩罚“复杂”模型。除了小样本情况外，这种差异应该可以忽略不计。