Neurociencia y Matemáticas: noviembre 2015

Introducción

Tras muchas semanas de retraso debido a una serie de contratiempos ¡Por fin tenemos la segunda parte de este tema!
He encontrado varios problemas para publicar este artículo:

La gripe, tres días en cama no ayudan a escribir esto.
Github desktop no está funcionando en mi laptop. Esto es un problema porque quería intentar hacer un repositorio útil para este artículo, en cambio sigo usando Github de manera pedestre. Bueno, poco a poco.
Mi obsesión por que los script se puedan ejecutar copiando y pegando. Odio cuando alguien cuelga un script que no funciona. (Igual dentro de unos meses este script deja de funcionar. Si eso pasase avísadme)

La parte buena de mi retraso es que este blog iba a ser dividido en varias partes, por lo que este post será especialmente largo para compensar la sequía.

Una vez más me gustaría agradecer a Álex Estudillo por los datos proporcionados (medias y desviaciones) y las imágenes.

Los datos los he generado basándome en unos datos de Alex tomando medias y desviaciones. Luego he simulado la matriz por sujetos con distribuciones normales.

La imagen sobre la que trabajaremos es la siguiente. Que por cierto, es propiedad intelectual de Álex, así que está aquí solo con fines de reproducción del ejemplo. Aun así por si acaso he modificado la imagen con GIMP.

Tratamiento de datos y descarga

Lo primero será leer los datos que es lo que hacemos en este script:

	needed <- list("dplyr", "data.table", # Libraries for data mining
	"jpeg", "ggplot2", "png", "grid", "geoR", "rgl", # Libraries for
	"repmis" # Libraries for web scrapping [3]
	)

	if(FALSE){ # If you want to install change for TRUE
	lapply(X = needed ,FUN =
	function(x){
	if(!require(x, character.only = T)){
	install.packages(x)
	}
	library(x,character.only = T)
	})
	}else{
	lapply(X = needed ,FUN =
	function(x){
	library(x,character.only = T)
	})
	}

	dir <- # "C:/Users/Javier Villacampa/Dropbox/Chabi/Blog/Eye track/"
	"C:/Users/usuario/Dropbox/Chabi/Blog/Eye track/" # Put your dir
	setwd(dir)
	dir.create(path = "Results")
	dir.create(path = "Data")
	rm(list = ls()); gc()


	# S read data -------------------------------------------------------------
	# Download the image
	url <- "http://2.bp.blogspot.com/-Qo8JF_ux808/VlEAanXvltI/AAAAAAAACds/TkzNkifUK4M/s1600/Imagen1.jpeg"
	download.file(url = url, destfile = "Data/Imagen1.jpeg", mode = "wb")


	image_1 <- readJPEG( source = "Data/Imagen1.jpeg")
	image_1 <- rasterGrob(image = image_1, interpolate=TRUE)
	# Download the csv from dropbox
	d <- repmis::source_DropboxData(file = "sample-1.txt", key = "1s4uo7u340snqtv", sep ="\t") %>% data.table()
	# E read data -------------------------------------------------------------



	# S Clean data ------------------------------------------------------------
	d[ , TRIAL_INDEX := factor(TRIAL_INDEX)]
	d[ , Subject := factor(Subject)]
	d <- d[ d$CURRENT_FIX_X <= 1000 & d$CURRENT_FIX_X >= 0 &
	d$CURRENT_FIX_Y <= 700 & d$CURRENT_FIX_Y >= 0, ]
	# S Clean data ------------------------------------------------------------


	# S Plot ------------------------------------------------------------------
	d <- d %>% data.frame()

	d2 <- copy(d) # Save a copy from the old data
	d <- d2[, c("CURRENT_FIX_X", "CURRENT_FIX_Y")]
	# E Plot ------------------------------------------------------------------

view raw Blog005 01 read data.R hosted with ❤ by GitHub

Solución 1

La primera solución es muy parecida a lo que hice en el post anterior, lo único que aporta
coord_fixed() es que nos fija las imágenes para que no tengamos problemas con el reescalamiento.
Las órdenes de dibujo que definen el objeto ggplot vienen de la librería ggplot2.

	# Solution 1
	# S plot density ----------------------------------------------------------
	p <-
	ggplot(d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y)) +
	annotation_custom(grob = image_1, # Image
	xmin= 0, # Coordinates to represent the image
	xmax= 1000,
	ymin= 0,
	ymax= 800
	# The image has 700 pixel. But the matrix has the median in 400 and if a put the image with 700px the eyes are in 350 (coord y).
	# The I changed the size to 800 px in order to adjust the eyes position to the data position
	) +
	coord_fixed(xlim = c(0, 800), ylim = c(0, 700)) + # fix the image in the coordinates that we deffined before
	geom_point(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.5) + # Point transaparecen
	stat_density2d(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y , fill = ..level.., alpha = ..level..),
	size= 100, bins= 50, geom='polygon') +
	theme_bw() + # Kind of theme. I strongly recomend theme_bw
	scale_fill_gradient( low = "green", # Lowest color value
	high = "red" # High color value
	) +
	scale_alpha_continuous(range=c(0.0, 1) , guide = FALSE) +
	# You can play with the range to show a better image. Range belongs to [0, 1] interval
	xlim(0, 1000) + # Control lim for x-axe
	ylim(0, 700) # Control lim for y-axe
	print(p)
	jpeg(filename = "Results/Solution1.jpeg", width = 1000, height = 700, units = "px")
	print(p)
	dev.off()
	# E plot density ----------------------------------------------------------

view raw 02 Solution1.R hosted with ❤ by GitHub

Solución 2

Aquí solo aporto una función para controlar los colores del heatmap, función que será utilizada durante todo el resto del ejercicio al quedar los colores bastante más aparentes. Ver lineas: 3 y 23.
colorRamp viene de la librería grDevices

	# Solution 2
	# S plot density with Palette ---------------------------------------------
	colfunc <- colorRampPalette(c("darkblue", "lightblue", "green", "yellow", "red")) # To create scale for representing the heatmap
	p <-
	ggplot(d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y)) +
	annotation_custom(grob = image_1, # Image
	xmin= 0, # Coordinates to represent the image
	xmax= 1000,
	ymin= 0,
	ymax= 800
	# The image has 700 pixel. But the matrix has the median in 400 and if a put the image with 700px the eyes are in 350 (coord y).
	# The I changed the size to 800 px in order to adjust the eyes position to the data position
	) +
	coord_fixed(xlim = c(0, 1000), ylim = c(0, 700)) + # fix the image in the coordinates that we deffined before
	geom_point(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.5) + # Point transaparecen

	stat_density2d(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y , fill = ..level.., alpha = ..level..),
	size= 100, bins= 50, geom='polygon') +
	theme_bw() + # Kind of theme. I strongly recomend theme_bw
	scale_fill_gradientn(colours=colfunc(100)) + # Create a color for filling the heatmap
	guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10)) +
	theme(legend.title=element_blank()) +
	scale_alpha_continuous(range=c(0.0, 1) , guide = FALSE) + # You can play with the range to show a better image. Range belongs to [0, 1] interval
	xlim(0, 1000) + # Control lim for x-axe
	ylim(0, 700) # Control lim for y-axe
	print(p)
	jpeg(filename = "Results/Solution2.jpeg", width = 1000, height = 700, units = "px")
	print(p)
	dev.off()
	# E plot density with Palette ---------------------------------------------

view raw 03 Solution2.R hosted with ❤ by GitHub

Solución 3

Los heatmaps con distribuciones están muy bien de cara a representar, pero lo que buscamos es mas bien un histograma sobre el plano que nos cuente el número de ocurrencias por punto. Para este fin podemos definir zonas de interés mediante hexágonos y de esta manera generamos el heatmap-histograma. Lo malo de esta solución es que es difícil generar los hexágonos de manera que queden exactamente en las ROIs que queremos y por tanto es una solución poco útil de cara al análisis.

stat_binhex nos permite crear tantos hexágonos como queramos para definir ROIs. El parámetro de bins nos da el número de hexágonos que caben en el eje x.

	# Solution 3
	# S wit hexagons ----------------------------------------------------------
	p <-
	ggplot(d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y)) +
	annotation_custom(grob = image_1, # Image
	xmin= 0, # Coordinates to represent the image
	xmax= 1000,
	ymin= 0,
	ymax= 800
	) +
	coord_fixed(xlim = c(0, 1000), ylim = c(0, 700)) +
	geom_point(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.5) + # Point transaparecen

	stat_binhex(bins = 50, colour = "gray", alpha = 0.5) +

	theme_bw() + # Kind of theme. I strongly recomend theme_bw
	scale_fill_gradientn(colours=colfunc(400)) +
	guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10)) +
	theme(legend.title=element_blank()) +
	scale_alpha_continuous(range=c(0.0, 1) , guide = FALSE) + # You can play with the range to show a better image. Range belongs to [0, 1] interval
	xlim(0, 1000) + # Control lim for x-axe
	ylim(0, 700) # Control lim for y-axe
	print(p)
	jpeg(filename = "Results/Solution3.jpeg", width = 1000, height = 700, units = "px")
	print(p)
	dev.off()
	# E wit hexagons ----------------------------------------------------------

view raw 04 Solution3.R hosted with ❤ by GitHub

Como veis el resultado es muy similar, pero ahora ya hemos definido algo aproximado a una zona de interés.

Solución 4

Ahora vamos a definir unas zonas de interés con cuadrados. Es la misma lógica que tienen los hexágonos sólo que los cuadrados tendrán mucho más juego al poderse elegir más fácilmente donde los colocamos.

stat_bin2d La orden utilizada será ésta y tiene un uso similar a la que tiene la orden stat_binhex, aunque como veremos en el último apartado es mucho más flexible.

	# Solution 4
	# S with boxes ------------------------------------------------------------
	p <-
	ggplot(d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y)) +
	annotation_custom(grob = image_1, # Image
	xmin= 0, # Coordinates to represent the image
	xmax= 1000,
	ymin= 0,
	ymax= 800
	) +
	coord_fixed(xlim = c(0, 1000), ylim = c(0, 700)) +
	geom_point(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.5) + # Point transaparecen

	stat_bin2d(bins = 10, colour = "gray", alpha = 0.5) +

	theme_bw() + # Kind of theme. I strongly recomend theme_bw
	scale_fill_gradientn(colours=colfunc(400)) +
	guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10)) +
	theme(legend.title=element_blank()) +
	scale_alpha_continuous(range=c(0.0, 1) , guide = FALSE) + # You can play with the range to show a better image. Range belongs to [0, 1] interval
	xlim(0, 1000) + # Control lim for x-axe
	ylim(0, 700) # Control lim for y-axe
	print(p)
	jpeg(filename = "Results/Solution4.jpeg", width = 1000, height = 700, units = "px")
	print(p)
	dev.off()
	# E with boxes ------------------------------------------------------------

view raw 05 Solucion5.R hosted with ❤ by GitHub

Solución 5

Esta solución consiste en crear la matriz de frecuencias a mano y representárla con la orden surface3d del paquete rgl. Lo bueno de esta solución es que ya tenemos la matriz de frecuencias lista para el análisis y que además nos da un dibujo en 3D que se puede mover. Esta solución es bonita pero poco útil.

	# Solution 5
	# S Surface 3d ------------------------------------------------------------
	nbins <- 20
	x.bin <- seq(0, 1000, length= nbins)
	y.bin <- seq(0, 700, length= nbins)

	freq <- as.data.frame(table(findInterval(x = d[,1], vec = x.bin),
	findInterval(x = d[,2], vec = y.bin))) # generates a frequency matrix
	freq[,1] <- as.numeric(freq[,1])
	freq[,2] <- as.numeric(freq[,2])

	freq2D <- diag(nbins)*0 # create a square matrix with nbin rows.
	freq2D[cbind(freq[,1], freq[,2])] <- freq[,3] # Generates frequence mastrix

	surface3d(x.bin,y.bin,freq2D/5, col="steelblue")
	# E Surface 3d ------------------------------------------------------------

view raw 06 Solution 5.R hosted with ❤ by GitHub

Solución 6

Basándonos en la matriz de frecuencias de la solución 5 podemos intentar hacer otro tipo de dibujo. Sólo que ahora aproximaremos por polinomios a una forma que se ajuste a los valores dados. El heatnap queda aparente y nos da una visión más general de a dónde han ido las miradas, pero no nos permite hacer un análisis lo suficientemente bueno.

	# Solution 7
	# S plot heatmap ----------------------------------------------------------
	x <- interaction(x.bin, y.bin, sep = "_") %>% levels %>% strsplit(split = "_") %>%
	unlist %>% as.numeric() %>% matrix(ncol= 2, byrow = T) %>% data.frame
	names( x) <- c("x.bin" , "y.bin" )
	elevation.df <- data.frame(freq2D = as.vector(freq2D), x.bin = x$x.bin, y.bin = x$y.bin )


	elevation.loess <- loess(freq2D ~ x.bin*y.bin, data = elevation.df, degree = 2, span = 0.25) # Predict the gradients
	elevation.fit <- expand.grid(list(x.bin = seq(0, 1000, 10), y.bin = seq(0, 700, 10))) # Fit the values in the curve
	z <- predict(elevation.loess, newdata = elevation.fit)
	elevation.fit$Height <- as.numeric(z)
	elevation.fit$Height[ is.na( elevation.fit$Height )] <- 0

	colfunc <- colorRampPalette(c("darkblue", "lightblue", "green", "yellow", "red"))
	p <- ggplot(elevation.fit, aes(x.bin, y.bin, fill = Height)) + geom_tile() +
	geom_point(data= elevation.fit, aes(x = x.bin, y =y.bin), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.05) +
	xlab("X Coordinate") + ylab("Y Coordinate") +
	scale_fill_gradientn(colours=colfunc(100)) +
	scale_x_continuous(expand = c(0,0)) +
	scale_y_continuous(expand = c(0,0))
	print(p)
	jpeg(filename = "Results/Solution6.jpeg", width = 1000, height = 700, units = "px")
	print(p)
	dev.off()
	# E plot heatmap ----------------------------------------------------------

view raw 07 Solucion6.R hosted with ❤ by GitHub

Analisis final

Ahora vamos a utilizar la solución 4, que es la que más me gustó, para terminar el análisis. Pero vamos a elegir las ROI de manera manual y no de manera automática. Para ello definiremos los cuadrados con unos vectores para generar las distintos ROIs.
stat_bin2d(breaks = list( x= seq(from = 0, to = 1000, by = 50),
y = seq(from = 0, to = 800, by = 50) )

	# S First deffine de plots ------------------------------------------------
	p <-
	ggplot(d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y)) +
	annotation_custom(grob = image_1, # Image
	xmin= 0, # Coordinates to represent the image
	xmax= 1000,
	ymin= 0,
	ymax= 800
	) +
	coord_fixed(xlim = c(0, 1000), ylim = c(0, 700)) +
	geom_point(data= d, aes(x = CURRENT_FIX_X, y =CURRENT_FIX_Y), # coordinates
	color="black", # Color point
	# position=position_jitter(w=0.01,h=0.01), # Point plot desviation
	alpha=0.5) + # Point transaparecen
	stat_bin2d(breaks = list( x= seq(from = 0, to = 1000, by = 50),
	y = seq(from = 0, to = 800, by = 50)) ,
	colour = "gray", alpha = 0.5) +
	# stat_bin2d(bins = 10, colour = "gray", alpha = 0.5) +
	theme_bw() + # Kind of theme. I strongly recomend theme_bw
	scale_fill_gradientn(colours=colfunc(10)) +
	guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10)) +
	theme(legend.title=element_blank()) +
	scale_alpha_continuous(range=c(0.0, 1) , guide = FALSE) + # You can play with the range to show a better image. Range belongs to [0, 1] interval
	xlim(0, 1000) + # Control lim for x-axe
	ylim(0, 700) # Control lim for y-axe
	print(p)
	# E First deffine de plots ------------------------------------------------

view raw 08.1 Create plot for the analysis hosted with ❤ by GitHub

Por último analizaremos las frecuencias con un test de proporciones. Los pasos son los siguientes:

Lo primero con ggplot_build hemos sacado la estructura de la matriz que genera el plot. Mucho más sencillo que como lo hemos hecho en las soluciones 5 y 6 (línea4)

Luego nos hemos fijado en los cuadrantes donde los sujetos se fijaron más de 500 veces.

Se observa el cuadrante (250,300] (350,400] Ojo izquierdo de la cara izquierda (línea 10)
También podemos ver el cuadrante (700,750] (350,400] que es el ojo derecho de la cara de la derecha. Ambos cuadrantes son los más calientes.

Por último vamos a comparar si los sujetos en general se han fijado más en un cuadrante o en el otro. Para ello hacemos un test de proporciones en el que nos da un p valor menor que 0.05 (línea 24), con lo que podemos concluir que efectivamente se han fijado más en un cuadrante que en otro.

	# S Analysis --------------------------------------------------------------
	plotData <- ggplot_build(plot = p)
	plotData <- plotData$data
	plotData <- plotData[[3]] # Extract the frecuency matrix
	head(plotData)
	summary(plotData)
	plotData[ plotData$count > 500,]
	# fill xbin ybin count density ymax ymin yint xmax xmin xint PANEL group
	# 69 #A0FF00 (200,250] (350,400] 617 0.05875071 400 350 8 250 200 5 1 1
	# 70 #E5ED00 (250,300] (350,400] 697 0.06636831 400 350 8 300 250 6 1 1
	# 79 #FF0000 (700,750] (350,400] 942 0.08969720 400 350 8 750 700 15 1 1
	# 80 #50FF00 (750,800] (350,400] 537 0.05113312 400 350 8 800 750 16 1 1
	# 88 #C5F800 (250,300] (400,450] 658 0.06265473 450 400 9 300 250 6 1 1
	plotDataDensity <- data.frame(Square = paste0(plotData$xbin ,"-" , plotData$ybin), Count = plotData$count)
	plotDataDensity$Density <- plotDataDensity$Count / sum(plotDataDensity$Count) *100
	head(plotDataDensity)
	summary(plotDataDensity)
	plotDataDensity[ plotDataDensity$Density > 1,]
	prop.test(x = c(plotDataDensity$Count[70], plotDataDensity$Count[79]),
	n = rep(sum(plotDataDensity$Count),2) )
	# 2-sample test for equality of proportions with continuity correction
	#
	# data: c(plotDataDensity$Count[70], plotDataDensity$Count[79]) out of rep(sum(plotDataDensity$Count), 2)
	# X-squared = 39.399, df = 1, p-value = 3.455e-10
	# alternative hypothesis: two.sided
	# 95 percent confidence interval:
	# -0.03067201 -0.01598577
	# sample estimates:
	# prop 1 prop 2
	# 0.06636831 0.08969720
	# E Analysis --------------------------------------------------------------

view raw 08.2 Analysis.R hosted with ❤ by GitHub

Reflexiones

No me termina de hacer chiste cómo están representados los datos. Ese ajuste que tuve que hacer en la solución uno (líneas 10 y 11 de la solución 1) y que he arrastrado todo el post me hace sospechar que algo no del todo bueno a pasado entre los datos y la representación.
Por si acaso advierto que en un estudio real no miraríamos lo que he mirado. En un análisis haríamos comparaciones entre cuadrante y sujetos (grupos, condiciones) con un glm y así ver qué diferencias ha habido según una hipótesis.
Los cuadrantes han de ser definidos con cuidado para definir las áreas de interés en zonas claves según la literatura. Aún así, como definimos los cuadrantes de manera manual, siempre tiene un punto arbitrario y peligroso. Mi recomendación es definir las ROIs antes de hacer los análisis para evitar sesgos o tentaciones.
Tengo total desconocimiento sobre las ROIs en este tipo de estudios y por tanto cualquier análisis más complejo que hiciese no podría ser tomado en serio.
Todo esto ha costado de hacer, pero estoy bastante contento con los resultados. Espero que la currada haya merecido la espera. Gracias por seguirme.

Bibliografía

Ideas 1[1]
Ideas 2[2]
Guia para descargase csv de Dropbox [3]
Crear lineas en html [4]

Neurociencia y Matemáticas

miércoles, 25 de noviembre de 2015

Pequeño arreglo en la entrada anterior.

lunes, 23 de noviembre de 2015

Representar datos de Eyetracker Parte 2