R: Линейные регрессии (более 1) на одном кадре данных. Результаты в виде таблицы - PullRequest
Новая
       0

R: Линейные регрессии (более 1) на одном кадре данных. Результаты в виде таблицы

1 голос
/ 23 марта 2020

Мой набор данных похож на: 

data <- tibble( "DATE_FIRE"= c("1989-07-31", "1989-07-31", "1989-07-31", "1989-07-31","1989-07-31","1989-08-31", "1989-08-31", "1989-08-31", "1989-08-31","1989-08-31"), 
       "FID" = c(1,1,1,1,1,2,2,2,2,2),
       "Date" = c(1988, 1989, 1990, 1991, 1992, 1988, 1989, 1990, 1991, 1992),
       "NDVI" = c( 0.9, 0.8, 0.1, 0.2, 0.3, 0.8, 0.85, 0.15, 0.30, 0.50))

data$DATE_FIRE <- as.Date(data$DATE_FIRE, format= "%Y-%m-%d")
data$FID <- as.factor(data$FID)

    > data
# A tibble: 10 x 4
   DATE_FIRE  FID    Date  NDVI
   <date>     <fct> <dbl> <dbl>
 1 1989-07-31 1      1988  0.9 
 2 1989-07-31 1      1989  0.8 
 3 1989-07-31 1      1990  0.1 
 4 1989-07-31 1      1991  0.2 
 5 1989-07-31 1      1992  0.3 
 6 1989-08-31 2      1988  0.8 
 7 1989-08-31 2      1989  0.85
 8 1989-08-31 2      1990  0.15
 9 1989-08-31 2      1991  0.3 
10 1989-08-31 2      1992  0.5

Он касается лесного пожара и его восстановления по значениям NDVI. По мере восстановления леса значение NDVI возрастает.

  1. DATE_FIRE: год, когда произошел пожар для каждого участка
  2. FID: идентификатор каждого участка
  3. Date: дата измерения NDVI
  4. NDVI: значение NDVI

Я хотел бы выполнить 2 линейные регрессии, одну для FID=1 и другую для FID=2, чтобы сравните их скорость восстановления. Однако я должен применить коэффициент восстановления ТОЛЬКО к NDVI значениям, соответствующим датам после пожара (определяемого DATE_FIRE). В случае FID = 1 я должен взять только строки 3, 4 и 5, поскольку строки 1 и 2 соответствуют измерениям до пожара.

Кроме того, я хотел бы иметь свои результаты в виде таблицы; что-то вроде: 

> desired_output
# A tibble: 2 x 4
    FID  beta    r2     p
  <dbl> <dbl> <dbl> <dbl>
1     1 0.1    1     0   
2     2 0.175  0.99  0.01

ЧТО Я ПОПРОБОВАЛ ТАК ДАЛЕЕ:

Установите DATE_FIRE в годах, чтобы быть сопоставимыми с Date: 

data$DATE_FIRE <- year(data$DATE_FIRE)

Тогда: 

data_d <- data %>%
  group_by(FID) %>%
  filter(Date > DATE_FIRE) %>%
  do(tidy(lm(NDVI ~ Date,data)))

Группировка работает, но не фильтр. Любая помощь будет приветствоваться!

1 Ответ

1 голос
/ 23 марта 2020

Один вариант, включающий dplyr, tidyr, lubridate, purrr и broom, может быть: 

data %>%
 group_by(DATE_FIRE, FID) %>%
 filter(Date > year(DATE_FIRE)) %>%
 nest() %>%
 mutate(model = map(data, ~ tidy(lm(NDVI ~ Date, data = .))),
        r2 = map_dbl(data, ~ summary(lm(NDVI ~ Date, data = .))$r.squared)) %>%
 unnest(model)

  DATE_FIRE  FID             data term        estimate std.error statistic  p.value    r2
  <date>     <fct> <list<df[,2]>> <chr>          <dbl>     <dbl>     <dbl>    <dbl> <dbl>
1 1989-07-31 1            [3 × 2] (Intercept) -199.     5.85e-11  -3.40e12 1.87e-13 1    
2 1989-07-31 1            [3 × 2] Date           0.1    2.94e-14   3.41e12 1.87e-13 1    
3 1989-08-31 2            [3 × 2] (Intercept) -348.     2.87e+ 1  -1.21e 1 5.24e- 2 0.993
4 1989-08-31 2            [3 × 2] Date           0.175  1.44e- 2   1.21e 1 5.24e- 2 0.993
...