用于鲁棒线性回归的哪个函数/包与glmulti一起工作(即,表现得像glm)？

Question

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

glmulti是一个R函数/包,用于一般线性模型的自动模型选择,它构造了给定因变量和一组预测变量的所有可能的一般线性模型,通过经典的 glm函数拟合它们,然后允许进行多模型推理(例如,使用从AICc,BIC得到的模型权重.glmulti在理论上也与任何其他函数一起工作,它返回系数,模型的对数似然和自由参数的数量(以及可能的其他信息？),与 glm的格式相同.

我的目标:具有强大错误的多模型推理

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

例如,我可以假设线性模型中的误差分布为t分布而不是正态分布.利用其峰度参数,t分布可以具有较重的尾部,因此对于异常值(与正态分布相比)更稳健.

但是,我不承诺使用t分配方法.我很满意任何回馈对数似然的方法,因此可以使用glmulti中的多模式方法.但是,这意味着,不幸的是我不能(例如,使用R中的知名强大的线性模型lmRob从稳健或lmrob从robustbase),因为他们没有对数似然框架下运作,因此无法一起工作glmulti.

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

RI的唯一强大的线性回归函数发现在对数似然框架下运行的是heavyLm(来自重包); 它用分布模拟错误.不幸的是,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)？

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

• 知道一个强大的线性回归函数,它(a)在对数似然框架下运行,(b)表现得像glm(因此可以与glmulti开箱即用).
• 有heavyLm已经有工作glmulti.

很感谢任何形式的帮助!

Answer 1

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

\n\n

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

\n\n

这是代码：

\n\n

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

\n\n

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

\n\n

Initialization...\nTASK: Exhaustive screening of candidate set.\nFitting...\nCompleted.\n\n> print(stackloss.glmulti)\nglmulti.analysis\nMethod: h / Fitting: glm / IC used: bic\nLevel: 1 / Marginality: FALSE\nFrom 8 models:\nBest IC: 117.892471265874\nBest model:\n[1] "stack.loss ~ 1 + Air.Flow + Water.Temp"\nEvidence weight: 0.709174196998897\nWorst IC: 162.083142797858\n2 models within 2 IC units.\n1 models to reach 95% of evidence weight.\n

\n\n

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

\n\n

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

\n