C:\JUAN\POST\2020-07-30-GAM-MODELS-FOR-NEUROSCIENCE-CODE-101
Modelos GAM para Neurociência - Código 101
2020-07-30
GAMS PARA DOENÇAS NEURODEGENERATIVAS
Bem, este é um resumo da primeira parte da tese ou, pelo menos, da ideia original que tínhamos pensado e que agora está ligeiramente posta de lado. Muito resumidamente, a ideia é criar a database inicial a partir dos PET e MRI já normalizados com Matlab e SPM. Quando esta for um data.frame estável com todos os nossos dados (como também será para os SCC), poderemos criar modelos de regressão não paramétricos nos quais representar funções que unam pontos com a mesma intensidade estimada de atividade cerebral.
O objetivo inicial é obter imagens (gam.vis e gam.check) semelhantes às apresentadas abaixo. O código que apresento aqui resume o processo para as obter, embora ainda falte realizar uma etapa extra: comparar estes resultados com os resultados nativos obtidos com SPM, um software implementado em Matlab que é o padrão de mercado no mundo do diagnóstico por imagem médica (pelo menos no caso da neuroimagem). Enquanto não conseguir demonstrar de forma fidedigna que a minha metodologia é melhor do que a do SPM, serão apenas imagens muito bonitas.
#install.packages(c("cowplot","magick"))
library(magick);library(cowplot);library(ggplot2)
p1 <- ggdraw() + draw_image("GAMexample1.png")
p2 <- ggdraw() + draw_image("GAMexample2.png")
p3 <- ggdraw() + draw_image("GAMexample3.png")
p4 <- ggdraw() + draw_image("GAMexample4.png")
plot_grid(
p1, p2,p3, p4,
ncol = 2)

Claro que este é o mesmo problema dos SCC, a comparabilidade a nível matemático. Claro que nos SCC estou um passo à frente, pois consigo obter imagens que, a nível visual, são praticamente idênticas no formato, e pelo menos demonstrar visualmente que os resultados são semelhantes já é um resultado positivo. Talvez não seja o resultado final desejado, mas pelo menos estaria a demonstrar que existe outra forma, mais simples, de obter resultados semelhantes e a partir de uma abordagem de dados funcionais.
Mas agora o que nos diz respeito é obter esses modelos sobre diferentes cortes e ver a sua utilidade a nível visual e o que nos indicam. Como compará-los entre si e como compará-los com o SPM é outra história bem diferente.
PREÂMBULO:
#install.packages(c("mgcv","gamair","oro.nifti","memsic"))
library(mgcv);library(gamair);library(oro.nifti);library(memisc)
VISUALIZAÇÕES BÁSICAS:
Para a visualização de imagens PET ou MRI, eu não recomendaria de forma alguma a utilização do R. Existem vários pacotes para trabalhar com imagens, mas é evidente que o objetivo deles é trabalhar com os dados subjacentes à imagem, e não tanto com a imagem em si. Para trabalhar com imagens a nível visual, recomendaria antes o uso do SPM (clássico), do MRIcro ou do Mango. Qualquer um deles se assemelha mais a um editor de imagens padrão, algo como um Photoshop de cérebros. De qualquer forma, às vezes pode ser necessário apenas verificar se está tudo correto, extrair o valor específico de um píxel, fazer uma visualização ortogonal…
Embora o número de funções disponíveis seja limitado, aqui ficam algumas das mais úteis:
LEITURA DE IMAGENS E DADOS GERAIS:
img_003_S_1059 <- readNIfTI("003_S_1059", verbose = FALSE, warn = -1, reorient = TRUE,
call = NULL, read_data = TRUE)
img_003_S_1059
## NIfTI-1 format
## Type : nifti
## Data Type : 16 (FLOAT32)
## Bits per Pixel : 32
## Slice Code : 0 (Unknown)
## Intent Code : 0 (None)
## Qform Code : 2 (Aligned_Anat)
## Sform Code : 2 (Aligned_Anat)
## Dimension : 79 x 95 x 79
## Pixel Dimension : 2 x 2 x 2
## Voxel Units : mm
## Time Units : sec
VISUALIZAÇÕES MAIS ÚTEIS
Uma lista das diferentes visualizações que podemos obter, embora, como já disse, sejam muito mais completas e interativas noutros programas.
Axial:
image(img_003_S_1059) #axial por defecto

Coronal:
image(img_003_S_1059,plane="coronal") #coronal

Sagital:
image(img_003_S_1059,plane="sagittal") #sagital

Plano ortográfico:
orthographic(img_003_S_1059,text="Plano Ortografico PET") #plano ortografico
## Warning in min(x, na.rm = na.rm): ningún argumento finito para min; retornando
## Inf
## Warning in max(x, na.rm = na.rm): ningun argumento finito para max; retornando -
## Inf

DADOS ASSOCIADOS:
Estas imagens, evidentemente, têm uma série de dados associados (caso contrário, que tipo de data science iria eu fazer aqui?). Em geral, estes dados só me são úteis em bulk para obter um data.frame útil. No entanto, potencialmente posso pedir o valor de um dado concreto numa coordenada concreta, e isto pode ser útil para verificar, em certas etapas do desenvolvimento do código, se tudo está a funcionar corretamente. Por via das dúvidas, aqui fica:
Dados associados à imagem:
img_data=img_data(img_003_S_1059)
qplot(as.vector(img_data),
geom="histogram",
main = "Histogram for PET values",
xlab = "PET levels",
fill=I("blue"),
col=I("red"),
alpha=I(.2))

negatives=img_data[img_data<0]
length(negatives) # de hecho hay muchos valores inferiores a 0 por lo que lo más rápido será, en la database, igualarlos a 0 (MRIcro los identifica como 0's d ehecho)
## [1] 7590
Aqui, por exemplo, já vemos no histograma que há muitos valores próximos de zero ou até negativos, pelo que mais adiante será necessário ter isso em conta. Isto deriva provavelmente de um processo de normalização que reduziu valores nulos para valores inferiores a 0.
Coordenadas específicas na imagem:
#Diferentes coordenadas de la imagen con valores negativos (comprobacion)
img_data[1,2,3]
## [1] 0
img_003_S_1059[45,45,16]
## [1] 0.5555722
img_003_S_1059[10,70,40];img_data[10,70,40] # valdría cualquiera de los dos comandos
## [1] 0.3873922
## [1] 0.3873922
img_data[66,34,24]
## [1] 0.8744069