Qué pasa si uso una regresión de poisson en vez de una regresión logística

Para un tema de mi trabajo voy a utilizar una regresión de poisson en vez de una regresión logística, el evento es si o no y no tiene nada que ver el tiempo, ni se puede contabilizar como un número, pero a efectos prácticos es mejor para mi usar una regresión de poisson. Entonces, ¿qué pasa si hago una poisson en vez de binomial? Como siempre si mi n es muy grande hay relación entre ambas distribuciones. Pero yo quiero saber si me puede clasificar mis registros igual una regresión de poisson y una binomial y se me ha ocurrido hacer un ejercicio teórico muy simple.

Construyo con SAS 10.000 datos aleatorios con las variables independientes x e y normalmente distribuidas y la variable dependiente z que es una función logística "perfecta" de x e y:

data logistica;
do i=1 to 10000;
x=rannor(8);
y=rannor(2);
prob=1/(1+exp(-(-10+5*x-5*y)));
z=ranbin(8,1,prob);
output;
end;
drop i;
run;

data entrenamiento test;
set logistica;
if ranuni(6)>0.8 then output test;
else output entrenamiento;
run;

proc freq data=entrenamiento;
tables z;
quit;

Separo los datos en entrenamiento y test y vemos que un 8% aproximadamente de mis registros tienen valor 1. Sobre estos datos hago una logística y una poisson y veo los parámetros Sigue leyendo Qué pasa si uso una regresión de poisson en vez de una regresión logística

Mapas municipales de Argentina con R

Municipios Buenos Aires

En respuesta a un lector del blog he elaborado de forma rápida una nueva entrada que nos permite realizar mapas por municipalidades para Argentina, ya hay entradas similares pero está bien que este mapa tenga su propia entrada para facilitar las búsquedas. El ejemplo es rápido y es probable que el código tenga algún fallo o error, si es así lo comentáis y lo solvento. Como es habitual nos dirigimos a la web del proyecto Global Administrative Areas (http://www.gadm.org/country) y nos descargamos el mapa de Argentina por municipios que es el nivel 2, una vez descargado pocas líneas de R:

library(sp)
library(RColorBrewer)

ub_argentina="C:\\mapas\\ARG_adm2.rds"

argentina = readRDS(ub_argentina)
b.aires = argentina[argentina$NAME_1=="Buenos Aires",]

aleatorio = rpois(length(unique(b.aires$NAME_2)),3)
b.aires@data=data.frame(aleatorio)
spplot(b.aires,col.regions=brewer.pal(5, "Blues"))

Leemos el mapa entero y hacemos un subconjunto de datos sólo con los municipios de Buenos Aires, pintamos unos datos aleatorios con la función spplot y con la librería RColorBrewer podemos dar un formato más elegante a nuestro mapa. Saludos.

Con los vehículos autónomos me quedo sin trabajo (o no)

"Los actuarios y los analistas de pricing os vaís a quedar sin trabajo con la llegada de los vehículos autónomos", yo ya he escuchado esta frase. Como podéis ver en mi perfil de Linkedin (un poco obsoleto) mi familia y yo comemos gracias a mi trabajo en equipos de actuarios encargados de dar precios tanto de nuevo negocio como en el momento de la renovación para compañías que operan fundamentalmente con el ramo de Automóviles. Llevo muchos años separando el "AZAHAAAAR" de lo estadísticamente explicable mediante relaciones lineales y, aunque esté feo decirlo, se me da muy bien [anda que no ha escuchado gente este discurso] Pero qué pasa si los coches ya no dependieran de las personas, fueran todos prácticamente iguales, nos diera igual la zona de circulación y no existiera la posibilidad de explicar nada estadísticamente porque no hubiera margen de error. Qué pasaría si los coches fueran autónomos.

Lo primero y más urgente, me quedan unos 25-30 años para jubilarme y tengo que repasar los principales hitos en la industria del vehículo autónomo para ver su evolución:

  • En 2013 National Highway Traffic Safety Administration's (como siempre en USA nos sacan ventaja) empieza a regular o a clasificar este tipo de vehículos: http://www.nhtsa.gov/About+NHTSA/Press+Releases/U.S.+Department+of+Transportation+Releases+Policy+on+Automated+Vehicle+Development
  • En 2014 algunos estados norteamericanos como California, Florida, Michigan, North Dakota y Tennessee ya crean legislación al respecto
  • En 2015 Tesla Motor Inc. incorpora a sus vehículos el Autopilot y parece que es el punto de inflexión.
  • En 2016 los que trabajamos en el sector asegurador nos empezamos a preocupar:
    • http://www.rand.org/pubs/research_reports/RR443-1.html
    •  https://www.adrianflux.co.uk/driverless-cars/
  • Ahora Tesla sale en las noticias debido a los accidentes y a su forma de financiarse.

Lo primero que se me pasa por la cabeza, en Europa a qué esperamos al respecto; cuándo van a empezar a legislar y lo segundo, parece que no se dan grandes avances o por lo menos no producen vértigo, a Google se le ha ocurrido crear su propia ciudad para ir probando pero como el camino sea construir ciudades vamos listos... Para un dinosaurio que sigue haciendo modelos lineales el vehículo autónomo le puede parecer ciencia ficción, sin embargo no podemos perder de vista que los coches afectan a la calidad de las ciudades y que los coches autónomos pueden mejorar mucho esa calidad.

En mi opinión los vehículos autónomos han venido para quedarse, necesitarán estar asegurados aunque se tratará de grandes flotas y creo que las empresas que los comercialicen deberán de tener sus propios seguros, así que es probable que tenga menos trabajo en 10-15 años pero también se me dan muy bien y tengo muchas ideas para los seguros multirriesgo y en salud y ahí hay mucho camino por andar, empezando por el seguro del Hogar que debería ser obligatorio así que no me preocupo mucho (de momento).  Saludos.

 

El #brexit con #rstats o como mover spatial data con R

bye_england

Animación con R para ilustrar el uso de la función de maptools elide de código "insultantemente" sencillo:

library(maptools)
library(animation)

#Mapa descargado de:
#http://www.arcgis.com/home/item.html?id=6d611f8d87d54227b494d4c3becef6a0

ub_shp = "/Users/raulvaquerizo/Desktop/R/mapas/world/MyEurope.shp"
europa = readShapeSpatial(ub_shp)
plot(europa)

europa_sin_uk = europa[europa$FIPS_CNTRY != "UK",]
uk = europa[europa$FIPS_CNTRY == "UK",]

saveGIF(
for (i in seq(0,5,by=0.1)){
plot(europa_sin_uk)
uk = elide(uk,shift=c(-i,1))
plot(uk,add=TRUE)},
interval=.3,
movie.name="/Users/raulvaquerizo/Desktop/R/animaciones/brexit/bye_england.gif"
)

Nos descargamos el mapa del link que os pongo y poco más que leer el shape con readShapeSpatial y crear dos objetos uno con Europa sin la isla y otro con la isla, elide nos permite desplazar un objeto de spatial data dentro del gráfico y lo metemos en un bucle y bye England. Ahora el que me mueva las Canarias en un shape file con más de 2 líneas de código me paga una cervecita. Saludos.

El paquete de R weatherData para la obtención de datos meteorológicos en España

Tenía pendiente un proyecto con modelos de Lee Carter y el paquete weatherData de R (¡toma!) pero como no lo voy a llevar a cabo nunca os traigo a estas líneas un paquete más que interesante de R que nos permite obtener datos de las estaciones meteorológicas de los aeropuertos del mundo (https://www.wunderground.com/history/airport/) y encima te lo pone como un objeto de R, qué más podemos pedir. En github tenéis una completa batería de ejemplos de uso. En el caso de que necesitemos descargar información meteorológica de España tenemos que irnos a http://weather.rap.ucar.edu/surface/stations.txt donde están listados todos los aeropuertos que recoge este sistema de información, buscamos SPAIN y nos interesa el "ICAO" que es el International Civil Aviation Organization, el código del aeropuerto vamos. Con estas premisas si quiero recoger las temperaturas de 2015 del aeropuerto de Albacete:

install.packages("weatherData")
library(weatherData)

anio = getWeatherForYear("LEAB",2015)

La información de los aeropuertos es bastante completa y están distribuidos por toda la geografía nacional además tenemos muchos datos a nuestra disposición y bien tabulados. Yo conozco más de uno que se está acordando de mi por no haber escrito sobre este paquete unos meses antes. Saludos.

Test de bondad de ajuste con SAS

Pregunta que me han hecho hoy. Cómo hacer un test de bondad de ajuste con SAS y la respuesta que he dado:

data datos_aleatorios;
do i=1 to 200000;
*GENERAMOS UNAS VARIABLES ALEATORIAS;
variable_gamma = rangam(89,450); 
variable_exponencial = ranexp(23)*100+0.17045;
output;
end;
run;

*ods select ParameterEstimates GoodnessOfFit ;
proc univariate data=datos_aleatorios;
   var var:;
   histogram /   gamma;
run;

Mucho cuidado con estos test de hipótesis. Yo suelo conformarme con ver la tabla de cuantiles. Saludos.

 

Truco SAS. Como leer PC Axis con SAS

Estoy leyendo información del INE que tiene que terminar cargándose en SAS y estos datos están en formato PC Axis. Existen macros en SAS para generar datasets a partir de PC Axis pero la verdad es que no he llegado a entender muy bien como funcionan y tras varios errores la mejor opción que he encontrado es emplear R y el paquete pxR que han creado algunos miembros de la Comunidad de R-Hispano. Como realizo esta tarea es más que sencillo:

En R realizamos la importación del archivo *.px:

nacionalidad = read.px("ubicacion\\seccion_censal_nacionalidad.px")
nacionalidad = data.frame(nacionalidad)
write.csv( nacionalidad, file = "ubicacion\\nacionalidad.csv" )

Hemos generado un csv que importamos desde SAS:

proc import datafile="ubicacion\nacionalidad.csv"
     out=nacionalidades
     dbms=csv
     replace;
     getnames=yes;
run;

También quería aprovechar esta entrada para comentaros que es preferible usar los viejos csv para mover archivos entre R  y SAS que usar librerías como SASxport que generan ficheros "portables" de SAS, aunque los ficheros "portables" garantizan que se puedan leer con distintas versiones de SAS este paquete tarda mucho (demasiado) tiempo en crear los archivos. Y si alguien tiene una versión más sencilla de la macro de SAS que mande el link. Saludos.

Cartografía digitalizada de España por sección censal

Por si no lo sabéis tenemos disponible en la web del INE un mapa de España por sección censal que podéis descargaros y realizar mapas con R de una forma que es más que conocida para los lectores del blog:

#mapas con secciones censales
library(maptools)
ub_shp = "/Users/raulvaquerizo/Desktop/R/mapas/cartografia_censo2011_nacional/SECC_CPV_E_20111101_01_R_INE.shp"
seccion_censal = readShapeSpatial(ub_shp)
barcelona = seccion_censal[seccion_censal$NMUN=="Barcelona",]
plot(barcelona)

Barcelona_mapa_seccion_censal

A ver si me animo y preparo una BBDD para que podáis acceder desde QGIS a una serie de mapas como este, además de los mapas por código postal. Aunque necesitaría un poco de ayuda técnica (ahí lo dejo). Saludos.
 

Nuevo y muy mejorado mapa de España por provincias con Excel

Nuevo_mapa_españa1

Hacía tiempo que no publicaba un mapa de España de Excel, aquí tenéis una nueva versión que mejora mucho  a las anteriores. La primera mejora y la que más destaca es que nos permite incluir datos, además ponemos los nombres de las provincias para todos aquellos que dominen poco la geografía española. Podemos pintar hasta 4 datos distintos que se pueden seleccionar en el desplegable que tenéis arriba. Ahora los colores van en dos escalas que podéis seleccionar vosotros:

Nuevo_mapa_españa_excel2

A la hora de meter los datos a nivel provincial es necesario ir a la hoja datos_mapa en ella tenéis los 4 datos que podéis pintar, estos datos irán en un ranking que  a la postre asigna colores a los shapes que componen el conjunto de imágenes que hace el mapa de Excel Sigue leyendo Nuevo y muy mejorado mapa de España por provincias con Excel

Adyacencia de polígonos con el paquete spdep de R

Cuando trabajamos con zonificación o geolocalización la adyacencia entre los elementos del estudio es relevante. En este caso quería trabajar con la adyacencia entre los polígonos que componen un archivo de datos espaciales shapefile y para entender mejor como podemos obtener la adyancecia entre polígonos creo que lo mejor es hacer un ejemplo con un mapa, en este caso un mapa de municipios de Barcelona. El primer paso es disponer del objeto con los datos espaciales, de esto ya he escrito mucho en el blog y por eso no me detengo mucho:

ub="./Desktop/R/mapas/ESP_adm4.rds"

#Creamos los objetos de R
espania = readRDS(ub)
barcelona = espania[espania$NAME_2=="Barcelona",]
plot(barcelona)
#Marcamos el centro de cada poligono
points(coordinates(barcelona))

adyacencia poligonos con R 1

Leemos el objeto con los datos municipales de España y hacemos un subset para quedarnos sólo con Barcelona y realizamos un mapa municipal de la provincia de Barcelona sencillamente usando plot. Podemos identificar todos los centroides de los polígonos que componen este objeto con la función coordinates, ahora lo que necesitamos identificar es la adyacencia entre estos puntos, la adyacencia entre los municipos de Barcelona. En mi caso localicé el paquete spdep de R, muy adecuado para trabajar con ponderaciones.

Os pongo paso por paso el código de R y luego comento como voy buscando REFERENCIAS para crear las adyacencias Sigue leyendo Adyacencia de polígonos con el paquete spdep de R