使用 R 绘制描绘 R2（Rsquared）值的 100x100 等高线图

Question

我有一个由 106 列和 28 行组成的遥感数据集。这些行与我的实例中的单个观察或单个图有关。第一列存储可以识别每个图的唯一ID。接下来的 100 列存储连续光谱带（band_x、band_x2、band_x3 等）中每个图的平均测量反射率值。其余 5 列存储在每个地块的田间测量的各种植物参数（例如叶绿素、氮、生物量等）的值。数据集或多或少如下所示：

PlotID    b1     b2    ....    b99    b100    biomass    nitrogen
1         0.11   0.16          0.40   0.41    10         52
2         0.09   0.11          0.41   0.40    19         35
3         0.10   0.19          0.43   0.49    18         72
4         0.13   0.10          0.44   0.39    16         46
...    

我希望创建等高线图，描绘与单个植物参数（例如生物量）相关的两个波段的所有可能组合的所有可能相关性的 R2（Rsquared）值。例如，等高线图需要呈现所有可能的简单比率组合（band_x1/band_x2）与单个特征之间相关性的 R2 值。此外，我希望将其复制到另外两种类型的索引中，即归一化差异索引 ((band_x2+band_x1)/(band_x2-band_x1)) 和简单差异索引 (band_x2-band_x1)。

我一直在研究 R 中的 contour.plot 语法和各种实际示例，但是，无论如何，没有一个与我所追求的有关。我以前看过这些图表，所以一定有办法生成它们。谁能帮帮我？

提前致谢！

编辑：为了澄清一些事情，这是我正在寻找重新创建的图表示例： http://image.slidesharecdn.com/2269e63a-1825-41b1-8d58-6901fd5b56ba-150102021118-conversion-gate01/95/thenkabailuavgermanyfinal1b-46-638.jpg?cb=1420186425

在 Herka 的帮助下，我现在已经根据以下代码重新创建了大部分情节（但大部分代码主要与图形相关）：

n_band=101
dat <- read.table("C:\\data.txt", header=TRUE)
res <- expand.grid(paste0("b", seq(from = 450, to = 950, by =   5)),paste0("b",seq(from = 450, to = 950, by = 5)),outcome=c("nitrogen"))

res$R2 <- apply(res, MARGIN=1,FUN=function(x){
return(cor(dat[,x[1]]/dat[,x[2]],dat[,x[3]])^2)
})

library(scales)
library(ggplot2)
p1 <- ggplot(res, aes(x=Var1, y=Var2, fill=R2)) +
geom_tile() +
facet_grid(~outcome)
p1 +
theme(axis.text.x=element_text(angle=+90)) +
geom_vline(xintercept=c(seq(from = 1, to = 101, by = 5)),color="#8C8C8C") +
geom_hline(yintercept=c(seq(from = 1, to = 101, by = 5)),color="#8C8C8C") +
labs(list(title = "Contour plot of R^2 values for all possible correlations  between Simple Ratio indices & Nitrogen Content", x = "Wavelength 1 (nm)", y = "Wavelength 2 (nm)")) +
scale_x_discrete(breaks = c("b450","b475","b500","b525","b550","b575","b600","b625","b650","b675","b700","b725","b750","b775","b800","b825","b850","b875","b900","b925","b950")) +
scale_y_discrete(breaks = c("b450","b475","b500","b525","b550","b575","b600","b625","b650","b675","b700","b725","b750","b775","b800","b825","b850","b875","b900","b925","b950")) +
scale_fill_continuous(low = "black", high = "green")

ContourPlot 我在接近我的最终目标时越来越安静，但仍有一些事情我想改变： - 有一个离散颜色的比例尺，最好依靠一个非常多样化但渐变的配色方案，以更好地识别具有最高 R2 值的波段组合。理想情况下，我希望对所有地块使用标准数量的类 (8)，每个类都包含相同数量的观察值。从而允许软件本身根据相关的每个参数的最小和最大 R2 值来确定中断值。 - 此外，我希望能够识别每个图中的最高值，或者更具体地说是它们的 (x,y) 坐标，以便我可以判断哪些波段产生最高的相关性。我使用了 which.min 和 which.max，但它们没有产生合理的结果，也没有 (x,y) 坐标。

Answer 1

这是一个示例，您可以如何解决此类问题。我已经对如何计算 R2 做了一个假设，但如果它是错误的，这很容易解决。

首先，我们模拟一些数据

set.seed(123)
n_band=100
dat <- data.frame(matrix(runif(28*n_band),ncol=n_band))
colnames(dat) <- paste0("b",1:n_band)
dat$biomass <- rpois(28,10)
dat$nitrogen <- rpois(28,10)
dat$ID <- 1:28

然后，我们观察到对于 band1、band2 和结果的每个组合，我们只需要存储一个数字 (R2)。因此，首先我们生成一个包含所有列名组合的数据框作为字符串：

res <- expand.grid(paste0("b",1:n_band),paste0("b",1:n_band),outcome=c("biomass","nitrogen"))

然后我们使用 apply 来获取每行 res 的 R2（因此每个组合）。由于 res 的每一行包含三个列名，我们可以使用它们来访问原始数据。

#ignore warnings; correlation between similar variables is missing
res$R2 <- apply(res, MARGIN=1,FUN=function(x){
  return(cor(dat[,x[1]]/dat[,x[2]],dat[,x[3]])^2)
})

那么绘图就很简单了：

library(ggplot2)
p1 <- ggplot(res, aes(x=Var1, y=Var2, fill=R2))+
  geom_tile() +
  facet_grid(~outcome)
p1