сравнить коэффициенты пошагово

Question

Я использую mtcars данные, чтобы объяснить мою проблему.Например, я пытаюсь оценить коэффициент регрессии для переменной cyl с mpg в качестве зависимой переменной и оценить изменения коэффициента путем включения других других переменных.

Шаг 1: lm(mpg ~ cyl, data = df), чтобы получить общий коэффициент для cyl

Шаг 2: добавьте каждую из других переменных по одной в модель на шаге 1, выберите одну с помощьюНаибольшее изменение (%) в коэффициенте cyl, и добавьте эту переменную в приведенной выше модели.

Шаг 3: повторите шаг 2, добавив каждую переменную из оставшихся переменных в приведенную выше модель, и снова найдите переменную с наибольшим изменением коэффициента «цил»;

Шаги ....: повторять, пока все переменные не будут включены в модель.

library(dplyr)
df <- mtcars %>% select(mpg, cyl, disp, hp, wt)

my_fun1 <- function(df=data) {
  out_df <- data.frame(matrix(ncol = 0, nrow = (length(df) - 1)))
  md_1 <- lm(mpg ~ cyl, data = df)
  out_df$Models[1] <- "Crude"
  out_df$Estimate[1] <- md_1$coefficients[2]
  pre_change <- 0
  to_rm <- 0
  for (k in 2:(length(df)-1)) {
    for (i in 3:length(df)) {
      if (!i %in% to_rm) {
        md_tmp <- update(md_1, . ~ . + df[[i]])
        change <- abs(100*(md_tmp$coefficients[2] - md_1$coefficients[2])/md_1$coefficients[2])
        dif <- md_tmp$coefficients[2] - md_1$coefficients[2]
        if (change >= pre_change) {
          out_df$Estimate[k] <- md_tmp$coefficients[2]
          out_df$Models[k] <- paste("+", names(df)[[i]])
          out_df$Diff[k] <- md_tmp$coefficients[2] - md_1$coefficients[2]
          picked <- names(df)[[i]]
          picked_i = i
          pre_change <- out_df$`Change (%)`[k] <- change
        }
      }
    }
    to_rm <- c(to_rm, picked_i)
    md_1 <- update(md_1, .~. + eval(as.name(paste(picked))))
    pre_change = 0
   }
  out_df
}

my_fun1(df = df)

После запуска выше я ожидал получить коэффициенты регрессии cyl на каждом шаге в следующем формате.

  Models  Estimate     Diff Change (%)
1  Crude -2.875790       NA         NA
2   + wt -1.507795 1.367995   47.56937
3   + hp -1.227420 0.280375   18.59504
4 + disp -1.227420 1.037274   45.80194

Однако шаги 1 и 2 дают правильные результаты, шаги 3 и 4 - неправильные.Мы ценим любые предложения.

Answer 1

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

my_fun2 <- function(dat, i) {
  fit <- lm(mpg ~ cyl, data=dat)
  crude <- fit$coef[2]
  # vectorized evaluation function
  # fits model and calculates coef and change
  evav <- Vectorize(function(i) {
    # create extension string from the "i"s
    cf.ext <- paste(names(dat)[i], collapse="+")
    # update formula with extensions
    beta <- update(fit, as.formula(
      paste0("mpg~cyl", 
             # paste already accepted coefs case they exist
             if (length(bests) != 0) {
               paste("", names(dat)[bests], sep="+", collapse="")
             },
             "+", cf.ext)
      ))$coe[2]
    # calculate Diff
    beta.d <- abs(crude - beta)
    # calculate Change %
    beta.d.perc <- 100 / crude*beta.d
    # set an attribute "cf.name" to be able to identify coef later
    return(`attr<-`(c(beta=beta, beta.d=beta.d, 
                      beta.d.perc=beta.d.perc), 
                    "cf.name", cf.ext))
  }, SIMPLIFY=FALSE)  # simplifying would strip off attributes
  # create empty vector bests
  bests <- c()
  # lapply evav() over the "i"s
  res <- lapply(i, function(...) {
    # run evav()
    i.res <- evav(i)
    # find largest change
    largest <- which.max(mapply(`[`, i.res, 2))
    # update "bests" vector within function environment with `<<-`
    bests <<- c(bests, i[largest])
    # same with the "i"s
    i <<- i[-largest]
    return(i.res[[largest]])
  })
  # summarize everything into matrix and give dimnames
  res <- `dimnames<-`(rbind(c(crude, NA, NA), 
                            do.call(rbind, res)), 
                      list(
                        c("crude", 
                          paste0("+ ", mapply(attr, res, "cf.name"))),
                        c("Estimate", "Diff", "Change (%)")))
  return(res)
}

Использование

my_fun2(mtcars[c("mpg", "cyl", "disp", "hp", "wt")], i=3:5)
#          Estimate     Diff Change (%)
# crude  -2.8757901       NA         NA
# + wt   -1.5077950 1.367995  -47.56937
# + hp   -0.9416168 1.934173  -67.25711
# + disp -1.2933197 1.582470  -55.02733

Проверка

Проверка Diff s:

fit <- lm(mpg ~ cyl, data=mtcars[c("mpg", "cyl", "disp", "hp", "wt")])
sapply(c("disp", "hp", "wt"), function(x) 
  unname(abs(fit$coe[2] - update(fit, as.formula(paste("mpg~cyl+", x)))$coe[2])))
#      disp        hp        wt 
# 1.2885133 0.6110965 1.3679952 
sapply(c("disp", "hp"), function(x) 
  unname(abs(fit$coe[2] - update(fit, as.formula(paste("mpg~cyl+wt+", x)))$coe[2])))
#     disp       hp 
# 1.090847 1.934173 
sapply(c("disp"), function(x) 
  unname(abs(fit$coe[2] - update(fit, as.formula(paste("mpg~cyl+wt+hp+", x)))$coe[2])))
#    disp 
# 1.58247 

Должно выглядеть хорошо.