User Tools

Site Tools


neuroimagen:pcor

Como hacer correlaciones parciales en R

La correlacion parcial es una medida de asociacion entre dos variables (Lee aquí)

Esto lo hice usando el paquete pcor

library(ppcor)
read.csv("full_paper_data.csv", header = TRUE, sep = ",")-> x
y.data <- data.frame(x["AB140"],x["edad"],x["Sexo"],x["APOE_mas"],x["Estud"],x["CREA"],x["LeftHippocampus"])
y.data <- y.data[complete.cases(y.data),]
pcor.test(y.data$AB140,y.data$LeftHippocampus, y.data[,c("edad","Sexo","APOE_mas","Estud","CREA")])

Escribiendo la matriz de covarianza

w <- pcor(y.data)
write.csv(w$estimate, file="ab140_cov.csv")

Hasta ahi lo básico. Si se quiere hacer lo mismo para varias variables hay que hacer su programilla en R. Necesitariamos tener en dos archivos separados las listas de variables. Por ejemplo, el codigo siguiente calcula las correlaciones parciales entre el SUVR medido en las regiones de interes y el amiloide en plasma. La lista de ROIs esta en el archivo pib_fsl.rois y la lista de medidas de amiloide en plasma en amyloids.list. Los datos estan en el archivo recerc_ab_pib.csv. el resultado lo escribe en el archivo pcors.txt.

library(ppcor)
read.csv("recerc_ab_pib.csv", header = TRUE, sep = ",")-> x
scan("pib_fsl.rois", what = character())->rois
scan("amyloids.list", what = character())-> betas
sink(file = "pcors.txt", append = TRUE, type = "output", split = TRUE)
for(i in 1:length(rois)){
	for(j in 1:length(betas)){
		y.data <- data.frame(x[betas[j]],x["NEdad"],x["CSexo"],x[rois[i]])
		y.data <- y.data[complete.cases(y.data),]
		pcv <- pcor.test(y.data[betas[j]],y.data[rois[i]], y.data[,c("NEdad","CSexo")])
		writeLines(paste("ROI:", rois[i],"BETA:", betas[j], "PCOR:", pcv$estimate, "PVALUE:", pcv$p.value))
	}
}
sink()

Para construir una tabla con estos resultados mejor hacerse un parser.

parse_pcors.pl
#!/usr/bin/perl
 
# Copyright 2013 O. Sotolongo <osotolongo@fundacioace.com>
 
use strict; use warnings;
use Data::Dump qw(dump);
my $ifile = "pcors.txt";
#my $ofile = "ab_cortical_volume_pcor.csv";
my $ofile = "ab_pib_recerc.csv";
my %covtable;
my $itsaroi; my $itsab;
open DIF, "<$ifile" or die "Error:\n$!";
while(<DIF>){
	#ROI: lh_frontalpole_GrayVol BETA: total_amyloid PCOR: 0.0326186447676574 PVALUE: 0.77890554238218
	if(my ($roi, $beta, $pcor, $pvalue) = /^ROI:\s(.*)\sBETA:\s(.*)\sPCOR:\s(.*)\sPVALUE:\s(.*)$/){
		$covtable{$roi}{$beta}{'pcor'} = $pcor;
		$covtable{$roi}{$beta}{'p-value'} = $pvalue;
		$itsaroi = $roi;
	}
}
close DIF;
 
open ODF, ">$ofile" or die "Could not open file\n";
foreach my $abeta (sort keys %{$covtable{$itsaroi}}){
	print ODF " ,".$abeta."_pcor, ".$abeta."_pvalue";
}
print ODF "\n";
foreach my $roi (sort keys %covtable){
	print ODF "$roi";
	foreach my $abeta (sort keys %{$covtable{$roi}}){
		print ODF ", $covtable{$roi}{$abeta}{'pcor'}, $covtable{$roi}{$abeta}{'p-value'}"
	}
	print ODF "\n";
}
close ODF;
neuroimagen/pcor.txt · Last modified: 2020/08/04 10:58 (external edit)