动态拟合两条直线并在 R 中获取这些直线之间的交叉点

Question

我有一个像

这样的光谱反射率数据

library(hsdar)
library(tidyverse)

##Create some data
parameter <- data.frame(N = seq(1, 1.5, 0.05), LAI = seq(1,6,0.5))
spec <- PROSAIL(parameterList=parameter)

然后我计算了数据的一阶导数

d1 <- derivative.speclib(spec)

我使用以下代码从 d1 object 中提取了数据帧

d1_df <- d1@spectra@spectra_ma
d1_wav <- d1@wavelength
colnames(d1_df) <- d1_wav
#Plotting of the data
matplot(d1_wav,t(d1_df[1:11,]),type='l', xlim = c(660, 800), ylim=c(-0.01,+0.01), xlab='Wavelength /nm',ylab='Reflectance')

然后我将 far-red（680 到 700 nm）和 NIR（725 到 760 nm）区域子集化为

d1_df %>% as.data.frame() %>% 
  setNames(paste0("WV_", names(.))) %>% 
  mutate(ID = seq.int(nrow(.))) %>%
  select(281:301, ID) %>% 
  pivot_longer(cols = -ID) %>% 
  separate(name, c("chr", "wv"), convert = T)

d1_df %>% as.data.frame() %>% 
  setNames(paste0("WV_", names(.))) %>% 
  mutate(ID = seq.int(nrow(.))) %>%
  select(326:361, ID) %>% 
  pivot_longer(cols = -ID) %>% 
  separate(name, c("chr", "wv"), convert = T)

现在如何对两个区域分别拟合两条直线，得到下图中每个ID对应的两条直线交点的x？

Answer 1

这个问题没有唯一的答案，因为没有唯一的反射线（每个 ID 都有自己的反射线，因此有自己唯一的交叉点）。如果我们像这样获取您的子集数据：

region_A <- d1_df %>% as.data.frame() %>% 
  setNames(paste0("WV_", names(.))) %>% 
  mutate(ID = seq.int(nrow(.))) %>%
  select(290:301, ID) %>% 
  pivot_longer(cols = -ID) %>% 
  separate(name, c("chr", "wv"), convert = T) %>%
  mutate(ID = factor(ID))

region_B <- d1_df %>% as.data.frame() %>% 
  setNames(paste0("WV_", names(.))) %>% 
  mutate(ID = seq.int(nrow(.))) %>%
  select(332:350, ID) %>% 
  pivot_longer(cols = -ID) %>% 
  separate(name, c("chr", "wv"), convert = T) %>%
  mutate(ID = factor(ID))

绘制出来，我们看到：

p <- ggplot(region_A, aes(x = wv, y = value, group = ID)) + 
      geom_line() +
      geom_line(data = region_B) 
p

如果我们推断这些线，我们可以看到它们以不同的波长交叉：

p <- p + geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE, 
                     aes(colour = factor(ID)), se = FALSE) +
         geom_smooth(method = "lm", formula = y ~ x, fullrange = TRUE, 
                     data = region_B, aes(colour = factor(ID)), se = FALSE) +
         coord_cartesian(ylim = c(0, 0.0125))
p

我们可以像这样对每条线进行线性回归：

modA <- lm(value ~ wv * ID, data = region_A)
modB <- lm(value ~ wv * ID, data = region_B)

我们可以定义一个函数，当两个模型在特定波长下的预测相同时，returns 0，如下所示：

meet_at <- function(X, ID)
{
  A <- predict(modA, newdata = list(wv = X, ID = ID))
  B <- predict(modB, newdata = list(wv = X, ID = ID))
  abs(A - B)
}

这使我们能够使用 optimise 函数为两条线中的每一条找到交叉点和 return 一个不错的结果数据框，如下所示：

df <- do.call(rbind, lapply(unique(region_A$ID), function(i) {
  wv <- optimize(meet_at, c(700, 740), ID = i)$minimum
  value <- predict(modA, newdata = list(wv = wv, ID = i))
  data.frame(wv, value, ID = as.character(i))
}))

df
#>           wv       value ID
#> 1   708.8861 0.004254394  1
#> 11  710.4923 0.005915650  2
#> 12  712.1372 0.007343448  3
#> 13  713.6095 0.008527553  4
#> 14  714.8414 0.009483770  5
#> 15  715.8220 0.010241372  6
#> 16  716.5676 0.010833544  7
#> 17  717.1078 0.011292029  8
#> 18  717.4764 0.011644701  9
#> 19  717.7071 0.011914912 10
#> 110 717.8309 0.012121712 11

如果我们在绘图上绘制这些点，我们就知道这些结果是正确的：

p + geom_vline(data = df, aes(xintercept = wv, colour = ID)) + 
    geom_point(data = df)

所以你的问题的答案是根据 ID 在 708 和 718 nm 之间的波长范围内发生交叉，具体细节根据 df

动态拟合两条直线并在 R 中获取这些直线之间的交叉点

Dynamically fitting two straight lines and getting the crossing point between those lines in R

r

lm

tidyverse