Метод предсказания для felm из пакета lfe

Ответ 1

В качестве обходного пути вы можете объединить felm, getfe и demeanlist следующим образом:

library(lfe)

lm.model <- lm(data=demeanlist(iris[, 1:2], list(iris$Species)), Sepal.Length ~ Sepal.Width)
fe <- getfe(felm(data = iris, Sepal.Length ~ Sepal.Width | Species))
predict(lm.model, newdata = data.frame(Sepal.Width = 3)) + fe$effect[fe$idx=="virginica"]

Идея состоит в том, что вы используете demeanlist для центрирования переменных, затем lm для оценки коэффициента на Sepal.Width с использованием центрированных переменных, предоставляя вам объект lm, над которым вы можете запустить predict. Затем запустите felm + getfe, чтобы получить условное среднее для фиксированного эффекта, и добавьте это к выходу predict.

Ответ 2

Это может быть не тот ответ, который вы ищете, но, похоже, автор не добавил никаких функций в пакет lfe, чтобы делать прогнозы по внешним данным, используя установленную модель felm. Основное внимание, по-видимому, уделяется анализу групповых фиксированных эффектов. Однако интересно отметить, что в документации пакета указано следующее:

Объект имеет некоторое сходство с объектом 'lm', а некоторые возможно, будут работать методы постпроцессинга, разработанные для lm. Это может однако необходимо, чтобы принудить объект к успеху с этим.

Следовательно, можно было бы принудить объект felm к объекту lm, чтобы получить дополнительную функциональность lm (если вся необходимая информация присутствует в объекте для выполнения необходимых вычислений).

Пакет lfe предназначен для работы с очень большими наборами данных, и было предпринято усилие для сохранения памяти: в результате этого объект felm не использует/не содержит qr-декомпозицию, в отличие от lm объект. К сожалению, процедура lm predict основана на этой информации для вычисления прогнозов. Следовательно, принудительное выполнение объекта felm и выполнение метода прогнозирования не удастся:

> model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species)
> class(model2) <- c("lm","felm") # coerce to lm object
> predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
Error in qr.lm(object) : lm object does not have a proper 'qr' component.
 Rank zero or should not have used lm(.., qr=FALSE).

Если вы действительно должны использовать этот пакет для выполнения прогнозов, вы могли бы написать свою собственную упрощенную версию этой функции, используя информацию, имеющуюся в объекте felm. Например, коэффициенты регрессии OLS доступны через model2$coefficients.

Ответ 3

Это должно работать в случаях, когда вы хотите игнорировать групповые эффекты в предсказании, предсказывают новые X и хотят только доверительные интервалы. Сначала он ищет атрибут clustervcv, затем robustvcv, затем vcv.

predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }

  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0

  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)

  if(se.fit | interval != "none"){
    if(!is.null(object$clustervcv)){
      vcov_mat <- object$clustervcv
    } else if (!is.null(object$robustvcv)) {
      vcov_mat <- object$robustvcv
    } else if (!is.null(object$vcv)){
      vcov_mat <- object$vcv
    } else {
      stop("No vcv attached to felm object.")
    }
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat)))
  }
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}

Ответ 4

Чтобы расширить ответ от pbaylis, я создал слегка длинную функцию, которая прекрасно расширяется, чтобы обеспечить более одного фиксированного эффекта. Обратите внимание, что вы должны вручную ввести исходный набор данных, используемый в модели felm. Функция возвращает список с двумя элементами: вектор предсказаний и фрейм данных на основе new_data, который включает предсказания и фиксированные эффекты в виде столбцов.

predict_felm <- function(model, data, new_data) {

  require(dplyr)

  # Get the names of all the variables
  y <- model$lhs
  x <- rownames(model$beta)
  fe <- names(model$fe)

  # Demean according to fixed effects
  data_demeaned <- demeanlist(data[c(y, x)],
                             as.list(data[fe]),
                             na.rm = T)

  # Create formula for LM and run prediction
  lm_formula <- as.formula(
    paste(y, "~", paste(x, collapse = "+"))
  )

  lm_model <- lm(lm_formula, data = data_demeaned)
  lm_predict <- predict(lm_model,
                        newdata = new_data)

  # Collect coefficients for fe
  fe_coeffs <- getfe(model) %>% 
    select(fixed_effect = effect, fe_type = fe, idx)

  # For each fixed effect, merge estimated fixed effect back into new_data
  new_data_merge <- new_data
  for (i in fe) {

    fe_i <- fe_coeffs %>% filter(fe_type == i)

    by_cols <- c("idx")
    names(by_cols) <- i

    new_data_merge <- left_join(new_data_merge, fe_i, by = by_cols) %>%
      select(-matches("^idx"))

  }

  if (length(lm_predict) != nrow(new_data_merge)) stop("unmatching number of rows")

  # Sum all the fixed effects
  all_fixed_effects <- base::rowSums(select(new_data_merge, matches("^fixed_effect")))

  # Create dataframe with predictions
  new_data_predict <- new_data_merge %>% 
    mutate(lm_predict = lm_predict, 
           felm_predict = all_fixed_effects + lm_predict)

  return(list(predict = new_data_predict$felm_predict,
              data = new_data_predict))

}

model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species)
predict_felm(model = model2, data = iris, new_data = data.frame(Sepal.Width = 3, Species = "virginica"))
# Returns prediction and data frame

Ответ 5

Я думаю, что то, что вы ищете, может быть пакетом lme4. Я смог получить прогноз для работы, используя это:

library(lme4)
data(iris)

model2 <- lmer(data = iris, Sepal.Length ~ (Sepal.Width | Species))
predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
       1 
6.610102 

Возможно, вам придется немного поиграть, чтобы указать конкретные эффекты, которые вы ищете, но пакет хорошо документирован, поэтому это не должно быть проблемой.