#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