Есть ли способ построить реальное уравнение регрессии, взяв параметры из моделей в R? - PullRequest
Новая
       15

Есть ли способ построить реальное уравнение регрессии, взяв параметры из моделей в R?

2 голосов
/ 02 июня 2019

данные: 

d <- data.frame(x = rnorm(100, 0, 1),
            y = rnorm(100, 0, 1),
            z = rnorm(100, 0, 1))

функция для 5 моделей 

library(splines)
func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)
  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d) 
mods[[1]]

stargazer(mods, type="text)

Я хочу построить реальные уравнения регрессии в реальном формате каждой из моделей, взяв параметры из подгонкиМодели и переменные ind автоматически внутри R, если это возможно.Например: для модели fit1, перехват = -0,20612, х = 0,17443, х = 0,03203.Тогда уравнение будет примерно таким: y = -0.206 + 0.174x + 0.032z и т. Д., И вы хотите перечислить эти уравнения всех моделей в таблице вместе с очень распространенными полезными статистическими данными, такими как R2, значение P, прил. R2, наблюдения и т. Д. Stargazerне показывает мне желаемый результат.Поэтому я хочу убедиться, что есть ли способ сделать это в R, не делая это вручную в Excel?

Заранее спасибо!

Ответы [ 3 ]

2 голосов
/ 02 июня 2019

Мы можем map - моды , используя функцию * J.R. здесь и broom::glance для модели R2, P-value и adj.R2. 

library(purrr)
library(broom)
map_dfr(mods, 
        function(x) data.frame('Eq'=regEq(lmObj = x, dig = 3), broom::glance(x), stringsAsFactors = FALSE), 
        .id='Model') 

  Model                                                                              Eq    r.squared adj.r.squared    sigma  statistic   p.value df    logLik      AIC      BIC
1  fit1                                                   y = 0.091 - 0.022*x - 0.027*z 0.0012601436   -0.01933243 1.028408 0.06119408 0.9406769  3 -143.1721 294.3441 304.7648
2  fit2                                              y = 0.093 - 0.022*x - 0.003*I(z^2) 0.0006154188   -0.01999045 1.028740 0.02986619 0.9705843  3 -143.2043 294.4087 304.8294
3  fit3 y = 0.093 - 0.248*poly(x, 3)1 - 0.186*poly(x, 3)2 - 0.581*poly(x, 3)3 - 0.031*z 0.0048717358   -0.03702840 1.037296 0.11627016 0.9764662  5 -142.9909 297.9819 313.6129
4  fit4        y = 0.201 + 0.08*ns(x, 3)1 - 0.385*ns(x, 3)2 - 0.281*ns(x, 3)3 - 0.031*z 0.0032813558   -0.03868575 1.038125 0.07818877 0.9887911  5 -143.0708 298.1416 313.7726
  deviance df.residual
1 102.5894          97
2 102.6556          97
3 102.2184          95
4 102.3818          95
1 голос
/ 02 июня 2019

Проблема в том, что каждая из ваших моделей не совсем идеально подходит для табличных данных, например, подгонка 3 возвращает 4 оценки, в то время как подгонка 1 возвращает только 3

Если вам удобны списки, я бы предположил, что они являютсяотличный способ хранения такой информации 

library(broom)
library(tidyverse)
library(splines)

d <- data.frame(x = rnorm(100, 0, 1),
                y = rnorm(100, 0, 1),
                z = rnorm(100, 0, 1))

func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)
  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d) 

list_representation<- map(mods,tidy)
1 голос
/ 02 июня 2019

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

Функция fit2text берет подобранный объект и выводит строку символов с текстовым представлением формулы.Аргумент round дает количество цифр, до которых коэффициенты округляются в результате.Аргумент rmI, если он равен TRUE, удаляет любое I (...) и просто оставляет ... внутри, предполагая, для простоты реализации, что выражение внутри не содержит скобок.Если FALSE, то I не удаляется.

Другая статистика может быть извлечена из summary(mods[[1]]) или broom::glance(mods[[1]])

fit2text <- function(fit, round = 2, rmI = TRUE) {
  fo <- formula(fit)  
  resp <- all.vars(fo)[1]
  co <- round(coef(fit), round)
  labs <- c(if (terms(fit, "intercept") == 1) "", labels(fit))
  p <- gsub("\\+ *-", "- ", paste(resp, "~ ", paste(paste(co, labs), collapse = " + ")))
  p2 <- if (rmI) gsub("I\\(([^)]+)\\)", "\\1", p) else p
  gsub(" +", " ", p2)
}
sapply(mods, fit2text)

, давая: 

                                                           fit1 
                                  "y ~ -0.11 - 0.05 x + 0.03 z" 
                                                           fit2 
                                "y ~ -0.07 - 0.05 x - 0.04 z^2" 
                                                           fit3 
"y ~ -0.11 - 0.43 poly(x, 3) - 1.05 z + 0.27 + 0.04 poly(x, 3)" 
                                                           fit4 
    "y ~ -0.55 + 0.23 ns(x, 3) + 0.79 z - 0.25 + 0.04 ns(x, 3)"

Примечание

Код вВопрос не воспроизводился, поскольку отсутствовали вызовы библиотеки, использовались случайные числа без set.seed и в коде были некоторые дальнейшие ошибки.Для ясности мы предоставляем следующий воспроизводимый код, который мы использовали для ввода данных для ответа выше. 

library(splines)
set.seed(123)

d <- data.frame(x = rnorm(100, 0, 1),
            y = rnorm(100, 0, 1),
            z = rnorm(100, 0, 1))

# function to fit 5 models
func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)

  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d)
...