Table of Contents
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;