Как запустить линейную регрессию по строкам и заполнить данные только элементами с NA в R?

Question

У меня есть фрейм данных, содержащий данные о плотности населения в городах, зависящей от расстояния до центра города («пространственные профили расстояния»).

Фрейм данных выглядит следующим образом (пример):

set.seed(1)
data <- data.frame(cities = c("city1","city2","city3"),
    km1 = runif(3,6,7),
    km2 = runif(3,5,6),
    km3 = runif(3,4,5),
    km4 = c(3.5,3.2,NA),
    km5 = c(NA,NA,NA)
)

«города» содержит название города или идентификатор, в то время как переменные «км1-км4» содержат журнал плотности населения на этом расстоянии.Обратите внимание, что наблюдение 3 в примере не имеет данных для км4;все города не имеют данных для км5.

То, что я пытаюсь достичь, - это экстраполировать, на сколько километров город расходится, когда плотность населения следует экспоненциальной функции.

Для этого я хочу сначала запустить линейную регрессию y ~ x для каждой строки таблицы, где y - переменные km1-kmX, а x - соответствующее расстояние до центра города (1,2,3...).

 lm(km1-kmX ~ distance) 

(Переменная «расстояние» в этом примере не определена, так как я не знаю, как включить ее в кадр данных. Но я надеюсь, что идея дойдет до конца)

Итак,для переменных city1 и city2 следует использовать km1-km4, в то время как для city3, очевидно, только km1-km3.

Полученные коэффициенты beta_0 и beta_1 должны быть затем сохранены как переменные в соответствующей строке.

Затем я хочу использовать коэффициенты для расчета логарифма плотности населения для отсутствующих переменных, вроде:

km4 = beta_0 * exp(beta_1*4) #for observation 3
km5 = beta_0 * exp(beta_1*5) #for all observations

Я знаю, что описание немного расплывчато;Я хотел быть максимально точным с усилением всех деталей идеи.Спасибо за любую помощь.

Хотя, вероятно, в этом нет необходимости, это платформа, на которой я использую R:
R версия 3.4.2 (2017-09-28)
Платформа: x86_64-w64-mingw32 / x64(64-разрядная версия)
Работает под: Windows> = 8 x64 (сборка 9200)

Answer 1

Во-первых, я думаю, что нам нужно изменить ваши данные из «широкого» в «высокий» формат.Это удовлетворит комментарий Райана о том, что вы не можете выполнять линейную регрессию с одной строкой - он технически прав, но я думаю, что он упускает из виду тот факт, что у вас на самом деле 4-5 наблюдений в строке, а не 1. ( Комментарий с момента удаления.)

(Второе: никогда не называйте переменную data. Если вы забудете создать ее в новом сеансе R, все функции, которые зависят от нее, будут сбои любопытным и часто неинтуитивным способом вместоожидаемое более простое сообщение об ошибке Error: object 'data' not found. Я буду использовать dat с вашим кодом создания.)

Это продемонстрировано на примере нескольких пакетов из tidyverse:

library(dplyr)
library(tidyr)
library(purrr)

Изменение формыВо-первых, вы перечисляете как km1, km2 и т. д., но это категориальные переменные, а не числа, и я предполагаю, что вы хотите, чтобы числа хранились в них.Так что то, что у вас есть в качестве имени столбца (km1), действительно должно быть данными (km = 1).(О, и я убираю NA, так как они не помогают кормить модель. Мы вернем их позже.)

datlong <- dat %>%
  gather(km, dens, -cities) %>%
  mutate(km = as.numeric(gsub("km", "", km))) %>%
  rename(city = cities) %>%
  filter(complete.cases(.))
datlong
#     city km     dens
# 1  city1  1 6.265509
# 2  city2  1 6.372124
# 3  city3  1 6.572853
# 4  city1  2 5.908208
# 5  city2  2 5.201682
# 6  city3  2 5.898390
# 7  city1  3 4.944675
# 8  city2  3 4.660798
# 9  city3  3 4.629114
# 10 city1  4 3.500000
# 11 city2  4 3.200000

Теперь проблема в том, как сделать регрессию для каждогогород.Во-первых, давайте немного «наведем порядок», поместив все городские данные в одну «ячейку» кадра.

datnested <- datlong %>%
  group_by(city) %>%
  nest(.key = "citydat")
datnested
# # A tibble: 3 x 2
#   city  citydat         
#   <fct> <list>          
# 1 city1 <tibble [4 x 2]>
# 2 city2 <tibble [4 x 2]>
# 3 city3 <tibble [3 x 2]>

Теперь мы можем запустить регрессию для каждого набора данных:

datmodel <- datnested %>%
  mutate(model = map(citydat, ~ lm(dens ~ km, data = .x)))
datmodel
# # A tibble: 3 x 3
#   city  citydat          model   
#   <fct> <list>           <list>  
# 1 city1 <tibble [4 x 2]> <S3: lm>
# 2 city2 <tibble [4 x 2]> <S3: lm>
# 3 city3 <tibble [3 x 2]> <S3: lm>

Заметили встроенные модели в раме?Каждый выглядит примерно так:

datmodel$model[[1]]
# Call:
# lm(formula = dens ~ km, data = .x)
# Coefficients:
# (Intercept)           km  
#       7.470       -0.926  

Теперь , что можно использовать в другом месте.Давайте запустим прогноз:

predkm <- 1:5
datpred <- datmodel %>%
  mutate(pred = map(model, ~ data_frame(km = predkm, preddens = predict(.x, newdata = data.frame(km=predkm)))))
datpred
# # A tibble: 3 x 4
#   city  citydat          model    pred            
#   <fct> <list>           <list>   <list>          
# 1 city1 <tibble [4 x 2]> <S3: lm> <tibble [5 x 2]>
# 2 city2 <tibble [4 x 2]> <S3: lm> <tibble [5 x 2]>
# 3 city3 <tibble [3 x 2]> <S3: lm> <tibble [5 x 2]>

Аналогично:

datpred$pred[[1]]
# # A tibble: 5 x 2
#      km preddens
#   <int>    <dbl>
# 1     1     6.54
# 2     2     5.62
# 3     3     4.69
# 4     4     3.77
# 5     5     2.84

Хорошо, так как мы можем получить один результирующий кадр?

datpredonly <- datpred %>%
  select(city, pred) %>%
  unnest()
datpredonly
# # A tibble: 15 x 3
#    city     km preddens
#    <fct> <int>    <dbl>
#  1 city1     1     6.54
#  2 city1     2     5.62
#  3 city1     3     4.69
#  4 city1     4     3.77
#  5 city1     5     2.84
#  6 city2     1     6.37
#  7 city2     2     5.36
#  8 city2     3     4.36
#  9 city2     4     3.35
# 10 city2     5     2.34
# 11 city3     1     6.67
# 12 city3     2     5.70
# 13 city3     3     4.73
# 14 city3     4     3.76
# 15 city3     5     2.78

Если вы хотитесравните с оригиналом (для ошибок и т. д.), попробуйте:

full_join(datlong, datpredonly, by = c("city", "km")) %>%
  arrange(city, km)
#     city km     dens preddens
# 1  city1  1 6.265509 6.543607
# 2  city1  2 5.908208 5.617601
# 3  city1  3 4.944675 4.691595
# 4  city1  4 3.500000 3.765589
# 5  city1  5       NA 2.839583
# 6  city2  1 6.372124 6.367239
# 7  city2  2 5.201682 5.361514
# 8  city2  3 4.660798 4.355788
# 9  city2  4 3.200000 3.350063
# 10 city2  5       NA 2.344337
# 11 city3  1 6.572853 6.671989
# 12 city3  2 5.898390 5.700119
# 13 city3  3 4.629114 4.728249
# 14 city3  4       NA 3.756380
# 15 city3  5       NA 2.784510

Итак, вы обсуждали использование экспоненциальной регрессии: это обрабатывается в одном вызове lm ранее в ходе выполнения.Не стесняйтесь переходить с dens ~ km на конкретные экспоненциальные формулы.

Я разбил все это на компоненты.Вот длинная цепочка.

predkm <- 1:5
datnestedmodels <- datlong %>%
  group_by(city) %>%
  nest(.key = "citydat") %>%
  mutate(
    model = map(citydat, ~ lm(dens ~ km, data = .x)),
    pred = map(model, ~ data_frame(km = predkm,
                                   preddens = predict(.x, newdata = data.frame(km=predkm))))
  )
datnestedmodels %>%
  select(city, pred) %>%
  unnest()

Если вы предпочитаете (или нуждаетесь) в «широком» формате:

datnestedmodels %>%
  select(city, pred) %>%
  unnest() %>%
  spread(km, preddens, sep = "")
# # A tibble: 3 x 6
#   city    km1   km2   km3   km4   km5
#   <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 city1  6.54  5.62  4.69  3.77  2.84
# 2 city2  6.37  5.36  4.36  3.35  2.34
# 3 city3  6.67  5.70  4.73  3.76  2.78

Answer 2

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

x <- as.numeric(gsub("\\D", "", names(data)[-1]))  # c(1, 2, 3, 4, 5)
na.lm <- function(r, x) ifelse(is.na(r), predict(lm(r ~ x), list(x = x)), r)
cbind(data[1], t(apply(data[-1], 1, na.lm, x = x)))

давая:

  cities      km1      km2      km3     km4      km5
1  city1 6.265509 5.908208 4.944675 3.50000 2.839583
2  city2 6.372124 5.201682 4.660798 3.20000 2.344337
3  city3 6.572853 5.898390 4.629114 3.75638 2.784510