Índices de similitud

#Aquí se definen parámetros generales y el espacio de trabajo
options(max.print = 100000, width=300)                                          # dimensiones salida
getwd()                                                                         # donde estamos ahora
## [1] "/home/sergio/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_5/R"
#setwd("~/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_5 /R")                                                                # donde vamos

Esta presentación intende no solo mostrar como programar a los índices de similitud, sino también comentar sobre una programación más organizada y sencilla. Básicamente, se sugiere de construir funciones que solo hagan una cosa bien y no muchas mal.

Es útil, cuando se tienen datos dicotómicos formando una tabla grande, como por ejemplo las tablas de vegetación, donde en cada muestra se anota la presencia de todas las especies de plantas que se encuentra, de sintetizar las relaciones entre especies (pero también entre sitios) a través de una matríz de similitud, donde con cero se indica que dos especies nunca están juntas y con 1 que si siempre están.

Hay muchísimos índices, así que puede ser interesante de construir una función que calcule uno cualquier de dichos índices, solo indicando su nombre. Por esto procedemos a lo siguiente. Si se tienen dos variables dicotómicas, se construye una tabla como la siguiente, que cruza el Infarto con el uso de aspirina o de placebo (estudio longitudinal). Si ya existe, se puede volver a los datos brutos como sigue. Como ejemplo, sacamos los datos de Ellenberg, que tienen la presencia/ausencia de 76 plantas en 25 sitios de campos europeos.

data                    <- read.csv("Ellenberg.csv",row.names = 1)              #  datos de vegetación 
head(data)                                                                      #  de Ellenberg
##           X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21 X22 X23 X24 X25
## Arre_elat  0  1  0  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
## Dact_glom  1  1  1  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
## Heli_pube  0  1  0  0  1  1  1  0  1   1   0   1   1   0   0   1   1   1   0   1   1   1   0   1   0
## Brom_erec  1  0  1  1  0  0  0  0  1   1   0   0   0   0   1   0   0   0   0   0   0   0   0   1   0
## Fest_ovin  1  1  0  0  0  0  0  0  1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
## Poa_prat   1  1  1  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
nr                      <- dim(data)[1]; nr                                     # número de filas/especies
## [1] 76
nc                      <- dim(data)[2]; nc                                     # número de columnas/sitios
## [1] 25
namrows                 <- rownames(data)                                       # nombres de filas
namcols                 <- colnames(data)                                       # nombres de columnas

Se define una matríz de associación entre filas

ass                     <- matrix(0,nrow=nr,ncol=nr, dimnames = list(namrows,namrows))  # vacía
start_time <- Sys.time()
for (i in (1:nr)) {                                                             # dos iteraciones 
  for (j in (1:nr)) {                                                           # para llenar a la matriz
    a                   <- 0
    b                   <- 0
    c                   <- 0
    d                   <- 0
    for (k in 1:nc) {                                                           # iteración sobre columnas 
      if (data[i,k] == 1){
        if (data[j,k] == 1){
          a             <- a+1                                                  # a  (p/p) 
        } else {
          b             <- b+1                                                  # b  (p/a) 
        }  
      } else {
        if (data[j,k] == 1){
          c             <- c+1                                                  # c (a/p)
        } else {
          d             <- d+1                                                  # d (a/a)
        }  
      }
    }
    sor                 <- 2*a/(2*a+b+c)                                        # sorensen
    ass[i,j]            <- sor                                                  # valores en ass
  }   
}
end_time <- Sys.time()
round(ass[1:5,1:5],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin
## Arre_elat      1.00      0.96      0.79      0.33      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32
## Brom_erec      0.33      0.44      0.27      1.00      0.55
## Fest_ovin      0.22      0.28      0.32      0.55      1.00
end_time-start_time
## Time difference of 3.55771 secs

De hecho, este proceso se puede simplificar, copiando mitad matríz sobre otra mitad.

ass2                    <- matrix(0,nrow=nr,ncol=nr, dimnames = list(namrows,namrows)) # vacía
start_time <- Sys.time()
for (i in (1:nr)) {                                                             # dos iteraciones 
  for (j in (1:i)) {                                                            # para llenar mitad matriz
    a                   <- 0
    b                   <- 0
    c                   <- 0
    d                   <- 0
   for (k in 1:nc) {                                                           # iteración sobre columnas 
      if (data[i,k] == 1){
        if (data[j,k] == 1){
          a             <- a+1                                                  # a  (p/p) 
        } else {
          b             <- b+1                                                  # b  (p/a) 
        }  
      } else {
        if (data[j,k] == 1){
          c             <- c+1                                                  # c (a/p)
        } else {
          d             <- d+1                                                  # d (a/a)
        }  
      }
    }
    sor                 <- 2*a/(2*a+b+c)                                        # sorensen
    ass2[i,j]           <- sor                                                  # valores en ass
    if (i!=j) ass2[j,i] <- sor                                                  # también otra mitad
  }                       
}
end_time <- Sys.time()
round(ass2[1:5,1:5],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin
## Arre_elat      1.00      0.96      0.79      0.33      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32
## Brom_erec      0.33      0.44      0.27      1.00      0.55
## Fest_ovin      0.22      0.28      0.32      0.55      1.00
end_time-start_time
## Time difference of 1.723321 secs

El ideal sería de derivar a,b,c,d de una tabla de contingencia, cruzando dos variables, o sea con el comando table(data[i,],data[j,]). Infelizmente esto no funciona, porque si una de las columnas tiene todos ceros o todos unos, la tabla resulta con solo un fila o una columna. Una programación más eficiente podría ser

  ass3                   <- matrix(0,nrow=nr,ncol=nr, dimnames = list(namrows,namrows)) # vacía
start_time <- Sys.time()
  for (i in (1:nr)) {                                                           # dos iteraciones 
    for (j in (1:i)) {                                                          # para llenar mitad matriz
      a                   <- length(which(data[i,]==1&data[j,]==1))             # a (p/p)
      b                   <- length(which(data[i,]==1&data[j,]==0))             # b (p/a)
      c                   <- length(which(data[i,]==0&data[j,]==1))             # c (a/p)
      d                   <- length(which(data[i,]==0&data[j,]==0))             # d (a/a)
      sor                 <- 2*a/(2*a+b+c)                                      # sorensen
      ass3[i,j]           <- sor                                                # valores en ass
      if (i!=j) ass3[j,i] <- sor                                                # también otra mitad
  }                       
}
end_time <- Sys.time()
round(ass3[1:5,1:5],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin
## Arre_elat      1.00      0.96      0.79      0.33      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32
## Brom_erec      0.33      0.44      0.27      1.00      0.55
## Fest_ovin      0.22      0.28      0.32      0.55      1.00
end_time-start_time
## Time difference of 10.32377 secs

Aún este tercero algoritmo es más compacto, sin embargo toma mucho más tiempo. Por tanto, podemos construir una función con el segundo.

sorensen                <- function(data) {
  nr                    <- dim(data)[1]
  nc                    <- dim(data)[2]
  namrows               <- rownames(data)
  ass                   <- matrix(0,nrow=nr,ncol=nr, dimnames = list(namrows,namrows)) # vacía
  for (i in (1:nr)) {                                                           # dos iteraciones 
    for (j in (1:i)) {                                                          # para llenar mitad matriz
      a                   <- 0
      b                   <- 0
      c                   <- 0
      d                   <- 0
      for (k in 1:nc) {                                                         # iteración sobre columnas
        if (data[i,k] == 1){
          if (data[j,k] == 1){
            a             <- a+1                                                # a  (p/p) 
          } else {
            b             <- b+1                                                # b  (p/a) 
          }  
        } else {
          if (data[j,k] == 1){
            c             <- c+1                                                # c (a/p)
          } else {
            d             <- d+1                                                # d (a/a)
          }  
        }
      }
      sor                 <- 2*a/(2*a+b+c)                                      # sorensen
      ass[i,j]            <- sor                                                # valores en ass
      if (i!=j) ass [j,i] <- sor                                                # también otra mitad
    }
  }                       
  ass
}

así que el resultado será

round(sorensen(data)[1:8,1:8],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin Poa_prat Briz_medi Koel_pyra
## Arre_elat      1.00      0.96      0.79      0.33      0.22     0.96      0.29      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32     0.75      0.40      0.21
## Brom_erec      0.33      0.44      0.27      1.00      0.55     0.44      0.50      0.73
## Fest_ovin      0.22      0.28      0.32      0.55      1.00     0.28      0.89      0.75
## Poa_prat       0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Briz_medi      0.29      0.33      0.40      0.50      0.89     0.33      1.00      0.67
## Koel_pyra      0.22      0.28      0.21      0.73      0.75     0.28      0.67      1.00

Nótese que el índice se puede calcular para las columnas solo transponendo a la matriz:

datb = t(data)
round(sorensen(datb)[1:8,1:8],2)
##      X1   X2   X3   X4   X5   X6   X7   X8
## X1 1.00 0.58 0.63 0.80 0.58 0.44 0.42 0.47
## X2 0.58 1.00 0.58 0.49 0.60 0.53 0.57 0.45
## X3 0.63 0.58 1.00 0.62 0.54 0.44 0.48 0.43
## X4 0.80 0.49 0.62 1.00 0.57 0.40 0.44 0.49
## X5 0.58 0.60 0.54 0.57 1.00 0.75 0.69 0.71
## X6 0.44 0.53 0.44 0.40 0.75 1.00 0.72 0.74
## X7 0.42 0.57 0.48 0.44 0.69 0.72 1.00 0.65
## X8 0.47 0.45 0.43 0.49 0.71 0.74 0.65 1.00

Si se quiere cambiar índice, se necesita de reescribir la función, solo cambiando la linia con el cálculo de la función. Sino, conviene la siguiente función:

similitud               <- function(data,simil=NULL) {
  nr                    <- dim(data)[1]
  nc                    <- dim(data)[2]
  namrows               <- rownames(data)
  ass                   <- matrix(0,nrow=nr,ncol=nr, dimnames = list(namrows,namrows)) # vacía
  for (i in (1:nr)) {                                                           # dos iteraciones 
    for (j in (1:i)) {                                                          # para llenar mitad matriz
      a                   <- 0
      b                   <- 0
      c                   <- 0
      d                   <- 0
      for (k in 1:nc) {                                                         # iteración sobre columnas
        if (data[i,k] == 1){
          if (data[j,k] == 1){
            a             <- a+1                                                # a  (p/p) 
          } else {
            b             <- b+1                                                # b  (p/a) 
          }  
        } else {
          if (data[j,k] == 1){
            c             <- c+1                                                # c (a/p)
          } else {
            d             <- d+1                                                # d (a/a)
          }  
        }
      }
      if        (simil=="sorensen") {
            sim           <- 2*a/(2*a+b+c)                                      # Sorensen
      } else if (simil=="vandermaarel") {
            sim         <- (2*a-(b+c))/(2*a+(b+c))                            # Van der Maarel
      }      
      ass[i,j]            <- sim                                                # valores en ass
      if (i!=j) ass [j,i] <- sim                                                # también otra mitad
    }
  }                       
  round(ass[1:8,1:8],2)
}

Como resultado se consigue

round(similitud(data,simil="sorensen")[1:8,1:8],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin Poa_prat Briz_medi Koel_pyra
## Arre_elat      1.00      0.96      0.79      0.33      0.22     0.96      0.29      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32     0.75      0.40      0.21
## Brom_erec      0.33      0.44      0.27      1.00      0.55     0.44      0.50      0.73
## Fest_ovin      0.22      0.28      0.32      0.55      1.00     0.28      0.89      0.75
## Poa_prat       0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Briz_medi      0.29      0.33      0.40      0.50      0.89     0.33      1.00      0.67
## Koel_pyra      0.22      0.28      0.21      0.73      0.75     0.28      0.67      1.00
round(similitud(data,simil="vandermaarel")[1:8,1:8],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin Poa_prat Briz_medi Koel_pyra
## Arre_elat      1.00      0.92      0.58     -0.33     -0.56     0.92     -0.43     -0.56
## Dact_glom      0.92      1.00      0.50     -0.12     -0.45     1.00     -0.33     -0.45
## Heli_pube      0.58      0.50      1.00     -0.45     -0.37     0.50     -0.20     -0.58
## Brom_erec     -0.33     -0.12     -0.45      1.00      0.09    -0.12      0.00      0.45
## Fest_ovin     -0.56     -0.45     -0.37      0.09      1.00    -0.45      0.78      0.50
## Poa_prat       0.92      1.00      0.50     -0.12     -0.45     1.00     -0.33     -0.45
## Briz_medi     -0.43     -0.33     -0.20      0.00      0.78    -0.33      1.00      0.33
## Koel_pyra     -0.56     -0.45     -0.58      0.45      0.50    -0.45      0.33      1.00

Aún mejor es de compartir la función similitud en una que calcula a,b,c,d y una que calcula el índice que se quiere. De esta manera, para ajuntar índices, solo hay que cambiar esta última.

La función matsim() calcula la matriz de similitud, según el tipo de similitud eligido. El cálculo se lo hace por columnas, a menos que no se indique rows=TRUE. El cálculo de la similitud lo hace la función similitud, que se puede programar con todos los índices que se quieren, en base a los valores a,b,c,d de la tabla de contingencia construida.

matsim                  <- function(data,simil=NULL,rows=FALSE) {
  if(rows) {
    data                <- t(as.matrix(data))                                   # si por filas, transpone 
  }
  nr                    <- dim(data)[1]                                         # número filas
  nc                    <- dim(data)[2]                                         # número columnas
  namcols               <- colnames(data)                                       # nombres columnas
  ass                   <- matrix(0,nrow=nc,ncol=nc, dimnames = list(namcols,namcols)) # vacía
  for (i in (1:nc)) {                                                           # dos iteraciones 
    for (j in (1:i)) {                                                          # para llenar mitad matriz
      a                   <- 0
      b                   <- 0
      c                   <- 0
      d                   <- 0
      for (k in 1:nr) {                                                         # iteración sobre filas
        if (data[k,i] == 1){
          if (data[k,j] == 1){
            a             <- a+1                                                # a  (p/p) 
          } else {
            b             <- b+1                                                # b  (p/a) 
          }  
        } else {
          if (data[k,j] == 1){
            c             <- c+1                                                # c (a/p)
          } else {
            d             <- d+1                                                # d (a/a)
          }  
        }
      }
      ass[i,j]            <- similitud(a,b,c,d,simil)                           # valores en ass
      if (i!=j) ass [j,i] <- ass[i,j]                                           # también otra mitad
    }
  }                       
  ass
}

La función similitud cálcula las similitudes eligidas.

similitud                 <- function(a,b,c,d,simil) {
  if        (simil=="sorensen") {
    sim                   <- 2*a/(2*a+b+c)                                      # Sorensen
  } else if (simil=="vandermaarel") {
    sim         <- (2*a-(b+c))/(2*a+(b+c))                                      # Van der Maarel
  }  
  sim
}

Los resultados son los siguientes:

round(matsim(data,simil="sorensen",row=TRUE)[1:8,1:8],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin Poa_prat Briz_medi Koel_pyra
## Arre_elat      1.00      0.96      0.79      0.33      0.22     0.96      0.29      0.22
## Dact_glom      0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Heli_pube      0.79      0.75      1.00      0.27      0.32     0.75      0.40      0.21
## Brom_erec      0.33      0.44      0.27      1.00      0.55     0.44      0.50      0.73
## Fest_ovin      0.22      0.28      0.32      0.55      1.00     0.28      0.89      0.75
## Poa_prat       0.96      1.00      0.75      0.44      0.28     1.00      0.33      0.28
## Briz_medi      0.29      0.33      0.40      0.50      0.89     0.33      1.00      0.67
## Koel_pyra      0.22      0.28      0.21      0.73      0.75     0.28      0.67      1.00
round(matsim(data,simil="sorensen")[1:8,1:8],2)
##      X1   X2   X3   X4   X5   X6   X7   X8
## X1 1.00 0.58 0.63 0.80 0.58 0.44 0.42 0.47
## X2 0.58 1.00 0.58 0.49 0.60 0.53 0.57 0.45
## X3 0.63 0.58 1.00 0.62 0.54 0.44 0.48 0.43
## X4 0.80 0.49 0.62 1.00 0.57 0.40 0.44 0.49
## X5 0.58 0.60 0.54 0.57 1.00 0.75 0.69 0.71
## X6 0.44 0.53 0.44 0.40 0.75 1.00 0.72 0.74
## X7 0.42 0.57 0.48 0.44 0.69 0.72 1.00 0.65
## X8 0.47 0.45 0.43 0.49 0.71 0.74 0.65 1.00
round(matsim(data,simil="vandermaarel",row=TRUE)[1:8,1:8],2)
##           Arre_elat Dact_glom Heli_pube Brom_erec Fest_ovin Poa_prat Briz_medi Koel_pyra
## Arre_elat      1.00      0.92      0.58     -0.33     -0.56     0.92     -0.43     -0.56
## Dact_glom      0.92      1.00      0.50     -0.12     -0.45     1.00     -0.33     -0.45
## Heli_pube      0.58      0.50      1.00     -0.45     -0.37     0.50     -0.20     -0.58
## Brom_erec     -0.33     -0.12     -0.45      1.00      0.09    -0.12      0.00      0.45
## Fest_ovin     -0.56     -0.45     -0.37      0.09      1.00    -0.45      0.78      0.50
## Poa_prat       0.92      1.00      0.50     -0.12     -0.45     1.00     -0.33     -0.45
## Briz_medi     -0.43     -0.33     -0.20      0.00      0.78    -0.33      1.00      0.33
## Koel_pyra     -0.56     -0.45     -0.58      0.45      0.50    -0.45      0.33      1.00