The Little SAS Book contiene un excelente ejemplo para ilustrar las diferencias entre SAS como lenguaje de programación y C++ mostrando lo complicado que puede resultar procesar conjuntos de datos con un lenguaje de propósito general. Son 28 líneas de código C++ y 5 de SAS para leer un fichero delimitado e imprimirlo por pantalla. Es un ejemplo perfecto de cómo SAS es un lenguaje de cuarta generación con un alto nivel de abstracción y expresividad.
Queremos revisar esta comparación bajo otra perspectiva mostrando cómo SAS es el lenguaje perfecto para manejar estructuras de datos complejas y lo fácil que resulta implementar con él algoritmos complejos.
(Nota: con R basta una línea para leer e imprimir un conjunto de datos delimitado: print( read.table( "fichero.txt", sep = ";" )).)
En particular, muestra un pedazo de código para resolver el problema de las ocho reinas. El problema se reduce a encontrar una permutación de los números 1:n tales que
El código de SAS con el que resuelven este problema es así de estético, expresivo y comprensible:
No hay más color que el del resaltador de sintaxis, creo. Y en cuanto a la introducción del artículo, serán mis lectores los que habrán de decidir si tiene más que ver con las nueve que con las ocho reinas.
Publiqué hace un tiempo una entrada en esta bitácora sobre el problema que representa la desigualdad de los tamaños muestrales a la hora de comprender cierto tipo de datos, como por ejemplo, los que trata de representar el gráfico
Los gráficos de embudo son una alternativa pensada para evitar este tipo de sesgos. Por ejemplo
relaciona la proporción de casos de cáncer con el tamaño de la población añadiendo, si se me permite el término, curvas de isosignificancia para facilitar la comparación entre entidades desiguales en tamaño.
Iba a responder a los comentarios de mi entrada sobre las Jornadas de R y, muy en particular a los de Fernando Fernández, uno de los más fieles lectores de esta bitácora, y me he extendido tanto que he acabado convirtiéndola en una nueva. Pido excusas por haber tal vez abusado de mis prerrogativas para auparme de esta manera.
Tanto a él como a otros les chirrió que escribiese comenzamos una nueva época que en el plazo de tres o cuatro años nos va a conducir, con casi total seguridad, a un escenario en el que [...] R se use de manera casi exclusiva en la enseñanza de la estadística en los niveles universitarios.
Conste en primer lugar que es una frase cosecha mía y no sé cuántos, si algunos, estarán de acuerdo conmigo. Y se trata de una mezcla de visión y de deseo. Es cierto que llamadas a la uniformidad —y a la inmersión lingüística— pueden provocar sentimientos de rechazo. También es cierto que se trata de una sola frase que compendia muy sucintamente horas de reflexiones, conversaciones y experiencias y que me deja, por lo tanto, expuesto a malinterpretaciones extremas. Yo mismo los sentí al escribir la frase y me obligó a releerla. ¿Me atrevería a dejarla tal cual?
Me concedí finalmente la venia para no cambiarle una coma obligándome en contrapartida a fundamentarla debidamente. Y esta es la ocasión propicia para hacerlo.
Fernando dice que consider[a] importante un esfuer[z]o por trabajar al menos con dos herramientas distintas y estoy de acuerdo con él. Creo que, de hecho, uno de los objetivos estratégicos (transcenciendo los contenidos de cada una de las asignaturas de la carrera de estadística) debería ser que los alumnos aprendiesen a programar.
He estudiado estadística en dos universidades. En la una, la estadística era esencialmente teórica y jamás usamos ordenador para ningún fin. En la otra sólo había SAS. Hace mucho que dejé la universidad y temo equivocarme. Pero sí que conozco a muchos recién licenciados, he tenido ocasión de haber hablado con algunos profesores, etc. Y me da la impresión de que el recién licenciado medio no sabe programar. En el mejor de los casos, lo han expuesto a R en dos asignaturas, a Gretl en la de series temporales, a WinBUGS si hizo un máster especializado, a SPSS en quién sabe qué y, a lo más, hizo un curso de SAS de 20 horas. Pero le pides, aquí y ahora, un gráfico y abre Excel; y le pides una regresión y pulsa la tecla del hiperespacio (la del juego Asteroids) para evadirse a un dimensión remota.
No hablo desde una perspectiva, si se me permite, mini-max —un mínimo entre las exigencias máximas— sino max-mini: a la vista de las circunstancias actuales, presento una perspectiva más deseable.
Si se me concede que un objetivo de mínimos es lograr que todo recién licenciado en estadística (y áreas afines) domine un, al menos un, lenguaje de programación, —e insisto, sin perjuicio de que algunos sean también expertos en otros— la siguiente pregunta que se plantea es ¿cuál?
Y creo que ahora, en este año, aunque las cosas pudieran cambiar dentro de cinco o diez más, debiera ser R. Por muchos motivos:
Porque R puede instalarse en cualquier plataforma, usarse en cualquier momento, para cualquier fin y gratuitamente: se es estadístico 24×7, no 8×5.
Porque expone a los estudiantes a herramientas y métodos de primera línea de investigación.
Porque se usa para resolver problemas reales.
Porque es un lenguaje (relativamente) moderno: está orientado a objetos, contiene aspectos funcionales, etc. que facilitan el aprendizaje de otros lenguajes tales como Python, Java, etc.
Porque —en contraposición a herramientas propietarias— no obliga a universidades (muchas de ellas públicas) a hacer el trabajo sucio de formación por cuenta de unos señores que tienen un negocio de venta de software.
Porque introduce a los estudiantes en el mundo del software libre y les da la oportunidad de conocer un modo muy valioso de hacer cosas —ciencia y software, entre otras— y divulgar conocimiento.
Porque su uso desincentiva la piratería de software.
Y, seguro, me dejo muchos otros.
Si se me pregunta por alternativas y por SAS en concreto, resulta que:
SAS no se usa apenas en estadística: basta con revisar la agenda del último SAS Forum España para darse cuenta de que la empresa está en otra cosa.
Gran parte de los programadores en SAS no han llamado jamás a PROC MIXED, PROC PHREG, etc. Son puros trasladadores de datos. Y yo doy fe: he pagado la hipoteca de un apartamento traslandando datos de aquí para allá con SAS.
SAS no enseña nada que valga la pena ser aprendido hoy en día en cuanto a técnicas de programación: pudo haber resultado innovador en los años setenta; hoy ya no: es un callejón sin salida que crea más vicios que abre puertas a otros lenguajes.
Quiero acabar esta entrada bastante larga analizando la transcendencia de esa descripción de R que es la de lengua vehicular —en inglés suele utilizarse el término lingua franca— de la estadística y que tiene que ver con el hecho de que el lenguaje da forma al pensamiento.
Y la estadística clásica —que constitutye el grueso de cuanto se enseña en la universidad— está lastrada por los condicionantes de quienes, sin la ayuda de los recursos computacionales existentes hoy en día, la desarrollaron hace un siglo. Y que el ámbito de los problemas a los que se aplica está igualmente restringido.
Hoy, en el siglo XXI, existen problemas y técnicas mucho más diversos que la nueva lengua vehicular —que puede ser la puerta de entrada a métodos basados en remuestreos, a la simulación, a la predicción, etc.— puede ayudar a poner encima de la mesa para renovar de arriba a abajo la disciplina.
Esta vez toca compararlos a la hora de aplicar el método del jackknife.
Primero, el código SAS que recomienda el autor del artículo, que calcula la curtosis de un conjunto de datos trivial (una muestra de 10k valores que siguen una distribución uniforme):
data test;
do i = 1 to 10000;
x = ranuni(1234);
output;
end;
keep x;
run;
data outb;
do replicate = 1 to numrecs;
do rec = 1 to numrecs;
set test nobs=numrecs point=rec;
if replicate ^= rec then output;
end;
end;
stop;
run;
ods listing close;
proc univariate data=outb;
var x;
by replicate;
output out=outall kurtosis=curt;
run;
ods listing;
proc univariate data=outall;
var curt;
output out=final mean=jmean std=jstd;
run;
Tarda en ejecutarse 110 segundos en mi máquina. He probado también a aplicar el procedimiento a conjuntos de datos con 50k muestras, pero el sistema no ha dado de sí: se ha quedado sin espacio en disco.
(Nota: véanse los comentarios a algunas de las entradas anteriores de la serie. Algunos de los lectores de esta bitácora, de alguna manera, han logrado evitar el cuello de botella que para SAS supone depender del acceso a disco manteniendo un uso de CPU decente que yo, dicho sea de paso, no he visto casi nunca: una CPU al 100% en un ordenador que corre SAS es una rareza digna de fotografiarse.)
En R tenemos varias alternativas. Tal vez la más simple sea
que utiliza la función jackknife del paquete bootstrap y tarda 13.12 segundos. Y con 50k muestras, unos 10 minutos (hummm… falta de linealidad). Pero merece la pena probar en este tipo de contextos el paralelismo en R.
Para ello, utilizaremos el paquete doSMP:
library(doSMP)
w <- startWorkers(workerCount = 2)
registerDoSMP(w)
La primera línea carga el paquete. Las otras dos registran los workers dicho de otra manera, especifican cuántos hilos abrirá R para ejecutar las tareas que se paralelicen y los activan. Obviamente, no tiene sentido asignar más hilos que CPUs tenga el sistema.
Para ejecutar las tareas se utiliza foreach. Esta función es, en cierto modo, similar a for. Puede cualificarse con %dopar% o con %do%: la primera opción ejecuta el código en paralelo; la segunda, secuencialmente.
El código (con las dos opciones) es:
# en paralelo
res <-foreach( i = 1:N,.packages = "moments") %dopar%
kurtosis( x[-i])# iterativo
res <-foreach( i = 1:N,.packages = "moments") %do%
kurtosis( x[-i])
Es necesario pasarle a foreach el nombre de los paquetes necesarios para ejecutar el código subsiguiente: si ejecutamos las líneas anteriores omitiendo el parámetro .packages, R se quejará por no poder encontrar la función kurtosis.
Los resultados no han resultado particularmente satisfactorios. Con la segunda opción, la secuencial, ha tardado 23 segundos. Con la primera, en paralelo, 14.58. Parece que el uso de foreach implica una sobrecarga computacional sustancial que consume los beneficios de la paralelización.
Construir un cojunto de datos simples (dos vectores, x e y).
Hacer una regresión de y sobre x y capturar los residuos.
Crear 1000 vectores y' distintos añadiendo a (la predicción de y) en el modelo anterior una reordenación de los residuos.
Crear los correspondientes 1000 modelos haciendo la regresión de cada sobre x.
Obtener el histograma del coeficiente de la regresión.
Es un caso de bootstrap en el que no se muestrean directamente los valores iniciales sino los residuos del modelo.
El código en SAS es el siguiente:
/* 0: data creation */
data temp1;
x=1; y=45; output;
do x = 2 to 29;
y = 3*x + 6*rannor(1234);
output;
end;
x=30; y=45; output;
run;
%let regressors = x;
%let indata = temp1;
/* 1: perform the regression and get the predicted and residual values */
proc reg data= &INDATA;
model y=®ressors;
output out=out1 p=yhat r=res;
run;
/* 2: split the data: only the residuals will require URS */
data fit(keep=yhat ®RESSORS order) resid(keep=res);
set out1;
order+1;
run;
/* 3: this doesn't do any sampling - it copies the FIT data set repeatedly */
proc surveyselect data=fit out=outfit method=srs samprate=1 rep=1000;
run;
/* 4: this does the WR sampling of residuals for each replicate */
data outres2;
do replicate = 1 to 1000;
do order = 1 to numrecs;
p = ceil(numrecs * ranuni(394747373));
set resid nobs=numrecs point=p;
output;
end;
end;
stop;
run;
/* 5: then the randomized residuals are merged with the unrandomized records */
data prepped;
merge outfit outres2;
by replicate order;
new_y=yhat+res;
run;
/* 6: the bootstrap process runs on each replicate */
proc reg data=prepped outest=est1(drop=_:);
model new_y=®ressors;
by replicate;
run;
/* 7: and the sampling distribution is aggregated */
proc univariate data=est1;
var x;
output out=final pctlpts=2.5, 97.5 pctlpre=ci;
run;
proc print;
run;
que corre en mi máquina en 37.37 segundos.
Ofrezco dos alternativas sustancialmente más sucintas en R. La primera es una reinterpretación literal en R del código anterior,
x <-2:29
y <-3* x +6*rnorm(length(x))
x <-c(1, x,30)
y <-c(45, y,45)
m0 <-lm( y ~ x )
yhat <- m0$fitted.values
res <- m0$residuals
resultados <-replicate(1000,lm( yhat +sample( res )~ x )$coefficients[["x"]])hist( resultados )
que corre en 2.46 segundos.
El segundo utiliza el paquete boot,
library(boot)
datos <-data.frame( y = yhat, x = x, res = res )
foo.boot <-function( datos, ind ){lm(y + res[ind]~ x,data = datos)$coefficients[["x"]]}
res <-boot( datos, foo.boot, R = 1000)plot( res )
que corre en 2.71 segundos.
El paquete boot es un candidato ideal para reimplementarlo utilizando paralelismo. Sin embargo, su mantenedor no parece estar por la labor. ¡Qué pena!
Corre el rumor de una posible venta de SAS. Pueden ser un simple rumor pero se non è vero, è ben trovato: el máximo responsable y accionista mayoritario de SAS, Jim Goodnight, tiene ya 68 años y la empresa está sufriendo el acoso de la competencia en muchos frentes.
SAS quiso pasar de ser una compañía que especializada en herramientas de estadística a otra que proporcionase un entorno completo de herramientas del tipo de las denominadas de business intelligence. Y en ese esfuerzo topó con los grandes. Y con los (cada vez menos) pequeños, como R, en su nicho originario. De ahí que, según datos de la consultora Gartner, sea el único de los grandes proveedores de este tipo de soluciones en perder cuota de mercado:
Existe una conjunción de factores. Existen rumores. Incluso apuestas.
Don’t be loopy! es el título de una presentación realizada en el SAS Global Forum de 2007. Tiene que ver con el motivo que me hizo en mi día abandonar SAS y buscar —entonces aún no lo conocía— el cobijo de R: sus limitaciones para todo lo que tiene que ver con simulaciones, remuestreos, jackknifes, bootstraps y similares.
El artículo muestra lo que debería ser el estado del arte para realizar este tipo de programas con SAS. En el primero de los problemas que estudia, que denomina bootstrap simple, muestrea 1.000 veces un conjunto de datos de 50.000 observaciones y calcula el valor de la curtosis para cada una de ellas. Finalmente, proporciona un intervalo de confianza para dicho valor.
El artículo recomienda no utilizar bucles. Y mucho menos, bucles usando las llamadas macros. Propone usar el PROC SURVEYSELECT para muestrear los datos. Incluso, dado que SURVEYSELECT lee el fichero de entrada del disco una vez por muestra —es decir, 1000 veces en total en el ejemplo—, propone el comando sasfile para copiar los datos en RAM. La sintaxis es la siguiente:
El código recomendado para resolver en SAS este problema es
data YourData;
do i = 1 to 50000;
x = ranuni(1234);
output;
end;
keep x;
run;
sasfile YourData load;
proc surveyselect data=YourData out=outboot
seed=30459584
method=urs samprate=1 outhits
rep=1000;
run;
sasfile YourData close;
ods listing close;
proc univariate data=outboot;
var x;
by Replicate;
output out=outall kurtosis=curt;
run;
ods listing;
proc univariate data=outall;
var curt;
output out=final pctlpts=2.5, 97.5 pctlpre=ci;
run;
En mi ordenador necesita 70 segundos para ejecutarse y crea un fichero, outboot.sas7bdat, que ocupa 1,2 GB en mi disco duro. El código equivalente en R es
que no genera ficheros en disco, no incrementa de manera perceptible la memoria utilizada y corre en 8 segundos sobre la misma máquina.
¿Será por esto que SAS no quiere venir a las jornadas de usuarios de R? ¿No será más bien que tendremos que ir nosotros a las suyas para enseñarles cómo se han de hacer las cosas?
Por ser viernes, traigo a estas páginas un vídeo tan pedagógico como ameno. Es la conferencia de Dick De Veaux dentro la M2010 Data Mining Conference auspiciada por SAS.
El autor repasa los siete pecados capitales de la minería de datos, a saber
No realizar las preguntas adecuadas
No entender el problema correctamente
No prestar suficiente atención a la preparación de los datos
Ignorar lo que no está ahí
Enamorarse de los modelos
Trabajar en solitario
Usar datos malos
Frente a ellas, propone las siguientes virtudes:
Define el problema
Prepara los datos usando conocimiento sobre el campo del que proceden
Mantente dispuesto y preparado para aplicar nuevas ideas y modelos
Ten en cuenta los valores no informados: crea variables derivadas
Trabaja en equipo
Asegúrate de la calidad de los datos
Usa modelos, no únicamente asociaciones
¡Ah! Y que nadie se pierda alrededor del minuto 7:30 el icono que aparece en la esquina inferior izquierda del escritorio del ordenador.
Acaba de llegarme la noticia de que la versión 9.3 de SAS (sí, el producto de esa empresa que no quiere saber nada de las III Jornadas de Usuarios de R a pesar de que las palabras de su director general en España nos hicieran creer a algunos lo contrario) que, como de costumbre, es lo mejor de lo mejor. Entre los cambios grandes y pequeños que aporta están:
No es necesario pasar de de SAS 9.1.3 a SAS 9.2 para instalar la versión 9.3
Para producir gráficos con ODS no es necesaria la licencia de SAS/GRAPH
Los procedimientos gráficos SGPANEL, SGPLOT, SGRENDER y SGSCATTER ha pasado de SAS/GRAPH a SAS Base.
El nuevo procedimiento FMM de SAS/STAT permite ajustar modelos de mezclas finitas
Se pueden leer tablas de JMP en SAS directamente.
Y algunas más que podrán consultar los interesados en la página de SAS 9.3.
Recientemente tuvo lugar la conferencia del nosequé de SAS en algún lugar de EE.UU. Alguien decidió rodar el siguiente vídeo:
En él aparecen algunos de los participantes en las conferencias realizando comentarios simpáticos. Pero conforme iba viendo desfilar rostros, no dejaba de pensar en que existía un patrón en la muestra.
Y me pregunto, ¿cuál es la pirámide de edad de la comunidad de usuarios de SAS? ¿Y la los de R? ¿Será SAS “el futuro”?
(la predicción de 
