我的数据集类似于:
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值上升。
DATE_FIRE
:每个地块发生火灾的年份FID
:每个图的IDDate
:NDVI的测量日期NDVI
:NDVI值[我想做的是执行2个线性回归,一个回归FID=1
,另一个回归FID=2
,以比较其恢复率。但是,我必须将恢复率ONLY应用于与发生火灾(由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)))
分组类型有效,但过滤器无效。任何帮助都将受到欢迎!
涉及dplyr
,tidyr
,lubridate
,purrr
和broom
的一个选项可以是:
data %>%
group_by(FID) %>%
filter(Date > year(DATE_FIRE)) %>%
nest(data = c(NDVI, Date)) %>%
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