Firma espectral y su desviación estándar utilizando R

Question

Estoy tratando de graficar una firma espectral usando R similar a la mostrada en la imagen. Esta gráfica fue obtenida utilizando el Plugin de Clasificación Semiautomática-QGIS. Las regiones semi-transparentes muestran el rango de valores, sin embargo me gustaría representar la desviación estándar.

Extraje los valores de píxel en cada banda y su clase (cobertura de tierra). Mi estructura de datos actual es un DataFrame como el siguiente:

rojo <- c(152, 28, 34, 36, 139, 144, 29, 153)
verde <- c(23, 78, 69, 68, 22, 18, 76, 26)
azul <- c(139, 18, 14, 29, 137, 162, 19, 149)
clase <- c(1, 2, 2, 2, 1, 1, 2, 1)

df <- data.frame(azul, verde, rojo, clase)
df$class <- as.factor(df$class)
names(df) <- c("430nm", "530nm", "620nm", "clase")

# mi estructura de datos
df

> df
    430   530   620 clase
1   139    23   152     1
2    18    78    28     2
3    14    69    34     2
4    29    68    36     2
5   137    22   139     1
6   162    18   144     1
7    19    76    29     2
8   149    26   153     1

Me gustaría mostrar la desviación estándar en lugar del rango y mantener la escala del eje X y la línea punteada basada en la longitud de onda de cada banda.

Answer 1

Creo que es posible que desees el Intervalo de Confianza del 95% y no la Desviación Estándar, la cual no proporciona información muy útil.

El Intervalo de Confianza es el error estándar * valor crítico para la distribución T P( t > cv ) = 0.025 = P( t < -cv )

Tus datos de ejemplo no eran lo suficientemente grandes como para demostrar efectivamente una solución estadística. Sin graficar múltiples clases, aquí tienes lo que buscas. El panel superior es el Intervalo de Confianza y el inferior es la +/- Desviación Estándar.

Esto crea un color de transparencia. Puedes obtener los valores RGB para un color con nombre usando col2rgb("blue").

p=75 # pct de transparencia
clr <- rgb(red=0, green=0, blue=255, max = 255, 
         alpha = (100 - p) * 255 / 100)

Aquí tienes algunos datos de ejemplo. Utilicé una simulación de series temporales porque produjo una distribución en serie agradable similar a una firma espectral.

f <- abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.7), n = 200)))
  f <- ((f - min(f)) * 255 / (max(f) - min(f)))

Ahora, grafica el error estándar y la desviación estándar. Estoy usando plot base pero puedes evitar tener que usar la función polygon si creas la(s) gráfica(s) usando ggplot2.

par(mfrow=c(2,1)) 

# Intervalo de confianza del 95%
se <- sd(f)/sqrt(length(f))
cv <- qt(0.95, df=(length(f)-1)) 
L = f -  cv * se
U = f + cv * se 
plot(f, type="l")
  polygon(c(1:length(f), rev(1:length(f))), c(L, rev(U)),  
          col = clr, border = FALSE)

# +/- 1 Desviación Estándar
L = f -  sd(f)
U = f + sd(f)
plot(f, type="l")
  polygon(c(1:length(f), rev(1:length(f))), c(L, rev(U)),  
          col = clr, border = FALSE)

Aquí tienes un ejemplo de múltiples clases de intervalos de confianza

d <- rbind(data.frame(class=rep(1,201), w=seq(0,3,0.01)[1:201], 
  dn=abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.9), n = 200)))), 
  data.frame(class=rep(2,201), w=seq(0,3,0.01)[1:201], 
  dn=abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.6), n = 200))))

ci <- function(x, p=0.95) qt(p, df=(length(x)-1)) * sd(x)/sqrt(length(x))   
clr <- c(rgb(0, 0, 255, max = 255, alpha = (100 - 80) * 255 / 100),
        rgb(255, 0, 0, max = 255, alpha = (100 - 80) * 255 / 100)) 

# Intervalo de confianza del 95%
plot(d[d$class == 1,]$w, d[d$class == 1,]$dn, type="n", 
     xlim=range(d$w), ylim=range(d$dn), 
     xlab="números digitales", ylab="longitud de onda")

f = d[d$class == 1,]$dn 
L = f - ci(f); U = f + ci(f)
  polygon(c(seq(0,3,0.01)[1:201], rev(seq(0,3,0.01)[1:201])),   
          c(L, rev(U)), col = clr[1], border = FALSE)
    lines(d[d$class == 1,]$w, d[d$class == 1,]$dn)    
f = d[d$class == 2,]$dn 
L = f -  ci(f); U = f + ci(f)
  polygon(c(seq(0,3,0.01)[1:201], rev(seq(0,3,0.01)[1:201])),   
          c(L, rev(U)), col = clr[2], border = FALSE)
    lines(d[d$class == 2,]$w, d[d$class == 2,]$dn)    
abline(v=c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2), col="grey")
  legend("topleft", legend=paste0("clase", c(1,2)),
         fill = clr, bg="white")