R: glmrob не может предсказать модели с отброшенными коллинеарными столбцами, в то время как glm может?

Question

Я учусь реализовывать надежные glms в R, но не могу понять, почему я не могу заставить glmrob прогнозировать значения из моих регрессионных моделей, когда у меня есть модель, в которой некоторые столбцы отбрасываются из-за коллинеарности. В частности, когда я использую функцию прогнозирования для прогнозирования значений из glmrob, она всегда дает NA для всех значений. Я не наблюдаю это, когда прогнозирую значения из тех же данных и модели, используя glm. Кажется, не имеет значения, какие данные я использую - до тех пор, пока в подобранной модели есть коэффициент NA (а NA не является последним коэффициентом в векторе коэффициентов), прогноз не работает.

Такое поведение верно для всех наборов данных и моделей, которые я пробовал, когда внутренний столбец отбрасывается из-за коллинеарности. Я включаю ложный набор данных, в котором два столбца удаляются из модели, что дает два NA в списке коэффициентов. И glm, и glmrob дают почти одинаковые коэффициенты, но прогнозирование работает только с моделью glm. Итак, мой вопрос: что я не понимаю о надежной регрессии, которая помешала бы моим моделям glmrob генерировать предсказанные значения?

library(robustbase)

#Make fake data with two categorial predictors
df <- data.frame("category" = rep(c("A","B","C"),each=6))
df$location <- rep(1:6,each=3)
val <- rep(c(500,50,5000),each=6)+rep(c(50,100,25,200,100,1),each=3)
df$value <- rpois(NROW(df),val)

#note that predict works if we omit the newdata parameter. However I need the newdata param
#so I use the original dataframe here as a stand-in.  
mod <- glm(val ~ category + as.factor(location), data=df, family=poisson)
predict(mod, newdata=df) # works fine

mod <- glmrob(val ~ category + as.factor(location), data=df, family=poisson)
predict(mod, newdata=df) #predicts NA for all values

Answer 1

Я копался в этом и пришел к выводу, что проблема не в моем понимании робастной регрессии, а скорее в ошибке в пакете robustbase. Функция pregnet.lmrob неправильно выбирает необходимые коэффициенты из модели перед прогнозированием. Нужно выбрать первые x не-коэффициентов NA (где x= ранг модельной матрицы). Вместо этого он просто выбирает первые x коэффициентов, не проверяя, являются ли они NA. Это объясняет, почему эта проблема появляется только для моделей, где NA не является последним коэффициентом в векторе коэффициентов.

Чтобы это исправить, я скопировал исходный код Предиката, используя:

getAnywhere(predict.lmrob)

и создал свою собственную функцию замены. В этой функции я сделал одну модификацию кода:

...
p <- object$rank
if (is.null(p)) {
    df <- Inf
    p <- sum(!is.na(coef(object)))
    #piv <- seq_len(p) # old code
    piv <- which(!is.na(coef(object))) # new code
}
else {
    p1 <- seq_len(p)
    piv <- if (p) 
        qr(object)$pivot[p1]
}
...

Я запустил несколько сотен наборов данных, используя это изменение, и оно сработало хорошо.