User Tools

Site Tools


neuroimagen:r_lmfit

Haciendo ajustes en R

Hay varias maneras de realizar regresiones lineales en R. Lo mas sencillo es calcular el coeficiente de correlación pero un enfoque mas avanzado puede ser el descrito en esta introduccion a R

Pearson, Spearman y demas pesca

La manera mas rapida de buscar correlaciones entre variables es utilizando la funcion cor. Aqui, por ejemplo, leo los datos de un archivo csv y escribo las correlaciones de pearson y spearman a otros dos archivos. Se supone que el archivo original es una tabla con los valores de las variables. Los archivos de salid son una matriz con las correlaciones.

read.csv("vol_ab_mri.csv", header = TRUE, sep = ";")-> x
cor(x, use="pairwise.complete.obs", method="pearson")-> y
y*y->z
write.csv(z,"pearson2.csv")
write.csv(y,"pearson.csv")
cor(x, use="pairwise.complete.obs", method="spearman")-> y
write.csv(y,"spearman.csv")

Modelos lineales

La funcion lm permite ajustar a un modelo lineal con covariables considerando o no la interacción. Mas despacio. Podemos a justar la variable $x$ a $f(x) = a_0 x + a_1$, o a $f(x) = a_0 x + a_1 y + a_2$, o incluso a $f(x) = a_0 x + a_1 y + a_2 x y + a_3$. ( de verdad, mirar aqui)

;-) funny, isn't it?

El siguiente codigo, lee los datos de una archivo csv y hace 668 ajustes covariados por edad.

Primeramente se seleccionan de la tabla importada las variables,

data.matrix(x[, 3])-> edad
data.matrix(x[, 16:684])-> brain
data.matrix(x[, 9:15])-> ab

Luego, dado que el numero de ajustes es elevado, se vuelca la salida de STDOUT a un archivo (ver estas notas sobre entrada y salida de datos).

sink(file = "covfit.txt", append = TRUE, type = "output", split = TRUE)

Ahora se crea un ciclo sobre cada variable a ajustar,

for(i in 16:684){
...
	}

se realiza el ajuste y se imprime al archivo,

	data.matrix(x[, i])-> b
	cf = lm(ab ~ b * edad)
	writeLines(paste("ROI:", names(x[i])))
	print(summary(cf))

:!: Aunque lo normal en R es que la salida se imprima en STDOUT dentro de los ciclos esto hay que hacerlo explicitamente con la funcion print().

Todo junto es algo asi:

covfit.R
read.table("full_ab_mri_norm_reorder.csv", header = TRUE, sep = ",")-> x
data.matrix(x[, 3])-> edad
data.matrix(x[, 16:684])-> brain
data.matrix(x[, 9:15])-> ab
sink(file = "covfit.txt", append = TRUE, type = "output", split = TRUE)
 
for(i in 16:684){
	data.matrix(x[, i])-> b
	cf = lm(ab ~ b * edad)
	writeLines(paste("ROI:", names(x[i])))
	print(summary(cf))
	}

Para construir una tabla con los resultados de cada ajuste hay que construir un parser. (O convertirse en chino? ;-))

parse_covfit.pl
#!/usr/bin/perl
 
# Copyright 2013 O. Sotolongo <osotolongo@fundacioace.com>
 
use strict; use warnings;
use Data::Dump qw(dump); 
my $ifile ="covfit_int.txt";
my $ofile = "covres_int.csv";
 
my %covtable;
my $itsaroi; my $itsab;
 
open DIF, "<$ifile" or die "Error:\n$!";
while(<DIF>){
	if((my $roi) = /ROI:\s(.*)$/){
		$itsaroi = $roi;
	}
	if((my $abeta) = /Response\s(.*)\s:/){
		$itsab = $abeta;
	}
	if(my ($mr2, $ar2) =/Multiple R-squared:\s+(\d+\.\d+e*\-*\d*),\s+Adjusted R-squared:\s+(-*\d+\.\d+e*\-*\d*)/){
		$covtable{$itsaroi}{$itsab}{'r2'} = $mr2;
		$covtable{$itsaroi}{$itsab}{'ar2'} = $ar2;
	}
	if(my ($pvalue) = /.*p-value:\s+(\d+\.\d+e*\-*\d*)$/){
		$covtable{$itsaroi}{$itsab}{'p-value'} = $pvalue;
	}
}
close DIF;
 
open ODF, ">$ofile" or die "Could not open file\n";
 
foreach my $abeta (sort keys %{$covtable{$itsaroi}}){
	print ODF " ,".$abeta."_r2, ".$abeta."_p";
}
print ODF "\n";
foreach my $roi (sort keys %covtable){
	unless(($roi =~ /.*NumVert$/) || ($roi =~ /.*FoldInd$/) || ($roi =~ /.*Std$/) || ($roi =~ /.*Curv.*/)){
		print ODF "$roi";
		foreach my $abeta (sort keys %{$covtable{$roi}}){
			print ODF ", $covtable{$roi}{$abeta}{'r2'}, $covtable{$roi}{$abeta}{'p-value'}"
		}
		print ODF "\n";
	}
}
close ODF;
neuroimagen/r_lmfit.txt · Last modified: 2020/08/04 10:58 (external edit)