Archivo de la categoría: Modelos

Lectura conjunta de “Analyzing Linguistic Data”

Mi compañero y amigo Carlos Gil junto con Victor Peinado van a iniciar un interesante proyecto. Un grupo de lectura del libro Analyzing Linguistic Data: A practical introduction to Statistics desde aquí nos hacemos eco de este proyecto que podéis seguir desde aquí.

El calendario del curso es:

  • An Introduction to R: (1 semana, del 6/05 al 12/05)
  • Graphical Data Exploration: (1 semana, del 13/05 al 19/05)
  • Probability Distributions: (1 semana, del 20/05 al 26/05)
  • Basic Statistical Methods: (2 semanas, del 27/05 al 9/06)
  • Clustering and Classification: (2 semanas, del 10/06 al 23/06)
  • Regression Modeling: (2 semanas, del 24/06 al 7/07)
  • Mixed Models: (2 semanas, del 8/07 al 21/07)

Será todo un éxito, como las anteriores convocatorias. Saludos.

Trucos SAS. Medir la importancia de las variables en nuestro modelo de regresión logística

Hoy quería proponeros una forma poco ortodoxa de medir la importancia de las variables en un modelo de regresión logística con SAS. La cuestión es: dado un modelo de regresión logística crear un ranking con las variables más importantes dentro del modelo. Para esta tarea recomiendo el uso de random forest, pero puede ser imposible disponer de un software que realice este tipo de modelos. Imaginemos un caso concreto, disponemos de SAS STAT y nos da reparo trabajar con R. Para este caso podemos emplear el siguiente truco. El AIC (criterio de información de Akaike) es un estadístico que relaciona el cociente de la verosimilitud con el número de parámetros del modelo que ajustamos. Cuanto menor sea este cociente mejor será nuestro modelo. Si eliminamos una variable del modelo ¿cuánto empeora este modelo? Esa será la filosofía que emplearemos para analizar la importancia de las variables presentes en nuestro modelo. En la línea habitual hacemos un ejemplo para que podáis copiar y pegar en vuestro SAS:

Vamos a crear un dataset preparado para hacer una regresión logística perfecta donde en un 10% de los casos sucede un evento:

*REGRESION LOGISTICA PERFECTA;
data logistica;
do i=1 to 10000;
  normal1=rannor(8);
  normal2=rannor(45);
  normal3=rannor(32);
  normal4=rannor(7);
  normal5=rannor(98);
  unif1=ranuni(2);
  unif2=ranuni(21);
  unif3=ranuni(22);
  unif4=ranuni(23);
  unif5=ranuni(24);
  prob=1/(1+exp(-(-3.16+0.1*normal1-0.2*normal2+0.3*normal3-0.4*normal4+0.5*normal5+
	   0.1*unif1+0.2*unif2+0.3*unif3+0.4*unif4+0.5*unif5)));
  sucede=ranbin(8,1,prob);
  *TRAMIFICAMOS LAS VARIABLES;
  normal1=round(rannor(8),0.1);
  normal2=round(rannor(45),0.2);
  normal3=round(rannor(32),0.3);
  normal4=round(rannor(7),0.4);
  normal5=round(rannor(98),0.5);
  unif1=round(ranuni(2),0.1);
  unif2=round(ranuni(21),0.2);
  unif3=round(ranuni(22),0.3);
  unif4=round(ranuni(23),0.4);
  unif5=round(ranuni(24),0.5);
  output;
end;
drop i;
run;

title "Logistica con un 10% aprox de casos positivos";
proc freq data=logistica;
tables sucede;
quit;

Incialmente necesitamos las variables presentes en el modelo y el ajuste inicial, también un conjunto de datos SAS con los nombres de las variables. Esto es un poco chapuza, pero si seguís el blog podéis hacer este código mucho más elegante [no os lo voy a dar todo hecho] Seguir leyendo Trucos SAS. Medir la importancia de las variables en nuestro modelo de regresión logística

Parámetro asociado a una Poisson con SAS

Mirad que he visto datos en mi vida. Y esos datos siguen muchas distribuciones. Y una de las distribuciones más habituales con las que me he encontrado es la distribución de poisson. Esta distribución tiene una característica muy interesante: la varianza es igual que la media. Y si la varianza no es igual a la media tenemos distribuciones de poisson sobredispersa o poisson infradispersa con propiedades muy interesantes y que se emplea mucho en el ámbito actuarial, aunque tendremos eventos con una distribución de poisson cuando estamos hablando de eventos independientes en intervalos de tiempo. No soy yo el más adecuado para escribir sobre el modelo matemático que tienen detrás estas distribuciones, pero si me gustaría mostraros como hacer mediante SAS con el PROC GENMOD algo tan básico como obtener el parámetro asociado a mi distribución de poisson y el intervalo de confianza al 95% para este parámetro. Vale que el parámetro es la media pero tengo que escribiros un código SAS importante ¿Y cómo lo calculamos?

 *SIMULAMOS EL NÚMERO DE VISISTAS A UNA WEB
 EN HORARIO DE OFICINA;
data poisson;
do hora=8 to 17 by 0.5;
visitas = ranpoi(2,15);
output;
end;
run;
*GRAFICAMOS LA DISTRIBUCIÓN;
proc gchart data=poisson;
vbar hora /freq=visitas discrete;
run;quit;
*EMPLEAMOS EL PROC GENMOD PARA LA OBTENCIÓN
 DEL PARÁMETRO Y EL INTERVALO;
proc genmod data=poisson;
   model visitas = / dist=poisson;
   ods output parameterestimates=param;
run;
*ESTÁ EN ESCALA LOGARÍTMICA;
data param;
set param;
if _n_=1;
lamda = exp(estimate);
sup=exp(uppercl);
min=exp(lowercl);
call symput ('lamda',lamda);
run;
*CON PROC MEANS;
proc means data=poisson mean clm ;
var visitas;
quit;

Una entrada que puede parecer una tontería, pero que es necesaria para desordenar mi conciencia. Espero que os sea de utilidad, un saludo.

Test de Levene con SAS

El Test de Levene se aplica para la igualdad de varianzas. Es un análisis de la varianza de las desviaciones de los valores muestrales respecto a una medida de tenedencia central. Parte de la hipótesis nula de igualdad de varianzas. Para realizar este test en SAS emplearemos el PROC GLM en combinación con la opción HOVtest. En la línea habitual vemos un ejemplo:


data datos;
 input presion @@;
 if _n_<=5 then grupo=1;
 else if _n_<=10 then grupo=2;
 else if _n_<=15 then grupo=3;
 else if _n_<=20 then grupo=4;
 else grupo=5;
cards;
180 172 163 158 147 173 158 170
146 152 175 167 158 160 143 182
160 162 171 155 181 175 170 155 160
;
run;
 Proc GLM Data=datos ;
  Class grupo;
  Model presion = grupo;
  Means grupo/HOVTest;
  ODS select HOVFTest;
 Quit;

Vemos que se trata de un código sencillo donde modelizamos con GLM la variable dependiente con la variable grupo y en MEANS indicamos con HOVtest que deseamos que se realice el test de Levene, con ODS seleccionamos sólo esa salida. Este código podemos parametrizarlo y crear una macro que nos permita replicar el código:

%Macro Levene(datos,grupo,variable);
Proc GLM Data=&datos ;
Class &grupo;
Model &variable=&grupo;
Means &grupo/HOVTest;
ODS select HOVFTest;
Quit;
%Mend Levene;
%levene(datos,grupo,presion);

Espero que sea de vuestra utilidad. Saludos.

La distribución tweedie

tweedie.png

Reconozco que hace muy poco tiempo que trabajo con las distribuciones tweedie. Un viejo dinosaurio que trabaja sobre todo con SAS se hace el sordo cuando le hablan de la distribución tweedie. Quizá sea el trabajo con SAS el que me ha nublado. Pero ahora que empiezo a trabajar con otras herramientas… Para comprender mejor la base teórica para este tipo de distribuciones os enlazo a la wikipedia. Pero despierta mi interés debido a que se puede considerar una gamma con punto de masa en el 0 ¡toma aberración matemática! Aspecto interesante.

Este tipo de distribución necesita 3 parámetros; p que nos indica el tipo de distribución en R se denomina power, mu que es la media y phi que es la desviación típica. En R disponemos del paquete tweedie para trabajar con este tipo de distribuciones, la función rtweedie es la que genera números aleatorios según una tweedie con parámetros p, mu y phi:

#install.packages(“tweedie”)
library(tweedie)
#Números aleatorios distribuidos según una tweedie
y1 <- rtweedie( 10000, p=1, mu=600, phi=1000)
summary(y1)y1.5 <- rtweedie( 10000, p=1.5, mu=600, phi=1000)
summary(y1.5)
y2 <- rtweedie( 10000, p=2, mu=600, phi=1000)
summary(y2)
par(mfrow=c(2,2))
plot(density(y1),main="Densidad con parámetro p=1")
plot(density(y1.5),main=" Densidad con parámetro p=1.5")
plot(density(y2),main=" Densidad con parámetro p=2")

Más que interesante el resultado obtenido. Con p=1 fijaos que forma tiene la función de densidad con distintos puntos de masa, con p=1.5 estamos ante una distribución que se parece mucho a una gamma salvo por un “pequeño detalle” tenemos muchos ceros. Con p=2 deberíamos tener una gamma. Pero es entre el 1 y el 2 donde esta distribución es más interesante sobre todo si se aplica a modelos de riesgo en estudios actuariales ya que nos permite modelizar directamente la prima de riesgo sin ser necesario modelizar el coste y la frecuencia por separado.

 

Da comienzo la lectura de “The Elements of Statistical Learning”

Interesante iniciativa de Juanjo Gibaja y Carlos Gil Bellosta consistente en leer el libro “The Elements of Statistical Learning”. La idea es crear un grupo de trabajo que lea de forma coordinada el libro, pero además van a implementar en R los ejercicios y ejemplos. Desde esta bitácora no vamos a colaborar en la lectura pero si vamos a colaborar en la divulgación de resultados. Seguiremos muy de cerca este proyecto.

El sobremuestreo ¿mejora mi estimación?

El sobremuestreo (oversampling) es una técnica de muestreo que se emplea habitualmente cuando tenemos una baja proporción de casos positivos en clasificaciones binomiales. Los modelos pueden “despreciar” los casos positivos por ser muy pocos y nuestro modelo no funcionaría. Para incrementar el número de casos positivos se emplea el sobremuestreo. Ejemplos habituales pueden ser los modelos de fraude, un 99% de las compras son correctas, un 1% son fraudulentas. Si realizo un modelo puedo estar seguro al 99% de que todas mis compras son correctas, en este caso hemos de realizar un sobremuestreo para incrementar nuestros casos de fraude y poder detectar los patrones.

Personalmente no sabría deciros el porcentaje de casos positivos a partir del cual sería necesario llevar a cabo un proceso de remuestreo. A mi particularmente me gusta hacerlo siempre. Por lo menos realizar algunas pruebas para identificar aquellas variables que son más influyentes y comenzar a eliminar aquellas que no van a funcionar. Busco exagerar. Tampoco me quiero mojar mucho sobre la proporción de casos positivos y negativos, pero si estamos realizando un nuevo muestreo podemos emplear perfectamente un 50% para ambos, aquí si que dependemos del número de registros con el que estemos trabajando ya que al final el sobremuestreo será la repetición de los casos positivos sobre la tabla de entrada del modelo.

Sin embargo, cuando ya tengo decidido como va a ser mi modelo no me gusta realizar sobremuestreo. Lo considero un paso previo (algún lector del blog considerará estas palabras incoherentes). Después de toda esta exposición teórico-práctica de malos usos de un dinosaurio  en realidad lo que cabepreguntarse es ¿mejora la estimación un modelo con sobremuestreo?

Abrimos R y Tinn-R y manos a la obr Seguir leyendo El sobremuestreo ¿mejora mi estimación?

Árboles de decisión con SAS Base (con R por supuesto)

Con SAS Base podemos hacer árboles de decisión porque tenemos R. Así de sencillo. Vamos a utilizar SAS para gestionar nuestros datos y R será la herramienta que utilicemos para la realización del modelo de árbol de decisión. Posteriormente emplearemos las reglas generadas por el modelo para etiquetar a nuestros clientes en SAS. Con esta entrada pretendo ilustrar una serie de ejemplos en los que comunico SAS con R. Una herramienta nos sirve para el tratamiento de datos y la otra la utilizaremos para realizar modelos que no están al alcance de SAS. Para realizar esta comunicación SAS-R os planteo la creación en SAS de ficheros de texto con las instrucciones en R y la ejecución en modo batch de R con ese código creado en SAS. Aquí tenéis punto por punto el ejemplo:

El primer paso, como es habitual, es crear un conjunto de datos SAS con datos aleatorios que nos sirva de ejemplo Seguir leyendo Árboles de decisión con SAS Base (con R por supuesto)

Trucos R. Llevar a SAS las reglas de un árbol de decisión

Vuelvo hoy con el uso de rpart para la creación de árboles de decisión con R. Pero hoy, además de realizar un modelo de árbol con R quiero presentaros una función que nos permite guardar las reglas generadas con nuestro modelo en un fichero de texto para su posterior utilización con SAS. Retomamos un ejemplo visto con anterioridad en la bitácora con ligeras modificaciones:
#Inventamos un objeto para realizar el modelo
#En una cartera de clientes nuestro modelo tiene que identificar
#cuales contratan un PVI
#
clientes=20000
saldo_vista=runif(clientes,0,1)*10000
saldo_ppi=(runif(clientes,0.1,0.2)*rpois(clientes,1))*100000
saldo_fondos=(runif(clientes,0.1,0.9)*(rpois(clientes,1)-1>0))*100000
edad=rpois(clientes,60)
datos_ini<-data.frame(cbind(saldo_vista,saldo_ppi,saldo_fondos,edad))
datos_ini$saldo_ppi=(edad<=68)*datos_ini$saldo_ppi
#Creamos la variable objetivo a partir de un potencial
datos_ini$potencial=runif(1,0,1)+
(log(edad)/(log(68))/100) +
runif(1,0,0.001)*(saldo_vista>5000)+
runif(1,0,0.001)*(saldo_fondos>10000)+
runif(1,0,0.007)*(saldo_ppi>10000)-
runif(1,0,0.2)
datos_ini$pvi=as.factor((datos_ini$potencial>=quantile(datos_ini$potencial,
0.90))*1)
#
#Empleamos rpart para la realización del modelo
#
library(rpart)
arbol=rpart(as.factor(pvi)~edad+saldo_ppi+saldo_fondos,
data=datos_ini,method="anova",
control=rpart.control(minsplit=30, cp=0.0008) )

Tenemos un objeto rpart llamado arbol. En este punto necesitamos disponer de las reglas generadas por el modelo para SAS, donde el módulo específico para poder realizar determinados modelos tiene un precio muy alto. Buscando en Google encontraremos este link. En él tenemos una genial función de R list.rules.rpart que nos permite identificar las reglas que ha generado el modelo Seguir leyendo Trucos R. Llevar a SAS las reglas de un árbol de decisión